Versteeg H., Malalasekra W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method

Подождите немного. Документ загружается.

13.5 ILLUSTRATIVE EXAMPLES 437

The radiative source term for the enthalpy equation may be calculated from

h,rad

=∇. q

(r) = 4

πκ

−

(r, s) (13.45)

Further details of the method and its derivation can be found in Modest

(2003), Fiveland (1982, 1988, 1991), Fiveland and Jessee (1994), Jamaluddin

and Smith (1988) and Hyde and Truelove (1977). The standard discrete

ordinates method is suitable for Cartesian and axisymmetric geometries,

but it is not directly applicable to non-orthogonal and unstructured grids.

However, a modiﬁed version of the discrete ordinates method that is suitable

for complex geometries has been demonstrated by Charette et al. (1997) and

Sakami et al. (1996, 1998).

13.4.5 The finite volume method

Raithby and co-workers (Raithby and Chui, 1990; Chui et al., 1992, 1993;

Chui and Raithby, 1993) have presented a control volume integration method

for the calculation of radiative heat transfer that shares several features with

the discrete ordinates method. Equations for intensity are solved for a set of

discrete directions, which span the total solid angle 4

. In addition to control

volume integration the ﬁnite volume method uses a control angle integration

as well (see Chai et al., 1994a, b). Further improvements to the original

method along with appropriate interpolation techniques for non-orthogonal

mesh systems are available for radiative heat transfer calculations in complex

geometries. The difﬁculty that arises in the application of this class of

methods to non-orthogonal cells is the evaluation of cell intensities and

handling of control angle overhanging. These are both caused by misalign-

ment of cell surfaces with discretised control angles. Several techniques have

been developed to overcome this problem: details can be found in Chai and

Modar (1996), Baek and Kim (1998), Baek et al. (1998) and Murthy and

Mathur (1998). An illustrative program that demonstrate the ﬁnite volume

method is available in the appendices of Modest (2003).

First we consider the one-dimensional radiative heating of two cold black

plates by a layer of emitting/absorbing hot gas at temperature T

and with

absorption coefﬁcient

. The problem considered is shown in Figure 13.10.

An analytical solution is available for this very simple problem: see e.g.

Shah (1979), Siegel and Howell (2002) or Modest (2003). We apply the MC,

DTM and discrete ordinates methods to calculate the heat ﬂux to the plate

surfaces and compare the solutions obtained with the analytical result. The

non-dimensional wall heat ﬂux is given by

= 1 − 2E

(

L) (13.46)

∑

i =1

Illustrative examples13.5

Example 13.1

One-

dimensional

emitting/

absorbing

problem with

prescribed

temperature

ANIN_C13.qxd 29/12/2006 04:49PM Page 437

438 CHAPTER 13 CALCULATION OF RADIATIVE HEAT TRANSFER

Table 13.2 MC statistics

Optical Total number Standard error

length of bundles S

CPU time

0.1 6.80×10

8.50×10

0.2 1.30×10

2.10×10

0.5 3.40×10

4.60×10

320

1.0 6.80×10

3.70×10

564

1.5 1.00×10

3.30×10

800

2.0 1.30×10

4.80×10

874

2.5 1.70×10

3.60×10

1047

3.0 2.00×10

3.30×10

1430

Discrete ordinates method (DOM)

The DOM solution to this simple problem can be obtained by applying

equations (13.37) and (13.38). The solution procedure and the S

and S

discrete ordinates solutions are available in Modest (2003).

Figure 13.10 Conﬁguration for

one-dimensional temperature

prescribed problem

where E

(

L) is an exponential integral, values of which can be found in

tabulated form in Siegel and Howell (2002), and L is the distance between the

two plates. In this expression the wall heat ﬂux has been non-dimensionalised

by means of the black-body emissive power of the gaseous medium. We have

computed solutions for optical path lengths

L = 0.1, 0.2, ..., 3.0; this

ranges from ‘almost transparent (optically thin)’ to ‘optically dense (optically

thick)’.

DTM and MC methods

Although the problem is essentially one-dimensional, the DTM and MC

methods solve the directional intensity distribution by ray tracing, which

