DOI: 10.1615/AtoZ.n.numerical_heat_transfer


Numerical heat transfer is a broad term denoting the procedures for the solution, on a computer, of a set of algebraic equations that approximate the differential (and, occasionally, integral) equations describing conduction, convection and/or radiation heat transfer.

The usual objective in any heat transfer calculation is the determination of the rate of heat transfer to or from some surface or object. In conduction problems, this requires finding the temperature gradient in the material at its surface. In convection problems, the temperature gradient in a fluid flowing over a surface is needed to find the heat flux at that surface. In both cases, the determination of the complete temperature distribution in the region of interest is needed as a first step, and in convection one must also find the velocity distribution. Thus, a full solution of the energy equation and perhaps also the equations of motion is required. These are partial differential equations, possibly coupled.

Radiation is somewhat different, involving surfaces separated (in general) by a fluid which may or may not participate in the radiation. If it is transparent, and if the temperatures of the surfaces are known, the radiation and convection phenomena are uncoupled and can be solved separately. If the surface temperatures are not specified, but are to be found as part of the solution, or if the fluid absorbs or emits radiant energy, then the two phenomena are coupled. In either case, the solution of the algebraic and/or integral equations for radiation is required in addition to that of the differential equations for convection and conduction.

The equations describing heat transfer are complex, having some or all of the following characteristics: they are nonlinear; they comprise algebraic, partial differential and/or integral equations; they constitute a coupled system; the properties of the substances involved are usually functions of temperature and may be functions of pressure; the solution region is usually not a simple square, circle or box; and it may (in problems involving solidification, melting, etc.) change in size and shape in a manner not known in advance. Thus analytical methods, leading to exact, closed form solutions, are almost always not available.

Two approaches are possible. In the first, the equations are simplified — for example, by linearization, or by the neglect of terms considered sufficiently small, or by the assumption of constant properties, or by some other technique until an equation or system of equations is obtained for which an analytical solution can be found. It could be said that an exact analytical solution will be obtained for an approximate problem. The solution will, to some extent, be in error, and it will not normally be possible to estimate the magnitude of this error without recourse to external information such as an experimental result.

The second approach is to use a numerical method. To do this, the continuous solution region is, in most methods, replaced by a net or grid of lines and elements. The solution variables—temperature, velocity, etc.—are not obtained at all of the infinite number of points in the solution region, but only at the finite number of nodes of the grid, or at points within the finite number of elements. The differential equations are replaced by set of linear (or, rarely, nonlinear) algebraic equations, which must and can be solved on a computer.

There is also an error in this approach: for example, if derivatives are replaced by finite differences, only an approximate value for the derivative will be obtained, in principle, if the problem is well-posed and if the solution method is well-designed, this error will approach zero as the grid is made increasingly fine. In practice, a fine but not infinitesimal grid must be used. The error can, in general, be estimated and it can also be reduced at the price of increased effort (meaning increased computer time). It could now be said that an approximate solution will be obtained for an exact problem.

In this entry, emphasis will be given to methods for conduction and convection problems, with only a brief mention of radiation. Attention will be limited to incompressible fluids except when buoyancy is important, in which case the Boussinesq approximation will be made. The ideas described can often be used in more general situations, but other methods are also available, especially for high speed flow, which are beyond the present scope.

The Governing Equations for Conduction and Convection

For an incompressible Newtonian fluid, the continuity, momentum and energy equations are


where ( )T denotes vector transpose, rg is a buoyancy (body force) term, and denotes the rate of internal generation of heat per unit volume. These equations represent a system of five equations in three dimensions, or four equations in two dimensions. For conduction in solids, DT/Dt in (3) is replaced by ∂T/∂t, and (1) and (2) are not relevant.

