A new radiative transfer scattering phase function discretisation approach with inherent energy conservation

International Journal of Heat and Mass Transfer 73 (2014) 789–803

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A new radiative transfer scattering phase function discretisation approach with inherent energy conservation Thomas H. Roos a,⇑, Thomas M. Harms b,1 a b

CSIR, PO Box 395, Pretoria 0001, South Africa Department of Mechanical and Mechatronic Engineering, University of Stellenbosch, Private Bag X1, Matieland 7602, South Africa

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 27 November 2013 Received in revised form 20 January 2014 Accepted 22 January 2014 Available online 19 March 2014 Keywords: Discrete Ordinates Phase function Discretisation Anisotropic scattering Packed bed

In the popular Discrete Ordinates Method (DOM) formulation of the Equation of Radiative Transfer (ERT), the 4p solid angle range of directions is divided into a finite number of discrete directions or ordinates. This requires that the continuous distribution of the scattering phase function of the medium under consideration must be discretised to suit the different number, weightings and directions of the SN ordinate set being used. This must be done such that the sum of scattered energy is conserved relative to the continuous distribution, and that the asymmetry factor g is similarly conserved. This paper introduces a discretisation technique with inherent energy conservation, suitable for any quadrature scheme. The technique was tested on two large sphere scattering phase function distributions of interest for packed bed radiative heat transfer: the analytic distribution for a diffusely reflecting sphere (a backscattering test case) and the distribution for a transparent sphere (n = 1.5) obtained by ray tracing (a test case with strong forward scatter and some back-scatter). In both cases the resultant discretised phase function distributions for the S4, S6 and S8 ordinate sets produced errors for the sum of scattered energy conservation of less than 0.035% and errors for g less than 1.3%. This demonstrates the inherent energy conservation of the method, as well as visible reductions in g errors. The phase function values for each case are tabulated in the paper. The major benefit of the method is the fact that computationally costly matrix calculations are avoided at run-time: the discretisation for a given scattering medium using a quadrature scheme of given order is performed only once beforehand, and the resultant distributions can be stored in an input file or look-up table for future computations with different boundary conditions, different meshes and even different geometries. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction 1.1. Fundamentals of the discrete ordinate method The solution of the Equation of Radiative Transfer (ERT) is required whenever a scattering, absorbing and/or emitting medium is present in a volume of radiative heat transfer interest:

dI r ¼ ^s  rI ¼ jIb  bI þ ds 4p

Z

Ið^sj ÞUð^sj ; ^sÞdXj

ð1Þ

4p

where Eq. (1) is valid for grey media (or non-grey media on a spectral basis) [1]. The final term on the right hand side of Eq. (1)

⇑ Corresponding author. Tel.: +27 82 33 278 33; fax: +27 12 349 1156. 1

E-mail addresses: [email protected] (T.H. Roos), [email protected] (T.M. Harms). Tel.: +27 21 808 4376; fax: +27 21 808 4958.

http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.01.052 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved.

is the in-scattering term. The scattering phase function Uð^sj ; ^sÞ represents the probability that a photon travelling in direction ^sj will be scattered in direction ^s (the direction represented on the left hand side by dI=ds ¼ ^s  rI). An alternative representation of Uð^sj ; ^sÞ is U(h), where h is the angle between the ^s direction of the beam being calculated in Eq. (1) and incoming direction ^sj of radiation being considered in the in-scatter term. The scattering phase function distribution U(h) is continuous with angle h, and is axisymmetric about the axis of the incoming beam. Of the different techniques available for simulating radiative heat transfer in participating media, the SN or Discrete Ordinates Method (DOM) is popular and widely used [2]. The development and applications of the DOM have been extensively described by many authors, so will not be repeated here. In the DOM, the continuous 4p solid angle range is divided into a finite number n of discrete directions or ordinates, so in Eq. (1), numerical quadrature (with quadrature weights wj for each direction ^sj ) replaces the integral over direction in the in-scattering term [3,4]:

790

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

Nomenclature DOM ERT f() g HG I k LA MDOM n N Nintervals Nw k ^s wi x, y, z xk

discrete ordinates method Equation of Radiative Transfer function of () asymmetry factor Henyey–Greenstein radiative intensity, (W m2 sr1) current interval number linear-anisotropic modified discrete ordinates method number of directions ^si in current SN quadrature set order of SN quadrature set number of intervals number of different weightings of ordinates in interval number k unit direction vector discrete direction weight for direction ^si directions in Cartesian coordinate system angle of lower boundary of kth interval

Greek symbols b extinction coefficient (j + r) (m1) g direction cosine with respect to the z-axis

1 4p

Z

f ð^sÞdX ¼

4p

n 1 X wj f ð^sj Þ 4p j¼1

polar angle between i and j directions absorption coefficient (m1) direction cosine with respect to the x-axis direction cosine with respect to the y-axis scattering coefficient (m1) azimuthal angle continuous scattering phase function discretised scattering phase function normalised scattering phase function step-wise discontinuous axisymmetric scattering phase function solid angle, sin h dh du

h

j l n

r u U Uij e ij U Ustep dX

Subscripts b blackbody cells pertaining to grid cells HG Henyey–Greenstein i direction of current beam being calculated j direction of incoming beam being considered in in-scattering term wall boundary surfaces of the computational volume x, y, z x, y or z directions

ð2Þ Si ¼ jIb þ

In the standard DOM, this leads to a system of n equations described (in Cartesian coordinates) by:

s:rIi ¼ ni

n @Ii @Ii @Ii rX þ gi þ li ¼ jIb  bIi þ wj Uij Ij ; @x @y @z 4p j¼1

i ¼ 1; 2; . . . ; n

Iðx; y; z; ^sÞ ¼ Id ðx; y; z; ^sÞ þ Is ðx; y; z; ^sÞ

ð4Þ

The direct (or ‘‘ballistic’’ [5]) intensity component Id originates from sources on the boundary surfaces of the computational volume, which may be either collimated or diffuse [7]. It is assumed in the transport equation that there is no in-scattering or emission contribution to the direct intensity component, so Eq. (3) for the direct component is simplified:

n

@Id @Id @Id þg þl þ bId ¼ 0 @x @y @z

ð5Þ

which therefore allows an analytical solution [6]:

Id ðx; y; z; ^sÞ ¼ Iwall ðxwall ; ywall ; zwall ; ^sÞebsðx;y;z;^sÞ

ð6Þ

