All Courses
All Courses
Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Success Stories
Abstract The given problem demonstrates the flow of heat along a unit square domain at a defined length. Using both explicit and implicit scheme, the steady state and transient state conduction is computed. Applying Jacobi, gauss seidel and SOR iterative solver with a defined CFL value, a comparison study is done to find…
Soudip Hazra
updated on 02 Jun 2020
Abstract
The given problem demonstrates the flow of heat along a unit square domain at a defined length. Using both explicit and implicit scheme, the steady state and transient state conduction is computed. Applying Jacobi, gauss seidel and SOR iterative solver with a defined CFL value, a comparison study is done to find out the optimal iterative solver in terms of no of iterations performed and computational time.
Governing equations:
The general formula for 2D heat conduction is given by
The first part of (1) represents transient or unsteady state and the second part represents steady state of heat conduction.
Now considering the steady state equation let us discretize using numerical schemes
From (1) we can represent the steady state as
Discretizing the above equation, we get
Denoting
Rearranging the equation, we can write as
or
Finally, the equation can be represented as
Using implicit solvers, we can write the equation with respect to the number of iterations m.
1.Jacobi solver
2.Gauss-Siedel solver
3.Succesive over relaxation solver
Where w is the relaxation factor taken as 1.1
Now considering the unsteady state equation let us discretize using two numerical schemes:
1.Explicit schemes
Discretizing the equation
we get
Where n denote the time step
Rearranging it we can write as
Here k1 and k2 represents CFL no or diffusion probability
If CFL no>0.5, the solution becomes unstable.
If CFL<0.5, the solution is stable.
2.Implicit schemes
Discretizing the equation for (n+1) time steps
Using iterative solvers on (3) with respect to iteration no m.
1.Jacobi solver
2.Gauss-Siedel solver
3.Succesive over relaxation solver
Objective:
To observe the computational time, no of iterations and thermal gradient effect by different iterative solvers both for steady and transient case.
Initial conditions and boundary conditions:
Left 400k
Right 800k
Top 600k
Bottom 900k
Workflow of the problem:
Steady state conduction
Transient state conduction
Algorithm
Algorithm for steady state conduction
Algorithm for transient state conduction explicit solvers
Algorithm for transient state conduction implicit solvers
Algorithm for computing the optimal value of omega for SOR solvers
MATLAB CODE
Steady state solver code
clear all
clc
close all
%%% define the number of grid points required in x and y direction
nx=10;%no of grids in x axis
ny=nx;%no of grids in y axis
%define the temperature domain matrix
T=ones(nx,ny)*300;
%boundary conditions of temperature domain
TL=400;
TR=800;
TT=600;
TB=900;
%defining the temperatures at the boundaries
T(2:ny-1,1)=TL;% temperature at the left boundary
T(2:ny-1,nx)=TR;% temperature at the right boundary
T(1,2:nx-1)=TT;% temperature at the top boundary
T(nx,2:ny-1)=TB;% temperature at the bottom boundary
%define the temperature at the corner points
T(1,1)=(TT+TL)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% inputs for the three function call are [T,nx,ny]
%%%%%%%%%%% function call for Jacobi iteration solver%%%%%%%%%%%
[jacobi_computation_time, TOLD, X, Y,jacobi_iteration]=steady_state_jacobi_solver(T,nx,ny);
figure(1)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
titletext1=sprintf('Thermal gradient plot of steady state heat conduction by Jacobi iterative solver\n');
titletext2=sprintf('Jacobi iteration no=%d\n',jacobi_iteration);
titletext3=sprintf('Jacobi computation time=%d',jacobi_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%% function call for Gauss-Siedel iteration solver%%%%%%%%%%%
[gs_computation_time, TOLD, X, Y,gs_iteration]=steady_state_gauss_siedel_solver(T,nx,ny);
figure(2)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
titletext1=sprintf('Thermal gradient plot of steady state heat conduction by Gauss-Siedel iterative solver\n');
titletext2=sprintf('Gauss-Siedel iteration no=%d\n',gs_iteration);
titletext3=sprintf('Gauss-Siedel computation time=%d',gs_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%% function call for Successive Over Relaxation iteration solver%%%%%%%%%%%
[sor_computation_time, TOLD, X, Y,sor_iteration]=successive_over_relaxation(T,nx,ny);
figure(3)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
titletext1=sprintf('Thermal gradient plot of steady state heat conduction by successive over relaxation iterative solver\n');
titletext2=sprintf('successive over relaxation iteration no=%d\n',sor_iteration);
titletext3=sprintf('successive over relaxation computation time=%d',sor_computation_time);
title([titletext1 titletext2 titletext3]);
function call for steady state implicit time integration jacobi solver
function [time1, TOLD, X, Y,jacobi_iteration]=steady_state_jacobi_solver(T,nx,ny)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
error=9e9;
tolerance=1e-4;
%define the constant k
k=2*(dx^2+dy^2)/(dx^2*dy^2);
[X,Y]= meshgrid(x,y);% generating mesh grid
TOLD=T;
jacobi_iteration=1;
tic
while(error>tolerance)
for i=2:nx-1%space marching step for inner loop in x direction
for j=2:ny-1%space marching step for inner loop in y direction
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=max(max(abs(TOLD-T)));
TOLD=T;
jacobi_iteration=jacobi_iteration+1;% computes the no of iteration performed
end
time1=toc;% computes the time to run the iteration
end
function call for steady state implicit time integration gauss siedel solver
function [time2, TOLD, X, Y,gs_iteration]=steady_state_gauss_siedel_solver(T,nx,ny)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
%define the constant k
k=2*(dx^2+dy^2)/(dx^2*dy^2);
[X,Y]= meshgrid(x,y);% generating mesh grid
TOLD=T;
error=9e9;
gs_iteration=1;
tic
while(error>tolerance)
for i=2:nx-1%space marching step for inner loop in x direction
for j=2:ny-1%space marching step for inner loop in y direction
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=max(max(abs(TOLD-T)));
TOLD=T;
gs_iteration=gs_iteration+1;% computes the no of iteration performed
end
time2=toc;% computes the time to run the iteration
end
function call for steady state implicit time integration SOR solver
function [time4, TOLD, X, Y,sor_iteration]=successive_over_relaxation(T,nx,ny)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
%define the constant k
k=2*(dx^2+dy^2)/(dx^2*dy^2);
[X,Y]= meshgrid(x,y);% mesh grid generation
TOLD=T;
error=9e9;
w=1.1;% relaxation factor
sor_iteration=1;
tic
while(error>tolerance)
for i=2:nx-1%space marching step for inner loop in x direction
for j=2:ny-1%space marching step for inner loop in y direction
T(i,j)=(1-w)*TOLD(i,j)+w*[(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=max(max(abs(TOLD-T)));
TOLD=T;
sor_iteration=sor_iteration+1;% computes the no of iteration performed
end
time4=toc;% computes the time to run the iteration
end
Transient state solver code
clear all
clc
close all
%%% define the number of grid points required in x and y direction
nx=10;%no of grids in x axis
ny=nx;%no of grids in y axis
alpha=0.1;%thermal diffusivity factor
dt=1e-2;
nt=1400;%no of timesteps
%define the temperature domain matrix
T=ones(nx,ny)*300;
%boundary conditions of temperature domain
TL=400;
TR=800;
TT=600;
TB=900;
%defining the temperatures at the boundaries
T(2:ny-1,1)=TL;% temperature at the left boundary
T(2:ny-1,nx)=TR;% temperature at the right boundary
T(1,2:nx-1)=TT;% temperature at the top boundary
T(nx,2:ny-1)=TB;% temperature at the bottom boundary
%define the temperature at the corner points
T(1,1)=(TT+TL)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% inputs for the function call are [T,nx,ny,alpha]
%%%%%%%%%%%%%%%%%%%%%%%%%EXPLICIT SOLVER%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% function call for explicit iteration solver%%%%%%%%%%%
[computation_time, TOLD, X, Y,explicit_iteration,CFL_explicit]=transient_conduction_explicit_solver(T,nx,ny,alpha,dt,nt);
figure(1)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Explicit iterative solver\n');
titletext2=sprintf('Explicit iteration no=%d\n',explicit_iteration);
titletext3=sprintf('Explicit solver computation time=%d',computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%IMPLICIT SOLVER%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% function call for implicit Jacobi iteration solver%%%%%%%%%%%
[jacobi_computation_time, TOLD, X, Y,jacobi_iteration,CFL_jacobi]=transient_conduction_implicit_jacobi_solver(T,nx,ny,alpha,dt,nt);
figure(2)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Jacobi solver\n');
titletext2=sprintf('Jacobi iteration no=%d\n',jacobi_iteration);
titletext3=sprintf('Jacobi solver computation time=%d',jacobi_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%% function call for implicit Gauss-Siedel iteration solver%%%%%%%%%%%
[gs_computation_time, TOLD, X, Y,gs_iteration,CFL_gs]=transient_conduction_implicit_gauss_siedel_solver(T,nx,ny,alpha,dt,nt);
figure(3)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Gauss-Siedel solver\n');
titletext2=sprintf('Gauss-Siedel iteration no=%d\n',gs_iteration);
titletext3=sprintf('Gauss-Siedelsolver computation time=%d',gs_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%SUCCESSIVE OVER RELAXATION%%%%%%%%%%%%%
%%%%%%%%%%% function call for SOR Jacobi iteration solver%%%%%%%%%%%
[sor_jacobi_computation_time, TOLD, X, Y,sor_jacobi_iteration,CFL_jacobi_sor]=transient_conduction_implicit_sor_jacobi_solver(T,nx,ny,alpha,dt,nt);
figure(4)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Jacobi SOR solver\n');
titletext2=sprintf('SOR jacobi iteration no=%d\n',sor_jacobi_iteration);
titletext3=sprintf('SOR jacobi solver computation time=%d',sor_jacobi_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%% function call for SOR Gauss-Siedel iteration solver%%%%%%%%%%%
[sor_gs_computation_time, TOLD, X, Y,sor_gs_iteration,CFL_gs_sor]=transient_conduction_implicit_sor_gauss_siedel_solver(T,nx,ny,alpha,dt,nt);
figure(5)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Gauss Siedel SOR solver\n');
titletext2=sprintf('SOR gauss siedel iteration no=%d\n',sor_gs_iteration);
titletext3=sprintf('SOR gauss siedel solver computation time=%d',sor_gs_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%plot of simulation time and no of iterations%%%%%%%%%%%%%%
iterations=[explicit_iteration jacobi_iteration gs_iteration sor_jacobi_iteration sor_gs_iteration];
simulationtime=[computation_time jacobi_computation_time gs_computation_time sor_jacobi_computation_time sor_gs_computation_time];
figure(6)
c=categorical({'explicit','jacobi','Gauss Siedel','SOR jacobi','SOR Gauss Siedel'});
bar(c,iterations,'facecolor',[1 0.5 0.5])
ylabel('No of iterations')
figure(7)
c=categorical({'explicit','jacobi','Gauss Siedel','SOR jacobi','SOR Gauss Siedel'});
bar(c,simulationtime,'facecolor',[1 0.1 0.8])
ylabel('simulation time')
function call for transient state explicit time integration solver
function [time, TOLD, X, Y,iteration,k]=transient_conduction_explicit_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);% generating meshgrid
TOLD=T;
error=9e9;
%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
iteration=0;
tic
for n=1:nt % time marching loop
for i=2:nx-1%space marching loop at inner nodes of x direction
for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=TOLD(i,j)+k1*(TOLD(i+1,j)-2*TOLD(i,j)+TOLD(i-1,j))+k2*(TOLD(i,j-1)-2*TOLD(i,j)+TOLD(i,j+1));
end
end
TOLD=T;% updates the old values of temperature
iteration=iteration+1;% computes the no of iteration performed
end
time=toc;% computes the time to run the simulation
end
function call for transient state implicit time integration jacobi solver
function [time, TOLD, X, Y,jacobi_iteration,k]=transient_conduction_implicit_jacobi_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating meshgrid
TOLD=T;
%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
jacobi_iteration=1;
tic
for n=1:nt % time marching loop
error=9e9;
while(error>tolerance)
for i=2:nx-1%space marching loop at inner nodes of x direction
for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=term1*TPREV(i,j)+term2*(TOLD(i+1,j)+TOLD(i-1,j))+term3*(TOLD(i,j-1)+TOLD(i,j+1));
end
end
error=max(max(abs(TOLD-T)));% computes the error at each iteration
TOLD=T;% updates the value of temperature
jacobi_iteration=jacobi_iteration+1;% computes the no of iteration performed
end
TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc;% computes the time to run the simulation
end
function call for transient state implicit time integration gauss siedel solver
function [time, TOLD, X, Y,gs_iteration,k]=transient_conduction_implicit_gauss_siedel_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);% generating meshgrid
TOLD=T;
%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%define the constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
gs_iteration=1;
tic
for n=1:nt % time marching loop
error=9e9;
while(error>tolerance)%convergence loop
for i=2:nx-1%space marching loop at inner nodes of x direction
for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=term1*TPREV(i,j)+term2*(TOLD(i+1,j)+T(i-1,j))+term3*(T(i,j-1)+TOLD(i,j+1));
end
end
error=max(max(abs(TOLD-T)));% computes the error at each iteration
TOLD=T;% updates the old values of temperature
gs_iteration=gs_iteration+1;% computes the no of iteration performed
end
TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end
function call for transient state implicit time integration SOR jacobi solver
function [time, TOLD, X, Y,sor_iteration,k]=transient_conduction_implicit_sor_jacobi_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;
%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
sor_iteration=1;
w=1.1;% relaxation factor
tic
for n=1:nt % time marching loop
error=9e9;
while(error>tolerance)% convergence loop
for i=2:nx-1%space marching loop at inner nodes of x direction
for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+TOLD(i-1,j))+term3*(TOLD(i,j-1)+TOLD(i,j+1)));
end
end
error=max(max(abs(TOLD-T)));% computes the error at each iteration
TOLD=T;% updates the value of temperature
sor_iteration=sor_iteration+1;% computes the no of iteration performed
end
TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end
function call for transient state implicit time integration SOR gauss siedel solver
function [time, TOLD, X, Y,sor_iteration,k]=transient_conduction_implicit_sor_gauss_siedel_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;
%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
sor_iteration=1;
w =1.1;% relaxation factor
tic
for n=1:nt % time marching loop
error=9e9;
while(error>tolerance)% convergence loop
for i=2:nx-1%space marching loop at inner nodes of x direction
for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+T(i-1,j))+term3*(T(i,j-1)+TOLD(i,j+1)));
end
end
error=max(max(abs(TOLD-T)));% computes the error at each iteration
TOLD=T;% updates the value of temperature
sor_iteration=sor_iteration+1;% computes the no of iteration performed
end
TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end
Code for computing optimal omega
%%%%%%%%%%%% Code to compute the optimal value of omega for SOR iterative solver%%%%%
clear all
clc
close all
%%% define the number of grid points required in x and y direction
nx=10;%no of grids in x axis
ny=nx;%no of grids in y axis
alpha=0.1;%thermal diffusivity factor
dt=1e-2;
nt=1400;%no of timesteps
%define the temperature domain matrix
T=ones(nx,ny)*300;
%boundary conditions of temperature domain
TL=400;
TR=800;
TT=600;
TB=900;
%defining the temperatures at the boundaries
T(2:ny-1,1)=TL;% temperature at the left boundary
T(2:ny-1,nx)=TR;% temperature at the right boundary
T(1,2:nx-1)=TT;% temperature at the top boundary
T(nx,2:ny-1)=TB;% temperature at the bottom boundary
%define the temperature at the corner points
T(1,1)=(TT+TL)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%IMPLICIT SOLVER%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%SUCCESSIVE OVER RELAXATION%%%%%%%%%%%%%
omega_range = linspace(0.5, 1.6, 100);%define the range of omega
%%%%%%%%%%% function call for SOR Gauss-Siedel and SOR Jacobi solver%%%%%%%%%%%
for i = 1:length(omega_range)%marching loop for the entire range of omega
omega = omega_range(i);
%stores the no of iteration value for each omega value for jacobi and GS solver
iter_gs(i)=transient_implicit_sor_gauss_siedel_optimal_omega(T,nx,ny,alpha,dt,nt, omega);%
iter_jacobi(i)= transient_implicit_sor_jacobi_optimal_omega(T,nx,ny,alpha,dt,nt, omega);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%Computing the optimal node and omega value for SOR Gauss Siedel and SOR Jacobi %%%%%%%%%%%%%%%%%%%%%%%%
for j=1:length(iter_gs)% marching loop for the entire length of iter
minimum_iteration_gs=min(iter_gs);%stores the minimum value of iteration along the entire range of omega
minimum_iteration_jacobi=min(iter_jacobi);%stores the minimum value of iteration along the entire range of omega
end
position_of_optimal_omega_gs=find(iter_gs==minimum_iteration_gs);% stores the nodal position of the minimum iteration value
optimal_value_of_omega_gs=omega_range(position_of_optimal_omega_gs);%stores the optimum omega value for optimum nodal position of omega
position_of_optimal_omega_jacobi=find(iter_jacobi==minimum_iteration_jacobi);% stores the nodal position of the minimum iteration value
optimal_value_of_omega_jacobi=omega_range(position_of_optimal_omega_jacobi);%stores the optimum omega value for optimum nodal position of omega
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%plotting function%%%%%%%%%%%%%%%%%%
figure(1)%for SOR Gauss Siedel
plot(omega_range, iter_gs);
xlabel ('omega values');
ylabel ('No of iterations taken for convergence');
titletext1=sprintf('Computing optimal omega for SOR Gauss Siedel iterative solver\n');
titletext2=sprintf('Optimal value of omega=%d\n',optimal_value_of_omega_gs);
titletext3=sprintf('Optimal nodal position of omega=%d\n',position_of_optimal_omega_gs);
titletext4=sprintf('Optimal no of iteration for optimal omega value=%d\n',minimum_iteration_gs);
title([titletext1 titletext2 titletext3 titletext4]);
%%%%%%%%%%%%%%%%
figure(2)% for SOR Jacobi
plot(omega_range, iter_jacobi);
xlabel ('omega values');
ylabel ('No of iterations taken for convergence');
titletext1=sprintf('Computing optimal omega for SOR Jacobi iterative solver\n');
titletext2=sprintf('Optimal value of omega=%d\n',optimal_value_of_omega_jacobi);
titletext3=sprintf('Optimal nodal position of omega=%d\n',position_of_optimal_omega_jacobi);
titletext4=sprintf('Optimal no of iteration for optimal omega value=%d\n',minimum_iteration_jacobi);
title([titletext1 titletext2 titletext3 titletext4]);
function call for transient state implicit time integration SOR gauss siedel solver for optimal omega
function sor_iteration = transient_implicit_sor_gauss_siedel_optimal_omega(T,nx,ny,alpha,dt,nt, omega)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;
%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
sor_iteration=1;
w = omega;% relaxation factor
tic
for n=1:nt % time marching loop
error=9e9;
while(error>tolerance)% convergence loop
for i=2:nx-1%space marching loop at inner nodes of x direction
for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+T(i-1,j))+term3*(T(i,j-1)+TOLD(i,j+1)));
end
end
error=max(max(abs(TOLD-T)));% computes the error at each iteration
TOLD=T;% updates the value of temperature
sor_iteration=sor_iteration+1;% computes the no of iteration performed
end
TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end
function call for transient state implicit time integration SOR jacobi solver for optimal omega
function sor_iteration = transient_implicit_sor_jacobi_optimal_omega(T,nx,ny,alpha,dt,nt, omega)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;
%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
sor_iteration=1;
w=omega;% relaxation factor
tic
for n=1:nt % time marching loop
error=9e9;
while(error>tolerance)% convergence loop
for i=2:nx-1%space marching loop at inner nodes of x direction
for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+TOLD(i-1,j))+term3*(TOLD(i,j-1)+TOLD(i,j+1)));
end
end
error=max(max(abs(TOLD-T)));% computes the error at each iteration
TOLD=T;% updates the value of temperature
sor_iteration=sor_iteration+1;% computes the no of iteration performed
end
TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end
Results and analysis
1.Steady state solution
Thermal plot of Jacobi solver
Thermal plot of Gauss Seidel solver
Thermal plot of SOR solver
Bar plot with respect to no of iterations
Bar plot with respect to simulation time
Iteration and computational table
2.Transient state solution
Thermal plot of Jacobi solver
Thermal plot of Gauss Seidel solver
Thermal plot of SOR Jacobi solver
Thermal plot of SOR Gauss Seidel solver
Bar plot with respect to no of iterations
Bar plot with respect to simulation time
Iteration and computational table
3.Optimal omega plot for SOR Jacobi solver
4.Optimal omega plot for SOR Gauss siedel solver
5.Optimal omega computational table
Conclusion
Leave a comment
Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.
Other comments...
Week 12 - Final Project - Modelling and Analysis of a Datacenter
Objective: Use macros to create computer room air conditioning units (CRACs), server cabinets, power distribution units (PDUs), and perforated floor tiles in the data center. Organize the model using groups. Include effects of gravity and turbulence in the simulation. Define object-specific meshing parameters. Create contours,…
10 Sep 2022 07:27 PM IST
Week 10 - MRF project
Objective: The main objective of this project is to provide analysis and standards for modelling MRF fan using ICEPAK. Theory: MRF is a simple, robust CFD technique used for modelling rotating turbo machinery. In MRF technique a constant speed of rotation is assigned. The mesh model does not rotate, instead the navier…
02 Jul 2022 09:19 AM IST
Week 9 - PCB Thermal Simulation
Theory: Electronic components generate heat energy and dissipate a significant amount onto PCBs. Thermal management techniques such as heat sinks are often used to accelerate the heat dissipation rate from the electronic component to the ambient. Heat sinks offer low thermal resistance from the junction to the case and…
12 Jun 2022 10:49 AM IST
Week 8 - Natural Convection-II
Aim: Electric voltage control panel thermal analysis. Electric voltage control panel: it is a component of an electrical distribution system that divides an electric power to feed into branch circuits, while providing a circuit breaker. Components of voltage panel: Enclosure: an enclosure is a cabinet for electrical or…
24 Apr 2022 08:23 AM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
© 2025 Skill-Lync Inc. All Rights Reserved.