Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Post Graduate Program in Hybrid Electric Vehicle Design and AnalysisPost Graduate Program in Computational Fluid DynamicsPost Graduate Program in CADPost Graduate Program in CAEPost Graduate Program in Manufacturing DesignPost Graduate Program in Computational Design and Pre-processingPost Graduate Program in Complete Passenger Car Design & Product DevelopmentSuccess Stories
Numerical solution of 2D heat conduction equation using iterative techniques in Octave / MATLAB
Objective:-> To solve the 2D steady and unsteady (Transient) heat conduction equations by using the point iterative techniques. Following iterative technique can be used. Jacobi Gauss-seidel Successive over-relaxation Given data:- Top Boundary = 600 K Bottom Boundary = 900 K Left Boundary = 400 K Right Boundary…
Smitesh Badnapurkar
updated on 24 Dec 2020
Objective:-> To solve the 2D steady and unsteady (Transient) heat conduction equations by using the point iterative techniques.
Following iterative technique can be used.
- Jacobi
- Gauss-seidel
- Successive over-relaxation
Given data:-
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
Let us consider a square plate with boundary length of 1m. The number of grid points on X-axis & Y-axis 10.
α= Thermal diffusivity = 1.4 (m^2/s)
Convergence criteria Tolerance = 0.0001
We know heat always transfer from high temperature point to low temperature point.
1. Steady state condition:- Under this condition temperature within system remains unchanged with respect to change in time.
As we know equation for 2D steady state, ∂2T∂x2+∂2T∂y2=0
Discretization,
Let us discretize terms by central differencing method.
Hence we get,
∂2T∂x2=Ti-1,j-2Ti,j+Ti+1,jΔx2
∂2T∂y2=Ti,j-1-2Ti,j+Ti,j+1Δy2
Where, i & j are subscripts and denote nodes (grid) in rows and columns
Hence,
Ti-1,j-2Ti,j+Ti+1,jΔx2+Ti,j-1-2Ti,j+Ti,j+1Δy2=0
∴2Ti,jΔx2+2Ti,jΔy2=Ti-1,j+Ti+1,jΔx2+Ti,j-1+Ti,j+1Δy2
∴Ti,j(2Δx2+2Δy2)=Ti-1,j+Ti+1,jΔx2+Ti,j-1+Ti,j+1Δy2
∴Ti,j×k=Ti-1,j+Ti+1,jΔx2+Ti,j-1+Ti,j+1Δy2
Where, k=(2Δx2+2Δy2)
MATLAB Code for Steady state condition:-
Let us create code by considering above condition and equation.
clear all
close all
clc
L=1; %Length of square domain in meter
nx=ny=10; %Number of nodes (gridpoints)
%Mesh creation
x=linspace(0,L,nx);
dx=L/(nx-1); %Meshing in X direction
y=linspace(0,L,ny);
dy=L/(ny-1); %Meshing in Y direction
%2D matrix creation
[X,Y] = meshgrid(x,y);
%Initialization
T=300*ones(nx,ny);
T(:,1)=400; %Temperature at left wall
T(:,nx)=800; %Temperature at right wall
T(1,:)=600; %Temperature at top wall
T(ny,:)=900; %Temperature at bottom wall
T(1,1)=(400+600)/2; %Temperature at top left corner
T(nx,1)=(400+900)/2; %Temperature at bottom left corner
T(1,ny)=(600+800)/2; %Temperature at top right corner
T(nx,ny)=(800+900)/2; %Temperature at bottom right corner
Told=T;
k=(2/(dx)^2 + 2/(dy)^2);
tol=1e-4; %Convergence criteria
%Input method number
method= input('Enter method number 1.Jacobi 2.GS 3.SOR: ')
tic;
%Jacobian method
if method==1
itr_jaco = 1;
error=9e9; %error reinitialization
%Convergence loop
while error > tol
%space loop
for i=2:nx-1
for j=2:ny-1
T(i,j)=((Told(i-1,j)+Told(i+1,j))/(dx)^2 + (Told(i,j-1)+Told(i,j+1))/(dy)^2)/k;
endfor
endfor
error_n=abs(max(T-Told));
error = max(error_n);
Told=T; %Reinitialization of temprature
itr_jaco = itr_jaco + 1;
%Plotting
figure(1)
colormap(jet);
[c,h] = contourf(x,y,T,'showtext','on');
set(gca,'YDIR','reverse');
colorbar;
xlabel('X coordinate');
ylabel('Y coordinate');
title(['Jacobian iteration value = ',num2str(itr_jaco)]);
pause(0.01);
%Ellapsed time = 542.28
endwhile
%GS
elseif method==2
itr_gs = 1;
error=9e9; %error reinitialization
%Convergence loop
while error > tol
%space loop
for i=2:nx-1
for j=2:ny-1
T(i,j)=((T(i-1,j)+Told(i+1,j))/(dx)^2 + (T(i,j-1)+Told(i,j+1))/(dy)^2)/k;
endfor
endfor
error_n=abs(max(T-Told));
error = max(error_n);
Told=T; %Reinitialization temprature
itr_gs = itr_gs + 1;
%plotting
figure(2)
colormap(jet);
[c,h] = contourf(x,y,T,'showtext','on');
set(gca,'YDIR','reverse');
colorbar;
xlabel('X coordinate');
ylabel('Y coordinate');
title(['GS iteration value = ',num2str(itr_gs)]);
pause(0.01);
%Ellapsed time = 382.46
endwhile
%SOR
elseif method==3
itr_sor = 1;
error=9e9; %error reinitialization
%for w=0.8, No. of iterations = 164
%for w=0.9, No. of iterations = 136
%for w=1.0, No. of iterations = 112
%for w=1.1, No. of iterations = 93
%for w=1.2, No. of iterations = 76
%for w=1.3, No. of iterations = 60
%for w=1.4, No. of iterations = 46
%for w=1.5, No. of iterations = 28
%for w=1.6, No. of iterations = 35
%for w=1.55, No. of iterations = 35
w=1.51; %over relaxation
%convergence loop
while error > tol
%space loop
for i=2:nx-1
for j=2:ny-1
T(i,j)=(1-w)*Told(i,j) + (w*((T(i-1,j)+Told(i+1,j))/(dx)^2 + (T(i,j-1)+Told(i,j+1))/(dy)^2)/k);
endfor
endfor
error_n=abs(max(T-Told));
error = max(error_n);
Told=T; %Reinitialization temperature
itr_sor = itr_sor + 1;
%plotting
figure(3)
colormap(jet);
[c,h] = contourf(x,y,T,'showtext','on');
set(gca,'YDIR','reverse');
colorbar;
xlabel('X coordinate');
ylabel('Y coordinate');
title(['SOR iteration value = ',num2str(itr_sor)]);
pause(0.01);
%Ellapsed time = 102.00
endwhile
end
Ellapsed_time=tocAbove code consist of all 3 iterative solvers. When user enter 1 then code of Jacobian iterative solver starts running and display output on graph. Similarly, when user enter 2 then code for Gauss-Seidel iterative solver starts running and when user enter 3 then SOR code will run.
I. Jacobian solver
Above plot shows temperature profile by Jacobian solver. For solution convergence Jacobian solver required 208 iterations.
II. Gauss- Seidel solver
Above plot shows temperature profile by GS solver. For solution convergence this solver required 112 iterations.
III. Successive over relaxation
Above plot shows temperature profile by SOR solver. For solution convergence this solver required 24 iterations only. We get least number of iterations for ω=1.51.
2. Transient condition:- For transient heat transfer conditions the temperature within system varies with respect to time.
- Explicit solver
We have equation for 2D heat transient state is, ∂T∂t-α×[∂2T∂x2+∂2T∂y2]=0
∂T∂t=α×[∂2T∂x2+∂2T∂y2]
Discretization,
∂2T∂x2=Ti-1,j-2Ti,j+Ti+1,jΔx2
∂2T∂y2=Ti,j-1-2Ti,j+Ti,j+1Δy2
∂T∂t=Tn+1i,j-Tni,jΔt
Where, 'i' & 'j' are subscripts and denotes nodes in rows and columns
'n' is superscripts and denotes timestep
Hence,
Tn+1i,j-Tni,jΔt=α×[Ti-1,j-2Ti,j+Ti+1,jΔx2+Ti,j-1-2Ti,j+Ti,j+1Δy2]n
∴(Tn+1i,j-Tni,j)=α×(Δt)[Ti-1,j-2Ti,j+Ti+1,jΔx2+Ti,j-1-2Ti,j+Ti,j+1Δy2]n
∴Tn+1i,j=Tni,j+α×(Δt)[Ti-1,j-2Ti,j+Ti+1,jΔx2+Ti,j-1-2Ti,j+Ti,j+1Δy2]n
∴Tn+1i,j=Tni,j-α×(Δt)×2Tni,jΔx2-α×(Δt)×2Tni,jΔy2+α×(Δt)[Ti-1,j+Ti+1,jΔx2+Ti,j-1+Ti,j+1Δy2]n
∴Tn+1i,j=(1-2α×ΔtΔx2-2α×ΔtΔy2)Tni,j+α×ΔtΔx2(Ti-1,j+Ti+1,j)n+α×ΔtΔy2(Ti,j-1+Ti,j+1)n
∴Tn+1i,j=(1-2k1-2k2)Tni,j+k1(Ti-1,j+Ti+1,j)n+k2(Ti,j-1+Ti,j+1)n
Where, k1=(α×ΔtΔx2)
k2=(α×ΔtΔy2)
MATLAB code for unsteady explicit solver:-
Let us create code by considering above condition and equation.
clear all
close all
clc
L=1; %Length of square domain in meter
nx=ny=10; %Number of nodes (gridpoints)
%Mesh creation
x=linspace(0,L,nx);
dx=L/(nx-1); %Meshing in X direction
y=linspace(0,L,ny);
dy=L/(ny-1); %Meshing in Y direction
%2D matrix creation
[X,Y] = meshgrid(x,y);
%Initialization
T=300*ones(nx,ny);
T(:,1)=400; %Temperature at left wall
T(:,nx)=800; %Temperature at right wall
T(1,:)=600; %Temperature at top wall
T(ny,:)=900; %Temperature at bottom wall
T(1,1)=(400+600)/2; %Temperature at top left corner
T(nx,1)=(400+900)/2; %Temperature at bottom left corner
T(1,ny)=(600+800)/2; %Temperature at top right corner
T(nx,ny)=(800+900)/2; %Temperature at bottom right corner
alpha=1.4; %Diffusivity constant
dt=1e-3; %timestep value in seconds
t = 15; %Total time
nt = t/dt; %Number of timesteps
Told=T;
k1 = (alpha*dt)/(dx)^2 ;
k2 = (alpha*dt)/(dy)^2;
cfl = k1+k2; %Total cfl value, must be less than 0.5
tic;
%Time loop
for k=1:nt
%Space loop
for i=2:nx-1
for j=2:ny-1
T(i,j)=(1 - 2*k1 -2*k2)*Told(i,j) + k1*(Told(i-1,j)+Told(i+1,j)) + k2*(Told(i,j-1)+Told(i,j+1));
endfor
endfor
%Reinitializing temperature values
Told=T;
endfor
%Plotting
figure(1)
colormap(jet);
[c,h] = contourf(x,y,T,'ShowText','on');
set(gca,'YDIR','reverse');
colorbar;
title({'Transient state Explicit condition';['CFL ',num2str(cfl)]});
xlabel('X co-ordinate');
ylabel('Y co-ordinate');
Ellapsed_time = toc
%Ellapse time = 155.29 sec
Above plot shows temperature profile by explicit solver. For transient explicit condition combined cfl value should be less than 0.5. If cfl value exceeds 0.5 then the solution will blow up.
b. Implicit solver
We have equation for 2D heat transient state is, ∂T∂t-α×[∂2T∂x2+∂2T∂y2]=0
∂T∂t=α×[∂2T∂x2+∂2T∂y2]
Discretization,
∂2T∂x2=[Ti-1,j-2Ti,j+Ti+1,jΔx2]n+1
∂2T∂y2=[Ti,j-1-2Ti,j+Ti,j+1Δy2]n+1
∂T∂t=Tn+1i,j-Tni,jΔt
Where, 'i' & 'j' are subscripts and denotes nodes in rows and columns
'n' is superscripts and denotes timestep
Hence,
Tn+1i,j-Tni,jΔt=α×[Ti-1,j-2Ti,j+Ti+1,jΔx2+Ti,j-1-2Ti,j+Ti,j+1Δy2]n+1
∴(Tn+1i,j-Tni,j)=α×(Δt)[Ti-1,j-2Ti,j+Ti+1,jΔx2+Ti,j-1-2Ti,j+Ti,j+1Δy2]n+1
∴Tn+1i,j=Tni,j+α×(Δt)[Ti-1,j-2Ti,j+Ti+1,jΔx2+Ti,j-1-2Ti,j+Ti,j+1Δy2]n+1
∴Tn+1i,j=Tni,j-α×(Δt)×2Tn+1i,jΔx2-α×(Δt)×2Tn+1i,jΔy2+α×(Δt)[Ti-1,j+Ti+1,jΔx2+Ti,j-1+Ti,j+1Δy2]n+1
∴Tn+1i,j+α×(Δt)×2Tn+1i,jΔx2+α×(Δt)×2Tn+1i,jΔy2=Tni,j+α×(Δt)[Ti-1,j+Ti+1,jΔx2+Ti,j-1+Ti,j+1Δy2]n+1
∴Tn+1i,j(1+2α×ΔtΔx2+2α×ΔtΔy2)=Tni,j+α×ΔtΔx2(Ti-1,j+Ti+1,j)n+1+ α×ΔtΔy2(Ti,j-1+Ti,j+1)n+1
∴ (1+2k1+2k2)Tn+1i,j=Tni,j+k1(Ti-1,j+Ti+1,j)n+1+ k2(Ti,j-1+Ti,j+1)n+1
Where, k1=(α×ΔtΔx2)
k2=(α×ΔtΔy2)
∴ k3Tn+1i,j=Tni,j+k1(Ti-1,j+Ti+1,j)n+1+ k2(Ti,j-1+Ti,j+1)n+1
Where, k3=(1+2k1+2k2)
∴Tn+1i,j=Tni,j+k1(Ti-1,j+Ti+1,j)n+1+ k2(Ti,j-1+Ti,j+1)n+1k3
MATLAB code for unsteady implicit solver:-
Let us create code by considering above condition and equation.
clear all
close all
clc
L=1; %Length of square domain in meter
nx=ny=10; %Number of nodes (gridpoints)
%Mesh creation
x=linspace(0,L,nx);
dx=L/(nx-1); %Meshing in X direction
y=linspace(0,L,ny);
dy=L/(ny-1); %Meshing in Y direction
%2D matrix creation
[X,Y] = meshgrid(x,y);
%Initialization
T=300*ones(nx,ny);
T(:,1)=400; %Temperature at left wall
T(:,nx)=800; %Temperature at right wall
T(1,:)=600; %Temperature at top wall
T(ny,:)=900; %Temperature at bottom wall
T(1,1)=(400+600)/2; %Temperature at top left corner
T(nx,1)=(400+900)/2; %Temperature at bottom left corner
T(1,ny)=(600+800)/2; %Temperature at top right corner
T(nx,ny)=(800+900)/2; %Temperature at bottom right corner
alpha=1.4; %Diffusivity constant
dt=1; %timestep value in seconds
t = 15; %Total time
nt = t/dt; %Number of timesteps
Told=T;
Tprev=T;
k1 = (alpha*dt)/(dx)^2 ;
k2 = (alpha*dt)/(dy)^2;
k3 = (1+2*k1+2*k2)^(-1);
itr_jaco = 0;
itr_gs = 0;
itr_sor = 0;
tol=1e-4; %Convergence criteria
%Input method number
method= input('Enter method number 1.Jacobi 2.GS 3.SOR: ')
%Jacobian method
tic;
if method==1
%Time loop
for k=1:nt
error=9e9; %Error reinitialization
%Convergence loop
while error > tol
%Space loop
for i=2:nx-1
for j=2:ny-1
T(i,j)= k3*Tprev(i,j) + k3*k1*(Told(i-1,j)+Told(i+1,j)) + k3*k2*(Told(i,j-1)+Told(i,j+1));
endfor
endfor
error_n=abs(max(T-Told));
error = max(error_n);
Told=T; %Updating old values
itr_jaco = itr_jaco + 1;
endwhile
Tprev = T; %Updating previous timestep values
%Plotting
figure(1)
colormap(jet);
[c,h] = contourf(x,y,T,'ShowText','on');
set(gca,'YDIR','reverse');
colorbar;
xlabel('X coordinate');
ylabel('Y coordinate');
title(['Jacobian iteration value = ',num2str(itr_jaco)]);
pause(0.01);
%Ellapsed time = 46.163 sec
endfor
%GS method
elseif method==2
%Time loop
for k=1:nt
error=9e9; %Error reinitialization
%Convergence loop
while error > tol
%Space loop
for i=2:nx-1
for j=2:ny-1
T(i,j)= k3*Tprev(i,j) + k3*k1*(T(i-1,j)+Told(i+1,j)) + k3*k2*(T(i,j-1)+Told(i,j+1));
endfor
endfor
error_n=abs(max(T-Told));
error = max(error_n);
Told=T; %Updating old values
itr_gs = itr_gs + 1;
endwhile
Tprev = T; %Updating previous timestep values
%Plotting
figure(2)
colormap(jet);
[c,h] = contourf(x,y,T,'ShowText','on');
set(gca,'YDIR','reverse');
colorbar;
xlabel('X coordinate');
ylabel('Y coordinate');
title(['GS iteration value = ',num2str(itr_gs)]);
pause(0.01);
%Ellapsed time = 42.444 sec
endfor
%SOR
elseif method==3
%Time loop
for k=1:nt
error=9e9; %Error reinitialization
%For w=0.8 no of iterations = 406
%For w=0.9 no of iterations = 342
%For w=1.0 no of iterations = 288
%For w=1.1 no of iterations = 242
%For w=1.2 no of iterations = 202
%For w=1.3 no of iterations = 166
%For w=1.4 no of iterations = 231
%For w=1.5 no of iterations = 095
%For w=1.6 no of iterations = 114
%For w=1.55 no of iterations = 104
%For w=1.51 no of iterations = 096
%For w=1.49 no of iterations = 095
%For w=1.48 no of iterations = 100
w=1.495; %SOR over relaxation
%Convergence loop
while error > tol
%Space loop
for i=2:nx-1
for j=2:ny-1
Tgs= k3*Tprev(i,j) + k3*k1*(T(i-1,j)+Told(i+1,j)) + k3*k2*(T(i,j-1)+Told(i,j+1));
T(i,j) = (1-w)*Told(i,j) + w*Tgs;
endfor
endfor
error_n=max(abs(T-Told));
error = max(error_n);
Told=T; %Updating old values
itr_sor = itr_sor + 1;
endwhile
Tprev = T; %Updating previous timestep values
%Plotting
figure(3)
colormap(jet);
[c,h] = contourf(x,y,T,'ShowText','on');
set(gca,'YDIR','reverse');
colorbar;
xlabel('X coordinate');
ylabel('Y coordinate');
title(['SOR iteration value = ',num2str(itr_sor)]);
pause(0.01);
%Ellapsed time = 40.728 sec
endfor
endif
ellapsed_time=toc; Above code consist of all 3 iterative solvers. When user enter 1 then code of Jacobian iterative solver starts running and display output on graph. Similarly, when user enter 2 then code for Gauss-Seidel iterative solver starts running and when user enter 3 then SOR code will run.
I. Jacobian solver
Above plot shows temperature profile by Jacobian solver. For solution convergence Jacobian solver required 506 iterations.
II. Gauss Seidel solver
Above plot shows temperature profile by GS solver. For solution convergence this solver required 288 iterations.
III. SOR solver
Above plot shows temperature profile by SOR solver. For solution convergence this solver required 94 iterations only. We get least number of iterations for ω=1.495.
| Jacobian solver | Gauss Seidel solver | SOR solver | |
| Number of iterations required to converge solution for steady State equation | 208 | 112 | 24 |
| Number of iterations required to converge solution for transient State implicit equation | 506 | 288 | 94 |
| Time elapsed to converge solution for steady State equation (Second) | 542.28 | 382.46 | 102.00 |
| Time elapsed to converge solution for transient State implicit equation (Second) | 46.163 | 42.444 | 40.728 |
| Table: Comparison of iterations and time elapsed for converging solution for iterative solvers | |||
Conclusion:- From above plots and table it's clear that computational time and number of iteration required to converge solution for point iterative solver are as follows-
Jacobian>Gauss-Seidel>SOR
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...
Be the first to add a comment
Read more Projects by Smitesh Badnapurkar (41)
Discretization of 3D intake manifold using GEM-3D
Aim:-> To perform discretization of intake manifold with converting it in 1D circuit model from 3D model and observing effect of discretization length on result. Theory: Discretization is the splitting of large parts into smaller sections to improve a model's accuracy. There are two ways in which a fluid…
07 Sep 2021 03:50 PM IST
Conjugate heat transfer analysis of exhaust Port
Aim:-> Conjugate heat transfer (CHT) analysis of exhaust port. Theory:- As we know heat is transferred by 3 modes i.e. Conduction, convection & radiation. Conduction:- is an exchange of energy by direct interaction between molecules of a substance containing temperature differences. Convection:- is…
20 Aug 2021 02:07 PM IST
Turbocharger Modelling using GT Power
Aim:-> Study of turbocharger and its modeling using GT Suit. Theory:- A turbocharger (technically a turbosupercharger), colloquially known as turbo, is a turbine-driven, forced induction device that increases an internal combustion engine's power output by forcing…
20 Aug 2021 08:11 AM IST
Comparison of CI and SI engine and calibration of Single cylinder CI-Engine
Aim:-> To compare SI and CI engine and perform basic calibration of single cylinder CI engine. Comparison of SI V/S CI Engine model: SI Engine CI Engine SI engine is internal combustion engine that operates on the principle of spark ignition. CI engine is internal combustion engine that operates on the principle…
20 Aug 2021 08:09 AM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
Our Company
4th Floor, BLOCK-B, Velachery - Tambaram Main Rd, Ram Nagar South, Madipakkam, Chennai, Tamil Nadu 600042.
Top Individual Courses
Top PG Programs
Skill-Lync Plus
Trending Blogs
© 2025 Skill-Lync Inc. All Rights Reserved.