Solving the Steady and Unsteady 2-D heat conduction problem on MATLAB
Objective:
- Steady-state analysis & Transient State Analysis
Solve the 2D heat conduction equation by using the point iterative techniques.
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
implement the following methods for solving Steady-State equations:
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation
Absolute error criteria is 1e-4
- Write code to solve the Transient-State equation using Explicit and Implicit schemes.
- Steady-state vs Unsteady analysis
- Compare fastest steady-state simulation vs fastest transient simulation and explain the results
Introduction:
To solve the 2D heat conduction equation by using the point iterative or time marching techniques, the Boundary conditions for the problem is specified as follows:
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
The General governing partial differential equation for heat transfer is:
`(delT)/(delt)=alpha*{(del^2T)/(delx^2) + (del^2T)/(dely^2) + (del^2T)/(delz^2) }`
where T is the temperature, `alpha` is thermal diffusivity and x,y,z are space variables and t is time variable.
The governing partial differential equation for 2D steady-state heat conduction case is:
`(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`
Discretization of GDE:
The steady-state GDE is elliptical in nature as it depends on the information from all directions in the mesh therefore, Central Differencing Scheme is used to discretize this second-order partial differential equation.
Let `T_(i,j)^m text( be the temperature at ) i^(th) text( row & )j^(th) text( column in the assumed matrix field at )m^(th) text( iteration.)`
This equation can be discretized as:
`(T_(i+1,j) -2*T_(i,j) +T_(i-1,j))/(Deltax^2) + (T_(i,j+1) -2*T_(i,j) +T_(i,j-1))/(Deltay^2)=0`
This equation can be rewritten to simplify as:
`(T_R -2*T_P +T_L)/(Deltax^2) + (T_T -2*T_P +T_B)/(Deltay^2) =0`
Where:
`T_P= T_(i,j) text( (at point p))`
`T_R= T_(i+1,j)text( (at right point))`
`T_L= T_(i-1,j)text( (at left point))`
`T_T= T_(i,j+1)text( (at top point))`
`T_B= T_(i,j-1)text( (at bottom point))`
On further simplification, we get:
`T_P=1/K*{(T_R+T_L)/(Deltax^2) + (T_T +T_B)/(Deltay^2)}` --------------------------------> (1)
Where:
`K= 2*(Deltax^2 +Deltay^2)/(Deltax^2*Deltay^2)`
(NOTE: all temperatures are at same time step for all steady-state equations)
Equation 1 gives the temperature at a point as a function of temperatures at neighboring points.
This equation can be solved by using fixed-point iterative solving methods like Jacobi, Gauss-Siedal, Successive Over Relaxation on Gauss-Siedal. To study the comparison between these solvers, we first derive there solving equations and then compare the solution, the number of iterations taken to converge on MATLAB.
Jacobi Method for 2D steady-state heat conduction:
In the Jacobi method, the values of `T_p` at `(n+1)^(th)` iteration is calculated based on surrounding temperatures at `n^(th)` iteration until the solution converges. This method is also known as the method of Simultaneous displacements as all values at `n^(th)` iteration are updated with the values at `(n+1)^(th)` iteration for each iteration until the solution converges.
The formula for the Jacobi Method can be written as:
`T_P^(m+1)=1/K*{(T_R^(m)+T_L^(m))/(Deltax^2) + (T_T^(m) +T_B^(m))/(Deltay^2)}`
Gauss-Seidel Method for 2D steady-state heat conduction:
In the Gauss-Siedal method, the values of `T_p` at `(n+1)^(th)` iteration is calculated based on surrounding temperatures at `(n+1)^(th)` iteration if available else temperatures at `n^(th)` iteration until the solution converges. This method is also known as the method of Successive displacements as all values at `(n+1)^(th)` iteration are updated in the formula as soon as they are computed for each iteration until the solution converges. This technique helps in faster convergence of the solution
The formula for the Gauss-Siedal Method can be written as:
`T_P^(m+1)=1/K*{(T_R^(m)+T_L^(m+1))/(Deltax^2) + (T_T^(m) +T_B^(m+1))/(Deltay^2)}`
SOR on Gauss-Seidel Method for 2D steady-state heat conduction:
The SOR method is an advanced method that is used to scale the solution of other iterative solvers. If the relaxation factor "Omega" `omega` is greater than 1, then the solution will converge fast, and if it's less than 1 then the solution will take more iterations to reach the desired solution. The formula of SOR applied on GS method is:
`T_P^(m) =(1-omega)* T_P^(m-1) + omega*(T_(GS))`
or
`T_P^(m+1)=T_P^(m)*(1-alpha)+alpha/K*{(T_R^(m)+T_L^(m+1))/(Deltax^2) + (T_T^(m) +T_B^(m+1))/(Deltay^2)}`
The governing partial differential equation for the unsteady-state case is :
`(dT)/(dt) = alpha*((del^2T)/(delx^2) + (del^2T)/(dely^2))`
Discretization of GDE:
The unsteady-state GDE is hyperbolic in nature as it depends on the information from left to the right direction in the mesh therefore, Forward Differencing Scheme is used to discretize this second-order partial differential equation.
Let `T_(i,j)^(n,m)text( be the temperature at ) i^(th) text( row & )j^(th) text( column in the assumed matrix field at )n^(th) text( time-step and )m^(th)text( iteration.)`
This equation can be discretized as:
`(T_P^(n+1) -T_P^(n))/(Deltat)= alpha*{{T_L^n-2*T_P^n+T_R^n}/(Deltax^2) +{T_B^n-2*T_P^n+T_T^n}/(Deltax^2) }`
on rearranging the equation, we get:
`T_P^(n+1) = T_P^(n) + (alpha*Deltat)/(Deltax^2)*{T_L^n-2*T_P^n+T_R^n} + (alpha*Deltat)/(Deltay^2)*{T_B^n-2*T_P^n+T_T^n}`
Let us consider `K1=alpha*(Deltat)/(Deltax^2)` and `K2=alpha*(Deltat)/(Deltay^2)`
By further simplification, we get:
`T_P^(n+1) = (1-2*K1-2*K2)*T_P^(n) + K1*{T_L^n+T_R^n} + K2*{T_B^n+T_T^n}`
Here, K1 and K2 are X and Y components of CFL number. The CFL number stands for Courant–Friedrichs–Lewy number. For a stable solution, the sum of components of CFL number should be equal to or less than 0.5. Therefore,
`K1=CFL_x`
`K2=CFL_y`
`CFL_x+CFL_y <=0.5`
`:. K1+K2<=0.5`
From this condition, the size of time step can be calculated by following formula:
`dt= (Deltay^2*Deltax^2)/(2*alpha*(Deltax^2+Deltay^2))`
Where:
`dx=X/(nx) text( ,X is size of body in x direction in meters & nx is no. of divisions in x direction)`
`dy=Y/(ny) text( ,Y is size of body in y direction in meters & ny is no. of divisions in y direction)`
This equation can be solved by using fixed-point iterative solving methods like Jacobi, Gauss-Siedal, Successive Over Relaxation on Gauss-Siedal. To study the comparison between these solvers, we first derive there solving equations and then compare the solution, the number of iterations taken to converge on MATLAB.
Explicit Method for 2D Unsteady-state heat conduction:
`T_P^(n+1) = (1-2*K1-2*K2)*T_P^(n) + K1*{T_L^n+T_R^n} + K2*{T_B^n+T_T^n}`
Implicit (Jacobi) Method for 2D Unsteady-state heat conduction:
In the Jacobi method, the values of `T_p` at `(n+1)^(th)` iteration is calculated based on surrounding temperatures at `n^(th)` iteration until the solution converges. This method is also known as the method of Simultaneous displacements as all values at `n^(th)` iteration are updated with the values at `(n+1)^(th)` iteration for each iteration until the solution converges.
The formula for the Jacobi Method can be written as:
`T_P^(n+1,m) = T_P^(n) /(K3)+ (K1)/(K3)*{T_L^(n+1,m-1)+T_R^(n+1,m-1)} + (K2)/(K3)*{T_B^(n+1,m-1)+T_T^(n+1,m-1)}`
`text(Where, ) K3=1-2*K1-2*K2`
Implicit (Gauss-Seidel) Method for 2D Unsteady-state heat conduction:
In the Gauss-Siedal method, the values of `T_p` at `(n+1)^(th)` iteration is calculated based on surrounding temperatures at `(n+1)^(th)` iteration if available else temperatures at `n^(th)` iteration until the solution converges. This method is also known as the method of Successive displacements as all values at `(n+1)^(th)` iteration are updated in the formula as soon as they are computed for each iteration until the solution converges. This technique helps in faster convergence of the solution
The formula for the Gauss-Siedal Method can be written as:
`T_P^(n+1,m) = T_P^(n) /(K3)+ (K1)/(K3)*{T_L^(n+1,m)+T_R^(n+1,m-1)} + (K2)/(K3)*{T_B^(n+1,m)+T_T^(n+1,m-1)}`
Implicit (SOR on Gauss-Seidel) Method for 2D Unsteady-state heat conduction:
The SOR method is an advanced method that is used to scale the solution of other iterative solvers. If the relaxation factor "Omega" `omega` is greater than 1, then the solution will converge fast, and if it's less than 1 then the solution will take more iterations to reach the desired solution. The formula of SOR applied on GS method is:
`T_P^(n+1,m) =(1-omega)* T_P^(n+1,m-1) + omega*(T_(GS))`
or
`T_P^(n+1,m) =(1-omega)* T_P^(n+1,m-1) + omega*(T_P^(n) /(K3)+ (K1)/(K3)*{T_L^(n+1,m)+T_R^(n+1,m-1)} + (K2)/(K3)*{T_B^(n+1,m)+T_T^(n+1,m-1)})`
`text(where, )omega text( is relaxation factor, if )omega>1 text( solution is over relaxed or if )omega<1 text( solution is under relaxed)`
(NOTE: for temperature values, Superscript of T denotes time step and subscript of T denotes the position of the point)
Initial Condition:
Initial conditions are visualized by creating a 10*10 matrix and assigning the wall temperatures as mentioned, all inside points are initially at 320K. A filled contour is drawn from this array to visualize temps at all points.
This matrix is visualized using "contourf" command in MATLAB, the rainbow color scheme is used to show the gradient of temperature:
Explanation of the Main program:
The main program is used to assign all initial, variable, and constant values. These values are then passed on to the respective functions to calculate and plot the results, final answers are retrieved from every function to draw a comparison. The comparison between initial conditions and other types of solvers, a subplot is used.
The bar graphs are used to compare the number of iterations and errors at the last iteration.
- figure(8) will show a 2*2 subplot of filled contours of the initial condition and final conditions of Steady-state solvers i.e. Jacobi, Gauss-Seidel, SOR on GS.
- figure(9) will show a 2*2 subplot of filled contours of the final conditions of Unsteady-state solvers i.e. Explicit, Jacobi, Gauss-Seidel, SOR on GS.
- figure(10) will show a bar plot of no. of iterations taken by Steady-state solvers to achieve error less than tolerance specified.
- figure(11) will show a bar plot of no. of iterations taken by Unsteady-state solvers to achieve steady-state.
Global Variables specified in the main program:
- No. of divisions in X & Y direction respectively
- nx=10
- ny=10
- `Deltax" & "Deltay` (dx & dy) is calculated by taking a difference between two consecutive points on X & Y direction respectively
- Thermal diffusivity constant `(alpha)` is set to a=1.5
- `Deltat` (dt) is time step size, which is calculated using values of `Deltax" & "Deltay` according CFL stability criteria
- Max. no. of the time step is specified as nt=1400
- Value of CFL numbers is calculated as K,K1,K2 using values of `Deltax" & "Deltay`
- error tolerance value for space loop of Steady-state solvers, tol=1e-8
- error tolerance value for space loop of Un-Steady state solvers, tol2=1e-4
- defining SOR factor: sor_factor=1.2
Main Program:
clear all
close all
clc
%This program calls iterative solver functions to solve a given problem
%Function 1 is to solve steady state problem using Jacobi iterative method
%Function 2 is to solve steady state problem using Gauss-Siedal iterative method
%Function 3 is to solve steady state problem using Succesive Over Relaxation(SOR) on GS iterative method
%Function 4 is to solve unsteady state problem using Explicit method
%Function 5 is to solve unsteady state problem using Implicit-JACOBI iterative method
%Function 6 is to solve unsteady state problem using Implicit-GS iterative method
%Function 7 is to solve unsteady state problem using Implicit-SOR-GS iterative method
nx=10; %no. of grid points in X-direction
ny=10; %no. of grid points in Y-direction
t_initial=320*ones(ny,nx); %initial condition for inner nodes
t_initial(2:ny-1,1)=400;%initial condition for boundary nodes
t_initial(1,2:nx-1)=600;%initial condition for boundary nodes
t_initial(2:ny-1,nx)=800;%initial condition for boundary nodes
t_initial(ny,2:nx-1)=900;%initial condition for boundary nodes
t_initial(1,1)=500;%initial condition for corner nodes
t_initial(1,nx)=700;%initial condition for corner nodes
t_initial(ny,1)=650;%initial condition for corner nodes
t_initial(ny,nx)=850;%initial condition for corner nodes
% creating grid
x=linspace(0,1,nx); % x grid points
dx=abs(x(2)-x(1)); %distance between 2 points in x direction
y=linspace(0,1,ny); % y grid points
dy=abs(y(2)-y(1)); %distance between 2 points in y direction
a=1.5; %thermal diffusivity
%calulating time step
dt= (0.5*(dy^2)*(dx^2))*(a*((dx^2)+(dy^2)))^(-1);
%number of time steps
nt=1400;
%final time
ftime=dt*nt;
k= 2*(dx^2+dy^2)/(dx^2*dy^2); % calculating value of K for steady state
k1= a*dt/dx^2;% calculating value of K1 for transient state
k2= a*dt/dy^2;% calculating value of K2 for transient state
%difining tolerence values
tol=1e-8; % error tolerance value for space loop of Steady state solvers
tol2=1e-4; %error tolerance value for space loop of Un-Steady state solvers
%defining sor factor
sor_factor=1.2;
% calling functions to solve the problem
[steady_jacobi,counter1,error1]= jacobi_func(t_initial,nx,ny,dx,dy,tol,k) %Function 1
[steady_gs,counter2,error2]= gs_func(t_initial,nx,ny,dx,dy,tol,k) %Function 2
[steady_sorgs,counter3,error3]= sorgs_func(t_initial,nx,ny,dx,dy,tol,k,sor_factor) %Function 3
[unsteady_exp,counter4,error4,ts1]= explicit_func(t_initial,nx,ny,k1,k2,nt,dt) %Function 4
[unsteady_jacobi,counter5,error5,ts2]= jacobi_unsteady(t_initial,nx,ny,k1,k2,tol2,nt,dt) %Function 5
[unsteady_gs,counter6,error6,ts3]= gs_unsteady(t_initial,nx,ny,k1,k2,tol2,nt,dt) %Function 6
[unsteady_sorgs,counter7,error7,ts4]= sorgs_unsteady(t_initial,nx,ny,k1,k2,tol2,nt,dt,sor_factor) %Function 7
% plotting a subplot to compare the results using a filled contour
figure(8)
subplot(2,2,1)
contourf(t_initial) % plotting initial conditions
colorbar; colormap(jet);
axis ij %to flip the Y-axis
axis equal
grid on
xlabel('X-axis')
ylabel('Y-axis')
title_text1=sprintf('Initial conditions');
title_text2=sprintf('Initial condition: Top 600K, Left 400K, Right 800K, Bottom 900K');
title({[title_text1],[title_text2]})
subplot(2,2,2)
contourf(steady_jacobi) %plotting final conditions obtained from jacobi method
colorbar; colormap(jet);
axis ij %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
title_text1=sprintf('Steady State-Jacobi Method, iteration number = %d',counter1);
title_text2=sprintf('temp. diff. = %.10f K',error1);
title({[title_text1],[title_text2]})
subplot(2,2,3)
contourf(steady_gs) %plotting final conditions obtained from gauss-siedal method
colorbar; colormap(jet);
axis ij %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
title_text1=sprintf('Steady State-GS Method, iteration number = %d',counter2);
title_text2=sprintf('temp. diff. = %.10f K',error2);
title({[title_text1],[title_text2]})
subplot(2,2,4)
contourf(steady_sorgs) %plotting final conditions obtained from SOR on gauss-siedal method method
colorbar; colormap(jet);
axis ij %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
title_text1=sprintf('Steady State-SOR+GS Method, iteration number = %d',counter3);
title_text2=sprintf('temp. diff. = %.10f K',error3);
title({[title_text1],[title_text2]})
figure(9)
subplot(2,2,1)
contourf(unsteady_exp)
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Explicit Method, Iteration number = %d',counter4);
title_text2=sprintf('Thermal Diffusivity= %.1f, timestep size = %.12f sec, timestep = %d, ',a,dt,ts1);
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',ftime,error4);
title({[title_text1],[title_text2],[title_text3]})
axis ij %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
subplot(2,2,2)
contourf(unsteady_jacobi)
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Implicit Method, Jacobi, iteration number = %d',counter5);
title_text2=sprintf('Thermal Diffusivity= %.1f, timestep size = %.12f sec, timestep = %d',a,dt,ts2);
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',ftime,error5);
title({[title_text1],[title_text2],[title_text3]})
axis ij %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
subplot(2,2,3)
contourf(unsteady_gs)
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Implicit Method,G-S, Iteration number = %d',counter6);
title_text2=sprintf('Thermal Diffusivity= %.1f, timestep size = %.12f sec, timestep = %d, ',a,dt,ts1);
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',ftime,error6);
title({[title_text1],[title_text2],[title_text3]})
axis ij %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
subplot(2,2,4)
contourf(unsteady_sorgs)
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Implicit Method, SOR on G-S, Iteration number = %d',counter7);
title_text2=sprintf('Thermal Diffusivity= %.1f, timestep size = %.12f sec, timestep = %d, ',a,dt,ts1);
title_text3=sprintf('at time = %.4f sec, error = %.12f K',ftime,error7);
title({[title_text1],[title_text2],[title_text3]})
axis ij %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
figure(10) % comparison of number of iteration and error in diffrent methods to solve steady-state problem
X1 = categorical({'Jacobi','Gauss-Siedal','SOR on GS'});
X1 = reordercats(X1,{'Jacobi','Gauss-Siedal','SOR on GS'});
data1=[counter1;counter2;counter3]
bar(X1,data1);
xlabel('Iterative Methods')
ylabel('Number of iteration required to get error < 1e-4')
legend('No.of iterations')
title('Steady-State solving methods comparison')
figure(11) % comparison of number of iteration and error in diffrent methods to solve un-steady state problem
X2 = categorical({'Explicit','Implicit-Jacobi','Implicit-GS','Implicit-SOR-GS'});
X2 = reordercats(X2,{'Explicit','Implicit-Jacobi','Implicit-GS','Implicit-SOR-GS'});
data3=[counter4;counter5;counter6;counter7]
bar(X2,data3,'g');
xlabel('Time marching methods')
ylabel({['Number of iteration done till error<10^(-8) or max Timestep=1400'],['(whichever is achieved first)']})
legend('No.of iterations')
title('Unsteady-State solving methods comparison')Steady-State Jacobi Method:
This function solves the steady-state problem using the Jacobi iterative method.
"t_jacobi" is initialized by copying initial temps, then a copy of "t_jacobi" is defined as "told_jacobi" to store the values at the previous iteration.
value of error is initialized as `9^9`.
An outer while loop is used to check the tolerance criteria, if the value of error goes less than the tolerance specified in the main program, it will stop the solver.
A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.
"told_jacobi" is updated with "t_jacobi" values after the space loop so that the old values changes after each iteration.
Error is calculated inside while loop as a maximum of difference in temperatures of all points between consecutive iterations.
A filled contour plot of "t_jacobi" is plotted within the while loop to visualize the convergence.
Function:
function [steady_jacobi,counter1,error1] = jacobi_func(t_initial,nx,ny,dx,dy,tol,k)
%Function 1 is to solve steady state problem using Jacobi iterative method
% d^2T/dx^2 + d^2T/dy^2 = 0
%term1 = (1/k)*(tl+tr)/(dx)^2)
%term2 = (1/k)*(tt+tb)/(dy)^2)
%tp = term1 + term2
t_jacobi=t_initial;
%creating a copy of t
told_jacobi=t_jacobi;
error=9e9; %defining arbitrary error value, error>tollerance
%convergence loop
counter_jacobi=0; % loop counter
while (error>tol) %tolerance checking for every space loop
for i=2:nx-1 %spatial loop x direction
for j=2:ny-1 %spatial loop y direction
term1 = (1/k)*(told_jacobi(i-1,j) + told_jacobi(i+1,j))/(dx^2);
term2 = (1/k)*(told_jacobi(i,j-1) + told_jacobi(i,j+1))/(dy^2);
t_jacobi(i,j)= term1 + term2; % calculating new temperature
end
end
error= max(max(abs(told_jacobi-t_jacobi))) %calculating error
told_jacobi=t_jacobi; % updating old temperatures
counter_jacobi=counter_jacobi+1;
figure(1)
contourf(t_jacobi)
axis ij
grid on
colorbar; colormap(jet);
title_text1=sprintf('Steady State-Jacobi Method, iteration number = %d',counter_jacobi);
title_text2=sprintf('Max. Temp. diff. b/w consecutive iterations = %.10f',error);
title({[title_text1],[title_text2]})
pause(0.009)
if (error<1e-8)
break
end
end
steady_jacobi=t_jacobi;
counter1=counter_jacobi;
error1=error;
endSteady-State Gauss-Siedal Method:
This function solves the steady-state problem using the Gauss-Seidel iterative method.
"t_gs" is initialized by copying initial temps, then a copy of "t_gs" is defined as "told_gs" to store the values at the previous iteration.
value of error is initialized as `9^9`.
An outer while loop is used to check the tolerance criteria, if the value of error goes less than the tolerance specified in the main program, it will stop the solver.
A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.
"told_gs" is updated with "t_gs" values after the space loop so that the old values changes after each iteration.
Error is calculated inside while loop as a maximum of difference in temperatures of all points between consecutive iterations.
A filled contour plot of "t_gs" is plotted within the while loop to visualize the convergence.
Function:
function [steady_gs,counter2,error2] = gs_func(t_initial,nx,ny,dx,dy,tol,k)
%steady state
%G-S method
% d^2T/dx^2 + d^2T/dy^2 = 0
%term1 = (1/k)*(tl+tr)/(dx)^2)
%term2 = (1/k)*(tt+tb)/(dy)^2)
%tp = term1 + term2
t_gs=t_initial;
%creating a copy of t
told_gs=t_gs;
error=9e9; %defining arbitrary error value, error>tollerance
%convergence loop
counter_gs=0; % loop counter
while (error>tol) %tolerance checking for every space loop
for i=2:nx-1 %spatial loop x direction
for j=2:ny-1 %spatial loop y direction
term3 = (1/k)*(t_gs(i-1,j) + t_gs(i+1,j))/(dx^2);
term4 = (1/k)*(t_gs(i,j-1) + t_gs(i,j+1))/(dy^2);
t_gs(i,j)= term3 + term4; % calculating new temperature
error= max(abs(told_gs(i,j)-t_gs(i,j))) %calculating error
end
end
told_gs=t_gs; % updating old temperatures
counter_gs=counter_gs+1;
figure(2)
contourf(t_gs)
axis ij
grid on
colorbar; colormap(jet);
title_text1=sprintf('Steady State-GS Method, iteration number = %d',counter_gs);
title_text2=sprintf('Max. Temp. diff. b/w consecutive iterations = %.10f',error);
title({[title_text1],[title_text2]})
pause(0.03)
if (error<1e-8)
break
end
end
steady_gs=t_gs;
counter2=counter_gs;
error2=error;
endSteady-State Successive Over Relaxation on Gauss-Siedal Method:
This function solves the steady-state problem using the SOR on the Gauss-Seidel iterative method.
"t_sorgs" is initialized by copying initial temps, then a copy of "t_sorgs" is defined as "told_sorgs" to store the values at the previous iteration.
value of error is initialized as `9^9`.
An outer while loop is used to check the tolerance criteria, if the value of error goes less than the tolerance specified in the main program, it will stop the solver.
A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.
"told_sorgs" is updated with "t_sorgs" values after the space loop so that the old values changes after each iteration.
Error is calculated inside while loop as a maximum of difference in temperatures of all points between consecutive iterations.
A filled contour plot of "t_sorgs" is plotted within the while loop to visualize the convergence.
Function:
function [steady_sorgs,counter3,error3] = sorgs_func(t_initial,nx,ny,dx,dy,tol,k,sor_factor)
%steady state
%SOR on G-S method
% d^2T/dx^2 + d^2T/dy^2 = 0
%term1 = (1/k)*(tl+tr)/(dx)^2)
%term2 = (1/k)*(tt+tb)/(dy)^2)
%tp = term1 + term2
t_sorgs=t_initial;
%creating a copy of t
told_sorgs=t_sorgs;
error=9e9; %defining arbitrary error value, error>tollerance
%convergence loop
counter_sorgs=0; % loop counter
while (error>tol) %tolerance checking for every space loop
for i=2:nx-1 %spatial loop x direction
for j=2:ny-1 %spatial loop y direction
term5 = (1/k)*(t_sorgs(i-1,j) + t_sorgs(i+1,j))/(dx^2);
term6 = (1/k)*(t_sorgs(i,j-1) + t_sorgs(i,j+1))/(dy^2);
term7 = (1-sor_factor)*told_sorgs(i,j);
t_sorgs(i,j)= term7 + sor_factor*(term5 + term6); % calculating new temperature
error= max(abs(told_sorgs(i,j)-t_sorgs(i,j))) %calculating error
end
end
told_sorgs=t_sorgs; % updating old temperatures
counter_sorgs=counter_sorgs+1;
figure(3)
contourf(t_sorgs)
axis ij
grid on
colorbar; colormap(jet);
title_text1=sprintf('Steady State-SOR+GS Method, iteration number = %d',counter_sorgs);
title_text2=sprintf('Max. Temp. diff. b/w consecutive iterations = %.10f',error);
title({[title_text1],[title_text2]})
pause(0.003)
if(error<1e-8)
break
end
end
steady_sorgs=t_sorgs;
counter3=counter_sorgs;
error3=error;
endUnsteady-State Explicit Method:
This function solves the Unsteady-state problem using the Explicit method.
"t_un_ex" is initialized by copying initial temps, then a copy of "t_un_ex" is defined as "told_unex" to store the values at the previous iteration.
value of error is initialized as `9^9`.
An outer FOR loop is used to change the time-step of the solver, once it reaches maximum timestep defined in the main program or if it converges before that, it will stop the solver.
A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.
Error is calculated inside the Y-direction space loop as a maximum of difference in temperatures of all points between consecutive iterations.
"told_unex" is updated with "t_un_ex" values after the space loop so that the old values changes after each iteration.
A filled contour plot of "t_un_ex" is plotted within the while loop to visualize the convergence.
an IF loop is used to check for the convergence.
Function:
function [unsteady_exp,counter4,error4,ts1] = explicit_func(t_initial,nx,ny,k1,k2,nt,dt)
%unsteady state
%trnsient state solution
%explicit method
% dT/dt = a*((d^2T/dx^2)+(d^2T/dy^2))
% (tpnew - tpold)/dt = a*((tl+2*tp+tr)/dx^2 + (tt-2*tp+tb)/dy^2)
t_un_ex=t_initial;
%creating a copy of t
told_unex=t_un_ex;
%convergence loop
counter_unex=0; % loop counter
for n=1:nt %time loop
error=9e9; %error re-initialisation for every time step
current_time = dt*n;
for i=2:nx-1 %spatial loop x direction
for j=2:ny-1 %spatial loop y direction
term1 = k1*(told_unex(i-1,j) + told_unex(i+1,j) - 2*told_unex(i,j));
term2 = k2*(told_unex(i,j-1) + told_unex(i,j+1) - 2*told_unex(i,j));
t_un_ex(i,j)= told_unex(i,j) + term1 + term2; % calculating new temperature
error= max(abs(told_unex(i,j)-t_un_ex(i,j))) %calculating error
end
end
counter_unex=counter_unex+1;
told_unex=t_un_ex; % updating old temperatures
figure(4)
contourf(t_un_ex)
axis ij
grid on;
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Explicit Method, iteration number = %d',counter_unex);
title_text2=sprintf('Time Step = %d, current time = %.4f',n,current_time);
title_text3=sprintf('Temp. diff. b/w time steps = %.10f',error);
title({[title_text1],[title_text2],[title_text3]});
pause(0.003)
if(error<1e-8)
break
end
end
unsteady_exp=t_un_ex;
counter4=counter_unex;
error4=error;
ts1=n;
endUnsteady-State Jacobi Method:
This function solves the Unsteady-state problem using the Implicit method using the Jacobi iterative solver.
"t_un_imp" is initialized by copying initial temps, then a copy of "t_un_imp" is defined as "told_unimp" to store the values at the previous iteration.
value of error is initialized as `9^9`.
An outer FOR loop is used to change the time-step of the solver, once it reaches maximum timestep defined in the main program or if it converges before that, it will stop the solver.
An outer while loop is used to check the tolerance criteria, if the value of error goes less than the tolerance specified in the main program, it will stop the space solver
A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.
Error is calculated inside the Y-direction space loop as a maximum of difference in temperatures of all points between consecutive iterations.
"told_unimp" is updated with "t_un_imp" values after the space loop so that the old values changes after each iteration.
A filled contour plot of "t_un_imp" is plotted within the while loop to visualize the convergence.
an IF loop is used to check for the convergence.
function[unsteady_imp,counter5,error5,ts2]= jacobi_unsteady(t_initial,nx,ny,k1,k2,tol2,nt,dt)
%unsteady state
%trnsient state solution
%implicit method, Jacobi iterative solver
% dT/dt = a*((d^2T/dx^2)+(d^2T/dy^2))
t_un_imp=t_initial;
%creating a copy of t
told_unimp=t_un_imp;
tprev=t_un_imp;
k3= 1/(1+2*k1+2*k2);
count=0;
%convergence loop
for n=1:nt %time loop
error=9e9; %error re-initialisation for every time step
current_time = dt*n;
while(error>tol2) %tolerance checking for every space loop
for i=2:nx-1 %spatial loop x direction
for j=2:ny-1 %spatial loop y direction
term1 = k1*(told_unimp(i-1,j) + told_unimp(i+1,j));
term2 = k2*(told_unimp(i,j-1) + told_unimp(i,j+1));
term3 = tprev(i,j); %temperature at that point in previous timestep
t_un_imp(i,j)= k3*(term3 + term1 + term2); % calculating new temperature
error= max(max(abs(told_unimp-t_un_imp))) %calculating error
end
end
count=count+1;
told_unimp=t_un_imp; % updating old temperatures
end
figure(5)
contourf(t_un_imp)
axis ij
grid on;
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Implicit Method-JACOBI, iteration number = %d',count);
title_text2=sprintf('Time Step = %d, current time = %.4f',n,current_time);
title_text3=sprintf('Temp. diff. b/w time steps = %.10f',error);
title({[title_text1],[title_text2],[title_text3]});
pause(0.003)
tprev=t_un_imp;
if(error<1e-8)
break
end
end
unsteady_imp=t_un_imp;
counter5=count;
error5=error;
ts2=n;
endUnsteady-State Gauss-Seidel Method:
This function solves the Unsteady-state problem using the Implicit method using the Gauss-Seidel iterative solver.
"t_un_imp" is initialized by copying initial temps, then a copy of "t_un_imp" is defined as "told_unimp" to store the values at the previous iteration.
value of error is initialized as `9^9`.
An outer FOR loop is used to change the time-step of the solver, once it reaches maximum timestep defined in the main program or if it converges before that, it will stop the solver.
An outer while loop is used to check the tolerance criteria, if the value of error goes less than the tolerance specified in the main program, it will stop the space solver
A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.
Error is calculated inside the Y-direction space loop as a maximum of difference in temperatures of all points between consecutive iterations.
"told_unimp" is updated with "t_un_imp" values after the space loop so that the old values changes after each iteration.
A filled contour plot of "t_un_imp" is plotted within the while loop to visualize the convergence.
an IF loop is used to check for the convergence.
function[unsteady_imp,counter5,error5,ts3]= gs_unsteady(t_initial,nx,ny,k1,k2,tol2,nt,dt)
%unsteady state
%trnsient state solution
%implicit method, G-S iterative solver
% dT/dt = a*((d^2T/dx^2)+(d^2T/dy^2))
t_un_imp=t_initial;
%creating a copy of t
told_unimp=t_un_imp;
tprev=t_un_imp;
k3= 1/(1+2*k1+2*k2);
count=0;
%convergence loop
for n=1:nt %time loop
error=9e9; %error re-initialisation for every time step
current_time = dt*n;
while(error>tol2) %tolerance checking for every space loop
for i=2:nx-1 %spatial loop x direction
for j=2:ny-1 %spatial loop y direction
term1 = k1*(t_un_imp(i-1,j) + told_unimp(i+1,j));
term2 = k2*(t_un_imp(i,j-1) + told_unimp(i,j+1));
term3 = tprev(i,j); %temperature at that point in previous timestep
t_un_imp(i,j)= k3*(term3 + term1 + term2); % calculating new temperature
error= max(max(abs(told_unimp-t_un_imp))) %calculating error
end
end
count=count+1;
told_unimp=t_un_imp; % updating old temperatures
end
figure(6)
contourf(t_un_imp)
axis ij
grid on;
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Implicit Method-GS, iteration number = %d',count);
title_text2=sprintf('Time Step = %d, current time = %.4f sec',n,current_time);
title_text3=sprintf('Temp. diff. b/w time steps = %.10f',error);
title({[title_text1],[title_text2],[title_text3]});
pause(0.003)
tprev=t_un_imp;
if(error<1e-8)
break
end
end
unsteady_imp=t_un_imp;
counter5=count;
error5=error;
ts3=n;
endUnsteady-State Successive Over Relaxation on Gauss-Seidel Method:
This function solves the Unsteady-state problem using the Implicit method using SOR on the Gauss-Seidel iterative solver.
"t_un_imp" is initialized by copying initial temps, then a copy of "t_un_imp" is defined as "told_unimp" to store the values at the previous iteration.
value of error is initialized as `9^9`.
An outer FOR loop is used to change the time-step of the solver, once it reaches maximum timestep defined in the main program or if it converges before that, it will stop the solver.
An outer while loop is used to check the tolerance criteria, if the value of error goes less than the tolerance specified in the main program, it will stop the space solver
A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.
Error is calculated inside the Y-direction space loop as a maximum of difference in temperatures of all points between consecutive iterations.
"told_unimp" is updated with "t_un_imp" values after the space loop so that the old values changes after each iteration.
A filled contour plot of "t_un_imp" is plotted within the while loop to visualize the convergence.
an IF loop is used to check for the convergence.
function[unsteady_imp,counter5,error5,ts4]= sorgs_unsteady(t_initial,nx,ny,k1,k2,tol2,nt,dt,sor_factor)
%unsteady state
%transient state solution
%implicit method, SOR on G-S iterative solver
% dT/dt = a*((d^2T/dx^2)+(d^2T/dy^2))
t_un_imp=t_initial;
%creating a copy of t
told_unimp=t_un_imp;
tprev=t_un_imp;
%convergence loop
k3= 1/(1+2*k1+2*k2);
count=0;
error=9e9; %error initialisation
for n=1:nt %time loop
current_time = dt*n;
error=9e9; %error re-initialisation for every time step
while(error>tol2) %tolerance checking for every space loop
for i=2:nx-1 %spatial loop x direction
for j=2:ny-1 %spatial loop y direction
term1 = k1*(t_un_imp(i-1,j) + told_unimp(i+1,j));
term2 = k2*(t_un_imp(i,j-1) + told_unimp(i,j+1));
term3 = tprev(i,j); %temperature at that point in previous timestep
t_un_imp(i,j)= (1-sor_factor)*told_unimp(i,j) + sor_factor*k3*(term3 + term1 + term2); % calculating new temperature
error= max(max(abs(told_unimp-t_un_imp))) %calculating error
end
end
count=count+1;
told_unimp=t_un_imp; % updating old temperatures
end
figure(7)
contourf(t_un_imp)
axis ij
colorbar; colormap(jet);
grid on;
title_text1=sprintf('Unsteady State-Implicit Method-SOR-GS, iteration number = %d',count);
title_text2=sprintf('SOR factor= %.2f, Time Step = %d, current time = %.4f',sor_factor,n,current_time);
title_text3=sprintf('Temp. diff. b/w time steps = %.10f',error);
title({[title_text1],[title_text2],[title_text3]});
pause(0.003)
tprev=t_un_imp;
if(error<1e-8)
break
end
end
unsteady_imp=t_un_imp;
counter5=count;
error5=error;
ts4=n;
endResults:
Screen capture of contours of different solvers while they solving are shown below:
- Steady-State:
- Unsteady-State:
- Figure(8): will show a 2*2 subplot of filled contours of the initial condition and final conditions of Steady-state solvers i.e. Jacobi, Gauss-Seidel, SOR on GS.
- Figure(9): will show a 2*2 subplot of filled contours of the final conditions of Unsteady-state solvers i.e. Explicit, Jacobi, Gauss-Seidel, SOR on GS.
- Figure(10): will show a bar plot of no. of iterations taken by Steady-state solvers to achieve error less than tolerance specified.
- Figure(11): will show a bar plot of no. of iterations taken by Unsteady-state solvers to achieve steady-state.
- Conclusion:
- The Fastest Steady-State solution is achieved by using SOR on Gauss-Siedal iterative solver and it was approx 66.5% faster than Jacobi iterative method and 34.75% faster than simple G-S iterative solver.
- The fastest Unsteady-State solution is achieved by using Explicit time marching technique and it was 85.89% faster than Implicit-Jacobi method, 80.91% faster than Implicit-GS method and 75.93% faster than Implicit-SOR on GS method in terms of number of spatial iterations required to achieve steady-state
- When overall Steady-State and Unsteady-State solvers are compared, it appears that the steady solver is very fast in achieving convergence even when compared to the fastest Unsteady-State solver.
But from the above observation, the Implicit scheme is unconditionally stable, so we can directly reach on the converged Steady-State in just one iteration if we manipulate the value of "dt" accordingly. So, out of all types of solver, Implicit-GS will be the fastest solver among all tested here to achieve convergence.