must be applied in a three-dimensional fashion. Simple trigonometric

expressions are sufﬁcient to determine all ray intersections in this problem.

For the DTM the solid angle was discretised using 10 × 10 for the N

and

the N

directions; hence, 100 rays were used per surface in the calculation.

For the MC solution the number of bundles ﬁred was increased from 6 × 10

to 2 × 10

as the optical path length of the medium increased from

L = 0.1

to 3.0. The standard error, S

, was of the order of 10

−4

for each MC simula-

tion performed. Table 13.2 shows the level of error and CPU time used for

MC solutions. CPU time depends on the computer system, but all DTM

solutions used less than 1 CPU second, which highlights the cost of MC

solutions.

ANIN_C13.qxd 29/12/2006 04:49PM Page 438

13.5 ILLUSTRATIVE EXAMPLES 439

Using the S

approximation the radiative heat ﬂux to walls in non-

dimensional form is given by

= w′

{exp(−

) − exp[−(

−τ

]} (13.47)

where

y and

Here

is the absorption coefﬁcient, y is the co-ordinate direction perpendicular

to plates, and w′

are appropriately summed weights for this one-dimensional

problem. For example, for the non-symmetric S

approximation

= 0.5, the

associated weight is

/2 = 1.5707, and four

-values are involved, giving the

summed weights w′

= 4 ×

/2 = 2

Using the non-symmetric S

approximation the non-dimensional heat

ﬂux is given by

= {exp(−

/0.5) −exp[−(

−

)/0.5]} (13.48)

With the S

discrete ordinates approximation the non-dimensional heat ﬂux

is given by

= 0.3945012{exp(−

/0.2958759) − exp[−(

−

)/0.2958759]

+ 0.6054088[exp(−

/0.9082483) − exp[−(

−

)/0.9082483]} (13.49)

The lower plate solution ( y = 0) using the S

approximation is as follows:

= [1 − exp(−

L/0.5)] (13.50)

and with S

= 0.3945012[1 − exp(−

L/0.2958759)]

+ 0.6054088[1 − exp(−

L/0.9082483)] (13.51)

Results

The predictions of the wall heat ﬂux for all three solution methods are

compared with the analytical solution in Table 13.3 and Figure 13.11. It

can be seen that all solutions except the DOM S

appear to give good

agreement. The MC solution is the most accurate of all. The S

discrete

ordinates approximation is too simplistic to produce close agreement, but

it does predict the correct trend. In terms of calculation cost the resources

used by the DOM are negligible compared with the CPU times for the MC

method.

− T

)

− T

)

− T

)

− T

)

N/2

∑

i =1

− T

)

ANIN_C13.qxd 29/12/2006 04:49PM Page 439

Next we consider the much more challenging, but realistic problem of the

three-dimensional IFRF furnace geometry, which has been widely used

in the radiative heat transfer literature to evaluate calculation methods. The

conﬁguration is outlined in detail in the paper by Jamaluddin and Smith

(1988). In this illustration we apply the MC, DTM and S

DOM methods

for the calculation of radiative ﬂux to the ﬂoor and roof of the furnace,

described in the IFRF-M3 Trial (Flame 10). The furnace dimensions are

6.0 m × 2.0 m × 2.0 m in the x-, y- and z-directions and the measured tem-

peratures by Hyde and Truelove (1977) are used to deﬁne the radiation

problem. In the absence of an analytical solution, the MC results along with

the standard error estimate are used as the benchmark solution. The discrete

transfer calculation has been performed with 400 rays per surface element,

while S

quadrature is used for the discrete ordinates solution.

440 CHAPTER 13 CALCULATION OF RADIATIVE HEAT TRANSFER

Figure 13.11 Comparision of

numerical results for one-

dimensional temperature

prescribed problem

Table 13.3 Comparison of results from different solution methods for the one-

dimensional temperature prescribed problem

Optical Exact

thickness solution

MC DTM DOM S

DOM S

0.1 0.16742 0.167702 0.16923 0.18127 0.17626

0.2 0.29779 0.296686 0.29610 0.32968 0.31349