The thermodynamic and transport properties ρ, Cp, λ and η of most materials are in general dependent on temperature and (particularly in the case of ρ) also on pressure. It is often possible to neglect this dependency, with considerable simplification to the computational procedure. If buoyancy is important, density must be considered as a variable and an equation of state is required. A common assumption for low speed flow (small Mach number, together with other restrictions) is to invoke the Boussinesq approximation, in which the variation of density is neglected except in the body force term of the momentum equations, namely the last term of (2). For a Boussinesq fluid, it is assumed that ρ = f(T). The usual assumption is a linear relationship:


where β is the volumetric expansion coefficient, and the subscript denotes a reference state. However water near 4°C, and some other substances, may require a more complex equation.

Equations (1)–(3) are written in terms of what are called the primitive variables: ρ, u, T and p. Because of difficulties which can arise with the specification and implementation of pressure boundary conditions, they are sometimes transformed into equations for the derived variables vorticity ζ and vector potential ψ, or their two-dimensional counterparts ζ and the stream function ψ, together in either case with ρ and T. The relationships between the primitive and derived variables are:


Vorticity and vector potential are thus related by


and taking the curl of the momentum Equation (2) leads to the vorticity transport equation


As well as the elimination of pressure and the need to solve the continuity equation (which is automatically satisfied by the introduction of the vector potential or stream function), this formulation has, in two dimensions, the advantage of requiring the solution of just three equations (3), (6) and (7) instead of four. In three dimensions, however, the number is actually increased from five to seven, but the elimination of pressure is still sometimes seen as making the transformation worthwhile.

Equations (2) and (3) are written in what is called the flux or gradient formulation. With the help of Eq. (1), they can also be written in the conservation form


which is mathematically (and physically) equivalent, but which is sometimes felt to offer a computational advantage because this formulation explicitly recognizes the principles of conservation of momentum and energy.

Non-dimensional equations

It is almost always advantageous to transform the dimensional equations into a non-dimensional form, by scaling the dimensional variables with respect to known reference or characteristic quantities. First, if a parametric solution is needed—to study the effects of changes in density, velocity, etc. over a range of values—then the number of independent parameters is usually reduced by such a transformation. Second, useful physical insight can often be garnered in the process. If characteristic quantities are selected for length, velocity, pressure, time, etc., the equations can be rewritten in terms of dimensionless variables, and the coefficients in the new equations involve parameters such as Re, Pr and Ra (or Gr).

For example, let θ = (T –T0)/(T1 - T0), where TFor example, let and T1 are two specified temperatures (perhaps a wall temperature and a free stream temperature), and use U, L and L/U as the characteristic velocity, dimension and time scale of the problem. Then the energy and vorticity transport equations (3) and (7) become


in which, for simplicity, properties have been taken as constants (other than ρ in the body force term, for which (4) has been used), Q denotes the non-dimensional heat generation rate

Boundary conditions

These equations must be accompanied by boundary conditions appropriate to the particular problem. For transient problems, an initial condition is required, and in all problems boundary values for all variables are needed. In heat transfer, three types of condition on T are encountered:

  • specified temperature: T = Tw at a boundary, where the wall temperature Tw is known; this is called a Dirichlet condition;

  • specified heat flux: –λ(∂T/∂n)w = qw, where n is normal to the surface, and the wall heat flux qw is known; this is called a Neumann condition;

  • specified heat transfer coefficient α: –λ(∂T/∂n)w = α(Tw – Ti); this is called a mixed condition.

A variation of the second condition occurs when the heat flux is by radiation from a source of known temperature Ts: −λ(∂T/∂n)w = in which F is a view factor and α is the Stefan-Bolzmann constant. When either the second or third conditions apply, the wall temperature Tw must be found as part of the solution.

In convection problems, boundary values of u and/or its derivatives are also required. Typically, if the flow in the x-direction is parallel and fully developed, ∂u/dx = v = w = 0. At the free surface of a liquid, tangential stresses are often taken to be zero and ∂u/dx = 0. On any solid surface, the no-slip condition must be applied (the tangential velocity components are equal to those of the surface) and in the absence of suction or blowing, the normal velocity component is zero. If derived variables are used, boundary conditions for ζ and ψ may be found from those for u:

Turbulence modeling

