Chapter 4
Steady State Heat Conduction in 2D and 3D: Numerical Methods

1  Finite Difference Method

The finite difference method is implemented for multi-dimensional systems in very much the same way as for one-dimensional systems. However, since more than one space coordinate is involved, the finite difference formulae must be suitable approximations to the partial derivatives appearing in the heat equation.

1.1  Cartesian Coordinates

The 2D steady state heat equation with internal heat generation in rectangular Cartesian coordinates is

x
( k T

x
) +

y
( k T

y
) + g(x,y) = 0
If k and g are assumed constant this becomes
2 T

x2
+ 2 T

y2
+ g

k
= 0
This is called the Poisson's equation. Furthermore, if g = 0 this reduces to
2 T

x2
+ 2 T

y2
= 0
and this is Laplace's equation. In any case, the resulting equation must be solved subject to suitable boundary conditions and the result is the function T(x,y).
Finite difference methods easily and readily yield approximate solutions to steady state heat conduction problems. Let's start with Laplace's equation in a rectangular plate of width a and height b subject to Dirichlet conditions as follows
T(x,0) = 0

T(0,y) = 0

T(a,y) = 0
and
T(x,b) = f(x)
First we create a rectangular mesh of points xi, i=1,2,3...,N and yj, j=1,2,3,...,M. Nodes x1, xN, y1 and yM are boundary nodes while the rest are interior nodes. The value of temperature obtained by numerical approximation at a specific nodal location (xi,yj), will be denoted by T(xi,yj) = Ti,j.
For the interior nodes a second order accurate finite difference approximation of Laplace's equation is
Ti-1,j - 2 Ti,j + Ti+1,j

Dx2
+ Ti,j-1 - 2 Ti,j + Ti,j+1

Dy2
= 0
Equations for the boundary nodes are simply
T(xi,0) = Ti,1 = 0
for i = 1,2,...,N
T(0,yj) = T1,j = T(xN,yj) = TN,j = 0
and
T(xi,yM) = Ti,M = f(xi)
The above is a system of N ×M simultaneous linear algebraic equations which must be solved to obtain the approximate values of the temperatures at the interior nodes. The system of equations is sparse but the structure may not be very nice depending on how the nodes are numbered. Direct methods of solution (e.g. Gauss elimination) could be employed but there are not efficient. Iterative solution techniques are much more effective, particularly when N and M are large, and are widely employed. Methods for the solution of the resulting system of interlinked simultaneous algebraic equations are described next.
To simplify, consider the case of Laplace's equation in a square plate a = b with N = M and a uniform mesh with Dx = Dy = d. With this, the finite difference equation for the interior node (i,j) becomes
Ti-1,j + Ti,j-1 - 4 Ti,j + Ti+1,j + Ti,j+1 = 0
This is sometimes called the five point formula.
To obtain a banded matrix the mesh points must be relabeled sequentially from left to right and from top to bottom. As mentioned before, the resulting system can be solved by Gaussian elimination if N and M are small and by point iteration when they are large. Moreover, point iteration algorithms are very easy to implement. Further solution efficiency is possible by employing a line by line combined TDMA and line iteration solution method in exchange for a slight increase in algorithmic complexity.
All iterative methods require an initial guess for the temperature field but the final result must be independent of the guess. The iteration converges when the calculated temperature values at one iteration differ little from those obtained at the subsequent iteration. Convergence indicates that the approximate solution closely satisfies the finite difference analog of the original partial differential equation.
In the Jacobi point iteration method the interior nodes are visited sequentially and an improved guess of the value of Ti,j is calculated using the previous guess values of all neighboring nodes. If Ti,j represents the new value and Ti,jo the old one the algorithm is
Ti,j = Ti-1,jo + Ti,j-1o + Ti+1,jo + Ti,j+1o

