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 :-
`{delT}/{delt} + alpha ({del^2T}/{delx^2}+{del^2T}/{dely^2}) = 0`
where
`alpha` 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) + T_(i+1,j)}}/{Deltax^2} + {{T_(i,j-1) -2*T_(i,j) + T_(i,j+1)}}/{Deltay^2})^n`
Rearranging the terms:-
`T_(i,j)^(n+1) = T_(i,j)^(n) + alpha*Deltat({{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}}^n)`
For Explicit Scheme we have :-
`T_(i,j)^(n+1) = T_(i,j)^n + alpha*Deltat({{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}}^n)`
For Implicit Scheme we have :-
`T_(i,j)^(n+1) = T_(i,j)^n + alpha*Deltat({{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}}^(n+1))`
let `alpha{Deltat}/{Deltax^2}` be represented as k1
and `alpha{Deltat}/{Deltay^2}` be represented as k2
Substituting it in the above equation we get:-
`T_(i,j)^(n+1) = T_(i,j)^n + {k1*(T_(i-1,j) - 2*T_(i,j) + T_(i+1,j)) + k2*(T_(i,j-1) - 2*T_(i,j) + T_(i,j+1))}^(n+1)`
Rearranging the terms for Implicit Scheme:-
`(1 +2K1 + 2K2)*T_(i,j)^(n+1) = T_(i,j)^n + K1*(T_(i-1,j) + T_(i+1,j))^(n+1) + K2*(T_(i,j-1) + T_(i,j+1))^(n+1)`
Rearranging the terms for Explicit Scheme:-
`T_(i,j)^(n+1) = T_(i,j)^n + k1{T_(i-1,j) -2*T_(i,j) + T_(i+1,j)} + k2{T_(i,j-1) -2*T_(i,j) + T_(i,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 :-
`{del^2T}/{delx^2} + {del^2T}/{dely^2} = 0`
Discritizing the above eliptic equation using central difference scheme:
`{del^2T}/{delx^2} = {T_(i-1,j) -2*T_(i,j) + T_(i+1,j)}/{Deltax^2}`
and
`{del^2T}/{dely^2} = {T_(i,j-1) -2*T_(i,j) + T_(i,j+1)}/{Deltay^2}`
substituting it in the above governing equation for steady state we get,
`{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`
rearranging the above equation we get,
`{T_(i-1,j) + T_(i+1,j)}/{Deltax^2} + {T_(i,j-1) + T_(i,j+1)}/{Deltay^2} = 2*T_(i,j)*({1}/{Deltax^2} + {1}/{Deltay^2})`
`T_(i,j) = {1}/{k}*({T_(i+1,j) + T_(i-1,j)}/{Deltax^2} + {T_(i,j-1) + T_(i,j+1)}/{Deltay^2})`
where k is
`k = 2{Deltax^2 + Deltay^2}/{Deltax^2*Deltay^2}`
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.
`T_(i,j)^(n+1) = T_(i,j)^(n) + k1(T_(i+1,j) - 2*T_(i,j) + T_(i+1,j) )^n + k2(T_(i,j+1) -2T_(i,j) + T_(i,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.
`T_(i,j)^(n+1) = T_(i,j)^n + k1(T_(i+1,j) -2T_(i,j) +T_(i-1,j))^(n+1) + k2(T_(i,j+1), -2T_(i,j) + T_(i,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:-
`T_(i,j)^(it+1) = {1}/{1+2k1+2k2}*(T_(i,j)^(it) +k1(T_(i+1,j)^(it) + T_(i-1,j)^(it)) +k2 (T_(i,j+1)^(it) + T_(i,j-1)^(it)))`
For Steady State Calculation
`T_(i,j)^(it+1) = 1/k((T_(i+1,j)^(it)+T_(i-1,j)^(it))/(trianglex^2) + (T_(i,j+1)^(it)+T_(i,j-1)^(it))/(triangley^2))`
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
`T_(i,j)^(it+1)=1/(1+2k_1+k_2)(T_(i,j)^(it)+k_1(T_(i+1,j)^(it)+T_(i-1,j)^(it+1)) + k_2(T_(i,j+1)^(it)+T_(i,j-1)^(it+1)))`
For Seady State Calculation:-
`T_(i,j)^(it+1) = 1/k((T_(i+1,j)^(it)+T_(i-1,j)^(it+1))/(Deltax^2) + (T_(i,j+1)^(it)+T_(i,j-1)^(it+1))/(Deltay^2))`
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 `omega`.
slope = `{T^(it+1) - T^(it)}/{it + 1 - it}`
`T^(it+1) = T^(it)+omega(slope)`
The temperature values are updated as :-
`T^(it+1) = (1-omega)*T^(it) + omega*T^(it + 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")
endDriver 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.