Most flows of practical importance are turbulent, which adds a considerable burden to their computation. In principle, turbulent flow must obey the Navier-Stokes equations. However, their use to calculate it is impractical because of the huge demands which would be made on computer memory and speed. A common alternative is to use the two equation k — ε turbulence model [Rodi (1993)], which starts with the time-averaged equations of motion in which all variables are represented as the sum of a mean and a fluctuating component. Using tensor notation (summation over repeated suffixes), these equations are


where overbars denote time-averaged quantities and primes denote turbulent fluctuations from the mean. S is a volumetric source term. The terms such as , when multiplied by the density, are called Reynolds stresses, and are unknown quantities. The closure problem of the k – ε turbulence model is based upon relating the Reynolds stresses to the mean quantities by


where ηt is the turbulent or eddy viscosity and k = is the kinetic energy of the turbulence. In turn, ηt is given by


where ε is the rate of dissipation of turbulent energy

and is seen to be proportional to the fluctuating vorticity. Finally, transport equations are constructed for k and ε:


The quantity σt is the turbulent Prandtl number; its value lies in the range 0.5–0.9 depending upon the nature of the flow. The first terms on the right hand sides of (15) and (16) represent diffusive transport of k and ε, respectively. The term labeled P in (15) is usually designated "production": it represents the transfer of kinetic energy from the mean to the turbulent motion. The buoyancy term G represents an exchange between the turbulent kinetic energy k and potential energy. The last term in (15) is viscous dissipation. In (16), the second and third terms on the right denote generation of vorticity by vortex stretching and viscous destruction of vorticity, respectively. Rf is the flux Richardson number, defined as –G/P.

The various constants in these equations are commonly set to the values given in the following table.

*For vertical or horizontal layers, respectively.

Solution Methods for Conduction and Convection

The most common solution methods are finite difference and finite volume methods on the one hand, and weighted residual methods on the other. Among the latter, the most common are finite element and spectral methods. The approach used depends upon whether the flow is steady or unsteady, and on whether primitive or derived variables are being used. A necessarily brief description of some common methods follows.

Finite difference methods

To illustrate a finite difference solution, consider the equation


with suitable initial and boundary conditions, where α = λ/ρcp is the thermal diffusivity of the material. For a conduction problem, u and v would be zero. For a convection problem there would be the need also to solve the equations of motion for u and v; these equations might be coupled to (17) if buoyancy is significant. We will consider (17) as it stands, assume that u and v are known, and suppose at first that T is also known at all boundaries of the solution region.


A rectangular grid is superimposed over the solution region, a portion of which is shown in Figure 1 The intention now is to find approximations to the values of the solution variables at the nodes of the grid.

Portion of a computational grid, showing the numbering system.

Figure 1. Portion of a computational grid, showing the numbering system.

The vertical grid lines are numbered i = 1, 2, 3 ... M along the x-axis and the horizontal lines are numbered j = 1, 2, 3 ... N along the y-axis. A superscript n is used to indicate time: t = t0 + nk, where t0 is the initial time (usually set to zero) and k is the size of the time step. Thus the temperature at the point (i,j) and at time n is denoted by .

FTCS scheme

The equations to be solved contain terms such ∂T/∂t, u ∂T/∂x and ∂2T/∂y2. Using the simplest difference operators, namely forward in time and central in space, these terms may be represented by and


respectively, where hx and hy are the grid or mesh sizes in the x and y directions, respectively [e.g., de Vahl Davis (1986)]. It is assumed for simplicity that these quantities are constant (although not necessarily equal); often, this is not the case. Expressions for the truncation errors O(k), etc. can be found from a Taylor series expansion of the finite difference approximations.

With these and similar finite difference approximations, and neglecting the truncation errors, (16) may be written


Equations like (19) may be written for every one of the (M-2)(N-2) internal mesh points. This is known as the forward time-central space (FTCS) scheme. (19) is a finite difference approximation (FDA) to (17), the partial differential equation (PDE).

Explicit solution