where Iwall ðxwall ; ywall ; zwall ; ^sÞ is the value of intensity at the wall (collimated or diffuse) in the direction of interest ð^sÞ. The remaining term on the right hand side of Eq. (4) is the diffuse component Is, which originates by scattering and emission from the participating medium within the computational volume, and (as in the conventional DOM) is solved numerically:

ni

@Isi @Is @Is þ gi i þ li i ¼ ðj þ rmi ÞIsi þ Si ; @x @y @z

i ¼ 1; 2; . . . ; n

ð7aÞ

B

C ð7bÞ

j¼1 j–i



The modified DOM or MDOM is an augmented version of the standard DOM described above. It was developed to address the known ‘‘ray effect’’ limitation of the conventional DOM, by splitting the intensity I in Eq. (3) into two components [5,6]:

1

n C Xd rB C BX wj Uij Isj þ I Ui wall C B C 4p B wall A @

rmi ¼ r 1  ð3Þ

0

1 wi Uii 4p

 ð7cÞ

In the above formulation [5], a mechanism is included to remove the forward-scattering peak from the in-scattering term so that it is treated as transmission, to reduce the number of iterations before convergence due to scattering. 1.2. Scattering phase function discretisation In order to be used in a DOM or MDOM analysis, the phase function U(h) must be discretised into discrete scattering values Uij (where i represents the current ordinate direction being calculated and j represents the incoming beam). One concern with the discretisation of the continuous phase function U(h) into discrete scattering values Uij is that the total R scattered energy fraction (1=4p 4p Uð^sj ; ^sÞdXj ) and asymmetry facR tor (g = 1/4p 4pU(h)cos h dX) might not be conserved after numerical quadrature replaces the integrals, i.e.:

1 4p 1 4p

Z

Uð^sj ; ^sÞdXj ¼ 1–

4p

Z

n 1 X wj Uij ; 4p j¼1

UðhÞ cos hdXj ¼ g–

4p

i ¼ 1; 2; . . . ; n

n 1 X wj Uij cos hij ; 4p j¼1

i ¼ 1; 2; . . . ; n

ð8Þ

ð9Þ

For isotropic scattering (U = 1), this is not a problem, as the directions ^si and weights wi in the SN approximation are chosen to satisfy conservation for the ‘‘zeroeth’’ moment:

1 4p

Z 4p

dX ¼ 1 ¼

n 1 X wj 4p j¼1

ð10Þ

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

Unfortunately, the scattering phase function of a participating medium is typically not isotropic, and in some cases is strongly anisotropic. Linear-isentropic (LA) scattering phase function distributions, which take the form [8]:

U ¼ 1 þ A cos h ¼ 1 þ 3g cos h

ð11Þ

have been shown after discretisation to conserve both scattered energy and g for three orders of quadrature (S4, S8, and S12) and for all values of g less than or equal to 1/3 [5]: values of g greater than 1/3 lead to unphysical negative values of U for cos h values near 1 in Eq. (11). Non-linear scattering phase function distributions, however, are not generally conservative after discretisation. Using the Henyey– Greenstein (HG) phase-function approximation:

UHG ðhÞ ¼

1  g2

ð12Þ

1:5

½1 þ g 2  2g cosðhÞ

which gives a sharp forward scattering peak, deviations for calculated values of scattered energy (from unity) and for calculated values of g (from the prescribed value of g) were shown to reach and exceed 1% at quite low values of g [5]. Table 1 shows data from [5] where values in boldface are quoted in the text, while the remaining values had to be digitised from graphs: values with decimal fractions from linear plots and integer values from logarithmic plots. As the value of prescribed g increases, it can be seen that conservation is better achieved using a higher order SN quadrature set. The requirement for energy conservation has been known for some time, and one response has been to ‘‘normalise’’ discrete values obtained from a continuous distribution. One example of normalisation is the method of Kim and Lee [9]:

e ij ¼ Uij  U

n 1 X wj Uij 4p j¼1

!1

Another example is that of Wiscombe [10] who makes use of a correction factor for both the i and j directions

e ij ¼ ð1 þ ai þ aj ÞUij U

ð14Þ

which when combined with Eq. (8) results in n 1 X wi ð1 þ ai þ aj ÞUi;j ¼ 1 4p i¼1

ð15Þ

The set of correction factors is then obtained by the matrix solution of a set of equations. Farhat et al. [11] compared the methods of Kim and Lee [9] and of Wiscombe [10], as well as a least squares method, for a range of forward and backscattering Mie-phase function distributions in two-dimensional analyses. Energy conservation was excellent, with average errors of about 1014 reported, but no reference was made to asymmetry factor. Boulet et al. [12], however, found

Table 1 Error values of scattered energy and asymmetry factor for HG phase-function for 3 quadrature sets at different values of g [5]. g

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75

in a three-dimensional radiative transfer analysis that the quality of the results using the methods of Kim and Lee [9] and of Wiscombe [10] was strongly dependent on the value of asymmetry factor g: (1) Good results were be obtained for a low value of asymmetry factor (g = 0.2). (2) At a high value of asymmetry factor (g = 0.8), the discretised phase function was deformed, leading to a reduction in effective scattering and overprediction of radiative flux. (3) For anisotropy that is acute (g = 0.93), the discretised phase function was so strongly deformed that the results predicted were equal to those of a nonparticipating medium. This lack of simultaneous conservation of scattered energy and asymmetry factor led Hunter and Guo [2] to develop a matrixbased normalisation technique

e ij ¼ ð1 þ Aij ÞUij U

ð16Þ

where the normalisation factors Aij need to be found such that the e ij values result in the inequality signs normalised phase function U in Eqs. (8) and (9) being replaced with equality signs: n 1 X e ij ¼ 1; wj U 4p j¼1

i ¼ 1; 2; . . . ; n

n 1 X wj Uij cos hij ¼ g; 4p j¼1

Scattered energy error (%)

Asymmetry factor error (%)

S4

S8

S12

S4

S8

S12

0.97 1.77 3.4 6.2 10.9 19.1 33 62

0.6 1.07 2.3 4.6 8.7 15.8 28.1 54

<0.5 <0.5 <0.5 <0.5 <0.5 <0.5 0.95 2.74

3.29 4.9 7.5 11 18 29 49 86

<0.5 <0.5 <0.5 0.4 1.0 2.3 5.1 12

<0.5 <0.5 <0.5 <0.5 <0.5 <0.5 0.61 2

i ¼ 1; 2; . . . ; n