0.5 0.55732 0.557363 0.55680 0.63212 0.57800

1.0 0.78076 0.780648 0.78062 0.86467 0.78516

1.5 0.88551 0.886071 0.88653 0.95021 0.88134

2.0 0.93970 0.939207 0.93974 0.98168 0.93251

2.5 0.96737 0.966947 0.96740 0.99326 0.96122

3.0 0.98211 0.981435 0.98214 0.99752 0.97763

Example 13.2

Radiation in

a three-

dimensional

furnace

geometry

ANIN_C13.qxd 29/12/2006 04:49PM Page 440

13.5 ILLUSTRATIVE EXAMPLES 441

Figure 13.12 shows the calculated incident radiative ﬂux on the ﬂoor and the

roof of the furnace, and it is clear that the differences between the 400-ray

discrete transfer, the S

discrete ordinates and the MC solutions are very

small. A closer look at the actual numerical values (not presented here)

reveals that the discrete ordinate solution has the highest deviation. This

could be attributed to its spatial differencing practice and the accuracy of the

quadrature set. Malalasekera et al. (2002) have also shown that lower-order

discrete ordinates solutions do not provide good agreement with MC results.

This example illustrates that careful selection of numerical discretisation

parameters for the DTM and DOM can yield solutions that are very close to

MC solutions. It is worth noting that MC solutions require a vast amount of

computer resources (more than 20 times larger compared with DTM and

DOM) and therefore are not suitable for practical calculations involving CFD.

Figure 13.12 Calculated

incident radiative ﬂux on the

ﬂoor and the roof of the furnace.

DTM (400) is the discrete

transfer method solution with

400 rays per point, MC (163M)

is the Monte Carlo solution with

163 million bundles, and DOM

(S16) is the discrete ordinates

solution

In this example we present a further comparison of methods in the L-

shaped geometry shown in Figure 13.13. This more complex geometry has

also been used in the literature to evaluate different radiation calculation

techniques (Malalasekera and James, 1995, 1996; Henson and Malalasekera,

1997a; Sakami et al, 1998; Hsu and Tan, 1997). The problem speciﬁcation is:

(i) all walls are black at a temperature of 500 K and (ii) the enclosure is ﬁlled

with an emitting/absorbing and non-scattering medium at a temperature of

1000 K. The L-shape presents additional complexity for solution methods

because in a CFD computation the geometry cannot be discretised using

Cartesian co-ordinates without unnecessarily wasting storage by blocking

off a large number of cells. For the application of the DTM the geometry

has been modelled using the non-orthogonal surface mesh shown in Fig-

ure 13.13. The standard discrete ordinates method is not applicable when a

non-orthogonal mesh is employed. Here we have used a special version of

the discrete ordinates method by Sakami et al (1998). Figure 13.14 shows the

predicted distribution of net radiative heat ﬂux along line B–B (Figure 13.13)

for different values of the absorption coefﬁcient

= 0.5, 1.0, 2.0, 5.0 and 10.0

using the MC, DTM (16 × 16 rays per hemisphere), DOM method (with the

Example 13.3

Radiation in

three-

dimensional

complex

geometry

Solution

ANIN_C13.qxd 29/12/2006 04:49PM Page 441

approximation), and also a general-purpose integral equation technique

known as the YIX method (Hsu et al, 1993). The discrete transfer solution

is from Malalasekera and James (1995), the MC results are by Maltby (1996)

and the YIX solution by Hsu and Tan (1997). The MC calculation has

been performed using a mesh of 204 volumes using symmetry and releasing

1 024 000 bundles per volume or per surface. It can be seen that results of all

methods agree very well. The example is a good illustration to show the ﬂexi-

bility of the DTM and MC methods, which can both be applied without any

modiﬁcations to complex meshes. The DOM requires a special formulation.

When adequate resources are used all methods produce equally good results

that are in agreement with the MC solution. Further applications of the

DTM and MC method in complex geometry applications can be found in

Henson (1998), Henson and Malalasekera (1997a), Malalasekera and James

(1995, 1996) and Hsu and Tan (1997).

In many practical situations where radiative heat transfer is important,

