BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

PPrakash Kumar5 Jul 2023
HTMLMATLAB

AIM– SOLVING THE STEADY AND UNSTEADY 2D HEAT CONDUCTION PROBLEM

Steady-State Analysis and Transient State Analysis

OBJECTIVE

Solve the 2D heat conduction equation by using the point iterative techniques.

Implement the following methods for solving implicit equations;

Jacobi

Gauss – seidel

Successive over-relaxation

Write the transient analysis code for both implicit and explicit schemes.

GIVEN

The Boundary conditions for the problem are as follows;

Top Boundary = 600 k

Bottom Boundary = 900k

Left Boundary = 400 k

Right Boundary = 800 k

Absolute error criteria = 1e - 4

SOLUTION

Thermal Conduction

Heat is transferred by conduction when adjacent atoms or molecules collide, or as several electrons move backwards and forwards from atom to atom in a diagnosed way so as not to form a macroscopic electric current, or as photons collide and scatter. Conduction is the most significant means of heat transfer within a solid or between solid objects in thermal contact.

The 2D temperature heat equation is a partial differential equation.

dTdt–α⋅(∂2T∂x2+∂T2∂y2)=0">dTdt–α⋅(d2Tdx2+dT2dy2)=0

Here, alpha is a positive coefficient called the thermal diffusivity of the medium.

Steady-State Conduction

Steady-state conduction is the form of conduction that happens when the temperature difference driving the conduction is constant so that the spatial distribution of temperatures in the conducting object does not change any further. Thus, all partial derivatives of temperature concerning space may either be zero or have nonzero values, but all derivatives of temperature at any point concerning time are uniformly zero. In steady-state conduction, the amount of heat entering any region of an object is equal to the amount of heat coming out, i.e.

∂2T∂x2+∂T2∂y2=0">d2Tdx2+dT2dy2=0

Transient Conduction

During any period in which temperatures change in time at any place within an object, the mode of thermal energy flow is termed transient conduction. Another term is "non-steady-state" conduction, referring to the time-dependence of temperature fields in an object. Non-steady-state situations appear after an imposed change in temperature at a boundary of an object. The formula for the 2D transient heat conduction equation is:

∂T∂t=(α)⋅(∂2T∂x2+∂2T∂y2)">dTdt=(α)⋅(d2Tdx2+d2Tdy2)

EXPLICIT AND IMPLICIT METHODS are approaches used in numerical analysis for obtaining numerical approximations to the solutions of time-dependent ordinary and partial differential equations, as is required in computer simulations of physical processes.

Iterative Method

In computational mathematics, an iterative method is a mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the n-th approximation is derived from the previous ones. A specific implementation of an iterative method, including the termination criteria, is an algorithm of the iterative method. An iterative method is called convergent if the corresponding sequence converges for given initial approximations.

Iterative Solver

  1. Jacobi iterative solver

The Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges.

  1. Gauss-Seidel iterative solver

The Gauss-seidel method, or the method of successive displacement, is an iterative method used to solve a system of linear equations similar to the Jacobi method. The temperature at n+1 iterations uses values of n+1 iterations that have already been computed. In this path, old approximated values are overwritten with new approximate values for the iterations performed. Unlike Jacobi, Gauss-Seidel requires one storage vector which can be overwritten with each computation. This helps the solution to converge faster than the Jacobi method.

  1. Successive over-relaxation iterative solver

SOR is a variant of the Gauss-Seidel method for solving a linear system of equations, resulting in faster convergence. This exploration takes the form of a weighted average between the previous iterate and the computed Gauss-Seidel iterate successively for each component.

Formula: Tt+1 = Tts . (1-w) + w. T(GS)n+1 or T(Jacobi)n+1

Here, w is the relaxation factor.

Range: 0 < w < 2

0 < w < 1 (Under-Relaxed state)

This condition is used when the temperature gradient is high and due to this condition, the time step does not overshoot and helps to converge at a faster rate.

1 < w < 2 (Over-Relaxed state)

This condition is used when the temperature gradient is low and due to this condition more significant time steps are used to obtain the steady-state, i.e. for the values to converge at a faster rate.

Nomenclature

TP = Ti, j

TL = Ti-1, j (Left boundary)

TR = Ti+1, j (Right boundary)

TT = Ti, j+1 (Top boundary)