ð17Þ

ð18Þ

and in the symmetry condition being satisfied:

e ij ¼ / e ji / ð13Þ

791

ð19Þ

For this system of equations, there are infinitely many possible solutions. The solution chosen is that where the matrix has the minimum norm. This approach was applied to the problem which gave rise to the results in Table 1 [5], where conservation of both scattered energy and asymmetry factor was achieved in each case. A specific recommendation was that in order to accurately conserve both scattered energy and asymmetry factor in HG phase functions, normalisation should be applied for g P 0.20, 0.40, and 0.60 for the S4, S8, and S12 quadrature sets respectively. In this paper, a new discretisation approach with inherent conservation of scattered energy will be developed. Conservation of asymmetry factor will be investigated for spherical particle scattering phase functions for packed bed radiative transfer. 2. SN discretisation SN (level-symmetric) is a popular quadrature system used in the DOM, where the order N (an even integer) of the SN approximation determines the number of directions using the relationship n = N(N + 2). The SN discretisation approach leads to a triangular arrangement of directions per octant, as illustrated for the S4 and S6 cases in Fig. 1. The direction vectors are determined by connecting the origin and the points where planes perpendicular to the three Cartesian axes intersect the unit sphere (creating in effect ‘‘latitude’’ lines) and each other. The quadrature sets are designed to be symmetric relative to the three axes. Each increase by 2 in order N of the SN quadrature set adds another ‘‘latitude’’ plane perpendicular to each of the three axes, increasing the number of intersection points and therefore directions. The sum of weights of all the directions per SN quadrature set totals 4p, which is the solid angle of the sphere. For the S4 set, the weights wj for the three directions per octant (and therefore 3  8 = 24 directions over all 8 octants) have equal

792

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

Fig. 1. Principal octant of level symmetric quadrature set: S4 (left) and S6 (right).

value (wj = 4p/24 ffi 0.523599, j = 1, 2, . . . , 24). The incrementing of the order N from 4 to 6 (resulting in the S6 set) doubles the number of directions, but adds new weighting level: the three directions closest to the principal axes have equal weightings (wlevel1 ffi 0.160952), which differ in value from the three directions furthest from the principal axes (wlevel2 ffi 0.362647). Consequently, the S8 quadrature set has 3 weighting levels, the S10 set has 4 weighting levels, and so on. The SN is not the only family of quadrature sets used with DOM. Fig. 2 shows examples of other quadrature sets, namely the P16T16, P16-EW, T6 and SRAP7 sets [13]. 3. Scattering phase function axisymmetric interval discretisation approach 3.1. Axisymmetric step-wise discontinuous Ustep(h) concept The discrete scattering values Uij are used in the final term on the right hand side of Eq. (3) for the standard DOM, and in the diffuse 1st in-scattering term within the brackets of Eq. (7b) of the MDOM. In both cases, the number of angles of interest for Uij is small, being limited to n (as there are n incoming ordinate directions j for each current ordinate direction i) and is independent of the number of the computational cells in the computational domain. This is not true for the direct scattering values Uiwall in the 2nd in-scattering term of Eq. (7b) of the MDOM. For each cell volume centroid in the computational domain, the number of inbound ‘‘ballistic’’ rays is determined by (and equal to) the number of boundary cells. For a Nx  Ny  Nz Cartesian spatial grid, the number of inbound rays per computational cell is given by:

  Ninbound rays ¼ 2  ðNx  Ny Þ þ ðNy  Nz Þ þ ðNx  Nz Þ

ð20Þ

While the number of inbound rays is equal for each computational cell, the angle an inbound ray from the same boundary element makes with the current ordinate direction i will be different for separate computational cells, as the cell centroids coordinates are all unique. This means the number of angles is given by:

Nangles ¼ Ninbound rays  Ncomputational cells  n ¼ 2nNx Ny Nz ðNx Ny þ Ny Nz þ Nx Nz Þ

ð21Þ

For a cubic domain where Nx cells = Ny cells = Nz cells = Nside, Eq. (21) gives 6Nside5 as the number of possible angles, placing a considerable computational burden on the normalisation matrix solver. The approach followed in this paper can be described as follows: (1) The actual scattering phase function U(h) of the medium is transformed into a new axisymmetric distribution Ustep(h) which is stepwise-discontinuous with h. Because the outscattering U(h) is symmetrical about the incoming ordinate direction j, the same in-scattering phase function U(h) is symmetrical about the current ordinate direction i: U(hij) = U(hji) so Ustep(hij) = Ustep(hji). Therefore, separate unique Ustep(h) distributions are only required for each weighting in the quadrature set being applied (i.e.: one distribution is needed for the S4 quadrature set, two for S6, three for S8, etc.) and each U(h) distribution under consideration. The major benefit of the method is the fact that computationally costly matrix calculations are avoided at run-time: Ustep(h) can be pre-calculated, and may subsequently be read from a file or look-up table in the DOM software programme for subsequent calculations with different boundary conditions, different meshes and even different geometries.

Fig. 2. Principal octant directional distributions for different quadrature sets: from left to right P16-T16, P16-EW, T6 and SRAP7 [13].

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

(2) For in-scattering in the standard DOM and for the diffuse inscattering term of the MDOM, each Uij value is determined from the appropriate Ustep(h) distribution purely using the value of cos(hij) as an argument. The Ustep(h) distribution is chosen on the basis of the medium scattering distribution, the quadrature set and the weighting of the current ordinate direction i. (3) For the direct in-scattering term of the MDOM, each Uiwall value required is similarly obtained from the appropriate Ustep(h) distribution, using the value of cos(hiwall) as an argument. 3.2. Axisymmetric interval calculation procedure

(the ‘‘north pole’’ in Fig. 3 right). Firstly, all the ordinate directions j = 1, 2, . . . , n are ranked in terms of the angle hj they make with the current ordinate direction i, sorted from largest (p radians or 180°) to smallest (0). We start with hj = p: there will be only one ordinate direction j that makes an angle of p with the current ordinate direction i. We denote this ordinate direction j = 1. The 1st angular interval k = 1 extending from h = p to h = x1 is determined only by the weighting w1 of ordinate direction j = 1, which from symmetry is equal to the weighting wi of the current ordinate direction i (w ffi 0.362646957):

w1 ¼

Z p x1