One way to solve (19) is to rearrange it to express as a function of the values of T at time n. Since the initial condition provides values for T at all mesh points at time t = 0 (i.e., n = 0), the new or n = 1 values of T may be calculated explicitly at each internal mesh point. Once this is complete, the values of T at n = 2 may be found, and so on. This is the easiest way of solving (19) and, for equations as simple as this, it is adequate. It does not, however, always work.

Consistency, convergence and stability

There are three important characteristics of any finite difference scheme: consistency, convergence and stability. A scheme is said to be consistent if the truncation errors such as those appearing in (18) tend to zero as the grid and time steps tend to zero, so that the FDA approaches the PDE. Clearly, those in (18) satisfy this requirement, but there are some methods — notably, the DuFort-Frankel scheme — which are not unconditionally consistent. Such methods must be used with considerable care.

Convergence goes one step further: it requires the solution of the FDA to approach the solution of the PDE as the step sizes go to zero. Convergence does not always follow from satisfying consistency, and it is in fact not normally possible to determine directly whether a method will be convergent.

The third characteristic of a scheme is its stability. A scheme is stable if any errors such as round-off errors in the computer decay to zero (more strictly, remain bounded) as the calculations progress, i.e., as n → ∞. Stability can be examined by several methods including the von Neumann or Fourier Series method and the Matrix method. For the FTCS scheme (19) for solving (17), it can be shown that stability requires αk/h2 ≤ 0.5, where h = min(hx, hy). Thus if the mesh is fine (small h) or if α is large, a small value of k and therefore a long computation time will be required.

A further limitation on stability can exist; the Courant or CFL condition, which is that uk/h must be no greater than unity at any mesh point. The physical significance of this condition is that the fluid should not move a distance greater than one mesh length during one time step, and it applies to explicit schemes for hyperbolic flows, e.g., for flows with little or no diffusion.

Although a direct determination of convergence is not usually possible, the Lax equivalence theorem [Richtmeyer and Morton (1967), p. 45] provides some help: "Given a properly posed linear initial value problem and a finite difference approximation to it that satisfies the consistency condition, stability is the necessary and sufficient condition for convergence". Unfortunately, most convection heat transfer problems are not linear, and the Lax theorem at best provides guidance to necessary but not always sufficient conditions.

The word "convergence" is also used in a different sense from that discussed above. In an iterative procedure for the solution of the algebraic FDAs, the error in each estimate of the solution must decrease with each successive iteration. If it does, the procedure is said to be convergent; if it does not, the procedure is of no use.

Coupled systems for convection problems

The FTCS method can be immediately and easily extended to convection problems using the derived variable formulation, Eqs. (6) and (7). The main complication arises from the need to determine boundary conditions on vorticity. A common formula is that known as the Woods method [Roache (1972)].

Boundary conditions on the vector potential (or, in two dimensions, the stream function) must be obtained from those on velocity; see Hirasaki and Heliums (1968, 1970).

From some given initial state, the FDA of (7) may be used to obtain new values of vorticity. The FDA of (6) is then solved for the vector potential or stream function, and new values for the velocity components obtained from the FDA of (5). The calculations for the time step are completed by using (19) to obtain new values of temperature.

Implicit solution

The limitation on the time step k imposed by the stability criterion led to the development of implicit methods which require more work per time step, but allow the use of larger values of k and therefore fewer time steps. Implicit methods involve the construction of FDAs in which some or ail of the spatial derivatives are evaluated at time level n + 1, rather than at n.

For example, in the Crank-Nicolson method, the FDAs of ∂T/∂x and ∂2T/∂x2 are written



respectively. These averages of the FDAs at n and n + 1 can be interpreted as central difference approximations to the derivatives at an intermediate time n + 1/2.

The application of (20) and the corresponding FDAs in the y direction will introduce unknown values of T at time level n + 1 to the equations to be solved. In two dimensions there will be five unknowns in each equation, namely the values of T at the grid point (i,j) and at the four surrounding grid points. In three dimensions there will be seven unknowns.

