1

Two-Dimensional Steady Heat Conduction

we solved various heat conduction problems in various geometries in a systematic but highly
mathematical manner by
(1) deriving the governing differential equation by performing an energy balance on a
differential volume element,
(2) expressing the boundary conditions in the proper mathematical form, and
(3) solving the differential equation and applying the boundary conditions to determine the
integration constants.
• This resulted in a solution function for the temperature distribution in the medium, and the
solution obtained in this manner is called the analytical solution of the problem. The
analytical solution of a problem is also referred to as the exact solution since it satisfies the
differential equation and the boundary conditions.
• Analytical solution methods are limited to highly simplified problems in simple geometries
with simple boundary conditions. But many problems encountered in practice involve
complicated geometries with complex boundary conditions or variable properties, and cannot
be solved analytically. In such cases, sufficiently accurate approximate solutions can be
2
obtained by computers using a numerical method.
• For steady state with internal heat sources and constant thermal conductivity, the conduction
equation can be written for two-dimensional heat flow as
𝜕 2 𝑇 𝜕 2 𝑇 𝑒ሶ
2
+ 2+ =0
𝜕𝑥 𝜕𝑦 𝑘
• The solution to this equation may be obtained by analytical, numerical, or graphical
techniques. We shall consider one method (numerical method) for analyzing
conduction in two-dimensional systems.
• The heat flow in the x and y directions may be calculated from the Fourier equations

• The numerical methods for solving differential equations are based on replacing the
differential equations by algebraic equations. In the case of the popular finite difference
method, this is done by replacing the derivatives by differences.
• This method is based on finite difference techniques which are ideally suited for solution
by means of high-speed, digital computers 3
 • The conduction-partial-differential equations are written in a finite-difference form as

• Consider a rectangular region in which heat

conduction is significant in the x- and y- direction.
Now divide the x-y plane of the region into a
rectangular mesh of nodal points spaced ∆x and ∆ y
apart in the x and y-direction, respectively, as shown
in Figure , and consider a unit depth of ∆ z =1.
• Our goal is to determine the temperatures at the nodes,
and it is convenient to number the nodes and describe
their position by the numbers instead of actual
coordinates. The nodal network for the finite difference
formulation of two dimensional conduction in
rectangular coordinates

4
 • A logical numbering scheme for two-dimensional problems is the double subscript
notation(m, n) where m = 0, 1, 2, . . . , M is the node count in the x-direction and n = 0, 1,
2, . . . , N is the node count in the y-direction. The coordinates of the node (m, n) are
simply x = m ∆x and y = n ∆ y, and the temperature at the node (m, n) is denoted by Tm, n

• Each node represents a certain control volume or

volume element, and its temperature is a measure
of the average temperature of the control volume.
For example, the temperature of the node m,n
may be viewed as the average temperature of the
surrounding dashed area in Figure .
• The numerical accuracy of the calculations
depends strongly on the number of designated
nodal points. If this number is large (a fine
mesh), extremely accurate solutions can be
obtained .

5
 • Now consider a volume element of size ∆x × ∆ y ×1 centered about a general interior node
(m, n) in a region in which heat is generated at a rate of 𝑒ሶ and the thermal conductivity k is
constant, as shown in Figure . Assuming the direction of heat conduction to be toward the
node under consideration at all surfaces, the energy balance on the volume element can be
expressed as

(for the steady case)

• The energy balance equation for the volume element is based

on the analogous to Kirchhoff’s law for electric circuits (under
steady state, the summation of heat flows into the volume
element equals zero).
6
• Assume the temperatures between the adjacent nodes to vary linearly and noting that the
heat transfer area is Ax= ∆ y ×1 = ∆ y in the x-direction and Ay = ∆x ×1 = ∆x in the y-
direction. The energy balance relation above becomes

Dividing each term by ∆x × ∆ y and simplifying gives

• In finite difference analysis, usually a square mesh is used for simplicity, and thus ∆x and
∆ y are taken to be the same. Then ∆x = ∆ y= l, then final equation of the nodal point m,n
becomes

7
• That is, the finite difference formulation of an interior node is obtained by adding the
temperatures of the four nearest neighbors of the node, subtracting four times the
temperature of the node itself, and adding the heat generation term. It can also be expressed
in this form, which is easy to remember:

• When there is no heat generation in the medium, the finite difference equation for an
interior node further simplifies to

which has the interesting interpretation that the temperature of each interior node is the
arithmetic average of the temperatures of the four neighboring nodes. This statement is also
true for the three-dimensional problems except that the interior nodes in that case will have
six neighboring nodes instead of four.

8
 ➢ Boundary Nodes
• When boundary conditions such as the specified heat flux and convection and/or radiation are
specified at the boundary, the finite difference equation for the node at that boundary is
obtaining by written an energy balance on the volume element at that boundary. The energy
balance with convection boundary condition (Figure ) is again expressed as

Then

9
 With ∆x = ∆ y , the final form of the boundary node is
ℎ ∆𝑥 ℎ ∆𝑥 0.5 (∆𝑥)2 𝑒ሶ
0.5 𝑇𝑚−1,𝑛 +0.5 𝑇𝑚+1,𝑛 + 𝑇𝑚,𝑛−1 + 𝑇∞ -(2+ ) 𝑇𝑚,𝑛 + =0
𝑘 𝑘 𝑘

• The dimensionless group h ∆x /k in the above formulations is called the Biot number
(Bi = h ∆x /k).
• Note that equations for convection boundaries may be converted to insulated boundaries
by simply setting Bi = 0 in the respective formula.

10
11
12
13
 ➢ Solving the Finite-Difference Equations