The discretisation approach followed will now be described. The directional computation domain of the DOM and MDOM (from a three-dimensional trigonometric perspective) can be described as a unit sphere with the surface markings of an irregular polyhedron with polygon-shaped sides. Fig. 3 (left) shows the markings for an S6 quadrature set. Each black dot represents where an ordinate direction j passing through the centre of the sphere penetrates the sphere surface. Each side of each polygon on the sphere surface is a line (a portion of a great circle on the sphere surface) separating two adjacent dots. Where the adjacent dots represent ordinate directions of equal weighting, the line separating them is equidistant from either dot. Where the weightings of adjacent dots are unequal, the line is closer to the dot representing the smaller weighting; the net result is that the surface area of each polygon-shaped area of the sphere (each representing a solid angle) is equal to the weighting wj of the corresponding ordinate direction j. As the sum of the solid angles must equal 4p, so must the sum of weightings (wj). Even though Fig. 3 (left) illustrates this for the S6 quadrature of the level-symmetric approach, the foregoing discussion is independent of quadrature scheme. Since the in-scattering phase function is rotationally symmetric around the axis of the current ordinate direction i (as discussed in Section 3.1), this sphere with irregular polyhedron markings needs to be transformed into a rotationally symmetric representation. This is done by replacing the polygon lines with a series of ‘‘latitude’’ lines as surface markings, where the ‘‘north pole’’ and ‘‘south pole’’ are the dots of the current ordinate direction i where it enters and leaves the sphere (see Fig. 3 right). The surface area between adjacent ‘‘latitude’’ lines must equal the sum of the weightings of the ordinate directions whose black dots lie between the lines. The process or algorithm followed to determine the positions of the ‘‘latitude’’ lines will now be described by means of an example: we choose the S6 quadrature set with an ordinate direction with large weighting (w ffi 0.362647) as current ordinate direction i

793

2p sin h dh ¼ 2p½ cos hpx1 ; w

1

) x1 ¼ cos1

2p

  1 ¼ cos1 ð0:942282944Þ ¼ 160:44

ð22Þ

This is illustrated in Fig. 4 (right): the ‘‘latitude line’’ closest to the ‘‘pole’’ contains only one black dot (the current ordinate direction i), and the area of the sphere surface enclosed by this latitude line (shaded dark orange in Fig. 4, right) is equal to the weighting w1. For the second (k = 2) angular interval from x1 to x2, it is first necessary to determine the number n of ordinate directions j where the included angle hj with the current ordinate direction i is equal to the next largest included angle after h1 = p, such that h2 < h1 = p. This process is illustrated in Fig. 4 (left). It can be seen that there are 7 ordinate directions j adjacent to the current ordinate direction i, but only one is closest (ordinate direction with large weighting w ffi 0.362647 indicated by brown arrow), therefore n = 1. This then leads to:

n  wj¼2 ¼

Z

x1

x2

2p sin hdh ¼ 2p½ cos hxx12 ;

n  w  j¼2 þ cos x1 2p   1 1  0:362647  0:942282944 ¼ cos 2p

) x2 ¼ cos1

¼ cos1 ð0:884565888Þ ¼ 152:20

ð23Þ

In certain cases, two errors can ensue. Firstly, it is possible that angle hj of ordinate direction j lies outside the calculated angular range between x1 and x2. This is rectified by extending the angular range by including subsequent ordinate directions j and recalculating x2 until hj lies inside the angular range between x1 and x2:

ðnl¼1  wj1 Þ þ ðnl¼2  wj2 Þ ¼

Z

x1 x2

2p sin h dh ¼ 2p½ cos hxx21

Fig. 3. Ordinate weighting representations of S6 quadrature set: conventional (left) and proposed rotationally symmetric representation (right).

794

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

z

x

y

Fig. 4. S6 quadrature set with current direction i (an ordinate of large weighting): the 7 adjacent ordinate directions where arrows of the same colour denote ordinate directions equidistant from i (left); and 1st 3 latitude lines proposed of rotationally symmetric discretisation (right).

  ðnl¼1  wj1 Þ þ ðnl¼2  wj2 Þ þ cos x1 2p

) x2 ¼ cos1

ð24Þ

The second possible error is that a subsequent ordinate direction j may be associated with an angle hj that may lie in the previous angular range, in this case between x1 and x2. The remedy is analogous to that mentioned above, ordinate direction j is included in the previous range and x2 is recalculated. Returning to our test case, we check whether ordinate j = 2 lies between with x1 and x2. The included angle between j = 2 and i is h2 ¼ 158:81 , which does indeed lie between x1 and x2. This second range is indicated by the dark yellow band in Fig. 4 (right) containing the single black dot. Moving to j = 3, we find that two ordinate directions with small weightings (w = 0.160952. . .), indicated by the red arrows in Fig. 4(left), are equidistant from the current ordinate direction i, so n = 2.

n  wj¼3 ¼

Z

x2

x3

Since h4 ¼ 139:92 , both h3 and h4 lie between 152.20° and 135.88°, so an error of the first kind is avoided. Unfortunately, an error of the second kind occurs in the next interval. The next closest ordinates to i are j = 5, two further ordinate directions with small weightings and equidistant from the current ordinate direction i, indicated by the green arrows in Fig. 4(left), so nl=3 = 2. h5 ¼ 137:62 , which lies in the previous interval between with x2 and x3 and so constitutes an error of the second kind. The x3 boundary must again be moved to accommodate the j = 5 ordinates. Continuing in the philosophy of Eq. (24):

ðnl¼1  w3 Þ þ ðnl¼2  w4 Þ þ ðnl¼3  w5 Þ ¼ ¼ 2p½ cos hxx23 ) x3 ¼ cos1

ðnl¼1  w3 Þ þ ðnl¼2  w4 Þ ¼

x2

  ð4  0:160951818Þ þ ð2  0:362646957Þ þ cos x2 2p ð27Þ

ð25Þ

The included angle between j = 3 and i is h3 ¼ 146:38 , which does not lie between with x2 and x3 and so constitutes an error of the first kind. The next closest ordinate direction with its weighting needs to be included, to move the x3 boundary. Next closest are j = 4, two ordinate directions with large weightings (w = 0.362647. . .), equidistant from the current ordinate direction i, indicated by the blue arrows in Fig. 4(left). This is a second group (l = 2) of ordinates and weightings n in this interval, so nl=2 = 2. Repeating the calculation as per the first line of Eq. (24): x1

2p sin h dh

¼ cos1 ð0:666666667Þ ¼ 131:81

  ðn  wj¼3 Þ ) x3 ¼ cos1 þ cos x2 2p   1 ð2  0:160952Þ  0:884565888 ¼ cos 2p