The algebraic equations to be solved thus constitute a system. can no longer be written explicitly in terms of known quantities. All the new values of T at time level n + 1 must be found by a solution of this system of equations, i.e., implicitly.

The solution may be found iteratively, using, for example, the Gauss-Seidel or successive over relaxation (SOR) methods for large linear systems, or it may be found directly, using an elimination method. Gaussian elimination or a similar method is usually impractical for the very large systems which result from the use of a fine mesh, and even SOR can be extremely demanding. The alternating direction implicit (ADI) method is frequently employed to overcome this difficulty.

In ADI, the derivatives are written implicitly, i.e., at time n + 1, in just one of the space directions, while those in the other direction(s) are written explicitly, i.e., at time n. Thus there are only three unknowns in each equation; for example, in the x-direction, they are the time n + 1 values of T at the points i + 1, i and i – 1. The system of equations is therefore tridiagonal. For such a system, a highly efficient elimination algorithm known as the Thomas algorithm [see, e.g., Roache (1972)] is available. The values of T thus computed are regarded as intermediate values. The equations are now written implicitly in the next coordinate direction, and explicitly in the other direction(s), using the intermediate values just computed. The direction in which implicit differencing is applied cycles through the two (or three) coordinates in turn. One time step is regarded as having been completed when the two (or three) intermediate steps have been taken. It is interesting to note that each intermediate ADI step is unstable but that the complete sequence of two (or, in three dimensions, three) steps is unconditionally stable.

Upwind differences and false diffusion

A stability limitation on the FTCS method has been noted arising from the use of central differences for the spatial derivatives. In highly convective flows, an additional, very severe restriction also occurs, associated with the value of a Reynolds (or Peclet) number based on the local velocity and the mesh size. If, in (19), Reh = uh/α > 2, unphysical oscillations in time and/or space will develop, leading eventually to a total breakdown of the solution. To overcome this problem, the idea of upwind differences has been developed: one-sided differences in which the direction of differencing is determined by the sign of the corresponding velocity. Thus, for u > 0, ∂u/∂x is approximated by (ui – ui−1)/hx. This device is very effective in preventing the instability but introduces a truncation error which is proportional to uh(∂2u/∂x2) and is thus like a diffusion term. In fact, for small values of viscosity and large values of velocity, this "false" diffusion can swamp the true diffusion, preventing an accurate solution from being obtained. Second or higher order upwinding, and other differencing schemes, have been developed to overcome this deficiency. Among such methods is the QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics) and other related methods, originated by Leonard (1988).

Boundary conditions

The implementation of Dirichlet boundary conditions is straight forward; indeed, it happens automatically when an FDA is constructed at a mesh point adjacent to a boundary. Neumann and mixed boundary conditions, involving an unknown boundary value but a known boundary gradient, are more complex. Consider the boundary point (l,j) shown in Figure 1, and suppose that the condition


must be imposed there. There are several ways in which this may be achieved. The simplest is to replace the derivative ∂T/∂x by the first order forward FDA (T2,j – Ti,j)/hx, so that (21) may be written


After the values of T at internal mesh points have been computed, the boundary values can be found from expressions such as (22). An improved accuracy may be obtained by using a second order forward FDA of the derivative. Better still is to use a central FDA at the point (l,j). This requires the introduction of a hypothetical mesh point outside the solution region, the value of T at which is obtain from the boundary condition. The boundary value of T can then be computed as part of the solution procedure.

Weighted residual methods*

* The outline given in this section draws heavily on the text of Fletcher (1991).

The characteristic feature of finite difference methods is that values are sought and obtained for the unknown quantities (velocity, temperature, etc.) at the node points of a computational mesh. The variation of these quantities between node points is implied by the order of the approximations used, but in practice, no use is made of this knowledge.

In weighted residual methods, on the other hand, the assumption is made that the solution can be represented by an appropriate combination of analytical functions over the whole of the solution region. In one dimension, for example, a solution of the form


