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.

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

¶x^{2}

+

¶^{2} T

¶y^{2}

+

g

k

= 0

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

¶^{2} T

¶x^{2}

+

¶^{2} T

¶y^{2}

= 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
x_{i}, i=1,2,3...,N and y_{j}, j=1,2,3,...,M. Nodes
x_{1}, x_{N}, y_{1} and y_{M} are boundary nodes
while the rest are interior nodes. The value of temperature
obtained by numerical approximation at a specific nodal location
(x_{i},y_{j}), will be denoted by T(x_{i},y_{j}) = T_{i,j}.
For the interior nodes
a second order accurate finite difference approximation of Laplace's
equation is

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

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 T_{i,j} is calculated using
the previous guess values of all neighboring nodes. If
T_{i,j} represents the new value and T_{i,j}^{o} the old one
the algorithm is

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 T_{i-1,j} and T_{i,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

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

The second term on the right hand side is the correction or
displacement in the value of the temperature at node (i,j)
from T_{i,j}^{o} to T_{i,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

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

¶z^{2}

= 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 r_{i}, z_{j} with i=1,2,...,N and j=1,2,...,M, and T(r_{i},z_{j}) » T_{i,j},

1

r

¶

¶r

(r

¶T

¶r

) »

1

r_{i}

(r ¶T/¶r)|_{ri+Dr/2} - (r ¶T/¶r)|_{ri-Dr/2}

Dr

=

=

1

r_{i}

(r_{i}+Dr/2)(

T_{i+1,j}-T_{i,j}

Dr

)-(r_{i}-Dr/2)(

T_{i,j}-T_{i-1,j}

Dr

)

Dr

Since

¶^{2} T

¶z^{2}

»

T_{i,j-1} - 2 T_{i,j} + T_{i,j+1}

Dz^{2}

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

1

r_{i}

(r_{i}+Dr/2)(

T_{i+1,j}-T_{i,j}

Dr

)-(r_{i}-Dr/2)(

T_{i,j}-T_{i-1,j}

Dr

)

Dr

+

T_{i,j-1} - 2 T_{i,j} + T_{i,j+1}

Dz^{2}

= 0

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

r_{i}-

d

2

r_{i}

T_{i-1,j} + T_{i,j-1} - 4 T_{i,j} +

r_{i}+

d

2

r_{i}

T_{i+1,j} + T_{i,j+1} = 0

The Gauss-Seidel point iteration algorithm becomes in this case

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

»

T_{N,j} - T_{N-1,j}

d

= 0

i.e. T_{N,j} = T_{N-1,j}.
Alternatively, if the boundary condition
on a rectangular plate at point (x=a,y) is

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

Note that all the above expressions can all be generically written as

a_{p} T_{P} =

å

a_{nb} T_{nb} + b

where the summation contains just T_{E} and T_{W} in the case of
uni-dimensional systems, contains T_{E}, T_{W}, T_{N} and T_{S}
in the case of 2D systems and all T_{E}, T_{W}, T_{N}, T_{S}, T_{T}
and T_{B} 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

J_{w} A_{w} - J_{e} A_{e} + J_{s} A_{s} - J_{n} A_{n} + g DV = 0

where
the J_{i}'s are the heat fluxes and
A_{i}'s are the cross-sectional areas of the finite volume faces
through which heat enters/leaves the finite volume.
Therefore the products J_{i} A_{i} are the rates of transport of
energy (J/w) through the various control volume faces.
The transport rates are given as

J_{e} A_{e} = [(

dx_{e-}

G_{P}

) + (

dx_{e+}

G_{E}

)]^{-1}(T_{P}-T_{E}) A_{e}

J_{w} A_{w} = [(

dx_{w-}

G_{W}

) + (

dx_{w+}

G_{P}

)]^{-1}(T_{W}-T_{P}) A_{w}

J_{n} A_{n} = [(

dy_{n-}

G_{P}

) + (

dy_{n+}

G_{N}

)]^{-1}(T_{P}-T_{N}) A_{n}

J_{s} A_{s} = [(

dy_{s-}

G_{S}

) + (

dy_{s+}

G_{P}

)]^{-1}(T_{S}-T_{P})

where G_{i} = k_{i}.
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

T_{P} =

1

4

[T_{E} + T_{W} + T_{N} + T_{S}]

which is identical to the expression obtained before using the method of finite differences.

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 S_{1}, S_{2} and S_{3}. 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 = T_{b}

on S_{1},

k ÑT ·n + q = 0

on S_{2}, and

k ÑT ·n + h (T - T_{¥}) = 0

on S_{3}.
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+

ó õ

S_{2}

q T ds +

1

2

ó õ

S_{3}

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 T^{e} is expressed as a linear combination of
the nodal values T_{i}; i=1,2,...,r by means of finite element basis functions, i.e.

T^{e} =

r å
i=1

f_{i}T_{i} = f^{T}T

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 T_{i} that minimize
I(T) » å_{e=1}^{N} I^{e}(T^{e}), i.e.

¶I

¶T_{i}

=

N å
e=1

¶I^{e}

¶T_{i}

= 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,

KT = 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

B^{T}DB d W+

ó õ

S_{3}

h f^{T} f ds

where

B =

é ê ê ê ê
ê ê ê ë

¶f_{1}

¶x

¶f_{2}

¶x

...

¶f_{r}

¶x

¶f_{1}

¶y

¶f_{2}

¶y

...

¶f_{r}

¶y

¶f_{1}

¶z

¶f_{2}

¶z

...

¶f_{r}

¶z

ù ú ú ú ú
ú ú ú û

and

D =

é ê ê
ê ë

k

0

0

0

k

0

0

0

k

ù ú ú
ú û

Finally, the forcing vector f is given by

f =

ó õ

W

g f^{T} d W-

ó õ

S_{2}

q f^{T} ds +

ó õ

S_{3}

h T_{¥} f^{T} 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
T_{E}X
by
T_{T}H,
version 3.85. On 30 Sep 2009, 16:33.