Z

x2

x3

2p sin h dh ¼ 2p½ cos hxx23 ;

¼ cos1 ð0:8333333Þ ¼ 146:44

Z

2p sin h dh ¼ 2p½ cos hxx12

The values of h3, h4 and h5, representing six ordinate directions, now all lie between x2 ¼ 152:20 and x3 ¼ 131:81 . This third range is indicated by the light yellow band in Fig. 4 (right) containing six black dots. This process is then continued until all the ordinate directions with their weightings are allocated into band intervals. Rearranging Eq. (24) to represent all intervals:

XNinterv als XNw k k¼1

j¼1

nj  wj

 k

¼

XNinterv als Z k¼1

xkþ1

! 2p sin h dh

xk

Z p n X wj ¼ 4p ¼ 2p sin h dh j¼1

ð28Þ

0

where Nw k represents the number of different ordinate weightings in the boundary interval k. 3.3. Axisymmetric interval results for SN ordinate sets

  ðn1  w3 Þ þ ðn2  w4 Þ ) x3 ¼ cos1 þ cos x2 2p   ð2  0:160952Þ  ð2  0:362647Þ ¼ cos1  0:884565888 2p ¼ cos1 ð0:719899222Þ ¼ 135:88 ð26Þ

Angular discretisations were subsequently performed for several SN ordinate sets. Fig. 5 shows that 10 intervals could be found for the S4 ordinate set. For the S6 ordinate set, Fig. 6 shows that 20 intervals were found from an ordinate direction of small weighting (w = 0.160952. . .), and 14 intervals from an ordinate direction of large weighting (w = 0.362647. . .). Fig. 7 shows the results for the S8 ordinate set: the numbers of intervals found from ordinate

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

795

Fig. 5. Angular positions of directional ordinates and resultant angular interval boundaries for S4 quadrature set.

Fig. 6. Angular positions of directional ordinates and resultant angular interval boundaries for S6 quadrature set from ordinate of small weighting (top) and large weighting (bottom).

directions of small, medium and large weightings were 18, 28 and 12 respectively. The interval boundaries for all weightings are given as values of cos h in Table 2 for S4 and S6 and in Table 3 for S8. The precision of the input data (direction cosines and their weights) used from [14] is 15 significant digits. This allows the values in Tables 2 and 3 to be published to 15 significant digits as well for two reasons: firstly, by comparing symmetrical (positive and negative) direction pairs the resultant accuracy of the discretisation boundaries can be seen to have degraded to 13 significant digits; and secondly, to allow validation of software written using the technique described in this paper. 3.4. Scattering phase function discretisation test case 1: large diffusely reflecting sphere The analytical scattering phase function for a large diffusely reflecting sphere is given in [15] as:

UðhÞ ¼

8 ðsin h  h cos hÞ 3p

ð29Þ

and as can be seen in Fig. 8, has zero forward scattering leading to a backscattering peak. The analytical solution of the scattered energy sum and the asymmetry factor are therefore given by:

1 4p

Z 4p

Z p 1 UðhÞ2p sin h dh 4p 0 Z p 2p 8 ðsin h  h cos hÞ sin h dh ¼ 4p 0 3p  p 4 1 ð4h  3 sinð2hÞ þ 2h cosð2hÞÞ ¼ 3p 8 0  4 3p ¼1 ¼ 3p 4

Ug ð^si ; ^sÞdXi ¼

Z 1 g ¼ cos h ¼ UðhÞ cos hdX 4p 4p Z p 1 ¼ UðhÞ2p sin h cos h dh 4p 0 Z p 2p 8 ðsin h  h cos hÞ sin h cos h dh ¼ 4p 0 3p  p 4 1 ¼ ð4 sinð3hÞ þ 9h cos h þ 3h cosð3hÞÞ 3p 36 0 4 h pi 4 ¼  ¼  ¼ 0:4444444 3p 9 3

ð30Þ

ð31Þ

The above two sets of calculations were then repeated using numerical integration of the scattering phase function using distributions

796

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

Fig. 7. Angular positions of directional ordinates and resultant angular interval boundaries for S8 quadrature set from ordinate of small weighting (top), medium weighting (centre) and large weighting (bottom).

of 2°, 1° and 0.5° resolution respectively. This was done to determine the effect of scattering phase function angular resolution on the results, as the analytical solutions were known. The results Table 2 Boundaries of discretisation intervals as values of cos h for S4 and S6. Ordinate set

S4

Weighting Boundary No

–

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

S6 Small

are shown in Table 4. It can be seen from the sums of scattered energy that an accuracy of three and four significant figures can be expected for g for scattering phase function angular resolutions of 2° and 0.5° respectively. Constant values of scattering phase function valid for each angular interval for each of the ordinate weightings of the S4, S6 and S8 sets were calculated. The method followed for each interval i was:

R xkþ1

Large

Uk ¼ 1.0 11/12 3/4 1/2 1/3 0 1/3 1/2 3/4 11/12 1.0

1.0 0.974383722541931 0.923151167625792 0.782100778417195 0.666666666666668 0.551232554916141 0.435798443165614 0.326848832374212 0.275616277458073 0.160182165707546 0.000000000000004 0.160182165707537 0.275616277458064 0.326848832374203 0.435798443165606 0.551232554916133 0.666666666666660 0.782100778417187 0.923151167625783 0.974383722541922 1.0

1.0 0.942282944124736 0.884565888249473 0.666666666666668 0.307717055875266 0.192282944124740 0.141050389208601 0.000000000000004 0.141050389208592 0.192282944124731 0.307717055875258 0.666666666666660 0.884565888249465 0.942282944124728 1.0

xk

UðhÞ2psin h dh PNw k j¼1 nj  wj

ð32Þ

In the above, the numerator represents the numerical integral from the lower (k) to upper (k + 1) boundary, while the denominator is the sum (from 1 to Nw k, where Nw k represents the number of different ordinate weightings in interval k) of the products of the number of ordinates of each weighting present in the interval and the value of respective weightings. Calculations were then performed to determine the effective values of asymmetry factor for each discretised scattering phase function distribution:

g¼

 1 X Nintervals X Nw k j¼1 nj  wj  cos hj  Uk k¼1 4p

ð33Þ

The results are shown in Table 5. In Table 5, the values of scattered energy error and g are shown to four significant figures, the highest resolution of the numerical integration determined from Table 4. The errors for scattered