would be assumed, where ai(t) are coefficients which in general may be time-dependent and which must be evaluated during the solution process, and where φi(x) are trial functions whose form might be chosen, at least in past, to be compatible with some anticipated features of the solution. The challenge is to choose appropriate functions f and then to find the coefficients a such that (23) satisfies the given differential equation and its boundary and initial conditions as well as possible at all points in the solution domain.

Suppose that the differential equation to be solved is written


where L is a differential operator, and the overbar denotes the exact solution. If an approximate solution T is inserted into (24) it will not, in general, evaluate to zero but to some residual R:


The coefficients ai(t) are found by requiring that the integral of the weighted residual over the solution region is zero;


This represents a system of M equations, m = 1, 2, 3, ... M, from which the coefficients ai may be obtained. If these coefficients are (as implied here) functions of t, the system is comprised of ordinary differential equations. In steady problems, a system of algebraic equations is obtained.

With different choices of weighting functions Wm, different weighted residual methods are obtained. In the finite volume method, the solution region is divided into M subdomains Dm within which Dm = 1 and beyond which Dm = 0. In the collocation method, Wm = δ(x - xm), where δ is the Dirac delta function; substitution into (26) then leads to the conclusion that finite difference methods are a form of collocation method. In the least squares method, Wm = ∂R/∂am, from which follows the requirement that the integral of R2 over the solution region will be a minimum. Finally, in the Galerkin method, the weighting functions are chosen from the same family as the trial functions themselves, and these are usually linear or low order polynomial functions.

We will now examine the finite volume and Galerkin methods in a little more detail.

Finite volume method

Consider the problem of solving Laplace's equation


in a two-dimensional region which, as shown in Figure 2, is divided into a number of quadrilateral elements. The directions of the computational grid lines have been deliberately chosen not to be coincident with the directions of the axes to show the power and generality of the method.

V finite volume mesh.

Figure 2. V finite volume mesh.

We insert (26) into (25), with W = 1, and apply Green's theorem, replacing the integral over the area ABCD with an integral around the contour ABCD, to obtain




Each derivative in (28) must be approximated on each of the four edges of the subdomain ABCD. For example,



and in turn

If the mesh is not too distorted, etc., so that

After further algebraic manipulation, we obtain



etc., and φA, XA, yA, etc. are evaluated as the averages of the respective four surrounding nodal values. Equation (31) is in a form convenient for iterative solution by successive over-relaxation.

It is now apparent that, in implementation, the finite volume method is very similar to the finite difference method. However, the FVM, being derived from an integral form of the governing equations, has better conservation properties than the FDM. It is also capable of handling computational domains of a complex shape more easily. It should also be noted that derivative boundary conditions can be easily implemented by direct substitution into one of the approximations (29).

SIMPLE methods

An important and popular class of methods originated with the SIMPLE algorithm (Semi-Implicit Method for Pressure-Linked Equations) due to Patankar and Spalding (1972) and described in Patankar (1981). They are FVMs, with the discretization of the primitive variable equations performed on a staggered grid as illustrated in Figure 3. The control volume (CV) associated with the point (i,j) is shown shaded. The different variables are identified with different points in the control volume: pressure, temperature and any other scalar with the center of the CV, and the u and v velocity components with the centers of the upstream and downstream faces of the CV in the x and y directions, respectively.

Main grid (solid lines) and staggered grid (broken lines). The horizontal and vertical arrows indicate the grid points for the u and v velocity components.

Figure 3. Main grid (solid lines) and staggered grid (broken lines). The horizontal and vertical arrows indicate the grid points for the u and v velocity components.

The basis of the method is the use of discretized momentum equations to advance from a known (or initial) state at time or iteration number n to estimate new velocities at time or iteration number n + 1. These estimates will not, in general, be accurate, as they do not satisfy continuity. Equations for corrections to the values of the pressure and the velocities are then employed, and the corrected pressure used as the starting point for the next cycle of calculations.

Improvements and extensions to SIMPLE, known as SIMPLER, SIMPLEC, SIMPLEX etc. have been developed since the original publication of the method in 1972.