TB = Ti, j-1 (Bottom boundary)

Equations used

2D heat conduction equation:

dTdt–α⋅(∂2T∂x2+∂T2∂y2)=0">dTdt–α⋅(d2Tdx2+dT2dy2)=0

In discrete form:

Ti+1,j^n-2Ti,j^n+Ti-1,j^n/dx2+Ti,j+1-2Ti,j^n+Ti,j-1^n/dy2=1/<Ti,j^n+1-Ti,j^n/dt

Here,

n represents the current time step and n+1 represents the next time step.

i, i+1 are the neighbouring points in the x-direction.

j, j+1 are the neighbouring points in the y-direction.

∠x">∠x is grid spacing in the x-direction.

∠y">∠y

is grid spacing in the y-direction.

The time derivative term is hyperbolic in nature which is to be discretised using a forward differencing scheme and the space derivatives are elliptical in nature which is to be discretised using a central differencing scheme.

Steady State 2D Heat Conduction Equation

∂2T∂x2+∂T2∂y2=0">d2Tdx2+dT2dy2=0

>Steady State Discretization

(TL – 2.TP + TR)/ ∠x2">∠x2

+ (TT – 2.TP + TB)/ ∠y2">∠y2

= 0;

  1. TL-2.TP+TR/dx2+TT-2.TP+TB/dy2=0
  • Therefore,

TP= K1 (TT + TR) + K2 (TT + TB)

Where,

K1 = ∠y22.(∠x2+∠y2)">∠y22.(∠x2+∠y2)

and K2 = ∠x22.(∠x2+∠y2)">∠x22.(∠x2+∠y2)

>Using iterative solvers for solving the conduction equation steady state:

Jacobi iterative solver

(TP)n+1 = K1 (TL + TR)n + K2 (TT + TB)n

Gauss-Seidel iterative solver

(TB)n+1 = K1 (TL + TR)n+1 + K2 (TT + TB)n+1

Successive over relaxation iterative solver

(TP)n+1 = (TP)n.(1-w) + w.( K1(TL +TR)n+1 + K2(TT+ TB)n+1)>

Transient State

∂T∂t=α⋅(∂2T∂x2+∂2T∂y2)">dTdt=α⋅(d2Tdx2+d2Tdy2)

>Unsteady-State Discretization

>Solving Explicitly

(TP)n+1 = (1 - 2.K1 – 2.K2).(TP)n + K1(TL + TR)n + K2(TT + TB)n

Here,

K1 = α.∠t∠x2">α.∠t∠x2

K2 = α.∠t∠y2">α.∠t∠y2

>Stability Criteria

K1 + K2 ≤">≤

0.5

>Solving Implicitly

(1 - 2.K1 – 2.K2).(TP)n+1 = (TP)n + K1 (TL + TR)n+1 + K2 (TT + TB)n+1

Here,

K1 = α.∠t∠x2">α.∠t∠x2

K2 = α.∠t∠y2">α.∠t∠y2

Jacobi iterative solver