797

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803 Table 3 Boundaries of discretisation intervals as values of cos h for S8. Ordinate set

S8

Weighting Boundary No

Small

Medium

Large

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

1.0 0.984207302327561 0.793045738464625 0.531585395344876 0.456954261535373 0.398115825398309 0.323484691588806 0.202621906982682 0.031585395344874 0.000000000000003 0.031585395344880 0.202621906982688 0.323484691588812 0.398115825398315 0.456954261535379 0.531585395344882 0.793045738464631 0.984207302327567 1.0

1.0 0.972746959207813 0.890987836831252 0.859402441486375 0.827817046141497 0.722746959207812 0.691161563862935 0.627990773173180 0.481021389995564 0.449435994650687 0.417850599305810 0.344365907717002 0.312780512372125 0.203768349203377 0.000000000000003 0.2037683492033830 0.3127805123721310 0.344365907717008 0.417850599305816 0.449435994650693 0.481021389995571 0.627990773173187 0.691161563862941 0.722746959207818 0.827817046141504 0.859402441486381 0.890987836831258 0.972746959207819 1.0

1.0 0.926515308411192 0.831759122376560 0.655243813965367 0.491725569212245 0.094756186034629 0.000000000000003 0.094756186034635 0.491725569212251 0.655243813965373 0.831759122376566 0.926515308411198 1.0

Fig. 8. Large diffusely reflecting sphere scattering phase function: comparison of analytical and discretised distributions for S4 quadrature set.

Table 4 Effect of angular resolution of large diffusely reflecting sphere scattering phase function on calculated sum of scattered energy and g. Property

Scattering phase function angular resolution

Sum scattered energy Scatt. energy error % g Calc. g Error %

1.000324 0.0324 0.444760 0.0710

2°

1° 1.000076 0.0076 0.444519 0.0168

0.5° 1.000019 0.0019 0.444461 0.0036

energy and g were determined relative to the angular resolutiondependent numerical integration values from Table 4, not the ideal values from Eq. (30) and (31). As can be seen, absolute values of scattered energy error and g error are everywhere less than 0.03% and 1.3% respectively. The resultant distributions for each angular interval for each of the ordinate weightings of the S4, S6 and S8 sets are shown in Figs. 8–10 respectively, and are listed in the Appendix in Tables A1–A6 for the 0.5° resolution. The precision chosen of 6 digits after the decimal point arises from a comparison of the results for the 1° and 0.5° resolutions.

798

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

Table 5 Large diffusely reflecting sphere values of g calculated from different discretisations of the scattering phase function distribution: comparison of non-discretised and S4, S6 and S8 discretisations. Quad. set

Weight

Scattering phase function angular resolution 2°

S4 S6 S8

– S L S M L

1°

0.5°

Scatt. energy error%

g Calc.

g Error%

Scatt. energy error%

g Calc.

g Error%

Scatt. energy error%

g Calc.

g Error%

0.0160 0.0110 0.0196 0.0252 0.0192 0.0134

0.4397 0.4431 0.4403 0.4421 0.4432 0.4391

1.142 0.365 1.004 0.599 0.353 1.272

0.0010 0.0083 0.0031 0.0041 0.0093 0.0016

0.4394 0.4429 0.4403 0.4417 0.4428 0.4390

1.147 0.360 0.949 0.635 0.396 1.248

0.0009 0.0025 0.0013 0.0029 0.0018 0.0020

0.4393 0.4430 0.4403 0.4416 0.4427 0.4389

1.151 0.330 0.936 0.649 0.400 1.249

Fig. 9. Large diffusely reflecting sphere scattering phase function: comparison of analytical and discretised distributions for S6 quadrature set.

Fig. 10. Large diffusely reflecting sphere scattering phase function: comparison of analytical and discretised distributions for S8 quadrature set.

In order to prove the inherent energy conservation of the method, Eq. (32) is rearranged: Z xkþ1  X 1 X Nintervals  X Nw k Uk kj¼1 nj  wj ¼ UðhÞ2p sin h dh; ) Uk j¼1 nj  wj k¼1 4p ki ¼

)

1 X Nintervals k¼1 4p



Z

xkþ1

!

The apparent discrepancy between the small (less than 0.03%) but finite scattered energy errors in Table 5 and implied scattered errors in Eq. (34) of zero, is due to the use of numerical integration of a discrete scattering phase function distribution in the former and analytical integration of an algebraic phase function distribution in the latter.

UðhÞ2p sin h dh

xk

XN 1 X Nintervals wk Uk j¼1 nj  wj k¼1 4p



Z p 1 UðhÞ2p sin h dh ¼ 1 ¼ 4p 0 ð34Þ

3.5. Scattering phase function discretisation test case 2: large transparent sphere with n = 1.5 The second test case is for spherical glass beads with a refractive index (n) of 1.5. Two ray tracing analyses were available. The

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

799

Fig. 11. Large transparent sphere scattering phase function: comparison of distributions of [16,17] in linear (left) and log (right) scales.

Fig. 12. Large transparent sphere scattering phase function: comparison of analytical and discretised distributions for S4 quadrature set, as well as HG distribution for g = 0.65.

Fig. 13. Large transparent sphere scattering phase function: comparison of analytical and discretised distributions for S6 quadrature set.

800

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

Fig. 14. Large transparent sphere scattering phase function: comparison of analytical and discretised distributions for S8 quadrature set.

Table 6 Large transparent sphere error values of g and scattered energy for S4, S6 and S8 discretisations compared with non-discretised scattering phase function distribution. Quadrature set

Weighting

Scattered energy error (%)

g Calc.

Error (%)

Nondiscretised S4 S6

–

–

0.6600

–

– Small Large Small Medium Large

0.0265 0.0348 0.0114 0.0096 0.0066 0.0058

0.6572 0.6549 0.6526 0.6516 0.6579 0.6549

0.421 0.775 1.118 1.273 0.308 0.776

S8

first analysis [16] had to be digitised from a log polar plot obtained from a pdf of the original paper. The second analysis [17] had been performed using the Zemax software, and the resultant distribution was digitally available in 1° resolution. Fig. 11 shows a comparison of the two distributions in both linear (left) and logarithmic (right) scales. A strong forward-scattering peak can be seen, with the characteristic ‘‘fishtail’’ backscattering pattern of dielectric spheres: this therefore provides a good test case for simultaneous forward and back-scattering peaks. Good agreement can be found between the two distributions, except the distribution from [17] has additional features, most visible on the logarithmic scale:

yielded a scattered energy sum of 1.0. This was achieved at h ¼ 0:000495 . Table 4 suggests that a precision of four significant figures is appropriate for a smoothly-varying discrete distribution with 1° resolution. An asymmetry factor value of 0.65997803 was obtained, which when corrected to four significant figures is 0.6600. As with the large diffusely reflecting sphere, constant values of scattering phase function valid for each angular interval for each of the ordinate weightings of the S4, S6 and S8 sets were calculated, with the resultant distributions shown in Figs. 12–14 respectively. As before, calculations were then performed to determine the effective values of asymmetry factor for each discretised scattering phase function distribution. The results are shown in Table 6, again to four significant figures. As with the diffusely scattering sphere, the results are very good: scattered energy error is everywhere below 0.035% and g error does not exceed 1.3% for any instance. A qualitative comparison can be made with the HG distribution with g = 0.65 in Table 1 (shown plotted in Fig. 12 for comparison), for which values of scattered energy error and g error are available. For the S4 quadrature set, the scattered energy error and g error are 0.0265% and 0.421% respectively, compared with values of 19.1% and 29% in Table 1. For the S8 quadrature set, the corresponding maximum absolute values are 0.0096% and 1.273% respectively, compared with 15.8% and 2.3% in Table 1. 4. Conclusions

(1) A strong forward-scattering peak of 52540 for h ¼ 0 , dropping sharply to 9.234 for h ¼ 1 . (2) A back-scattering peak of 17.187 for h ¼ 180 , dropping sharply to 4.605 for h ¼ 179 . (3) Two additional small peaks between cos h values of 0.0 and 0.4. The good 1° resolution of the [17] distribution resulted in it being used to develop the discretised distributions: the fact that the two halves (left and right) of the [16] distribution did not lie on top of each other gave rise to accuracy uncertainty. The same process was followed as that for the large diffusely reflecting sphere. No exact analytical phase function distribution was available as an absolute reference, so only numerical integration could be used for this purpose. The scattered energy sum was found to greatly exceed unity at a value of 1.99996, due to the strong forward peak. Rather than truncating this peak, an additional point was added between h ¼ 0 and h ¼ 1 with a scattering phase function value of 9.234 (the same value as for h ¼ 1 ). The position of the point was varied until the numerical integration

A scattering phase function discretisation approach with inherent scattered energy conservation for the Discrete Ordinates Method has been developed. The method has been tested on two large sphere scattering phase function distributions of interest for packed bed radiative transfer, that of a diffusely reflecting sphere (a back-scattering test case) and for a transparent sphere with n = 1.5 (a test case with strong forward scatter and some back-scatter). In both cases the resultant discretised phase function distributions for the S4, S6 and S8 ordinate sets produced errors for the sum of scattered energy conservation of less than 0.035% and errors for g less than 1.3%. This demonstrates the inherent energy conservation of the method, as well as visible reductions in g errors. While the method was only tested for three orders of the SN quadrature scheme in this paper, it should be applicable to any quadrature scheme of any order with similar results. The major benefit of the method is the fact that computationally costly matrix calculations are avoided at run-time: the discretisation for a given scattering medium using a quadrature scheme of given order is performed only once beforehand, and the resultant distributions

801

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

can be stored in an input file or look-up table for future computations with different boundary conditions, different meshes and even different geometries.

forming the ray tracing of the transparent sphere and extracting the scattering phase function.

Acknowledgement

Appendix A. Discretised scattering phase function values for S4, S6 and S8 ordinate sets

The work described in this paper has been supported by the CSIR. The authors would like to thank Derek Griffith [17] for per-

Tables A1–A6.

Table A1 Discretised scattering phase function values for diffusely reflecting and transparent spheres for S4 set. Interval No

1 2 3 4 5 6 7 8 9 10

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere

1.0 11/12 3/4 1/2 1/3 0 1/3 1/2 3/4 11/12

11/12 3/4 1/2 1/3 0 1/3 1/2 3/4 11/12 1.0

2.563828 2.279010 1.856973 1.480196 1.086893 0.642379 0.369262 0.190243 0.056702 0.007762

1.335706 0.050466 0.041177 0.050028 0.073457 0.223179 0.670872 1.642374 3.861111 7.168509

Table A2 Discretised scattering phase function values for diffusely reflecting and transparent spheres for small weighting ordinate of S6 set. Interval No

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere

1.0 0.974384 0.923151 0.782101 0.666667 0.551233 0.435798 0.326849 0.275616 0.160182 0.0 0.160182 0.275616 0.326849 0.435798 0.551233 0.666667 0.782101 0.923151 0.974384

0.974384 0.923151 0.782101 0.666667 0.551233 0.435798 0.326849 0.275616 0.160182 0.0 0.160182 0.275616 0.326849 0.435798 0.551233 0.666667 0.782101 0.923151 0.974384 1.0

2.633978 2.539562 2.320763 2.049925 1.824238 1.612975 1.420246 1.289362 1.159951 0.959377 0.745546 0.579138 0.486107 0.403309 0.297077 0.200428 0.118191 0.046603 0.009606 0.001415

1.195622 1.317234 0.140063 0.041123 0.040881 0.043810 0.054801 0.049767 0.044985 0.101530 0.143327 0.251196 0.381257 0.564631 0.926297 1.489509 2.376509 4.186911 6.683891 8.500012

Table A3 Discretised scattering phase function values for diffusely reflecting and transparent spheres for large weighting ordinate of S6 set. Interval No

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere

1.0 0.942283 0.884566 0.666667 0.307717 0.192283 0.141050 0.0 0.141050 0.192283 0.307717 0.666667 0.884566 0.942283

0.942283 0.884566 0.666667 0.307717 0.192283 0.141050 0.0 0.141050 0.192283 0.307717 0.666667 0.884566 0.942283 1.0

2.594970 2.456856 2.157029 1.606530 1.209320 1.083232 0.945847 0.757441 0.638242 0.542644 0.307384 0.088740 0.020771 0.004483

1.121029 0.836749 0.044126 0.047070 0.044426 0.046014 0.109000 0.138558 0.193911 0.295675 0.970470 3.040492 5.522901 7.695449

802

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

Table A4 Discretised scattering phase function values for diffusely reflecting and transparent spheres for small weighting ordinate of S8 set. Interval No

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere

1.0 0.984207 0.793046 0.531585 0.456954 0.398116 0.323485 0.202622 0.031585 0.0 0.031585 0.202622 0.323485 0.398116 0.456954 0.531585 0.793046 0.984207

0.984207 0.793046 0.531585 0.456954 0.398116 0.323485 0.202622 0.031585 0.0 0.031585 0.202622 0.323485 0.398116 0.456954 0.531585 0.793046 0.984207 1.0

2.648932 2.402539 1.929315 1.614458 1.497796 1.385799 1.229546 1.011832 0.870521 0.827417 0.699559 0.528265 0.424047 0.357791 0.295930 0.163087 0.032845 0.000637

1.335592 0.504412 0.041196 0.043452 0.047416 0.059425 0.040430 0.078673 0.147234 0.127581 0.162060 0.317019 0.510883 0.694242 0.926450 1.920191 5.120885 8.795724

Table A5 Discretised scattering phase function values for diffusely reflecting and transparent spheres for medium weighting ordinate of S8 set. Interval No

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere

1.0 0.972747 0.890988 0.859402 0.827817 0.722747 0.691162 0.627991 0.481021 0.449436 0.417851 0.344366 0.312781 0.203768 0.0 0.203768 0.312781 0.344366 0.417851 0.449436 0.481021 0.627991 0.691162 0.722747 0.827817 0.859402 0.890988 0.972747

0.972747 0.890988 0.859402 0.827817 0.722747 0.691162 0.627991 0.481021 0.449436 0.417851 0.344366 0.312781 0.203768 0.0 0.203768 0.312781 0.344366 0.417851 0.449436 0.481021 0.627991 0.691162 0.722747 0.827817 0.859402 0.890988 0.972747 1.0

2.632463 2.499889 2.370600 2.299715 2.154349 2.014173 1.920871 1.723282 1.562829 1.508280 1.419682 1.333338 1.222102 0.990678 0.718747 0.533363 0.456871 0.403221 0.351705 0.322149 0.244637 0.162039 0.128960 0.086883 0.049981 0.035633 0.014949 0.001540

1.184433 0.983845 0.060785 0.045716 0.041926 0.040411 0.040704 0.042010 0.044480 0.046931 0.053112 0.068648 0.044350 0.089575 0.157085 0.308870 0.438058 0.562829 0.708831 0.813406 1.199882 1.810569 2.192509 2.936149 3.928800 4.554192 6.130847 8.464889

Table A6 Discretised scattering phase function values for diffusely reflecting and transparent spheres for large weighting ordinate of S8 set. Interval No

1 2 3 4 5 6 7 8 9 10 11 12

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere

1.0 0.926515 0.831759 0.655244 0.491726 0.094756 0.0 0.094756 0.491726 0.655244 0.831759 0.926515

0.926515 0.831759 0.655244 0.491726 0.094756 0.0 0.094756 0.491726 0.655244 0.831759 0.926515 1.0

2.575693 2.379028 2.089685 1.758470 1.282392 0.913334 0.786862 0.500487 0.229247 0.107023 0.034512 0.006436

1.233622 0.266019 0.041415 0.041698 0.048549 0.139315 0.130212 0.418667 1.300703 2.611233 4.689081 7.361471

T.H. Roos, T.M. Harms / International Journal of Heat and Mass Transfer 73 (2014) 789–803

References [1] M.F. Modest, Radiative Heat Transfer, second ed., Academic Press, San Diego, 2003. [2] B. Hunter, Z. Guo, Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media, Int. J. Heat Mass Transfer 55 (2012) 1544–1552. [3] B.G. Carlson, K.D. Lathrop, Transport theory—the method of discrete-ordinates, in: H. Greenspan, C.N. Kelber, D. Okrent (Eds.), Computing Methods in Reactor Physics, Gordon & Breach, New York, 1968. [4] W.A. Fiveland, Discrete-ordinates methods for radiative heat transfer in isotropically and anisotropically scattering media, J. Heat Transfer 109 (1987) 809–812. [5] B. Hunter, Z. Guo, Phase-function normalisation in the 3-D discrete-ordinates solution of radiative transfer – Part I: Conservation of scattered energy and asymmetry factor, Numer. Heat Transfer B (2012) 203–222. [6] M.A. Ramankutty, A.L. Crosbie, Modified discrete-ordinates solution of radiative transfer in three-dimensional rectangular enclosures, J. Quant. Spectrosc. Radiat. Transfer 60 (1) (1998) 103–134. [7] C. Caliot, Numerical Methods in Radiative Transfer: Introduction to DOM, FVM and MCM, Promes – CNRS, 11 June 2010. [Online]. Available: sfera.sollab.eu/ downloads/Schools/Caliot_DOM_FVM_MCM.pdf. [Accessed 30 January 2012]. [8] J.S. Truelove, Three-dimensional radiation in absorbing–emitting–scattering media using the discrete-ordinates approximation, J. Quant. Spectrosc. Radiat. Transfer 39 (1988) 27–31.

803

[9] T.-K. Kim, H. Lee, Effect of anisotropic scattering on radiative heat transfer in two-dimensional rectangular enclosures, Int. J. Heat Mass Transfer 31 (1988) 1711–1721. [10] W.J. Wiscombe, On initialisation, error and flux conservation in the doubling method, J. Quant. Spectrosc. Radiat. Transfer 18 (1976) 637–658. [11] H. Farhat, N. Abdallah, M.N. Borjini, R. Méchi, R. Saïd, Comparison between Kim, Wiscombe and least squares scattering phase function renormalization methods, Int. J. Therm. Sci. 45 (2006) 1127–1139. [12] P. Boulet, A. Collina, J.L. Consalvi, On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media, J. Quant. Spectrosc. Radiat. Transfer 104 (2007) 460–473. [13] B. Hunter, Z. Guo, Comparison of quadrature schemes in DOM for anisotropic scattering radiative transfer analysis, Numer. Heat Transfer B 63 (2013) 485– 507. [14] G. Colomer, Numerical methods for radiative heat transfer (PhD Thesis), Universitat Politècnica de Catalunya, 2006. [15] J.R. Howell, R. Siegel, M.P. Mengüç, Thermal Radiation Heat Transfer, fifth ed., CRC Press, Boca Raton, 2011. [16] B. Variot, T. Menigault, G. Flamant, Modelling and optimization of a two-slab selective volumetric solar receiver, Sol. Energy 53 (4) (1994) 359–368. [17] D. Griffith, Private Communication, CSIR, Pretoria, 2013.