Galerkin finite element method

The Galerkin FEM is probably the most widely used of the "true" FEMs (in contrast to the FVM), although related methods including collocation, least squares and boundary elements are also used. To illustrate it, consider the Poisson equation


with on the boundaries of a square region. In the Gaierkin FEM, the approximate solution of (32) is written directly in terms of the nodal unknowns, i.e.,


where φi(x,w) are interpolating or approximating functions of some assumed form. Normally, linear or quadratic forms are used.

Consider the portion of the solution region shown in Figure 4. In this two-dimensional example, the local solution is interpolated or approximated, separately, in each of the four elements A, B, C and D which surround the (j,k) node.

Global and local notation.

Figure 4. Global and local notation.

A coordinate system (ξ, η) such that −1 ≤ ξ, η ≤ 1 is introduced for each element, as illustrated for element C in Figure 4. Then a suitable bilinear form for each approximating function φi is


Substitution of (33) into (32), application of the Galerkin formulation (26) and some manipulation yields


where W is the vector of unknown nodal values wi, an element of B is given by


and an element of G is given by


The parameter m ranges over all the internal nodes, i.e., 2 ≤ j ≤ NX − 1, 2 ≤ k ≤ NY − 1. The integrals in (36) are best evaluated over each element, individually, using the element-based coordinates (ξ, η). Only the four elements around node m make nonzero contributions to the integrals in (36) and (37).

For a uniform grid, with hx = hy = h, the following is eventually obtained:


Equation (38) represents a system of linear equations which can be solved by, for example, successive overrelaxation.

As a matter of interest, the finite difference representation of (32), using central difference approximations, would be


which is a five point formula, compared with (38) which is a nine point formula.

Spectral method

The spectral method uses the same form of approximate solution as the traditional Galerkin method. Like the latter, the trial functions Wm and weighting functions Wm are nonzero throughout the computational domain. On the other hand, the unknown coefficients a; in (23) cannot be identified with nodal unknowns.

In the spectral method, the trial and weighting functions are orthogonal:

Fourier Series, Legendre Polynomials and Chebyshev Polynomials are examples of orthogonal functions often used in spectral methods.


Whereas conduction and convection heat transfer are described by partial differential equations, the basic governing equations for radiation are integral equations. The integral nature of the equations arises from the fact that it is necessary to sum the radiant energy exchanges over surfaces, the radiant intensity at which depends upon the relative orientation to, and separation from other surfaces. When combined modes of heat transfer occur, integral and differential equations are both relevant.

Under idealized conditions—uniform surface temperatures, constant properties, etc.—radiation calculations may be complex because of a particular geometry, but are not difficult and can be undertaken using classical methods. Under more realistic conditions, the challenge is somewhat greater.

The complexities which are present in heat transfer problems were mentioned at the start of this entry. For radiation, in particular, they are compounded by a geometry which is generally difficult and by the fact that the material properties relevant to radiation are dependent on the properties of the radiation itself (especially wave length) as well as on the nature and surface characteristics of the materials.

The integral nature of the equations has demanded different approaches to their solution. As outlined by Edwards (1985) and Chung (1988), these include the methods of discrete ordinates, Monte Carlo, invariant imbedding, zone analysis and finite element imbedding.

In the method of discrete ordinates, a set of fixed directions is chosen to approximate the solid angle integrals over all necessary directions by a sum over the chosen discrete set of directions.

In the Monte Carlo method, on the other hand, a sequence of random directions from any starting point is chosen and the paths and histories of "photons" taking those initial directions are followed. For a given photon, its scattering, transmission, reflection and/or absorption are determined based on the geometry of the problem and the properties of the materials involved. A statistical analysis of the behavior of tens or hundreds of thousands of photons leads to an accurate result. The method can be used for complete radiation calculations, or for the more limited task of determining radiation view factors in complex geometries which can then be used in conjunction with more conventional methods of calculation.