(TP)n+1 = ((TP)n + K1 ((TL)n + (TR)n) + K2 ((TT)n + (TB)n) / (1 - 2.K1 – 2.K2)

Gauss-Seidel iterative solver

(TP)n+1 = ((TP)n + K1 ((TL)n+1 + (TR)n) + K2 ((TT)n + (TB)n+1) / (1 - 2.K1 – 2.K2)

Successive over-relaxation iterative solver

(TP)n+1 = (TP)n . (1 - w) + w. [((TP)n + K1 ((TL)n+1 + (TR)n) + K2 ((TT)n + (TB)n+1) / (1 - 2.K1 – 2.K2)]

ASSUMPTIONS

>Consider computation area as a square.

>nx = ny; i.e., the number of nodes in the x-direction is equal to the number of nodes in the y-direction.

>nx = ny = 10

>Apart from the boundary conditions, the inner nodes have an initial temperature value of 300K.

Boundary Conditions

>Top Boundary = 600 k

>Bottom Boundary = 900k

>Left Boundary = 400 k

>Right Boundary = 800 k

code to study the STEADYand UMSTEADY state for 2D heat conduction

MATLAB CODE

clear all
close all
clc
% Defining domain
nx = 10;
ny = 10;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% Grid size
dx = x(2)-x(1);
dy = y(2)-y(1);
% Defining temperature
T = 300*ones(nx,ny);
% Boundary Conditions
T(1,1:end) = 600;
T(end,1:end) = 900;
T(1:end,1) = 400;
T(1:end,end) = 800;
% Temperature at edges
T(1,1) = (T(1,2)+T(2,1))/2;
T(end,1) = (T(end-1,1)+T(end,2))/2;
T(1,end) = (T(1,end-1)+T(2,end))/2;
T(end,end) = (T(end-1,end)+T(end,end-1))/2;
% Absolute error criteria
tol = 1e-4;
% Initial temperature
Told = T;
% Some constants
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
% Relaxation Factor
w = 1.045 ;
% alpha
% Thermal diffusivity of the medium
A = 1.4;
% time interval
dt = 1/(2*((A/dx^2)+(A/dy^2)));
% constants
k1 = dt*A/dx^2;
k2 = dt*A/dy^2;
T_initial = T;
% Relaxation factor for unsteady state
W = 1.15;
% Function for respective solution
[error_j,iter1] = jacobi_steady_function(nx,ny,T,Told,dx,dy,x,y,k,tol);
[error_g,iter2] = gauss_seidel_function(tol,nx,ny,T,Told,dx,dy,x,y,k);
[error_s,iter3] = successive_over_relaxation_function(tol,nx,ny,T,Told,dx,dy,x,y,k,w);
[error,iter] = unsteady_explicit_function(nx,ny,k1,k2,Told,T,tol);
[error_ja,iter_ja] = unsteady_jacobi_function(nx,ny,k1,k2,T,Told,tol,T_initial)
[error_gs,iter_gs] = unsteady_gauss_seidel_function(nx,ny,k1,k2,T,Told,tol,T_initial)
[error_sor,iter_sor] = unsteady_successive_over_relaxation_function(nx,ny,k1,k2,T,Told,tol,T_initial,W)
FUNCTION CODE
Jacobi iteration method for steady-state
function [error_j,iter1] = jacobi_steady_function(nx,ny,T,Told,dx,dy,x,y,k,tol)
% dummy error value
error_j = 9e9;
% Definig iteration variable to calculate no of iterations
iter1 = 0;
% convergence loop
while (error_j > tol)
% iteration loop
for i = 2:nx-1
for j = 2:ny-1
% Formula to calculate the temprature at each grid point
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 for each iteration
error_j = max(max(abs(Told-T)));
% updating the temperature values for the next iteration
Told = T;
% Counting the no of iterations
iter1 = iter1 + 1;
end
% Plotting
figure(1)
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar
colormap("jet")
set(gca,'ydir','reverse')
title(sprintf('No. of Iterations = %d',iter1));
end
OUTPUT

Gauss-Seidel iterative method for steady-state

function [error_g,iter2] = gauss_seidel_function(tol,nx,ny,T,Told,dx,dy,x,y,k)
% Definig iteration variable to calculate no of iterations
iter2 = 0;
% dummy error value
error_g = 9e9;
% convergence loop
while (error_g > tol)
% iteration loop
for i = 2:(nx-1)
for j = 2:(ny-1)
% Formula to calculate the temprature at each grid point
T(i,j) = (1/k)*(((T((i-1),j)+T((i+1),j))/dx^2)+((T(i,(j-1))+T(i,(j+1)))/dy^2));
end
end
% error for each iteration
error_g = max(max(abs(Told-T)));
% updating the temperature values for the next iteration
Told = T;
% Counting the no of iterations
iter2 = iter2 + 1;
end
% Plotting
figure(2)
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar
colormap("jet")
set(gca,'ydir','reverse')
title(sprintf('No. of Iterations = %d',iter2));
end
OUTPUT
Successive over-relaxation method for steady-state
function [error_s,iter3] = successive_over_relaxation_function(tol,nx,ny,T,Told,dx,dy,x,y,k,w)
% Definig iteration variable to calculate no of iterations
iter3 = 0;
% dummy error value
error_s = 9e9;
% convergence loop
while (error_s > tol)
% iteration loop
for i = 2:(nx-1)
for j = 2:(ny-1)
% Formula to calculate the temprature at each grid point
T(i,j) = Told(i,j)*(1-w)+w*((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 for each iteration
error_s = max(max(abs(Told-T)));
% updating the temperature values for the next iteration
Told = T;
% Counting the no of iterations
iter3 = iter3 + 1;
end
% Plotting
figure(3)
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar
colormap("jet")
set(gca,'ydir','reverse')
title(sprintf('No. of Iterations = %d',iter3));
end
OUTPUT
Solving for unsteady-state explicitly
function [error,iter] = unsteady_explicit_function(nx,ny,k1,k2,Told,T,tol)
% Dummy Error
error = 9e9;
% Iteration Value
iter = 0;
% Convergence loop
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
X = Told(i-1,j) - 2*Told(i,j) + Told(i+1,j);
Y = Told(i,j-1) - 2*Told(i,j) + Told(i,j+1);
T(i,j) = Told(i,j) + k1*X + k2*Y;
end
end
% Max Error for the iteration
error = max(max(abs(Told - T)));
% Updating Temperature Value
Told =T;
% No. of iteration
iter = iter + 1;
end
% Plotting
figure(4)
[C,H] = contourf(T);
clabel(C,H);
axis ij
title(sprintf('no. of iterations = %d', iter));
end
OUTPUT
Unsteady state using Jacobi iteration
function [error_ja,iter_ja] = unsteady_jacobi_function(nx,ny,k1,k2,T,Told,tol,T_initial)
% Dummy Error
error_ja = 9e9;
error_t = 9e9;
% Iteration Value
iter_ja = 0;
% Time loop
while (error_t > tol)
% Convergence loop
while (error_ja > tol)
% Space loop
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
H = Told(i-1,j) + Told(i+1,j);
V = Told(i,j-1) + Told(i,j+1);
T(i,j) = (T_initial(i,j)*term1) + (term1*k1*H) + (term1*k2*V);
end
end
% error for convergence loop
error_ja = max(abs(Told-T));
% Temperature value
Told = T;
% Calculating no. of iteration
iter_ja = iter_ja + 1;
end
% error for time loop
error_t = max(max(abs(T_initial - T)));
% Assigning value
T_initial = T;
error_ja = 9e9;
end
% Plotting
figure(5)
contourf(T);
clabel(contourf(T))
axis ij
title(sprintf('No. of iterations = %d',iter_ja))
end
OUTPUT
Unsteady state using Gauss-Seidel iteration
function [error_gs,iter_gs] = unsteady_gauss_seidel_function(nx,ny,k1,k2,T,Told,tol,T_initial)
% Dummy Error
error_gs = 9e9;
error_t = 9e9;
% iteration value
iter_gs = 0;
% time loop
while (error_t > tol)
% convergence loop
while (error_gs > tol)
% space loop
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
H = T(i-1,j) + T(i+1,j);
V = T(i,j-1) + T(i,j+1);
T(i,j) = (T_initial(i,j)*term1) + (term1*k1*H) + (term1*k2*V);
end
end
% max error
error_gs = max(abs(Told-T));
% temperature value
Told = T;
% no. of iterations
iter_gs = iter_gs + 1;
end
% max error for time loop
error_t = max(max(abs(T_initial - T)));
T_initial = T;
error_gs = 9e9;
end
% Plotting
figure(6)
contourf(T);
clabel(contourf(T))
axis ij
title(sprintf('No. of iterations = %d',iter_gs))
end
OUTPUT
Unsteady state using Successive over relaxation iteration
function [error_sor,iter_sor] = unsteady_successive_over_relaxation_function(nx,ny,k1,k2,T,Told,tol,T_initial,W)
% Dummy Error
error_sor = 9e9;
error_t = 9e9;
% iteration
iter_sor = 0;
% time loop
while (error_t > tol)
% convergence loop
while (error_sor > tol)
% space loop
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
H = T(i-1,j) + T(i+1,j);
V = T(i,j-1) + T(i,j+1);
T(i,j) = (T_initial(i,j)*term1) + (term1*k1*H) + (term1*k2*V);
T(i,j) = (Told(i,j))*(1-W) + W*T(i,j);
end
end
% max error
error_sor = max(abs(Told-T));
% temperature
Told = T;
% no. of iterations
iter_sor = iter_sor + 1;
end
error_t = max(max(abs(T_initial - T)));
T_initial = T;
error_sor = 9e9;
end
% Plotting
figure(7)
contourf(T);
clabel(contourf(T))
axis ij
title(sprintf('No. of iterations = %d',iter_sor))
end
OUTPUT
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.