e.g. furnaces and combustors, the participating medium consists of a gaseous

mixture of combustion products. Its radiative properties, such as the absorp-

tion coefﬁcient, required for radiation calculations have to be evaluated by

442 CHAPTER 13 CALCULATION OF RADIATIVE HEAT TRANSFER

Figure 13.13 L-shaped

geometry

Figure 13.14 Comparison

of numerical results for the

L-shaped geometry

Calculation of

radiative

properties in

gaseous mixtures

13.6

ANIN_C13.qxd 29/12/2006 04:49PM Page 442

13.7 SUMMARY 443

considering the composition, temperature and pressure of the mixture.

Accurate calculations have to take into account the spectral characteristics of

radiative properties but, for practical purposes, total emissivity values and

scattering coefﬁcients can be used. There are a number of well-established

techniques of estimating total properties, including spectral band models,

correlations for gas absorptivities and emissivities, k-distribution methods

and weighted-sum-of-grey-gases model. The most widely used approximate

radiative property calculation method for combustion problems is the weighted-

sum-of-grey-gases (WSGG) model (see above references). Another popular

model in the radiation literature is the RADCAL program of Groshandler

(1993), which takes into account detailed spectral properties. Detailed dis-

cussions of the techniques can be found in Siegel and Howell (2002),

Brewster (1992), Denison and Webb (1993, 1995) and Modest (2003).

When radiation is involved in combustion simulations the overall success

of predictions depends on many modelling aspects: ﬂow and the turbulence

model, the combustion model, degree of detailed chemistry included, the

radiative property calculation methods and the radiation model. All these

are coupled and, therefore, given the available resources, the best possible

submodelling practice should be used to achieve a good overall result. More

advanced property calculations and ﬁner radiation calculations will give the

beneﬁt of higher accuracy and better resolution in radiation source terms.

The level of detail required in the application and affordable cost of a calcu-

lation is the major factor which determines the choice. It is well known that

temperatures are overpredicted if radiation is neglected. As minor species

such as NO are very sensitive to temperature, CFD combustion predictions

need to include radiation in order to achieve accurate predictions of these

pollutant species. In order to save computing resources it is sometimes pos-

sible to carry out the radiation calculations on a much coarser mesh than the

one used for the ﬂuid ﬂow and combustion calculations. The source terms

for the ﬂuid ﬂow calculation can be appropriately interpolated from the

coarse mesh calculation. This practice is quite adequate in many situations,

saving considerable computational overheads.

Most commercial CFD packages provide the user with the option of one or

more of the radiation models described in this chapter. Their main features

can be summarised as follows:

• Monte Carlo method: the most general and versatile of all models.

It requires very large amounts of computing resources and is, therefore,

unsuitable for general-purpose CFD calculations. The probabilistic

basis of the method enables its users to generate error estimates. The

MC method is useful in situations where solutions with quantiﬁable

uncertainty are required, e.g. for benchmarking and validation of

other models.

• Discrete transfer method: an economical general-purpose algorithm

based on ray tracing. The DTM has been successfully applied to a wide

range of combusting ﬂows. It is suitable for all types of structured and

unstructured meshes. The DTM is limited to isotropic scattering, and

is non-conservative (non-zero sum of the radiative heat ﬂuxes incident

upon the bounding surfaces of an enclosure); Coelho and Carvalho

Summary13.7

ANIN_C13.qxd 29/12/2006 04:49PM Page 443

(1997) have proposed a conservative DTM formulation. Uncertainty

estimates are available for a limited number of simple cases.

• Discrete ordinates method: a method based on the numerical quadrature

of angular direction integration and control volume integration for

path integration. The original method was limited to Cartesian and

axisymmetric grids, but the Sakami et al. (1996, 1998) DOM enables

applications involving unstructured meshes and more complex

geometry. The method in general is efﬁcient and applicable to scattering

and non-scattering problems. Higher-order discrete ordinates are

required to achieve accuracy and one of the drawbacks of the method

is that it is non-conservative and can suffer from what are known as ‘ray

effects’ mainly due to the discretisation practice used in the formulation.

• Finite volume method: a method based on integration over the control