Invariant or differential imbedding has been used in astrophysics calculations by Chandraskhar (1960). The method is based on the principle of invariance of the emergent radiation from a semi-infinite plane-parallel layer or slab to the addition or subtraction of layers of arbitrary optical thickness. Discrete ordinates are used to convert integrals to summations. The method is ideal for layers or semi-infinite half-spaces, but has difficulties in more complicated geometries.

In zone analysis [Ozisik (1973); Siegel and Howell (1981)], the surfaces of an enclosure are subdivided into zones over which the temperatures, radiation properties and incident and emitted radiant intensities are assumed to be uniform. The accuracy of the results depends on the fineness of the subdivisions and the validity of the assumptions of uniformity.

For specularly transmitting slabs, and for reflecting or diffusely transmitting and reflecting slabs whose bidirectional properties are known from, for instance, differential imbedding, a multislab system may be built up by finite element imbedding. Interference effects may be included as well as polarization and directional effects. Finite element methods may also be used to find radiation view factors between surfaces, and to solve problems of combined conduction-convection-radiation heat transfer.


Anderson, J. D. (1994) Computational Fluid Dynamics: The Basics with Applications, McGraw-Hill, New York.

Chandraskhar, S. (i960) Radiative Heat Transfer, Dover, New York.

Chung, T. J. (1988) Integral and integro–differential systems, Chapter 14 of Handbook of Numerical Heat Transfer (Eds. W. J. Minkowycz, E. M. Sparrow, G. E, Schneider and R. H. Fletcher), Wiley, New York.

de Vahl Davis, G. (1986) Numerical Methods in Engineering and Science. Allen & Unwin, London.

Edwards, D. K. (1985) Gas radiation transfer, Part 6 of Chapter 14 of Handbook of Heat Transfer Fundamentals (Eds. W. M. Rohsenow, J. P. Hartnett and E. N. Ganic), 2nd edn., McGraw-Hill, New York.

Fletcher, C. A. J. (1991) Computational Techniques for Fluid Dynamics, Vols. I and II, 2nd edn., Springer-Verlag, Berlin, Heidelberg.

Hirasaki, G. S. and Hellums, J. D. (1968) A general formulation of the boundary conditions on the vector potential in three dimensional hydrodynamics, Quart. Appl. Math., 26, 331-342.

Hirasaki, G. J. and Hellums, J. D. (1970) Boundary conditions on the vector and scalar potentials in viscous three-dimensional hydrodynamics. Quart. Appl. Math., 2H, 293-296.

Leonard, B. P. (1988) Elliptic systems: finite-difference method IV, Chapter 9 of Handbook of Numerical Heat Transfer (Eds. W. J. Minkowycz, E. M. Sparrow, G. E. Schneider and R. H. Fletcher), Wiley. New York.

Ozisik, M. N. (1973) Radiative Transfer. Wiley-Interscience, New York.

Patankar, S. V. (1980) Numerical Heat Transfer and Fluid Flow, Hemisphere/McGraw-Hill, New York.

Patankar, S. V. and Spalding, D. B. (1972) A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer, 15, 1787–1806. DOI: 10.1016/0017-9310(72)90054-3

Reddy, J. N, and Gartling, D. K. (1994) The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press, Boca Raton, FL.

Richtmeyer, R. D. and Morton, K. W. (1967) Difference methods for initial-value problems, Interscience, New York.

Rodi, W. (1993) Turbulence Models and Their Application in Hydraulics—A State of the Art Review, 3rd edn., A. A. Balkema, Rotterdam/Brookfield.

Siegel, R. and Howell, J. R. (1981) Thermal Radiation Heat Transfer. 2nd edn., Hemisphere, Washington.

Torrance, K. K. (1985) Numerical methods in heat transfer, Chapter 5 of Handbook of Heat Transfer Fundamentals (Eds. W. M. Rohsenow, J. P. Hartnett and E. N. Ganic), McGraw-Hill, New York.

Number of views: 33344 Article added: 2 February 2011 Article last modified: 4 February 2011 © Copyright 2010-2019 Back to top