Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Problem Statement:
To solve for 2D heat conduction equation for steady state and unsteady state by explicit and implicit method with given boundary conditions.
Using implicit method
Solve using different iterative solvers i.e., Jacobi, Gauss seidel and Successive over relaxation for both steady and transient state heat conduction equation.
The accuracy of the solution obtained in both explicit and implicit method should be in accordance with the tolerance value mentioned.
Boundary conditions:
Top boundary =600K
Bottom boundary =900K
Left boundary =400K
Right boundary = 800K
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 statisfied leads to stable solution. Because of the stability criteria small time steps has to be taken to obtain stable solutions. With decrease in time step the number of computation increases and hence the computational time required for the solution to converge is less and accuracy of the solution is high.
Implicit Method:
When a computation is amde for dependent variable are defined by set of equtions and a matrix is involved to find the unknown variable. This computation is called implicit method. This method is unconditionally stable as a result, no stability criteria is required to obtain a stable results and the computational time required for the solution to converge is high and accuracy of the solution is less relative to explicit method because this method uses simultaneous equtions and uses guess values as the initial values to find the approximated values. If the guess values are closer to the desired value the computation time is less. As there is no stability criteria requirement larger time steps can be used.
In order to solve for dependent variable iterative solvers are used. Iterative solvers used are:
1. Jacobi iterative solver
2. Gauss seidel iterative solver
3. Successive over relaxation iterative solver.
2D Heat conduction Equation:
`(delT)/(delt) = alpha*((del^2T)/(delx^2) +(del^2T)/(dely^2))`
The time derivative term `(delT)/(delx)` is hyperbolic in nature which can be discretized using forward differencing scheme and the spacial derivative are elliptical in nature which can be discretized using central differencing scheme.
Discretization of time derivative using forward difference scheme since the time derivative is varying in time direction which can be solved by marching method.
`((delT)/(delt))_(i,j)^n = (T_(i,j)^(n+1) - T_(i,j)^n)/(Deltat) + O(Deltat)`
Discretization of spacial derivatives using central differencing scheme since the spacial derivatives are elliptical in nature, which uses information from neighbouring points to obtain a solution. Below are the discretization for x and y spacial derivatives.
`((del^2T)/(delx^2))_(i,j)^n = (T_(i+1,j)^n - 2*T_(i,j)^n + T_(i-1,j)^n)/(Deltax)^2 + O(Deltax^2)`
`((del^2T)/(dely^2))_(i,j)^n = (T_(i,j+1)^n - 2*T_(i,j)^n + T_(i,j-1)^n)/(Deltay)^2 + O(Deltay^2)`
From the PDE heat conduction equation,
`(T_(i+1,j)^n - 2*T_(i,j)^n + T_(i-1,j)^n)/(Deltax)^2 + (T_(i,j+1)^n - 2*T_(i,j)^n + T_(i,j-1)^n)/(Deltay)^2 = (1/alpha)*(T_(i,j)^(n+1) - T_(i,j)^n)/(Deltat)`
where n is current time step and n+1 is next time step.
i,j+1 are the neighbouring points in x direction, j, j+1 are the neighbouring points in y direction.
`Deltax` is grid spacing in x direction. `Deltay`is grid spacing in y direction.
Iterative Solver:
1. Jacobi iterative solver:
In this method the temperature at (n+1)th iteration is found using old values of the temperature defined at nth iterations. Using the approximation obtained iterations procedure it continued untill the desired accuracy of the solution is obtained.
2. Gauss Seidal iterative solver:
In this type iterative solver, 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, Gauss seidal require one storage vector which can be overwritten with each computation.
3. Successive over relaxation iterative solver.
In this type of iterative solver relaxation factor is used which helps in increasing the convergence rate. It goes by the formula
`T^(n+1) = T^n*(1-w) + w*T(GS) or T(jacobi)`
where w is between O-1 it is called under relaxed state.
This condition is used when the temperature gradient is low and due to this condition larger time steps are used to obtain the steady state i.e., for the values to converge at faster rate.
`T^n` is the value of temperature for previous iterations.
T(GS) or T(jacobi) is the temperature of value for (n+1)th iteration using gauss seidel or jacobi iterative solver.
Steady State 2D Heat Conduction Equation:
Steady state is where the fluid properties does not change with respect to time.
From the 2D heat conduction equation since the system is in steady state time derivative is zero and diffusivity is constant.
`(delT)/(delt) = 0`
Hence we the equation for steady state as:
`(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`
Discretization for steady state:
`(T_L - 2*T_P + T_R)/(Deltax^2) + (T_T - 2*T_P +T_B)/(Deltay^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 = (Deltay^2)/(2*(Deltax^2 + Deltay^2))`
`K_2 = (Deltax^2)/(2*(Deltax^2 + Deltay^2))`
Solving the 2D Steady heat conduction equation 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-w) + w*(K_1(T_L+T_R)^(n+1) + K_2(T_R + T_B)^(n+1))`
Transient state:
Unsteady state is where the fluid properties change with respect to time.
The formula for 2D transient heat conduction equation is:
`(delT)/(delt) = alpha*((del^2T)/(delx^2) +(del^2T)/(dely^2))`
Discretization for Unsteady state heat conduction equation:
`(T_P^(n+1) - T_P^n)/(Deltat) - alpha*((T_L -2*T_P + T_R)/(Deltax^2) + (T_T -2*T_P + T_B)/(Deltay^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 2D transient heat conduction equation explicitly we get,
`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 * (Deltat)/(Deltax^2)`
`K_2 =alpha * (Deltat)/(Deltay^2)`
Stability Criteria:
`K_1 + K_2 <= 0.5` Stability criteria for 2D heat conduction when solving explicitly.
Solving for 2D Transient Heat Conduction equation implicitly we get,
`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 * (Deltat)/(Deltax^2)`
`K_2 =alpha * (Deltat)/(Deltay^2)`
here n+1 is the next step and n is the preceeding 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
`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
`T_P^(n+1) = T_P^n(1-w) +w*[(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 Condition:
Boundary conditions are defined at the boundaries of independent variable of the function. In this problem mesh size of 10*10 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 10*10 shown above the corner boxes are average of the two boundaries intersecting at that point. The values other than the boundaries are takenas guess values which is in this case is 300K(initial condition for transient heat conduction equation).
Main Script:
%code to study the behaviour of steady state and transient using explicit
%and implicit method
clear all
close all
clc
%defining spacial domain
nx = 10;
ny = 10;
%specifying length of the 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);
%BC
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
T(1,1)=(T(1,2)+T(2,1))/2;
T(1,end)=(T(1,(end-1))+T(2,end))/2;
T(end,1)=(T(end,2)+T((end-1),1))/2;
T(end,end)=(T((end-1),end)+T(end,(end-1)))/2;
%Making a copy of the initial temperature matrix
Told = T;
T_pre=T;
%tolerance
tol=1e-4;
%Coefficient obtained from the derivation for steady state equation
k = 2*((dx^2+dy^2)/(dx^2*dy^2));
%diffusivity coefficient in the 2D heat conduction equation
alpha = 1.4;
%Calculation of time step
dt=(0.5)/((alpha/dx^2)+(alpha/dy^2));
%Co-efficient multiplied for x diffuse term
k1 = dt*alpha/dx^2;
%co-efficient multiplied for y diffuse term
k2 = dt*alpha/dy^2;
%Relaxation factor for successive over relaxation factor for steady state
%heat conduction
a=1.5;
%Relaxation factor for successive over relaxation factor for transient
%heat conduction
a1 = 1.15;
%Function to calculate the temperature os steady state using jacobi
%iterative solver
[error_j,iter1] = jacobi(nx,ny,T,Told,k,dx,dy,tol);
%Function to calculate the temperature os steady state using guess seidal
%iterative solver
[error_gs,iter2] = gs(nx,ny,T,Told,k,dx,dy,tol);
%Function to calculate the temperature os 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 temperature of transient heat conduction explicitly
[error_dt,iter_dt] = transient_heat_conduction(nx,ny,k1,k2,Told,T,tol);
%Function to calculate the temperature 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 temperature 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 temperature of unsteady state using successive
%over relaxation 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 the main Script:
- Specifying the number of grid points in the spacial domain x and y direction. In this case is 10 grid points in x and y direction.
- Specifying the length of the domain in x and y space domain
- Defining the spacial domain for the length and grid points specified.
- Calculating the grid spacing in the x and y space domain
- Defining temperature with the guess values 300K for the number of the grid points defined
- Applying the boundary condition as mentioned in the problem
- The accuracy of the solution is measured with computed approximated solution and compared with the tolerance value so that the converged solution is in the desired accuracy.
- A copy of the initial temperature matrix is made inTold and T_pre which is used while calculating the temperature values of n+1 time steps explicitly and implicitly
- Calculating `K, K_1,K_2`. here `alpha`is taken 1.4
- Function is called to calculate the temperature values at n+1 time step for both staedy and unsteady state using different iterative solvers and explicit method is also used for calculating temperature for unsteady state heat conduction equation.
Explanation of Functions defined in the main script:
In the functions dummy error is assigned to error variable so that the while loop condition would satisfy and iterations are calculated.
The error is calculated and updated for each iteration. The error calculated is coded to take only the absolute values since negative values will not statisfy the while loop condition and iterations will stopped and the results obtained will not be a converged results.
Number of iterations taken for the solution to converge is stored in the iter variable. This is done by increamenting the iter variable by 1 for each iteration.
As the iteration proceeds to calculate temperature at n+1 time step temperature at nth time step is required and these temperature values should be updated for each iteration. If the temperature values are not updated for each iteration the solution will not converge and while loop will run infinite times.
Contourf plot is used to check if the solutions are converged or not i.e., the plot does not change as the approximated temperature value changes in the order of tolerance value defined.
Told stores the temperature values of the preceeding time step.
T stores the temperature values of the n+1 time step.
This logic is used in the below functions to calculate the temperature of Steady state Heat Conduction equation for the desired accuracy.
Functions are as follows:
Function to calculate the temperature of Steady state using jacobi iterative solver
%Function to calculate the temperature 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)
%iterative loop
for i = 2:(nx-1)
for j = 2:(ny-1)
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
error1 = (abs(Told-T))
error_j = max(error1,[],'all')
Told=T
iter1=iter1+1
end
figure(1)
[C,h] = contourf(T)
clabel(C,h)
axis ij
title('number of iteration for jacobi method for steady state heat conduction',iter1)
end
This functions uses jacobi iterative solver, temperature at n+1 time step is calculated using temperature at nth time step.
Function to calculate the temperature of Steady state using Gauss Seidal iterative solver
%Function to calculate the temperature of Steady state using Gauss Seidal iterative solver
function [error_gs,iter2] = gs(nx,ny,t,told,k,dx,dy,tol)
% Dummy error value
error_gs = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter2 = 0;
%convergence loop
while (error_gs>tol)
%Iteration loop
for i = 2:(nx-1)
for j = 2:(ny-1)
%Formula to calculate the temperature at each grid
%point
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));
end
end
%error for each iteration
error1 = abs(told-t);
error_gs = max(error1,[],'all');
%updating the temperature 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 reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij;
title('Number of iterations for Gauss seidal method steady state heat conduction',iter2)
end This function uses Gauss seidal iterative solver in which temperature at n+1 time step is calculated using n+1 approximated temperature in the previous iterations.
Function to calculate the temperature of Steady state using Successive over relaxation iterative solver
%Function to calculate the temperature of Steady state using Successive over relaxation iterative solver
function [error_sor_gs,iter3] = sor_gs(nx,ny,t,told,k,dx,dy,tol,a)
%Dummy error value
error_sor_gs = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter3 = 0;
%Convergence loop
while (error_sor_gs>tol)
%iteration loop
for i = 2:(nx-1)
for j = 2:(ny-1)
%Formula to calculate the temperature 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
%temperature 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 temperature 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 reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij;
title('number of iteration for Successive over relaxation method for steady state heat conduction', iter3)
end
This function uses Successive over relaxation to calculate the temperature at n+1 time step. It uses temperature values calculated using Gauss seidal method and relaxation factor which helps in increasing the convergence rate.
Functions mentioned below uses Unsteady state head conduction equation formula to calculate the temperature values at n+1 time step.
Function to calculate the temperature of Unsteady state using explicit method
%Function to calculate the temperature of Unsteady state using explicit
%method
function [error_dt,iter_dt]=transient_heat_conduction(nx,ny,k1,k2,told,t,tol)
%dummy error
error_dt = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_dt =0;
%convergence loop
while (error_dt>tol)
%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 represent temperature values in x direction
V = told(i,j-1)-2*told(i,j)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature at each grid
%point
t(i,j)= told(i,j)+k1*H+k2*V;
end
end
%error for each iteration
error1 = abs(told-t);
error_dt = max(error1,[],'all');
%updating the temperature 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_dt = iter_dt+1;
end
%Plot representing the converged results.
figure(4)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij
title('Number of iterations using transient heat conduction explicitly',iter_dt);
endFunction to calculate the temperature in Unsteady state using jacobi iterative solver
As we are now solving for the unsteady state heat conduction equation, steady state is achieved after a certain number of time steps.
Hence 1st while loop is used to check for the steady state of the solution by calculating the error between temperature between successive time steps and comparing it with the tolerance values. At steady state the temperature values do not change with increases in number of iterations and error value reduces below the tolerance value . With this the converge point is found out.
2nd while loop functions in the similar way as 1st while loop. But in 2nd while loop it finds the converged value of temperature for the particular time step.
This logic is used in the below functions which are used to solve 2D heat conduction equation implicitly.
%Function to calculate the temperature 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)
% dummy error
error_jac=9e9;
error_f = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_jac=0;
%While loop to stop iterations when the temperature values converge
%whle 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
%Co-efficient obtained while deriving for transient heat conduction equation
term1= (1+(2*k1)+(2*k2))^-1;
%k1 and k2 represent the cfl value for x and y
%domain repectively.
H = told(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
V = told(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature 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(error1,[],'all');
%updating the temperature 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 the convergenc rate.
error_f = max(max(abs(t_pre-t)));
t_pre = t;
%Reseting the error value for the covergence loop to compute for the
%next time step.
error_jac = 9e9;
end
%Plot representing the converged results.
figure(5)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h)
axis ij
title('Number of iterations for implicit Jacobi',iter_jac);
end
Function to calculate the Temperature in Unsteady state using Gauss seidal iterative solver.
%Function to calculate the temperature 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)
%Definition of dummy error
error_f = 9e9;
error_gs_imp=9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_gs_imp=0;
%While loop to stop iterations when the temperature values converge
%whle 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
%Co-efficient obtained while deriving for transient heat conduction equation
term1= (1+(2*k1)+(2*k2))^-1;
%k1 and k2 represent the cfl value for x and y
%domain repectively.
H = t(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
V = t(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature 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 temperature 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 the convergenc rate.
error_f = max(max(abs(t_pre-t)));
%Updating the previous temperature values to computed values.
t_pre = t;
%Reseting the error value for the covergence loop to compute for the
%next time step.
error_gs_imp = 9e9;
end
%Plot representing the converged results.
figure(6)
[C,h]=contourf(t);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values 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 in Unsteady state using successive over relaxation iterative solver
%Function to calculate the temperature of Unsteady state using successive
%over relaxation 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)
%Definition of dummy error
error_f = 9e9;
error_sor_gs_imp = 9e9;
%defining iteration variable to store and display the number of
%iterations computed for the values to converge.
iter_sor_gs_imp =0;
%While loop to stop iterations when the temperature values converge
%whle 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
%Co-efficient obtained while deriving for transient heat conduction equation
term1= (1+(2*k1)+(2*k2))^-1;
H = t(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
V = t(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
%Formula to calculate the temperature at each grid
%point
t(i,j)= (t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
%Successive over relaxation Formula to calculate the
%temperature at each grid point.
t(i,j)= (told(i,j))*(1-a1)+ a1*t(i,j);
end
end
%error for each iteration
error1 = abs(told-t);
error_sor_gs_imp = max(error1,[],'all');
%updating the temperature 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 the convergenc rate.
error_f = max(max(abs(t_pre-t)));
%Updating the previous temperature values to computed values.
t_pre = t;
%Reseting the error value for the covergence loop to compute for the
%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 reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h)
axis ij
title('Number of iteration for implicit SOR',iter_sor_gs_imp);
end
Output Graph:
The above plots shows number of iterations required by the iterative solver to obtain the desired accurate solution.
Inference from the Plots:
Successive over relaxation iterative solver convergence rate is higher than the jacobi and gauss seidal method for both steady and unsteady state heat conduction equation. The number of iterations required by the successive over relaxation method is less compared to jacobi and gauss seidal method.
Below Table gives the number of iterations required by iterative solvers to converge.
| Steady State Heat Conduction Equation | |
| Iterative Solver | No. of iterations |
| Jacobi | 207 |
| Gauss Seidal | 111 |
| Successive over Relaxation | 28 |
| Unsteady State Heat Conduction Equation. | |
| Iterative Solver | No. of iterations |
| Jacobi | 1982 |
| Gauss Seidal | 1458 |
| Successive over Relaxation | 1148 |
Conclusion:
Implicit method can be used to problems when converged results are needed in less computation time as they are unconditionally stable. With increases in time step size the convergence rate is increased which leads to steady state with less no. of iterations.
Explicit method results time accurate solutions as it is governed by stability criteria when not satisfied it leads to unstable solution and requires more iterations to attain steady state.
But in the code provided above time step size remains constant which is calculated with respect to the stability criteria. Hence for this time step size the explicit scheme leads to faster convergent rate and as the time step size is small the number of iterations increases for implicit scheme which results slower convergence rate.