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
Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
The governing equation for 2-D heat conduction is : 2-D Transient :- ∂T∂t+α(∂2T∂x2+∂2T∂y2)=0∂T∂t+α(∂2T∂x2+∂2T∂y2)=0 where α is the thermal diffusivity of the medium. Descritizing the above equation we get:- `{T_(i,j)^(n+1) - T_(i,j)^(n)}/{Deltat} + alpha*({{T_(i-1,j) -2*T_(i,j)…
Vipul Anand
updated on 26 Mar 2021
The governing equation for 2-D heat conduction is :
2-D Transient :-
∂T∂t+α(∂2T∂x2+∂2T∂y2)=0
where
α is the thermal diffusivity of the medium.
Descritizing the above equation we get:-
Tn+1i,j-Tni,jΔt+α⋅({Ti-1,j-2⋅Ti,j+Ti+1,j}Δx2+{Ti,j-1-2⋅Ti,j+Ti,j+1}Δy2)n
Rearranging the terms:-
Tn+1i,j=Tni,j+α⋅Δt({Ti-1,j-2⋅Ti,j+Ti+1,jΔx2+Ti,j-1-2⋅Ti,j+Ti,j+1Δy2}n)
For Explicit Scheme we have :-
Tn+1i,j=Tni,j+α⋅Δt({Ti-1,j-2⋅Ti,j+Ti+1,jΔx2+Ti,j-1-2⋅Ti,j+Ti,j+1Δy2}n)
For Implicit Scheme we have :-
Tn+1i,j=Tni,j+α⋅Δt({Ti-1,j-2⋅Ti,j+Ti+1,jΔx2+Ti,j-1-2⋅Ti,j+Ti,j+1Δy2}n+1)
let αΔtΔx2 be represented as k1
and αΔtΔy2 be represented as k2
Substituting it in the above equation we get:-
Tn+1i,j=Tni,j+{k1⋅(Ti-1,j-2⋅Ti,j+Ti+1,j)+k2⋅(Ti,j-1-2⋅Ti,j+Ti,j+1)}n+1
Rearranging the terms for Implicit Scheme:-
(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
Rearranging the terms for Explicit Scheme:-
Tn+1i,j=Tni,j+k1{Ti-1,j-2⋅Ti,j+Ti+1,j}+k2{Ti,j-1-2⋅Ti,j+Ti,j+1})
While, in explicit scheme we come accross unstablity condition, so the value of k1 and k2 has been calculated based on von neumann stability method.
k1+k2 <= 0.5
2-D Steady State :-
∂2T∂x2+∂2T∂y2=0
Discritizing the above eliptic equation using central difference scheme:
∂2T∂x2=Ti-1,j-2⋅Ti,j+Ti+1,jΔx2
and
∂2T∂y2=Ti,j-1-2⋅Ti,j+Ti,j+1Δy2
substituting it in the above governing equation for steady state we get,
Ti-1,j-2⋅Ti,j+Ti+1,jΔx2+Ti,j-1-2⋅Ti,j+Ti,j+1Δy2=0
rearranging the above equation we get,
Ti-1,j+Ti+1,jΔx2+Ti,j-1+Ti,j+1Δy2=2⋅Ti,j⋅(1Δx2+1Δy2)
Ti,j=1k⋅(Ti+1,j+Ti-1,jΔx2+Ti,j-1+Ti,j+1Δy2)
where k is
k=2Δx2+Δy2Δx2⋅Δy2
Solution Schemes:-
Explicit Scheme:-
In this scheme temperature at node of at next time step, being the only unkown in the equation is calculated by using information of preceding time step only.
Tn+1i,j=Tni,j+k1(Ti+1,j-2⋅Ti,j+Ti+1,j)n+k2(Ti,j+1-2Ti,j+Ti,j-1)n
Implicit Scheme:-
The temperature values of next time step is calculated by solving simultaneous difference equation for the entire domain.
For optimization of time and space of computational problem, iterative technique is used.
Tn+1i,j=Tni,j+k1(Ti+1,j-2Ti,j+Ti-1,j)n+1+k2(Ti,j+1,-2Ti,j+Ti,j-1)n+1
Iterativite schemes:-
1. Jacobi Method
This method solves the difference equation for the entire domain interatively. I reduces the cost of matrix calculations, which gets heavier when the domain of interest is large.
For Transient Calculation:-
Tit+1i,j=11+2k1+2k2⋅(Titi,j+k1(Titi+1,j+Titi-1,j)+k2(Titi,j+1+Titi,j-1))
For Steady State Calculation
Tit+1i,j=1k(Titi+1,j+Titi-1,j△x2+Titi,j+1+Titi,j-1△y2)
it represents the iteration number.
2. Gauss-Seidel Method
A modified form of Jacobian, uses the updated value from the present iteration for known nodes. This results in fast convergence.
For transient State Calculation
Tit+1i,j=11+2k1+k2(Titi,j+k1(Titi+1,j+Tit+1i-1,j)+k2(Titi,j+1+Tit+1i,j-1))
For Seady State Calculation:-
Tit+1i,j=1k(Titi+1,j+Tit+1i-1,jΔx2+Titi,j+1+Tit+1i,j-1Δy2)
it represents the iteration no.
3. Successive over Relaxation Model
The method is modified form of Gauss-Seidal. The solution is updated with knowledge of slope and a multiplying factor ω.
slope = Tit+1-Titit+1-it
Tit+1=Tit+ω(slope)
The temperature values are updated as :-
Tit+1=(1-ω)⋅Tit+ω⋅Tit+1
it is the iteration no.
The above numerical formation has been used in the matlab code described below:-
A class named conductionmodel has been developed in order to facilitate numerical claculation using different schemes.
classdef conductionModel
%Contains solver for Steady as well as Unsteady State
% The class initializes geometry as well as temperature values set to
% default. Also the geometry and grid size can be controlled by input
% through parameters.
% The parameters of the class are as follows:-
% Nx = No. of node point in x direction
% Ny = No. of Node points in y direction
% L = Length of the geometry.
% W = Width of the geometry.
% dt = time step to be used in the computation
% t = time of the solution.
properties
%% Defining The Geometry and Initial Condiotions
% No. of Nodes x
Nx = 11;
% No. of Nodes y
Ny = 11;
% Length
L = 1;
% Width
W = 1;
% X - Axis Nodes
x
% Y - Axis
y
% Mesh
mesh_xx
mesh_yy
% Initial Error
error = 9e9;
% Allowable Error
tolerance = 1e-4;
% Time Step
dt = 0.00125;
%
alpha = 2;
% Initial Conditions
T
dx
dy
t = 2;
% As omega should be grearter than 1
omega = 1.2;
end
methods
%% Class Constructor
% Constructor initializes the class taking input from the script
% file, entered via parameters.
function obj = conductionModel(Nx,Ny,L,W,dt, t)
%CONDUCTIONMODEL Construct an instance of this class
% Detailed explanation goes here
obj.Nx = Nx;
obj.Ny = Ny;
obj.L = L;
obj.W = W;
obj.dt = dt;
obj.x = linspace(0,obj.L,obj.Nx);
obj.y = linspace(0,obj.W,obj.Ny);
[obj.mesh_xx, obj.mesh_yy] = meshgrid(obj.x, obj.y);
% Applying initial temperature to the entire space domain
obj.T = 300*ones(obj.Nx, obj.Ny);
obj.dx = obj.x(2)-obj.x(1);
obj.dy = obj.y(2)-obj.y(1);
obj.t = t;
end
%% Temperature Boundary Conditions
function outputArg = boundaryCondition(obj,Tl,Tr,Tt,Tb)
%Set the temperature boundary condition
% The function has 4+1 input parameters, The first parameter
% is decated to the super class, rest are described below:-
% Tl = Left boundary temperature value
% Tr = Right boundary temperature value
% Tf = Top boundary temperature value
% Tb = Bottom boundary temperature value
% Applying Left, Right, Top and Bottom
% Boundary temperature values
obj.T(:,1) = Tl;
obj.T(:,end) = Tr;
obj.T(1,:) = Tt;
obj.T(end,:) = Tb;
% Applying average temperature boundary condition at
% cornors as given in the question.
obj.T(1,1) = (Tl+Tt)/2;
obj.T(1,end) = (Tt+Tr)/2;
obj.T(end,1) = (Tb+Tl)/2;
obj.T(end,end) = (Tr+Tb)/2;
outputArg = obj;
end
%% Solvers Definition Unsteady State
% The section of Unsteady State contains the itterative solvers
% function for transient state calculations. Each function has an
% input variable of obj, it represts the class instance.
% Jacobian unsteady
function [T, it_counter] = jacobiUnsteady(obj)
% T, a copy of Temperature values for calculation
T = obj.T;
% T_old, a copy of T to store previous iteration data.
T_old = T;
% T_P, a copy of T to store previous time step data.
T_P = T;
% k1, k2, term1, term2, term3 constants calculated
k1 = obj.alpha*(obj.dt/(obj.dx^2));
k2 = obj.alpha*(obj.dt/(obj.dy^2));
term1 = 1/(1+2*(k1+k2));
term2 = term1*k1;
term3 = term1*k2;
% it_counter, total itteration counter
it_counter = 0;
% Time iteration
for it = 1:(obj.t/obj.dt)
% Error initialization for each time step
error = obj.error;
% Convergence Loop
while error > obj.tolerance
for i = 2:obj.Nx-1
for j = 2:obj.Ny-1
H = T_old(i-1,j) + T_old(i+1,j);
V = T_old(i,j+1) + T_old(i,j-1);
T(i,j) = term1*T_P(i,j) + term2*H + term3*V;
end
end
it_counter = it_counter + 1 ;
% Calculation of maximum error
error = max(max(abs(T_old - T)));
% updating iteration temperature data
T_old = T;
end
% Updating converged temperature data for the present time
% step
T_P = T;
end
end
% Gauss-Seidal Method
function [T, it_counter] = gaussSeidelUnsteady(obj)
% T, a copy of Temperature values for calculation
T = obj.T;
% T_old, a copy of T to store previous iteration data.
T_old = T;
% T_P, a copy of T to store previous time step data.
T_P = T;
% k1, k2, term1, term2, term3 constants calculated
k1 = obj.alpha*(obj.dt/(obj.dx^2));
k2 = obj.alpha*(obj.dt/(obj.dy^2));
term1 = 1/(1+2*(k1+k2));
term2 = term1*k1;
term3 = term1*k2;
% it_counter, total itteration counter
it_counter = 0;
% Time iteration
for it = 1:(obj.t/obj.dt)
% Error initialization for each time step
error = obj.error;
% Convergence Loop
while error > obj.tolerance
for i = 2:obj.Nx-1
for j = 2:obj.Ny-1
H = T(i-1,j) + T_old(i+1,j);
V = T_old(i,j+1) + T(i,j-1);
T(i,j) = term1*T_P(i,j) + term2*H + term3*V;
end
end
it_counter = it_counter + 1;
% Calculation of maximum error
error = max(max(abs(T_old - T)));
% updating iteration temperature data
T_old = T;
end
% Updating converged temperature data for the present time
% step
T_P = T;
end
end
% Successive Over Relaxation
function [T, it_counter] = sorUnsteady(obj)
% T, a copy of Temperature values for calculation.
T = obj.T;
% T_old, a copy of T to store previous iteration data.
T_old = T;
% T_P, a copy of T to store previous time step data.
T_P = T;
% k1, k2, term1, term2, term3 constants calculated.
k1 = obj.alpha*(obj.dt/(obj.dx^2));
k2 = obj.alpha*(obj.dt/(obj.dy^2));
term1 = 1/(1+2*(k1+k2));
term2 = term1*k1;
term3 = term1*k2;
obj.omega = 1.2;
% it_counter, total iteration counter
it_counter = 0;
% Time iteration
for it = 1:(obj.t/obj.dt)
error = obj.error;
% Convergence Loop
while error > obj.tolerance
for i = 2:obj.Nx-1
for j = 2:obj.Ny-1
H = T(i-1,j) + T_old(i+1,j);
V = T_old(i,j+1) + T(i,j-1);
temp = ((T_P(i,j)*term1)+(H*term2)+(V*term3));
T(i,j) = ((1-obj.omega)*T_old(i,j))+...
(obj.omega*temp);
end
end
it_counter = it_counter + 1;
% Calculation of maximum error
error = max(max(abs(T_old - T)));
% updating iteration temperature data
T_old = T;
end
% Updating converged temperature data for the present time
% step
T_P = T;
end
end
% Explicit unsteady
function [T, it_counter] = explicitUnsteady(obj)
% T, a copy of Temperature values for calculation.
T = obj.T;
% T_old, a copy of T to store previous time step data.
T_old = T;
% k1, k2 constants calculated.
% As per von Neumann stability method k1+k2 is <= 0.5 for
% explicit solver.
k1 = obj.alpha*(obj.dt/(obj.dx^2));
k2 = obj.alpha*(obj.dt/(obj.dy^2));
% it_counter, total iteration counter
it_counter = 0;
% Time iteration/marching
for it = 1:(obj.t/obj.dt)
for i = 2:obj.Nx-1
for j = 2:obj.Ny-1
H = T_old(i-1,j)-2*T_old(i,j)+T_old(i+1,j);
V = T_old(i,j-1)-2*T_old(i,j)+T_old(i,j+1);
T(i,j) = T_old(i,j) + k1*(H) + k2*(V);
end
end
it_counter = it_counter + 1;
% Here, T_old is storing previous time step data.
T_old = T;
end
end
%% Solvers Definition Steady State
% Jacobian Steady Method
function [T, it_counter] = jacobiSteady(obj)
% T, a copy of Temperature values for calculation
T = obj.T;
% T_old, a copy of T to store previous iteration data.
T_old = T;
% Constant calculated
k = 2*(obj.dx^2+obj.dy^2)/(obj.dx^2*obj.dy^2);
% Error is initialized
error = obj.error;
% it_counter, total iteration counter.
it_counter = 0;
% Convergence Loop
while error > obj.tolerance
for i = 2:obj.Nx-1
for j = 2:obj.Ny-1
H = (T_old(i-1,j)+T_old(i+1,j))/(obj.dx^2);
V = (T_old(i,j-1)+T_old(i,j+1))/(obj.dy^2);
T(i,j) = (H+V)/k;
end
end
it_counter = it_counter + 1;
% Calculation of maximum error
error = max(max(abs(T_old - T)));
% Updating previous iteration temperature data.
T_old = T;
end
end
% Gauss-Seidel Method
function [T, it_counter] = gaussSeidelSteady(obj)
% T, a copy of Temperature values for calculation
T = obj.T;
% T_old, a copy of T to store previous iteration data.
T_old = T;
% Constant calculated
k = 2*(obj.dx^2+obj.dy^2)/(obj.dx^2*obj.dy^2);
% Error is initialized
error = obj.error;
% it_counter, total iteration counter
it_counter = 0;
% Convergence Loop
while error > obj.tolerance
for i = 2:obj.Nx-1
for j = 2:obj.Ny-1
H = (T(i-1,j)+T_old(i+1,j))/(obj.dx^2);
V = (T(i,j-1)+T_old(i,j+1))/(obj.dy^2);
T(i,j) = (H+V)/k;
end
end
it_counter = it_counter + 1;
% Calculation of maximum error
error = max(max(abs(T_old - T)));
% Updating temperature values of previous iteration
T_old = T;
end
end
% Successive Over Relaxation
function [T, it_counter] = sorSteady(obj)
% T, a copy of Temperature values for calculation
T = obj.T;
% T_old, a copy of T to store previous iteration data.
T_old = T;
% Constant calculated
k = 2*(obj.dx^2+obj.dy^2)/(obj.dx^2*obj.dy^2);
% Initializing Error
error = obj.error;
% it_counter, total iteration counter
it_counter = 0;
% Convergence Loop
while error > obj.tolerance
for i = 2:obj.Nx-1
for j = 2:obj.Ny-1
H = (T(i-1,j)+T_old(i+1,j))/(obj.dx^2);
V = (T(i,j-1)+T_old(i,j+1))/(obj.dy^2);
temp = (H+V)/k;
T(i,j) = ((1-obj.omega)*T_old(i,j))+...
(obj.omega*(temp));
end
end
it_counter = it_counter + 1;
% Calculation of error
error = max(max(abs(T_old - T)));
% Updating previous iteration temperature values.
T_old = T;
end
end
end
end
Another class of name Result has been developed to storing and handling the resulting variables.
classdef Result
% Property:-
% T will hold the temperature values
% it will hold iteration no. for the solution to converge
% title will hold the string value of Method employed to solve the
% equation.
properties
T
it
title
end
methods
function obj = Result(T, it, title)
obj.T = T;
obj.it = it;
obj.title = title;
end
end
end
plotter function is defined for plotting data:-
function fig = plotter(xx,yy,T,it,ix, method)
% A function for platting contour plots
% xx,yy holds the x and y cordinated of mesh
% T = Templerature value to be plotted
% it = total iteration no to obtain solution
% method = which method was used to obtain the result
fig = figure(ix+1);
[l,m] = contourf(xx,yy,T);
clabel(l,m);
colorbar
title({method,"No of iteration = "+ it});
xlabel("X- Direction")
ylabel("Y - Direction")
end
Driver script for utilizing the developed formulation is described below:-
clear all
close all
clc
% Declaring No. of nodex is X - Direction.
Nx = 11;
% Declaring No. of Nodes in y - Direction.
Ny = 11;
% Length of the Domain
L = 1;
% Width of the Domain
W = 1;
% Time Step
dt = 0.00125;
% Time of Simulation
t = 2;
% a, conductionModel object is instantiated
a = conductionModel(Nx, Ny, L, W, dt, t);
% Mesh grid is copied from the object for plotting
xx = a.mesh_xx;
yy = a.mesh_yy;
% Boundary Condition is applied
a = boundaryCondition(a,400,800,600,900);
% A contour plot is presented for boundary condition discription
figure(1)
[l,m] = contourf(xx,yy,a.T);
clabel(l,m)
title("Initial Condition")
xlabel("X - Axis")
ylabel("Y - Axis")
colorbar
% The model is solved for Steady as well as unsteady cases as descibed
% below.
[T, it] = jacobiUnsteady(a);
[T1, it1] = sorUnsteady(a);
[T2, it2] = explicitUnsteady(a);
[T3, it3] = gaussSeidelUnsteady(a);
[T4, it4] = jacobiSteady(a);
[T5, it5] = gaussSeidelSteady(a);
[T6, it6] = sorSteady(a);
% Result array is declared for ease in plotting the output data.
% Result class is declared to hold individual results.
Result_array = [Result(T,it, "Jacobian Unsteady State");...
Result(T1, it1, "SOR Unsteady State");...
Result(T2, it2, "Explicit unsteady State");...
Result(T3, it3, "Gauss Seidel Unsteady State");...
Result(T4, it4, "Jacobian Steady State");...
Result(T5, it5, "Gauss Seidel Steady State");
Result(T6,it6, "SOR Steady State")];
for i = 1:length(Result_array)
% Plotter function has been written for routine plotting.
plotter(xx, yy, Result_array(i).T, Result_array(i).it, i, Result_array(i).title);
end
Output:-
Comparing the plots and iterations of unsteady state implicit analysis. No. of iteration has been reducing from jacobian to sor method with no perceivable change in result. It means that SOR is the best performer. Since the iteration in SOR method is less, so time consumed in order to converge to aproximate solution should also be less as compared to other two.
Same Pattern can be seen in Steady state also:-
No. of iteration in explicit analysis cant be directly compared with the implicit method as it has fixed no. of iteration, which is equal to the no. of time steps.
Also, Explicit analysis suffers a stability cause. The time step chosen in the above alaysis is such that the explicit soultion is stable. It can be easily shown that implicit solution can be obtained with less no. of iteration without compromising with the result.
As the time step has been choosen keeping in mind the stability of explicit solution, The absolute values on no. of iteration cant be compared between explicit and implicit methods.
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 Vipul Anand (37)
Senstivity Analysis - GRI Mechanism
Introduction: Sensitivity analysis given an overview of magnitude of variation of an independent variable over dependent variable. It is the study that helps in investigation of factors that has the potential to alter the given outcome. Since our study is focussed on the chemical reactions, the reactants and their initial…
04 Jul 2022 02:31 PM IST
FULL HYDRO case set up (PFI)
Objective:- Simulate PFI - engine in Converge CFD Geometry The prepared geometry from the previous section has been used in the present simulation. Simulation Type of Simulation Crank angle-based IC engine simulation The given parameters were, RPM of 3000 Connecting rod of length 0.18m Stroke of length 0.09m…
04 Jul 2022 02:25 PM IST
Simulation of Flat Plate Heat Sink
Detailed Objectives and Motivation Flat plate heats sinks are most common cooling solution for low to medium heat dissipating electronic equipment. Its ease in manufacturing attracts diverse range of material selection with required thermal properties. Selection of material is not the only part and parcel in manufacturing…
04 Jul 2022 02:24 PM IST
Simulation of Flow over a Cylinder and Explaining the phenomenon of Karman Vortex Street
Project Report - Flow over a Cylinder Introduction Flow past object, has been a topic of interest since long time. The phenomenon is very common in engineering designs. The complex nature of physics and behaviour of transport phenomenon around the object is still being actively studied. The nature of flow changes at critical…
04 Jul 2022 01:49 PM IST
Related Courses
0 Hours of Content
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.