angle for direction integration and the control volume method for path

integration developed by Raithby and co-workers. The main advantage

of this method is that it is fully conservative and applicable to arbitrary

geometries. One of the disadvantages is that the method can incur errors

due to what is known as ‘solid angle overhanging’.

In this chapter we have presented the key aspects of the more popular

methods, which should clarify the role played by the numerical parameters

associated with each model (ray number, order of discrete ordinates, grey gas

models etc.). The interested reader should follow up the relevant literature

to understand the ﬁner details of these methods and also to research the

many other methods that we were forced to omit for the sake of brevity.

444 CHAPTER 13 CALCULATION OF RADIATIVE HEAT TRANSFER

ANIN_C13.qxd 29/12/2006 04:49PM Page 444

In the derivations of Chapter 4 a linear proﬁle assumption was used to cal-

culate the gradients (

∂φ

∂

x), (

∂φ

∂

y) etc. at the faces of the control volume.

For simple diffusion problems the practice was shown to give reasonably

accurate results even for coarse grids. In Example 4.3 we also observed that

by reﬁning the grid the accuracy of the solution can be improved. Grid

reﬁnement is the main tool at the disposal of the CFD user for the improve-

ment of the accuracy of a simulation. The user would typically perform a

simulation on a coarse mesh ﬁrst to get an impression of the overall features

of the solution. Subsequently the grid is reﬁned in stages until no (signiﬁc-

ant) differences of results occur between successive grid reﬁnement stages.

Results are then called ‘grid independent’. Here we brieﬂy demonstrate the

theoretical basis of this method of accuracy improvement and compare the

order of discretisation schemes as a measure their efﬁcacy.

Consider the equally spaced one-dimensional grid (spacing ∆x) shown in

Figure A.1.

Appendix A Accuracy of a flow simulation

Figure A.1

For a function

(x) the Taylor series development of

(x +∆x) around

the point i at x is

(x +∆x) =

(x) +

∆x +

+ . . . (A.1)

In our notation we use discrete values

and

for

(x) and

(x +∆x)

respectively so that equation (A.1) can be written as

∆x +

+ . . . (A.2)

∆x

∂

∂φ

∂

∆x

∂

∂φ

∂

ANIN_Z01.qxd 29/12/2006 10:06 AM Page 445

446 APPENDIX A

This may be rearranged to give

∆x =

−

− . . . (A.3)

=−

− ...

=+Truncated terms (A.4)

By neglecting the truncated terms which involve the multiplying factor ∆x

we may write

≈ (A.5)

The error involved in the approximation (A.5) is due to neglect of the trun-

cated terms. Formula (A.3) suggests that the truncation error can be reduced

by decreasing ∆x. In general the truncated terms of a ﬁnite difference scheme

contain factors ∆x

. The power n of ∆x governs the rate at which the error

tends to zero as the grid is reﬁned and is called the order of the difference

approximation. Therefore equation (A.5) is said to be ﬁrst-order in ∆x and

we write

=+O(∆x) (A.6)

Since it uses values at points E and P (where x

> x

) to evaluate the gradi-

ent (

∂φ

∂

x) at P, formula (A.6) is called a forward difference formula with

respect to point P.

Similarly we may derive a backward difference formula for (

∂φ

∂

x) at P

from

(x −∆x) =

(x) −

∆x +

+ . . . (A.7)

After some algebra we ﬁnd the backward difference formula for (

∂φ

∂

x) at P:

=+O(∆x) (A.8)

Equations (A.7) and (A.8) are both ﬁrst-order accurate. The backward and

forward difference formulae described here involve values of

at two points

only.

By subtracting equation (A.7) from (A.1) we get

(x +∆x) −

(x −∆x) = 2

∆x +

+ . . . (A.9)

∆x

∂

∂φ

∂

−

∆x

∂φ

∂

∆x

∂

∂φ

∂

−

∆x

∂φ

∂

−

∆x

∂φ

∂

−

∆x

∂φ

∂

∆x

∂

−

∆x

∂φ

∂

∆x

∂

∂φ

∂

ANIN_Z01.qxd 29/12/2006 10:06 AM Page 446