4
Note that in the Jacobi method, the temperature field is only updated once all the nodes have been visited.
Typically, nodes are visited from left to right and from bottom to top. I.e. (2,2),(3,2),...,(N-1,2),(2,3),(3,3),...,(N-1,M-1). Note that in this scheme, when visiting node (i,j) nodes (i-1,j) and (i,j-1) would have already been visited (and improved guesses Ti-1,j and Ti,j-1 would be available. The Gauss-Seidel iteration takes advantage of this and uses calculated values as soon as they are available. The algorithm is
Ti,j = Ti-1,j + Ti,j-1 + Ti+1,jo + Ti,j+1o

4
As a result of updating as the calculation proceeds, the Gauss-Seidel method converges at a faster rate than the Jacobi method.
It is possible to increase even more the speed of convergence by using the Successive Over relaxation (SOR) method. First, note that the Gauss-Seidel algorithm can be rewritten as
Ti,j = Ti,jo + Ti-1,j + Ti,j-1 - 4 Ti,jo + Ti+1,jo + Ti,j+1o

4
The second term on the right hand side is the correction or displacement in the value of the temperature at node (i,j) from Ti,jo to Ti,j.
The SOR method converges at a faster rate because larger corrections are done at each iteration, i.e. an acceleration parameter w is introduced such that
Ti,j = Ti,jo + w Ti-1,j + Ti,j-1 - 4 Ti,jo + Ti+1,jo + Ti,j+1o

4
where the so called, over relaxation parameter w is 1 w 2.

1.2  Cylindrical and Spherical Coordinates

The procedures used to obtain finite difference formulae for cylindrical and spherical coordinates are very similar. Here, it will be illustrated for the case of cylindrical coordinates only. Consider steady state heat conduction inside a short, hollow cylinder (inner radius a, outer radius b, height h, constant thermal conductivity k) under Dirichlet boundary conditions and such that the temperature is independent of the angular location q (axisymmetric system) and without any internal heat generation. Laplace's equation in cylindrical coordinates (r,z) becomes
1

r

r
(r T

r
) + 2 T

z2
= 0
Let the boundary conditions be given by
T(r,0) = 0

T(a,y) = 0

T(b,y) = 0
and
T(r,h) = f(r)
An approximation by finite differences can be readily obtained using the same procedure described above. However, the resulting finite difference formulae are slightly different. Specifically, introducing a uniform mesh of nodes spaced by Dr, Dz and located at ri, zj with i=1,2,...,N and j=1,2,...,M, and T(ri,zj) Ti,j,
1

r

r
(r T

r
) 1

ri
(r T/r)|ri+Dr/2 - (r T/r)|ri-Dr/2

Dr
=
= 1

ri
(ri+Dr/2)( Ti+1,j-Ti,j

Dr
)-(ri-Dr/2)( Ti,j-Ti-1,j

Dr
)

Dr
Since
2 T

z2
Ti,j-1 - 2 Ti,j + Ti,j+1

Dz2
the finite difference analog of Laplace's equation in circular cylindrical coordinates for axisymmetric systems becomes
1

ri
(ri+Dr/2)( Ti+1,j-Ti,j

Dr
)-(ri-Dr/2)( Ti,j-Ti-1,j

Dr
)

Dr
+ Ti,j-1 - 2 Ti,j + Ti,j+1

Dz2
= 0
Assuming now, for simplicity that Dr = Dz = d, rearrangement yields the corresponding five point formula
ri- d

2

ri
Ti-1,j + Ti,j-1 - 4 Ti,j +
ri+ d

2

ri
Ti+1,j + Ti,j+1 = 0
The Gauss-Seidel point iteration algorithm becomes in this case
Ti,j =
ri- d

2

ri
Ti-1,j + Ti,j-1 +
ri- d

2

ri
Ti+1,jo + Ti,j+1o

4

1.3  Derivative Boundary Conditions

Newmann, Robin and the fourth power radiation boundary conditions involve derivatives of the temperature. The simplest approach is to use one sided, finite differences. The tradeoff is that the resulting formula are only first order accurate and a fine mesh is required for good results.
For instance, if the boundary condition on a rectangular plate at point (x=a,y) is
T

x
= 0
a finite difference approximation of this is
T

x
TN,j - TN-1,j

d
= 0
i.e. TN,j = TN-1,j.
Alternatively, if the boundary condition on a rectangular plate at point (x=a,y) is
-k T

x
= qa
a finite difference approximation is
-k T

x
-k TN,j - TN-1,j

d
= qa
If the boundary condition is instead
-k T

x
= h (T - T)
a finite difference approximation is
-k T

x
-k TN,j - TN-1,j

d
= h (TN,j - T)

2  Finite Volume Method

Consider the problem of estimating the temperature T(x,y,z) in a three-dimensional brick x [0,X], y [0,Y], z [0,Z] undergoing steady state heat conduction with constant internal heat generation. The heat equation for this case becomes

x
( k T

x
) +

y
( k T

y
) +

z
( k T

z
) + g = 0
subject to specified conditions at the boundaries of the domain.
The finite volume formulation regards the brick as composed of a set of adjoining volumes each containing a node. The dimensions of the typical volume are Dx ×Dy ×Dz and the typical node P has six neighbors E, W, N, S, T and B. The finite volume method proceeds by integrating the heat equation over the typical volume element and the result is
aP TP = aE TE + aW TW + aN TN + aS TS + aT TT + aB TB + b
where
aE = ke Dy Dz

dxe

aW = kw Dy Dz

dxw

aN = kn Dz Dx

dyn

aS = kw Dz Dx

dys

aT = kt Dx Dy

dzt

aB = kb Dx Dy

dzb

b = g Dx Dy Dz
and
aP = aE + aW + aN + aS + aT + aB
Note that all the above expressions can all be generically written as
ap TP =
anb Tnb + b
where the summation contains just TE and TW in the case of uni-dimensional systems, contains TE, TW, TN and TS in the case of 2D systems and all TE, TW, TN, TS, TT and TB in the case of three dimensional systems.
The discrete equation derived using the FV method can be rearranged to produce a very appealing and physically meaningful expression. Consider a two dimensional system. The generic discrete form of the conservation equation is
Jw Aw - Je Ae + Js As - Jn An + g DV = 0
where the Ji's are the heat fluxes and Ai's are the cross-sectional areas of the finite volume faces through which heat enters/leaves the finite volume. Therefore the products Ji Ai are the rates of transport of energy (J/w) through the various control volume faces. The transport rates are given as
Je Ae = [( dxe-

GP
) + ( dxe+

GE
)]-1(TP-TE) Ae

Jw Aw = [( dxw-

GW
) + ( dxw+

GP
)]-1(TW-TP) Aw

Jn An = [( dyn-

GP
) + ( dyn+

GN
)]-1(TP-TN) An

Js As = [( dys-

GS
) + ( dys+

GP
)]-1(TS-TP)
where Gi = ki.
Many special forms can be obtained by simplification of the above. For instance, for a rectangle with constant thermal properties, without internal heat generation at steady state and using Dx = Dy one obtains
TP = 1

4
[TE + TW + TN + TS]
which is identical to the expression obtained before using the method of finite differences.

3  Finite Element Method

Consider a three dimensional body occupying the volume W, inside which heat is generated at a rate g and whose outer surface G is subjected to three distinct boundary conditions on portions S1, S2 and S3. In rectangular Cartesian coordinates the governing equation is
·(k T ) + g =

x
(k T

x
) +

y
(k T

y
) +

z
(k T

z
) + g = 0
and the boundary conditions are
T = Tb
on S1,
k T ·n + q = 0
on S2, and
k T ·n + h (T - T) = 0
on S3.
The variational statement of the above problem involves the minimization of the functional
I(T) = 1

2



W 
[k ( T

x
)2 + k ( T

y
)2 + k ( T

z
)2 - 2 g T] dW+


S2 
q T ds + 1

2



S3 
h (T-T)2 ds
In the finite element method the domain W is first subdivided into a number of regions with simple polygonal shape and the temperature within each element Te is expressed as a linear combination of the nodal values Ti; i=1,2,...,r by means of finite element basis functions, i.e.
Te = r

i=1 
fiTi = fTT
where f is the basis function matrix and T is the vector of nodal temperatures. In the Ritz formulation of the finite element method one looks for the values of Ti that minimize I(T) e=1N Ie(Te), i.e.
I

Ti
= N

e=1 
Ie

Ti
= 0
for i=1,2,....,M where M is the total number of nodes in the finite element mesh and N is the total number of elements. The result is a set of M linear algebraic equations for the M called the finite element equation, namely,
K T = f
Solving the system yields the individual components of the solution vector T (i.e. the nodal values of the temperature).
The matrix K is sometimes called the global conductance matrix and is analogous to the stiffness matrix encountered when applying the finite element method in structural mechanics. The conductance matrix is given by
K =


W 
BT D B d W+


S3 
h fT f ds
where
B =







f1

x
f2

x
...
fr

x
f1

y
f2

y
...
fr

y
f1

z
f2

z
...
fr

z








and
D =



k
0
0
0
k
0
0
0
k




Finally, the forcing vector f is given by
f =


W 
g fT d W-


S2 
q fT ds +


S3 
h T fT ds
Note that for insulated boundaries (i.e. h = q = 0) the corresponding integrals vanish and no term appears reflecting this contribution in the finite element equation (natural boundary condition).
In order to assure convergence of the finite element model as the element size decreases, the basis functions must satisfy the standard requirements of compatibility and completeness, i.e. the basis functions selected must provide for continuity of the temperature at the interface between any two elements as well as for the continuity of temperature and heat flux inside each element. Standard polynomials of order n satisfy the requirement.
As with one-dimensional systems, the simplest global basis functions for an individual node are those that have the value of 1 at the node, decay in a linear fashion to zero as all neighboring nodes are approached and are identically zero for all nodes. Higher accuracy with fewer elements can be obtained by using higher order polynomials. Quadratic Lagrangian polynomials are commonly employed



File translated from TEX by TTH, version 3.85.
On 30 Sep 2009, 16:33.

Updated: 2009-09-30, 16:33