• Once the nodal network has been established and an appropriate finite-difference equation
has been written for each node, the temperature distribution may be determined. The
problem reduces to one of solving a system of linear, algebraic equations.
• If N nodal points are selected, a set of N linear, algebraic equations is obtained. It can be
solved by matrix inversion or by a Gauss-Seidel iteration method for the values of
unknown temperatures of N nodal points.
• For small number of unknown temperatures (small number of linear, algebraic equations),
the solution can be found by hand using the relaxation method.

• Formulation as a Relaxation method. Consider a system of N finite-difference

equations corresponding to N unknown temperatures. Identifying the nodes by a single
integer subscript, rather than by the double subscript (m, n). In this technique the nodal
equation is set equal to some residual 𝜀𝑚 as follows:

14
Then the following procedure is considered:
1- Values of the nodal temperatures are assumed.
2- The value of the residual 𝜀 for each node is calculated from the respective equation and the
assumed temperatures.
3- The larger residual is “relaxed” to (𝜀 < 0.1) and the corresponding temperature of the nodal
point is calculated.
4- The residuals related to the new value of the nodal temperature are calculated.
5- The procedure is continued until all residuals are sufficiently close to zero.

15
 • The relaxation method is explained by the following
very simple example, as shown in the figure.
• Identifying the nodes by a single integer subscript (1, 2,
3, and 4), rather than by the double subscript (m, n), then
the four equations for interiors nodes 1, 2, 3, and 4 T= 50 °C
would be

T= 0 °C

Rearrange the nodal equations and set these equations to some residual 𝜀 as follows:

For approximate solution 𝜀 1, 𝜀 2, 𝜀 3, and 𝜀 4 ≈ 0 (for example, 𝜀 max ≤ 1). 16

Let’s begin with a guess values of T1, T2, T3, and T4. The guess values may be to roughly
interpolate from the boundary temperatures. For example T1 may be between 100°C and 500°C.

1. Calculate the values of 𝜀 1, 𝜀 2, 𝜀 3, and 𝜀 4 .

2. The larger 𝜀 is relaxed to zero and the corresponding temperature is calculated. For
example if 𝜀 2 is the larger value, put 𝜀 2 = 0 and calculate T2 from Equation (2).
3. Since T2 is found in Equations (1) and (4), then the residuals 𝜀 1 and 𝜀 4 are calculated with
the new value of T2.
4. Repeat steps 3 and 4 until 𝜀 1, 𝜀 2, 𝜀 3, and 𝜀 4 ≤ 𝜀 max.
The following Table shows the results of each iteration.

17
Then, T1 = 231.15, T2 = 218.55, T3 = 106.05, and T4 = 93.65.

18
• It is good practice to verify that a numerical solution has been correctly formulated this can be
achieved by performing an energy balance on a control surface surrounding all nodal regions
whose temperature have been evaluated. The temperatures should be substituted into the
energy balance equation, and if the balance is not satisfied to a high degree of precision, the
finite difference equations should be checked for errors.
• For previous example, T1 = 231.15, T2 = 218.55, T3 = 106.05 and T4 = 93.65. The energy
balance equation, required that the summation of heat input to the system equals zero.

19
20
 ➢ One-Dimensional Steady Heat Conduction
• The energy balance method for one dimensional heat conduction through plane wall is
based on subdividing the medium into a sufficient number of volume elements and then
applying an energy balance on each element. This is done by first selecting the nodal points
(or nodes) at which the temperatures are to be determined and forming elements (or control
volume) over the nodes by drawing lines through the midpoints between the nodes.
• To demonstrate the approach, again consider steady one-
dimensional heat transfer in a plane wall of thickness L
with heat generation 𝑒ሶ (x) and constant conductivity k. The
wall is now subdivided into M equal regions of thickness
∆x =L/M in the x-direction, and the divisions between the
regions are selected as the nodes. Therefore, we have M
+1 nodes labeled 0, 1, 2, . . . , m -1, m, m +1, . . . , M, as
shown in Figure

21
 • Elements are formed by drawing vertical lines through the midpoints between the nodes.
Note that all interior elements represented by interior nodes are full-size elements (they
have a thickness of Dx), whereas the two elements at the boundaries are half-sized.
• Assuming the heat conduction to be into the element on all surfaces, an energy balance on
the element can be expressed as

which simplifies to

22
➢ Boundary Conditions
The finite difference form of various boundary conditions can be obtained by replacing
ሶ
𝑄𝑙𝑒𝑓𝑡 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 by suitable expression as follows:

1. Specified Heat Flux Boundary Condition

Special case: Insulated Boundary(𝑞ሶ 0 =0)

23
2. Convection Boundary Condition

3. Radiation Boundary Condition

4. Combined Convection and Radiation Boundary

Condition

24
 ➢ Verifying the Accuracy of the Solution
• It is good practice to verify that a numerical solution has been correctly formulated by
performing an energy balance on a control surface surrounding all nodal regions
whose temperature have been evaluated. The temperatures should be substituted into the
energy balance equation, and if the balance is not satisfied to a high degree of precision, the
finite difference equations should be checked for errors.
• Another option for validating a numerical solution involves comparing results with those
obtained from exact solutions for the physical problem for which there exist exact
solutions.
• Even when the finite-difference equations have been properly formulated and solved,
the results may still represent a coarse approximation to the actual temperature field. This
behavior is a consequence of the dimensions of the volume element (∆x and ∆y) and of
finite difference approximations, such as k(∆y.1)(Tm-1,n-Tm,n)/ ∆x , to Fourier’s law of
conduction, -k(dy.1)𝜕T/𝜕x. The finite-difference approximations become more accurate as
the nodal network refined (∆x and ∆y are reduced)

25
Example2

26
27
28
Exact Solution is required
The conduction equation is

29