Conservation of scattered energy and asymmetry factor in the new Rotationally Symmetric Spherical Discretisation scheme

International Journal of Heat and Mass Transfer 101 (2016) 205–225

Contents lists available at ScienceDirect

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

Conservation of scattered energy and asymmetry factor in the new Rotationally Symmetric Spherical Discretisation scheme T.H. Roos a,⇑, T.M. Harms b, C.G. du Toit c a

CSIR, PO Box 395, Pretoria 0001, South Africa Department of Mechanical and Mechatronic Engineering, University of Stellenbosch, Private Bag X1, Matieland 7602, South Africa c Department of Mechanical and Nuclear Engineering, North-West University, Private Bag X6001, Potchefstroom 2520, South Africa b

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 20 November 2015 Received in revised form 2 March 2016 Accepted 2 March 2016 Available online 28 May 2016 Keywords: Discrete Ordinates Phase function Discretisation Anisotropic scattering Packed bed

A Rotationally Symmetric Spherical Discretisation (RSSD) technique has recently been developed for the discretisation of scattering phase functions in the Discrete Ordinates Method (DOM) in a consistent and rational manner. RSSD has inherent energy conservation, is suitable for any quadrature scheme used, minimising run-time matrix calculation, and to date has been used for the S4, S6 and S8 ordinate sets. In this follow-up paper, the RSSD boundary intervals for the four weightings of the S10 ordinate set are presented. The effect of asymmetry factor g and angular resolution of the discrete input scattering phase function distribution ‘‘grid” on final RSSD errors (in scattered energy conservation and calculated asymmetry factor) is explored, making use of the Henyey–Greenstein (HG) family of distributions. It is shown that the RSSD scattered energy error ne RSSD is a strong function of the grid resolution, indicating that ne RSSD can be controlled down to user-required levels by judicious choice of grid resolution and by local grid refinement in regions of steep gradients and peaks in the scattering phase function distribution (in this paper error levels as low as 0.002% have been achieved). RSSD asymmetry factor error ng RSSD behaves differently to ne RSSD , however, displaying a strong grid resolution dependence in a ‘‘grid sensitivity zone” of combinations of coarse grid resolution and high values of prescribed g (for the HG family of scattering phase functions, above g values of 0.9, 0.95 and 0.98 for grid angular resolutions of 1°, 0.5° and 0.25° respectively). Outside of the ‘‘grid sensitivity zone” g error is insensitive to grid resolution and achieves maximum (positive and negative) values of between
2.7% and 1% between prescribed g values of 0.7 and 0.9. This band decreases to a range from
1% to 1% for the S4, S6 and S10 ordinate sets. The grid resolution conclusions are supported by examining RSSD discretisation for two real scattering phase function distributions: for a diffusely reflecting large sphere with an asymmetry factor of
0.4444, grid angular resolutions of 2°, 1°, 0.5° and 0.25° resulted in the decreasing maximum ne RSSD values of 0.0449%, 0.0136%,
0.0050% and
0.0015% respectively, but near-constant maximum ng RSSD values of
1.225%,
1.239%,
1.247% and
1.250% respectively, while for a transparent large sphere with an asymmetry factor of 0.660, grid angular resolutions of 1°, 0.5° and composite 0.5°/0.25° resolution resulted in the decreasing maximum ne RSSD values of 0.0685%,
0.0402% and
0.0046% respectively and near-constant maximum ng RSSD values of
1.286%,
1.310% and
1.313% respectively. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Whenever a scattering, absorbing and/or emitting medium is present in a volume of radiative heat transfer interest, the Equation of Radiative Transfer (ERT) needs to be solved:

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

Z

4p

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

⇑ Corresponding author. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.03.005 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.

ð1Þ

Eq. (1) is valid for grey media, or non-grey media on a spectral basis [1]. The term generally presenting the most difficulty is the in-scattering term, the final term on the right hand side of Eq. (1). Here Uð^sj ; ^sÞ is the scattering phase function, which describes 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). The scattering phase function can also be represented as Uð hÞ, where h is the angle between the ^s direction of the beam being calculated in Eq. (1) and incoming direction ^sj of radi-

206

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

Nomenclature DOM ERT fðÞ g HG I k LA MDOM n N Ninterv als Nwk s ^s wi x, y, z xk

Discrete Ordinates Method Equation of Radiative Transfer function of () RU asymmetry factor g ¼ 1=4p 4p ðhÞ cos h dX Henyey–Greenstein radiative intensity (W m
2 sr
1) 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 distance unit direction vector discrete direction weighting for direction ^si directions in Cartesian coordinate system angle of lower boundary of kth interval

g

direction cosine with respect to the y-axis polar angle between i and j directions absorption coefficient (m
1) direction cosine with respect to the z-axis direction cosine with respect to the x-axis, error scattering coefficient (m
1) 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 ~ ij U Ustep dX

Subscripts b blackbody HG Henyey–Greenstein i direction of current beam being calculated j direction of incoming beam being considered in inscattering term

Greek symbols b extinction coefficient (j + r) (m
1)

ation being considered in the in-scatter term. The scattering phase function distribution Uð hÞ is axisymmetric about the axis of the incoming beam ^sj , and is continuous with angle h. The ERT can be solved analytically for only a small fraction of radiative heat transfer problems. For more general radiative heat transfer problems involving participating media, the Discrete Ordinates Method (DOM) is widely used [2]. In the DOM, the continuous 4p solid angle range is divided into a finite number n of discrete directions or ordinates, so numerical quadrature (with quadrature weights wj for each direction ^sj ) replaces the integral over direction in Eq. (1) in the in-scattering term [3,4]:

1 4p

Z 4p

f ð^sÞdX ¼

1 Xn wj f ð^sj Þ j¼1 4p

ð2Þ

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

@Ii @Ii @Ii r Xn s  rIi ¼ ni wU I; þ gi þ li ¼ jIb
bIi þ j¼1 j ij j @x @y @z 4p i ¼ 1; 2; . . . ; n

ð3Þ

While many quadrature ordinate sets exist for use in the DOM, such as the Legendre–Chebyshev (PN-TN), Legendre-equal weight (PN-EW), triangle-tessellation (TN) and spherical ring approximation (SRAPN) sets [5], the most well-known is the level-symmetric (SN) family of ordinate sets, where the order N (an even integer) in the subscript determines the number of ordinate directions using the relationship n ¼ NðN þ 2Þ. The modified DOM or MDOM is an augmented version of the standard DOM described above. It was developed to address a known limitation of the conventional DOM (the so-called ‘‘ray effect”), by splitting the intensity I in Eq. (3) into direct and scattered components, described in [6–8] and [9]:

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

ð4Þ

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).

A concern with the discretisation of the continuous phase function Uð hÞ into discrete scattering values Uij is that the total scatR tered energy fraction (1=ð4pÞ 4p Uð^sj ; ^sÞdXj ) and asymmetry R factor (g ¼ 1=ð4pÞ 4p Uð hÞ cos h dX) might not be conserved after numerical quadrature replaces the integrals [2], i.e.:

1 4p 1 4p

Z 4p

Z 4p

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

n 1 X wj Uij ; 4p j¼1

Uð hÞ cos h dXj ¼ g–

i ¼ 1; 2; . . . ; n

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

ð5Þ

i ¼ 1; 2; . . . ; n

ð6Þ

As discussed in [9] and reiterated here for convenience, conservation of total scattered energy fraction and asymmetry factor g occurs for two classes of scattering distributions:  Isotropic scattering (U ¼ 1)  Linear-anisotropic (LA) scattering [1,10] (U ¼ 1 þ A cos h ¼ 1 þ 3g cos h) for all values of g less than or equal to 1/3 [6] Values of g greater than 1/3 lead to unphysical negative values of U. Table 1 Error values of scattered energy and asymmetry factor for HG phase-function for 3 quadrature sets at different values of g.[6] Scattered energy error (%)

Asymmetry factor error (%)

g

S4

S8

S12

S4

S8

S12

0.4

0.97

0.6

<0.5

<0.5

<0.5

0.45

1.77 3.4 6.2 10.9 19.1 33 62 124 270 711 3161 20,609 335,015

1.07 2.3 4.6 8.7 15.8 28.1 54 107 235 620 2781 18,298 299,296

<0.5

3.29 4.9

<0.5

<0.5

<0.5 <0.5 <0.5 <0.5 0.95 2.74 11 32 110 610 4298 71,449

7.5 11 18 29 49 86 157 320 787 3318 21,032 337,759

<0.5 0.4 1.0 2.3 5.1 12 27 65 195 938 6244 100,768

<0.5 <0.5 <0.5 <0.5 0.61 2 5 18 70 414 2938 47,921

0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 0.98 0.995

207

Number of ordinates

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

5

S10 smallest

4

boundaries

3 2 1 0 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

Number of ordinates

cos θ 5

S10 2nd smallest boundaries

4 3 2 1 0 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

cos θ Fig. 1. Angular positions of directional ordinates and resultant angular interval boundaries for S10 quadrature set from ordinate of smallest weighting (top) and 2nd smallest weighting (bottom).

Unfortunately, the scattering phase function of a participating medium is typically neither isotropic nor linear-anisotropic. Such non-linear scattering phase function distributions are not generally conservative after simple discretization. Using the Henyey–Greenstein (HG) phase-function approximation:

UHG ðhÞ ¼

1
g2

ð7Þ

1:5

½1 þ g 2
2g cosðhÞ

Number of ordinates

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) reached and exceeded 1% at quite low values of g [6] (see Table 1). The data in Table 1 was obtained from different sources in [6]. The underlined values (with a precision of two decimal places)

~ ij ¼ f ðUij Þ U

ð8Þ

While two-dimensional analysis resulted in excellent energy conservation using these techniques [13], in three-dimensional analysis [14] the quality of the results depended strongly on the value of asymmetry factor g. While results were good for

5

S10 2nd largest boundaries

4 3 2 1 0 -1

Number of ordinates

are all quoted in the text from [6], while the single decimal precision values and the integer values were obtained from digitising graph plots from [6] with linear and logarithmic y-axes respectively. As the value of prescribed g increases, 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 [11,12]

-0.8

-0.6

-0.4

-0.2

0

cos θ

0.2

0.4

0.6

0.8

1

3

S10 largest boundaries 2

1

0 -1

-0.8

-0.6

-0.4

-0.2

0

cos θ

0.2

0.4

0.6

0.8

1

Fig. 2. Angular positions of directional ordinates and resultant angular interval boundaries for S10 quadrature set from ordinate, 2nd largest weighting (centre) and largest weighting (bottom).

208

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

Table 2 Boundaries of discretisation intervals as values of cos h for S10 quadrature set. Boundary No.

1 2 3 4 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

Weighting values Smallest

2nd smallest

2nd largest

Largest


1.0
0.997235444097142
0.947235255560311
0.863394888224938
0.863394888224938
0.600000377073665
0.464409333302683
0.417173700668710
0.369938068034737
0.307671657954962
0.254906913515273
0.207671280881300
0.171066546179900
0.000000000000006 0.171066546179889 0.207671280881289 0.254906913515262 0.307671657954951 0.369938068034726 0.417173700668699 0.464409333302672 0.600000377073653 0.701129043755646 0.863394888224927 0.947235255560300 0.997235444097131 1.0


1.0
0.984969222554198
0.954907667662594
0.892641257582819
0.892641257582819
0.774338513158457
0.721573768718769
0.684969034017369
0.648364299315969
0.595599554876280
0.524532443085891
0.487927708384491
0.451322973683091
0.356851708415145
0.210123675393706
0.162888042759733
0.000000000000006 0.162888042759722 0.210123675393695 0.356851708415133 0.451322973683079 0.487927708384480 0.524532443085880 0.595599554876268 0.648364299315957 0.684969034017357 0.721573768718757 0.774338513158446 0.845405624948835 0.892641257582808 0.954907667662583 0.984969222554187 1.0


1.0
0.981697632649300
0.945092897947900
0.892328153508211
0.892328153508211
0.727296810451916
0.619938068034736
0.569937879497905
0.488240812459099
0.458179257567495
0.421574522866095
0.374338890232122
0.171067488864056
0.123831856230083
0.071067111790395
0.065537999984679
0.000000000000006 0.065537999984667 0.071067111790383 0.123831856230072 0.171067488864045 0.374338890232111 0.421574522866084 0.458179257567484 0.488240812459088 0.569937879497893 0.619938068034724 0.727296810451905 0.845092520874227 0.892328153508200 0.945092897947889 0.981697632649289 1.0


1.0
0.976382183683014
0.949999811463169
0.913902110763809
0.913902110763809
0.809715517016348
0.755315448966288
0.740284671520486
0.674746671535813
0.609715705553180
0.547235632633976
0.426381995146186
0.408079627795486
0.393048850349684
0.325253328464193
0.310222551018391
0.158586661797527
0.137519738543969
0.018302367350706
0.000000000000006 0.018302367350694 0.137519738543958 0.158586661797515 0.310222551018379 0.325253328464181 0.393048850349672 0.408079627795474 0.426381995146174 0.547235632633964 0.609715705553169 0.674746671535801 0.740284671520475 0.755315448966276 0.809715517016336 0.913902110763798 0.949999811463158 1.0

g ¼ 0:2., progressive deformation of the phase function with g led to reduction in effective scattering and over-prediction of radiative flux for g ¼ 0:8, and by g ¼ 0:93 the discretised phase function was so strongly deformed that the results predicted were equal to those of a non-participating medium. A matrix-based normalisation technique has been developed [2] to address the lack of simultaneous conservation of scattered energy and asymmetry factor

~ ij ¼ ð1 þ Aij ÞUij U

ð9Þ

where the normalisation factors Aij need to be found such that the ~ ij values result in the inequality signs normalised phase function U in Eqs. (5) and (6) being replaced with equality signs:

1 X ~ ij ¼ 1; j ¼ 1nj
1 wj U 4p n 1 X ~ ij cos hij ¼ g; wj U 4p j¼1

i ¼ 1; 2; . . . ; n

i ¼ 1; 2; . . . ; n

ð10Þ

ð11Þ

and in the symmetry condition being satisfied:

~ ji ~ ij ¼ U U

ð12Þ

Since there are infinitely many possible solutions for this system of equations, the solution chosen is that where the matrix has the minimum norm. When this technique was applied to the problem which gave rise to the results in Table 1 [6], 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. The Rotationally Symmetric Spherical Discretisation (RSSD) technique has recently been developed [9] for the rational discretisation of scattering phase functions in the Discrete Ordinates Method (DOM). The RSSD technique has inherent energy conservation, and is suitable for any quadrature scheme. For each weighting w of a quadrature set (SN or otherwise), a set of intervals in angle h is calculated where a constant discretised phase function value U applies over each interval. The first interval (k ¼ 1) extends from h ¼ p to h ¼ x1 , where Dh ¼ p
x1 is determined only by the weighting w1 of the current ordinate direction:

w1 ¼

Z p x1

2p sin h dh ¼ 2p½
cos hpx1

ð13Þ

For each subsequent angular interval k, extending from xk to xkþ1 , it is necessary to determine the number nj of ordinate directions j of weighting wj which make an included angle hj with the current ordinate direction i such that xk > hj > xkþ1 . When the order N of the quadrature set SN is higher than 4, more than one weighting value is possible. Therefore, in each subsequent interval k, a number N wi of sets of amounts nj of incident ordinates of weighting wj pass through the angular range xk to xkþ1 , such that the angular interval may be defined by single or multiple weightings under conditions described in [9]:

Z N wk X nj  wj ¼ j¼1

xk

xkþ1

2p sin h dh ¼ 2p½
cos hxxkkþ1

ð14Þ

209

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

100 10

% error

1

1° energy

0.5° energy

0.25° energy

0.5/1.0 en

0.25/0.5 en

0.1°/0.25° en

1° g

0.5° g

0.25° g

0.5/1.0 g

0.25°/0.5° g

0.1°/0.25° g

0.1 0.01

0.001 0.0001 0.4

0.5

0.6

0.7 Prescribed g value

0.8

0.9

1

100

% error

10

1° energy

0.5° energy

0.25° energy

0.5/1.0 en

0.25/0.5 en

0.1°/0.25° en

1° g

0.5° g

0.25° g

0.5/1.0 g

0.25°/0.5° g

0.1°/0.25° g

1

0.1

0.01 0.94 Prescribed g value Fig. 3. Grid errors ne grid and ng grid for discrete HG phase function distributions (top), close-up showing local maxima (bottom).

The intervals and interval boundary values are dependent only upon the weightings of the quadrature set used, and are independent of the scattering phase function applied. RSSD intervals have been calculated and published for the different weightings of the S4, S6 and S8 ordinate sets [9]. Once the range of applicable boundary intervals have been determined, and a continuous or discrete scattering phase function distribution UðhÞ has been chosen, the constant phase

function values applicable in each interval may be determined from:

intxkþ1 UðhÞ2p sin h dh k PNwk j¼1 nj  wj x

Uk ¼

ð15Þ

Rearranging Eq. (15), the inherent energy conservation of the method can be proven:

100

10

S6 s

S6 l

S8 s

S8 m

S8 l

S10 smallest

S10 2nd smallest

S10 2nd largest

S10 largest

1 % error

S4

0.1

0.01

0.001

0.0001 0.4

0.5

0.6

0.7 Prescribed g value

0.8

Fig. 4. ne RSSD for discrete HG phase function distributions of 1° resolution.

0.9

1

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

0.03

0.15

0.02

0.1

0.01 % relave energy error

% relave energy error

210

0 -0.01 -0.02 -0.03 -0.04 -0.05 -0.06

S4

S6 s

S6 l

S8 s

S8 m

0 -0.05 S10 smallest

-0.1

S10 2nd smallest -0.15 S10 2nd largest -0.2

S8 l

-0.07

S10 largest

-0.25 0.4

0.5

0.6 0.7 0.8 Prescribed g value

0.9

1

0.4

0.5

0.6 0.7 0.8 Prescribed g value

0.9

1

0.9

1

0.9

1

0.035

0.02

S6 l 0.01

0.03

S6 s

S8 m

S8 s

0.025

S8 l

0.02

% relave energy error

S4

0.015 % relave energy error

0.05

0.005 0 -0.005

S10 smallest S10 2nd smallest S10 2nd largest

0.015

S10 largest

0.01 0.005 0

-0.005 -0.01

-0.01 0.4

0.5

0.6 0.7 0.8 Prescribed g value

0.9

0.4

1

0.6 0.7 0.8 Prescribed g value

0.01

0.005 S4

S6 s

0.005

S6 l

S8 s

0

S8 m

S8 l

0.003

% relave energy error

0.004 % relave energy error

0.5

0.002 0.001 0 -0.001

-0.005 -0.01 -0.015

S10 smallest

-0.02 S10 2nd smallest -0.025

S10 2nd largest

-0.03

S10 largest

-0.035 -0.002

-0.04 0.4

0.5

0.6 0.7 0.8 Prescribed g value

0.9

1

0.4

0.5

0.6 0.7 0.8 Prescribed g value

Fig. 5. ne RSSDrelativ e for discrete HG phase function distributions of uniform angular resolution: 1.0° (top), 0.5° (middle) and 0.25° (bottom).

Uk

N wk X

Z nj  wj ¼

xk

j¼1

¼

)

xkþ1

UðhÞ2p sin h dh)

NX Nwk v als X 1 inter Ui nj  wj 4p k¼1 j¼1 !

!

NX v als Z xkþ1 1 inter UðhÞ2p sin h dh 4p k¼1 xk

Z p NX Nwk v als X 1 inter 1 ðUk nj  wj Þ ¼ UðhÞ2p sin h dh ¼ 1 4p k¼1 4p 0 j¼1

ð16Þ

A major benefit of the RSSD 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 lookup table for future computations with different boundary conditions, different meshes and even different geometries. The RSSD technique was applied to two large (geometric optics) sphere scattering phase function distributions: diffusely reflecting spheres and transparent spheres (with refraction index n of 1.5) [9]. These distributions are of interest as test cases for packed bed radiative heat transfer [15]. The analytical scattering phase function distribution for diffusely reflecting (large) spheres [16]

UðhÞ ¼

8 ðsin h
h cos hÞ 3p

ð17Þ

211

0.03

0.15

0.02

0.1

0.01

0.05

0

% relave energy error

% relave energy error

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

-0.01 -0.02 -0.03 -0.04

S4

S6 s

-0.05

S6 l

S8 s

-0.06

S8 m

-0.05 S10 smallest

-0.1

S10 2nd smallest -0.15

S10 2nd largest

-0.2

S8 l

-0.07

S10 largest

-0.25 0.4

0.5

0.6 0.7 0.8 Prescribed g value

0.9

1

0.4

0.6 0.8 Prescribed g value

1

0.01

0.005 S4

S6 s

S6 l

S8 s

0.005

0.004 0.003

S8 m

0 -0.005

S8 l

% energy error

% relave energy error

0

0.002 0.001 0

-0.01 -0.015

S10 smallest

-0.02 S10 2nd smallest -0.025 S10 2nd largest

-0.03

-0.001

S10 largest

-0.035

-0.002

-0.04 0.4

0.5

0.6 0.7 0.8 Prescribed g value

0.9

1

0.4

0.5

0.6 0.7 0.8 Prescribed g value

0.9

1

Fig. 6. ne RSSDrelativ e for discrete HG phase function distributions with twin angular resolutions of 0.5/1.0° (top) and 0.1/0.25° (bottom).

gives an asymmetry factor g value of
0.44444, and when plotted shows zero forward scattering leading to a backscattering peak. Discrete scattering phase function distributions of 2°, 1° and 0.5° angular resolution were generated from the analytic distribution. These were in turn used to generate RSSD discretised phase function distributions for the S4, S6 and S8 ordinate sets, and produced errors for the sum of scattered energy conservation equal to or less than 0.0252 %, 0.0093% and 0.0029% for the 2°, 1° and 0.5° angular resolutions respectively, and errors in calculated asymmetry factor g of 1.272%, 1.248% and 1.249% respectively. The energy error decreased as grid resolution improved, but the g error, while small, appeared insensitive to grid resolution [9]. A ray tracing analysis performed on a transparent sphere (n = 1.5) resulted in a discrete scattering phase function distribution with an angular resolution of 2° [17]. The distribution exhibited strong forward scatter with the characteristic ‘‘fishtail” back-scatter distribution of dielectric spheres. Trapezoidal numerical integration yielded a value of asymmetry factor g of 0.6600. The RSSD discretisations of the phase function distribution for the S4, S6 and S8 ordinate sets produced maximum errors for the sum of scattered energy conservation and calculated asymmetry factor g of 0.035% and 1.273% respectively. The order of energy and g error was similar to the 2° angular resolution of the previous diffusely reflecting sphere test case, despite the large difference in asymmetry factor [9]. The purpose of this paper is to further explore, in a systematic manner when using RSSD, the combined effect on energy conservation error and calculated asymmetry factor error of:  the angular resolution of the input discrete scattering phase functions (such as may be obtained by ray tracing), and

 the asymmetry factor of scattering phase functions. The data for errors of energy and asymmetry factor conservation from [6] summarised in Table 1 were obtained using the HG family of phase functions, and therefore these phase functions serve as a good set of test cases for this purpose. 2. Interval boundaries for the S4, S6, S8 and S10 ordinate sets The boundaries of the discretisation intervals for the S4, S6 and S8 ordinate sets are given in [9], while for the S10 set were calculated (using the four S10 weightings obtained from [18]: wsmallest ffi 0:017370217029714, w2nd smallest ffi 0:094441160002949, w2nd largest ffi 0:114997165624522, wlargest ffi 0:148395116470556) using the method described in [9]. The numbers of intervals for the four weightings of the S10 ordinate set (26, 32, 32 and 36 respectively for the smallest to largest weightings, see Figs. 1 and 2) are larger than those applicable for SN ordinate sets of lower order N (10 for S4; 20 and 14 for the small and large weightings of S6; and 18, 28 and 12 for the small, medium and large weightings of S8 [9]), which allow for higher resolution of the scattering phase function being modelled. The interval boundaries for S10 are given (as values of cos h) to 15 significant digits in Table 2, for three reasons:  The precision of the input data (direction cosines and weights, [18]) is 15 significant digits.  By comparison of the positive and negative direction pairs, the resultant accuracy of the discretisation boundaries has degraded to 13 significant digits; and

212

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

100

2.5

S4

S4

2 S6 s S6 l S8 s

1

S8 m 0.1

S8 l S10 s'est

0.01

S6 l

0.5

S8 s

0 S8 m

-0.5

S8 l

-1 -1.5

S10 s'est

0.4

0.5

0.6 0.7 0.8 0.9 Prescribed g value

1

2.5

S10 2 l'est

-3 0.4

S10 l'est

0.5

0.6 0.7 0.8 0.9 Prescribed g value

1

1.2

S4

2

S10 2 s'est

-2.5

S10 2 l'est

0.001

S6 s

S10 l'est S4

0.8

S6 s

0.4

1

S6 l

0.5

S8 s

0

S8 m

-0.5 S8 l

-1 -1.5

S10 s'est

% absolute g error

1.5 % absolute g error

1

-2 S10 2 s'est

S6 l

0 -0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6

S10 s'est

-2

-2

S10 2 s'est S10 2 l'est

-3 0.4

0.5

0.6 0.7 0.8 0.9 Prescribed g value

1

1.2

S10 2 l'est

-2.8 0.4

S10 l'est

0.5

0.6 0.7 0.8 0.9 Prescribed g value

1

1.2

S4

0.8

S10 2 s'est

-2.4

-2.5

S6 s

S10 l'est S4

0.8

S6 s

0.4

0.4 S6 l

0 -0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6 S10 s'est

-2 -2.4 -2.8 0.4

0.5

0.6 0.7 0.8 0.9 Prescribed g value

% absolute g error

% absolute g error

S6 s

1.5 % absolute g error

% absolute g error

10

0

S6 l

-0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6

S10 s'est

-2

S10 2 s'est

-2.4

S10 2 l'est

-2.8

S10 2 s'est S10 2 l'est 0.4

1 S10 l'est

0.5

0.6

0.7

0.8

0.9

Prescribed g value

1 S10 l'est

Fig. 7. ng RSSD for discrete HG phase function distributions: Left – Uniform angular resolutions of 1.0° (top), 0.5° (centre), 0.25° (bottom), Right – Composite angular resolutions of 0.5°/1.0° (top), 0.25°/0.5° (centre) and 0.1°/0.25° grid (bottom).

 To allow validation of software written using the technique described in [9]. 3. Error calculations in HG distributions 3.1. ‘‘Grid” error in g and energy conservation (before RSSD discretization) Discrete distributions of HG scattering phase functions were constructed for prescribed values of g varying from 0.4 to 0.995 at different angular resolutions of the computational ‘‘grid”: three uniform angular resolutions (1°, 0.5°, 0.25°); and three composite angular resolutions (0.5°/1.0°, 0.25°/0.5° and 0.1°/0.25°). In the

composite grid, the finer resolution was applied to a small interval containing the forward scattering peak extending from 0° to 4.261°, and the coarser resolution from 4.261° to 180°, in a manner analogous to grid refinement in computational fluid dynamics and finite element analysis. The angular value of 4.261° corresponds to cos h ¼ 0:997235, the boundary of the smallest 1st interval of the SN ordinate sets under consideration (1st boundary from the smallest weighting of the for S10 ordinate set). Using trapezoidal numerical integration, the resultant sum of scattered energy

Esum grid ¼

Ndiscrete 1 X Uj 2p sinðhj ÞDhj 4p j¼1

ð18Þ

213

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

1.2

S4

1.2

S4

0.8

S6 s

0.8

S6 s

0.4 S6 l

0 -0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6

% relave g error

% relave g error

0.4

S10 s'est

-2.8 0.4

0.5

0.6 0.7 0.8 0.9 Prescribed g value

1

-0.8

S8 m

-1.2

S8 l

-1.6 -2.4

S10 2 l'est

-2.8

S10 s'est S10 2 s'est S10 2 l'est 0.4

0.5

S10 l'est

0.6 0.7 0.8 0.9 Prescribed g value

1 S10 l'est

1.2

S4

1.2

S4

0.8

S6 s

0.8

S6 s

0.4 S6 l

0 -0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6 S10 s'est

-2 -2.4 -2.8 0.4

0.5

0.6 0.7 0.8 0.9 Prescribed g value

S6 l

0 % relave g error

% relave g error

S8 s

S10 2 s'est

0.4

-0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6 S10 s'est

-2

S10 2 s'est

-2.4

S10 2 l'est

-2.8

S10 2 s'est S10 2 l'est 0.4

1

0.5

S10 l'est

1.2 0.8

0.6 0.7 0.8 0.9 Prescribed g value

1 S10 l'est

S4

1.2

S4

S6 s

0.8

S6 s

0.4

0.4 S6 l

0 -0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6 S10 s'est

-2 -2.4 -2.8 0.4

0.5

0.6 0.7 0.8 0.9 Prescribed g value

% relave g error

% relave g error

-0.4

-2

-2 -2.4

S6 l

0

-0.4

S8 s

-0.8

S8 m

-1.2

S8 l

-1.6 S10 s'est

-2

S10 2 s'est

-2.4

S10 2 l'est

-2.8

1

S6 l

0

S10 2 s'est S10 2 l'est 0.4

S10 l'est

0.5

0.6 0.7 0.8 0.9 Prescribed g value

1 S10 l'est

Fig. 8. ng RSSDrelativ e for discrete HG phase function distributions: Left – Uniform angular resolutions of 1.0° (top), 0.5° (centre), 0.25° (bottom), Right – Composite angular resolutions of 0.5°/1.0° (top), 0.25°/0.5° (centre) and 0.1°/0.25° grid (bottom).

and calculated asymmetry factor

g calc grid ¼

Ndiscrete 1 X Uj 2p sinðhj Þ cosðhj ÞDhj 4p j¼1

ng grid ¼ ð19Þ

were calculated. The grid errors before RSSD discretization, therefore, are:  Grid energy error:

ne grid ¼

Ndiscrete 1 X Uj 2p sinðhj ÞDhj
1 ¼ Esum grid
1 4p j¼1

 Grid g error:

ð20Þ

N discrete 1 X Uj 2p sinðhj Þ cosðhj ÞDhj
g prescribed 4p j¼1

¼ g calc grid
g prescribed

ð21Þ

For all cases except the composite resolutions for prescribed g values below 0.9, the grid error values (ne grid ; ng grid ) are positive in sign, i.e. the sum of scattered energy and g value is overpredicted. The error values are negative, however, for prescribed g values below 0.89 for the composite resolutions. To allow plotting on a logarithmic scale, the absolute values of the errors (jne grid j; jng grid j) are shown (Fig. 3, top), and where the error value changes sign for the composite resolutions (at g values between 0.875 and 0.9), a break is shown in the distributions.

214

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

5

Scaering phase funcon Φ

4.5

image: diff. refl. sphere

4

HG: g = 0.4444

3.5 3 2.5 2 1.5 1 0.5 0 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

cos θ Fig. 9. Comparison of two scattering phase function distributions where g ¼ 0:4444: HG distribution and mirror image of large diffusely reflecting sphere.

3

analycal distribuon

Scaering phase funcon Φ

2.5

S10 smallest S10 2nd smallest

2

S10 2nd largest 1.5

S10 largest

1

0.5

0 -1

-0.8

-0.6

-0.4

-0.2

0

cos θ

0.2

0.4

0.6

0.8

1

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

Table 3 Effect of angular resolution on grid errors for large diffusely reflecting sphere scattering phase function. Property

Esum grid ne grid % g calc grid ng grid %

Scattering phase function angular resolution 2°

1°

0.5°

0.25°

1.0001553 0.01553
0.444593 0.03340

1.0000318 0.00318
0.444476 0.00718

1.0000063 0.00063
0.444450 0.00133

1.0000013 0.00013
0.444446 0.00029

In each case, the grid g error is greater than the grid energy error ðng grid > ne grid Þ. Further, for both energy error ne grid and g error ng grid :  Error values increase for increasing values of prescribed g: ngrid ðg small Þ < ngrid ðg large Þ, before reaching a local maximum of about 1% error at g values between 0.97 and 1.0 (Fig. 3, bottom), and then increase again rapidly.  Error values increase for decreasing angular resolution of the discrete phase functions energy error:

ngrid ðDhlarge Þ > ngrid ðDhsmall Þ  Use of a coarse/fine composite grid reduces the level of error below that obtainable by a coarse-only uniform grid, to a level between that obtainable by uniform grids with resolutions of the coarse and fine regions of the composite grid respectively:

ngrid ðDhsmall Þuniform < ngrid ðDhsmall ; Dhlarge Þcomposite < ngrid ðDhlarge Þuniform Error is therefore driven to a degree by the resolution of the tiny 0°–4.261° part of the phase function distribution, pointing to the backward-difference ‘‘endpoint” of g ¼ 1:0 in the trapezoidal numerical integration being the driver of this error.

3.2. RSSD scattered energy error The discrete HG phase function distributions discussed earlier were discretized using the RSSD approach for the different weight-

215

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

Table 4 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 and S10 discretisations. Quad. Set

Weight

Scattering phase function angular resolution 2°

1°

Scatt. energy error %

gcalc

RSSD

0.5°

g error %

Scatt. energy error %

gcalc

RSSD

0.25°

g error %

Scatt. energy error %

gcalc

RSSD

g error %

Scatt. energy error %

gcalc

RSSD

g error %

S4

–

0.0077


0.4391


1.200

0.0019


0.4394


1.140

0.0012


0.4393


1.150

0.0003


0.4393


1.152

S6

S L

0.0224 0.0132


0.4431
0.4404


0.298
0.914

0.0012 0.0004


0.4429
0.4403


0.337
0.937

0.0007 0.0019


0.4430
0.4403


0.334
0.934


0.0002 0.0005


0.4430
0.4403


0.333
0.938

S8

S M L

0.0097 0.0449 0.0056


0.4417
0.4430
0.4390


0.615
0.323
1.225

0.0074 0.0136 0.0031


0.4416
0.4427
0.4389


0.632
0.382
1.239

0.0021 0.0016 0.0025


0.4416
0.4427
0.4389


0.649
0.400
1.247

0.0001 0.0005 0.0001


0.4415
0.4427
0.4389


0.653
0.401
1.250

S10

S 2S 2L L

0.0394
0.0305 0.0256 0.0088


0.4435
0.4431
0.4437
0.4436


0.206
0.301
0.158
0.181

0.0084
0.0039 0.0024
0.0079


0.4432
0.4434
0.4435
0.4436


0.282
0.242
0.219
0.192


0.0002
0.0023 0.0003
0.0050


0.4432
0.4434
0.4435
0.4436


0.286
0.236
0.222
0.194

0.0002 0.0008
0.0003
0.0015


0.4432
0.4434
0.4435
0.4436


0.286
0.231
0.223
0.191

Table 5 Values of forward
and backscatter points in large transparent sphere scattering phase function distributions of different resolutions. Grid resolution

1.0° 0.5° Twin 0.5°/0.25°

Forward scattering peak points

Backscatter peak points

1st

2nd

Inserted

Value (0°)

Angle

Value

52,540 210,110 840,439

1.0° 0.5° 0.25°

9.2344 8.2720 8.2647

ings of the S4, S6, S8 and S10 ordinate sets. The RSSD scattered energy error ne RSSD values

ne RSSD ¼ Esum RSSD
1 ¼

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

i ¼ 1; 2; . . . ; n

ð22Þ

for a uniform angular resolution of 1° (Fig. 4) form a distribution similar to the grid energy error ne grid distribution for the uniform 1° grid resolution (Fig. 3, top). Since the ne RSSD values in Fig. 4 span values from below 0.01 to about 50 and differ visibly from each other on a logarithmic scale, comparison with ne grid is difficult. For this reason, it is useful to compare the differences between Esum grid and Esum RSSD values expressed as percentages of the Esum grid values (ne RSSDrelativ e in Fig. 5 for uniform resolutions and Fig. 6 for composite resolutions):

ne RSSDrelativ e

Esum RSSD
Esum grid ¼  100% Esum grid

ð23Þ

RSSD energy error ne RSSD therefore comprises the grid energy error ne grid (Fig. 3, top) with an additive envelope of positive and negative corrections ne RSSDrelativ e (Figs. 5 and 6), nondimensionalised with respect to Esum grid :

ne RSSDrelativ e % þ Esum grid Þ
1 100 ne RSSDrelativ e % ¼ ðEsum grid
1Þ þ Esum grid  100

1st

2nd

Angle

Value (180°)

Angle

Value

0.028343° 0.014533° 0.006506°

17.187 34.042 67.967

179° 179.5° 179.75°

4.605 8.830 17.226

 ne RSSDrelativ e forms smoothly continuous functions of g, with the shape of a damped oscillation. Values begin at zero for g ¼ 1:0, peak at an absolute maximum (of positive or negative value) between g values of 0.75 and 0.98, then oscillate between local maxima and minima as g decreases in value (an exception to this trend is the medium weighting of the S8 ordinate set at an angular grid resolution of 0.25°: the maximum closest to the g = 1.0 point at g = 0.95 is not a global maximum since the minimum at g = 0.75 has a larger absolute value).  The extent of the envelope of ne RSSDrelativ e is small, everywhere smaller than 0.23% of Esum grid for the grid angular resolutions investigated in this paper. This means that ne RSSD is a strong function of grid energy error ne grid :

ne RSSD ¼ ne grid þ Esum grid 

)ne grid


ne RSSDrelativ e % 100

0:23%  ðEsum grid Þ2 < ne RSSD 100 < ne grid þ

Esum grid  1:0;

0:23%  ðEsum grid Þ2 100

ðEsum grid Þ2  1:0

)ne grid
0:0023 < ne RSSD < ne grid þ 0:0023

 The extent of the envelope of ‘‘relative error” itself decreases as the grid resolution improves: the maximum ‘‘relative error” drops from about 0.23% for 1° grid resolution to 0.036% for 0.25° grid resolution:

ne RSSD ¼ Esum RSSD
1 ¼ ðEsum grid 

ne grid
0:00036 < ne RSSD 0:25 < ne grid þ 0:00036 )ne RSSD

ne RSSDrelativ e % ¼ ne grid þ Esum grid  100

The following trends are observed:

ð25Þ

ð24Þ

ð26Þ

As the magnitude of grid error decreases with increasing grid resolution, the RSSD errors decrease faster, with the grid error being the limiting value. This suggests RSSD energy error can

216

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

1000000

1 deg resoluon

100000

Scaering phase funcon Φ

0.25 deg resoluon 10000 1000

100 10 1 0.1 0.01 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

cos θ

Scaering phase funcon Φ

10 9

1 deg resoluon

8

0.25 deg resoluon

7

6 5 4 3 2 1 0 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

cos θ Fig. 11. Large transparent sphere scattering phase function: comparison of distributions of 1° and 0.25° angular resolution in log (top) and linear (bottom) scales.

be controlled to required levels by judicious choice of grid resolution.  The ne RSSDrelativ e distributions for composite resolutions (Fig. 6) are practically identical to those of uniform distributions (Fig. 5) comprising the same angular resolution as the coarse region of the composite distribution. While the choice of composite resolutions over uniform resolutions makes no appreciable change to ne RSSDrelativ e , it does improve grid error ne grid (Fig. 3). This means that final RSSD energy accuracy (dependent on both grid error ne grid and relative error ne RSSDrelativ e ) can be improved by increasing grid resolution in regions of steep gradients and peaks in the scattering phase function distribution, without needing increased resolution elsewhere. 3.3. RSSD error in calculated g value The error in calculated g values after RSSD discretization, ng RSSD , is given by:

ng RSSD ¼ g calc RSSD
g prescribed ¼

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

i ¼ 1; 2; . . . ; n

ð27Þ

ng RSSD was determined for the discrete HG phase function distributions mentioned earlier with uniform (1.0°, 0.5° and 0.25°) and composite (0.5°/1.0°, 0.25°/0.5° and 0.1°/0.25°) angular resolutions

(Fig. 7). For the uniform angular resolution of 1.0° (Fig. 7, top left), the contrast between the high values of ng RSSD (about 59%) for a g value of 0.995 and the low ng RSSD values (between 1% and
2.7%) at lower g values necessitated plotting absolute values against a logarithmic scale. Between prescribed g values of 0.95 and 0.995, ng RSSD differs significantly for the uniform resolutions, implying strong grid resolution dependence. The distributions for the composite graphs are practically identical to those of uniform distributions comprising the same angular resolution as the fine region of the composite distribution. For prescribed g values below 0.95, however the distributions of ng RSSD for both uniform resolutions and composite resolutions all appear remarkably similar, indicating that g error is insensitive to grid resolution. To gain more insight into the systematic effect of grid resolution on g error, the differences between the values of calculated g obtained from the grid (g calc grid ) and from RSSD discretization (g calc RSSD ) respectively were calculated (as done in the previous section for sum of scattered energy) as percentages of g calc grid (Fig. 8):

ng RSSDrelativ e ¼

g calc RSSD
g calc grid  100% g calc grid

ð28Þ

RSSD g error ng RSSD therefore comprises the g energy error ng grid (Fig. 3, top) with an additive envelope of positive and negative corrections (Fig. 8), non-dimensionalised with respect to g calc grid :

217

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

ng RSSD ¼ g calc RSSD
g prescribed   ng RSSDrelativ e % þ g calc grid
g prescribed ¼ g calc grid  100 ¼ ðg calc grid
g prescribed Þ þ g calc grid 

)ng RSSD ¼ ng grid þ g calc grid 

ng RSSDrelativ e % 100

ng RSSDrelativ e % 100

ð29Þ

The ng RSSDrelativ e distributions display the following trends:  The ng RSSDrelativ e distributions are practically identical for all resolutions, with very minor variations. This indicates that the ng RSSDrelativ e envelope, unlike the ne RSSDrelativ e envelope, is insensitive to grid resolution.  ng RSSDrelativ e , like ne RSSDrelativ e , forms smoothly continuous functions of g, with the shape of a damped oscillation. ng RSSDrelativ e values begin at zero for g ¼ 1:0, peak at the 1st local maximum (of positive or negative value) between g values of 0.8 and 0.97, then oscillate between local maxima and minima as g decreases in value.  All ng RSSDrelativ e values lie in a band between
2.7% and 1%:

ng RSSD

4. Error calculations in real scattering phase function distributions

ng RSSDrelativ e % ¼ ng grid þ g calc grid  100

 2  2 2:7%  g calc grid 1:0%  g calc grid < ng RSSD < ng grid þ )ng grid
100 100 0 < g calc grid < 1:0; 0 < ðg calc grid Þ2 < 1:0; )ng grid
0:027 < ng RSSD < ng grid þ 0:01

ð30Þ

 This band decreases to a range from
1% to 1% for the S4, S6 and S10 ordinate sets:

ng grid
0:01 < ng RSSD < ng grid þ 0:01

 The effect of the underprediction by the small weightings of the S6, S8 and S10 ordinate sets in the maximum error zone (between prescribed g values of 0.7 and 0.9) appears to be offset to some degree by the overprediction of the large weightings of the same ordinate sets. Comparing Figs. 7 and 8, it can be deduced that the effect of grid resolution is most visible at combinations of high values of prescribed g and course grid resolution, i.e. it is expected that grid resolutions coarser than 1° would be visible at prescribed g values below 0.95. This is illustrated in Fig. 3 (top): grid energy error (and therefore grid g error) will exceed 1% at progressively lower values of prescribed g for decreasing grid resolution. Outside this ‘‘grid sensitivity zone” (i.e. in the case of sufficiently fine grid resolution and/or lower values of prescribed g), g error is insensitive to grid resolution. The implication of this is that the magnitude of the large g errors of Table 1 can be avoided: when using the RSSD approach, the magnitude of the g error can be bounded to between
2.7% and 1% (even to
1% to 1% for the S4, S6 and S10 ordinate sets) even at high prescribed values of g, as long as grid error is managed by choosing an appropriately fine scattering phase function resolution for the g value of the scattering phase function concerned.

ð31Þ

The next step is to test these findings against ‘‘real” scattering phase function distributions:  Large diffusely reflecting sphere  Large transparent sphere Since the publication of [9] (the results of which are quoted in the introduction to this paper), it was discovered that for the numerical integration of the aforementioned scattering phase function distributions, the numerical treatment of the backwarddifference endpoints (cos h ¼
1 and cos h ¼ 1) of the grid were incorrect, which would lead to inaccurate error results. In addition, the RSSD interval boundaries for S10 ordinate set are now available (see Table 2).

10

Griffith

9

S10 smallest

8

Phase funcon Φ

S10 2nd smallest 7

S10 2nd largest

6

S10 largest

5

4 3 2 1 0 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

cos θ Fig. 12. Large transparent sphere scattering phase function: comparison of analytical and discretised distributions for S10 quadrature set.

218

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

Table 6 Large transparent sphere values of g, scattered energy error and g error for the three scattering phase function grid resolutions: comparison of non-discretised and RSSD S4, S6, S8 and S10 discretisations. Quadrature set

Weight

Scattering phase function angular resolution 1°

0.5°

Scattered energy error %

Calculated g

Twin

g error %

Scattered energy error %

0.66025

Calculated g

g error %

Scattered energy error %

Calculated g

0.66034

g error %

Grid

–

S4

–

0.0220

0.65752


0.4140

0.00296

0.65771


0.3984

0.00176

0.66011 0.65748


0.3984

S6

S L


0.0161
0.0047

0.65528 0.65289


0.7536
1.1155


0.04015 0.00458

0.65553 0.65296


0.7279
1.1173


0.00014 0.00424

0.65502 0.65274


0.7724
1.1167

S8

S M L

0.0144
0.0018 0.0024

0.65176 0.65824 0.65519


1.2861
0.3044
0.7673


0.00433
0.00247 0.03698

0.65169 0.65838 0.65498


1.3096
0.2962
0.8124


0.00278

0.00455
0.00097

0.65145 0.65818 0.65510


1.3125
0.2929
0.7590

S10

S 2S 2L L


0.0014
0.0227 0.0185 0.0685

0.65752 0.65910 0.65847 0.65919


0.4141
0.1744
0.2703
0.1613

0.00179 0.00458
0.00081 0.00261

0.65755 0.65935 0.65864 0.65946


0.4230
0.1497
0.2572
0.1334

0.00198
0.00174 0.00296 0.00199

0.65732 0.65919 0.65846 0.65924


0.4228
0.1403
0.2507
0.1325

14 transparent sphere

Scaering phase funcon Φ

12

HG: g = 0.660 10

8 6 4

2 0 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

cos θ Fig. 13. Comparison of two scattering phase function distributions where g ¼ 0:660: HG distribution and large transparent sphere.

The correct numerical treatment was applied to the grid endpoints, and an updated set of RSSD scattering phase function values, including for the four weightings of S10 ordinate set, for the large diffusely reflecting sphere and large transparent sphere is given in Tables A1–A10. 4.1. Scattering phase function discretisation test case 1: large diffusely reflecting sphere Since the theoretical values for sum of scattered energy and asymmetry factor g are unity and -0.44444 respectively [9], the determination of grid error is straightforward. The analytical scattering phase function distribution was used to generate discrete scattering phase function distributions of 2°, 1°, 0.5° and 0.25° angular resolution. Using trapezoidal numerical integration, the sum of scattered energy Esum grid and asymmetry factor g calc grid were calculated (Table 3). Comparing the results for the 2°, 1° and 0.5° angular resolutions with those given in [9], the percentage error of both scattered energy (ne grid Þ and g (ng grid Þhas approximately halved since correct numerical treatment of the backwarddifference grid endpoints has been applied. There is not a g ¼
0:4444 point among the HG grid error distributions in Fig. 3 (top) for comparison with the values in Table 3.

The RSSD boundaries for the different ordinate sets are all symmetrical about cos h ¼ 0, so appropriate insight will be derived from an HG distribution where g ¼ þ0:4444. Fig. 9 shows a comparison of an HG distribution (g ¼ þ0:4444) and the diffusely reflecting sphere scattering phase function, mirrored about the cos h ¼ 0 axis. HG distribution ne grid values for g ¼ þ0:4444 are 0.00520%, 0.00105% and 0.00023% for angular resolutions of 1°, 0.5° and 0.25° respectively. These values are different from, but remain a largely consistent fraction (164%–183%) of, the values in Table 3. The HG distribution ng grid values for g ¼ þ0:4444 (0.00982 %, 0.00198% and 0.00055% for angular resolutions of 1°, 0.5° and 0.25° respectively) are more variable fractions of the values in Table 3 (130%–190%) than for the ne grid values, but the grid error trends appear consistent. Despite having the same value of g, the HG errors are larger than those for the diffusely scattering sphere because the distribution is more concave with cos h (see Fig. 9), leading to more averaging error and greater backward-difference grid endpoint error. The discrete scattering phase function distributions of 2°, 1°, 0.5° and 0.25° angular resolution were discretized using the RSSD approach for the different weightings of the S4, S6, S8 and S10 ordinate sets. Since a finer grid resolution (0.25°) was available for the large diffusely reflecting sphere than in [9], it was used to generate

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

the large diffusely reflecting sphere RSSD discretised scattering phase function values in Tables A1–A10. The RSSD discretisation for the S10 ordinate set is shown in Fig. 10. The sum of scattered energy Esum RSSD and asymmetry factor g calc RSSD for the different ordinate sets were then calculated, and the errors ne RSSD and ng RSSD determined (Table 4). The same trends emerge as determined from the HG RSSD distributions:  The sum of scattered energy error continually improves with increasing grid resolution  The g error, while small, is largely insensitive to grid resolution. Differences can also be identified. For the large diffusely reflecting sphere scattering phase function, the maximum g error occurs for the large weighting of the S8 ordinate set (with values of
1.225 to
1.250) followed by the S4 ordinate set, irrespective of grid resolution. For the HG distribution with a value of prescribed g of 0.4444, however, the maximum g error occurs for the small weighting of the S8 ordinate set (with values of
1.197 to
1.230) followed by the large weighting of the S6 ordinate set (Fig. 7). Clearly, g error depends not only upon the value of prescribed g but also the distribution itself.

4.2. Scattering phase function discretisation test case 2: Large transparent sphere with n=1.5 The second test case is for spherical transparent spheres with a refractive index (n) of 1.5. A scattering phase function distribution in 1° resolution had been used in [9], obtained from ray tracing analysis using the Zemax software [17]. The distribution is characteristic of dielectric spheres: a strong forward-scattering peak with a ‘‘fishtail” backscattering pattern. Subsequent to the work published in [9], two more ray tracing analyses were performed [19], yielding scattering phase function distributions of 0.5° and 0.25° resolution. The 0.25° resolution distribution was not used in its entirety, but was used to create a twin 0.5°/0.25° resolution grid, with the higher resolution in regions of steep gradients and peaks. Fig. 11 shows a comparison of the 1° resolution grid and the 0.25° resolution inserts, (in both logarithmic and linear scales) but without the 0.5° resolution grid, in order to highlight the contrast. In all three distributions, the difference in scattering phase function value between the first two points (h ¼ 0 and the adjacent point: h ¼ 1 ; h ¼ 0:5 or h ¼ 0:25 ) is so great that numerical integration fails to achieve a unity result for sum of scattered energy. While the difference in scattering phase function value between the final two points (h ¼ 180 and the adjacent point: h ¼ 179 ; h ¼ 179:5 or h ¼ 179:75 , describing the backscatter peak) are also significant, it is dwarfed by the forward scatter difference. Rather than truncating the forward peak, an additional point was inserted in each case between the first two points, where the scattering phase function value of the inserted point is the same as that of the second point, but the angular position of the inserted point was varied until the numerical integration yielded a scattered energy sum of 1.0. All the relevant values are given in Table 5. Once the angular position had been determined in each case for the inserted point that resulted in a sum of scattered energy of unity (again with the correct numerical treatment at the endpoints for numerical integration), the value of g was numerically determined for each scattering phase function grid resolution. The grid values of g obtained are 0.66025, 0.66034 and 0.66011 for the 1.0°, 0.5° and twin resolution respectively. Ideal values of g for this test case are unavailable for comparison purposes, so instead the values of energy error and g error were calculated for an HG distri-

219

bution (with prescribed g value of 0.6601) for indicative purposes: energy error values of 0.01258%, 0.00269% and 0.00064% were obtained for angular resolutions of 1.0°, 0.5° and twin respectively, while the corresponding g errors were 0.01826%, 0.00393% and 0.00095%. The discrete scattering phase function distributions of 1.0°, 0.5° and twin angular resolution were discretized using the RSSD approach for the different weightings of the S4, S6, S8 and S10 ordinate sets. As with the diffusely scattering sphere test case, an updated set of RSSD scattering phase function values for the large diffusely reflecting sphere is given in Tables A1–A10. The RSSD discretisation for the S10 ordinate set is shown in Fig. 12. The sum of scattered energy and asymmetry factor for the different ordinate sets were then calculated, and the errors determined (Table 6). As is the case with the diffusely scattering sphere and the HG distributions, the sum of scattered energy error continually improves with increasing grid resolution, while the g error is insensitive to grid resolution. For the large transparent sphere scattering phase function, the maximum g error occurs for the small weighting of the S8 ordinate set (with values of
1.2861 to
1.3125) followed by the large weighting of the S6 ordinate set, irrespective of grid resolution. For an HG distribution with a value of prescribed g of 0.6601, however, the maximum g error also occurs for the small weighting of the S8 ordinate set, but for larger error values (with values of
2.157 to
2.200) followed by the small weighting of the S6 ordinate set (Fig. 7). It is suspected that the maximum g error occurs for the same weighting of the same ordinate set for both the HG and transparent sphere distributions due to the similarity of ‘‘concavities” of the two distributions (Fig. 13). As with the large diffusely reflecting sphere scattering phase function, the difference in value in error values for the HG distribution and the large transparent sphere for the same value of g reflects the influence of the distribution itself

5. Conclusions In addition to the discretisation boundaries for the different weightings of the S4, S6 and S8 ordinate sets given in [9], the set of boundary values for the four weightings of the S10 ordinate set have been calculated and presented. The effect of angular resolution and asymmetry factor (g) of discrete scattering phase functions on energy conservation error and calculated asymmetry factor error when using RSSD discretization has been explored, making use of the Henyey–Greenstein (HG) family of distributions. Distributions of RSSD energy error ne RSSD and RSSD g error ng RSSD show close similarity to those of the ‘‘grid” energy and g errors error values (ne grid ; ng grid ) which arise from the resolution of the discrete scattering phase functions before RSSD discretization, and have the following attributes:  The error in calculated asymmetry factor g is always greater than scattered energy error.  Error values increase for increasing values of prescribed g  Error values increase for coarsening angular resolution of the discrete phase function grid  Use of a composite grid with coarser and finer regions reduces the level of error below that obtainable by a coarse-only uniform grid, to a level between that obtainable by uniform grids with resolutions of the coarse and fine regions of the composite grid respectively. When the RSSD energy and g errors are expressed as differences between the values obtained from the grid and from RSSD

220

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

discretization as relative errors, further systematic trends are observed:  The ‘‘relative errors” of both energy and g are smoothly continuous functions of g, with the shape of a damped oscillation. ‘‘Relative error” values begin at zero for g ¼ 1:0, peak at an absolute maximum (of positive or negative value) between g values of 0.66 and 0.98, then oscillate between local maxima and minima as g decreases in value. RSSD error (both energy and g) therefore comprises the grid error with a superimposed envelope of positive and negative RSSD corrections.  The extent of the envelope of ‘‘relative error” is small. RSSD relative energy error ne RSSDrelativ e is everywhere smaller than 0.23% (for 1° grid resolution), and RSSD relative g error ng RSSDrelativ e lies in a band between
2.7% and 1% (dropping to a range from
1% to 1% for the S4, S6 and S10 ordinate sets)  The extent of the envelope of RSSD relative energy error ne RSSDrelativ e itself decreases as the grid resolution improves: the maximum relative error drops from about 0.23% for 1° grid resolution to 0.036% for 0.25° grid resolution. The use of composite grids does not improve relative energy error.  The distributions of RSSD relative g error ng RSSDrelativ e , by contrast, are practically identical for all grid resolutions, both uniform and composite, with very minor variations. This indicates that ng RSSDrelativ e is insensitive to grid resolution. As the magnitude of grid energy error decreases with increasing grid resolution, the RSSD energy relative error decreases faster, leaving the grid error as the limiting value. This suggests scattered energy error using RSSD can be controlled to required levels by judicious choice of grid resolution. While the use of composite grids does not improve relative energy error, it improves grid energy error. The final RSSD energy accuracy, which is dependent on both grid error and relative error, can therefore be improved by use of composite grids, by increasing grid resolution in regions of steep gradients and peaks in the scattering phase function distribution, without needing increased resolution elsewhere. Calculated asymmetry factor error resulting from RSSD discretisation behaves differently to RSSD scattered energy error. A ‘‘grid sensitivity zone” exists at combinations of high values of prescribed g and course grid resolution, where g error displays strong grid resolution dependence. For the HG family of scattering phase functions, the ‘‘grid sensitivity zone” comprises regions where g > 0:9, g > 0:95 and g > 0:98 for scattering phase function angular resolutions of 1°, 0.5° and 0.25° respectively. Outside of the ‘‘grid sensitivity zone, g error:  Is insensitive to grid resolution  Exhibits maximum (positive and negative) values between prescribed g values of 0.7 and 0.9.

 Lies in a band between
2.7% and 1% (decreasing to a range from
1% to 1% for the S4, S6 and S10 ordinate sets). The effect of the underprediction by the small weightings of the S6, S8 and S10 ordinate sets in the maximum error zone (prescribed g values of 0.7 and 0.9) appears to be offset to some degree by the overprediction of the large weightings of the same ordinate sets. The implication of this is that the magnitude of the large g errors of Table 1 can be avoided, and that the magnitude of the g error can be bounded to between
2.7% and 1% when using the RSSD approach (even at high prescribed values of g), as long as grid error is managed by choosing an appropriately fine scattering phase function resolution for the g value of the scattering phase function concerned. The grid resolution conclusions above are supported by two real scattering phase function distributions: for an asymmetry factor of
0.4444, grid angular resolutions of 2°, 1°, 0.5° and 0.25° resulted in maximum RSSD g errors of
1.225%,
1.239%,
1.247% and
1.250% respectively, while for an asymmetry factor of 0.660, grid angular resolutions of 1° and 0.5° resulted in maximum RSSD g errors of
1.2861% and
1.3096%. The previously reported benefit of RSSD discretization of scattering phase function distributions for use in DOM [9] is 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 lookup table for future computations with different boundary conditions, different meshes and even different geometries. It can now be concluded that by judicious choice of discrete scattering phase function grid resolution:  Errors in sum of scattered energy can be controlled to userrequired levels (in this paper error levels as low as 0.002% have been achieved).  Errors in calculated g can be bounded to between
2.7% and 1%, and to between
1% to 1% for the S4, S6 and S10 ordinate sets.

Acknowledgements The work described in this paper has been supported by the CSIR.

Appendix A. Discretised scattering phase function values for S4, S6 and S8 ordinate sets See Tables A1–A10.

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.563337 2.278949 1.856963 1.480254 1.086864 0.642407 0.369186 0.190281 0.056712 0.007725

1.337173 0.045977 0.041186 0.049990 0.073435 0.223106 0.671181 1.642555 3.862763 7.166036

221

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225 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.633622 2.539600 2.320655 2.049853 1.824247 1.613028 1.420222 1.289367 1.160027 0.959283 0.745639 0.579057 0.486103 0.403329 0.296993 0.200390 0.118152 0.046979 0.009601 0.001312

1.188399 1.344834 0.131579 0.041162 0.040882 0.043759 0.055201 0.048822 0.044983 0.101346 0.143290 0.251565 0.380966 0.565955 0.927348 1.490317 2.376719 4.188614 6.688911 8.491919

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.594244 2.456646 2.157055 1.606541 1.209333 1.082983 0.945713 0.757574 0.638492 0.542626 0.307367 0.088733 0.020770 0.004449

1.116974 0.840326 0.043822 0.047083 0.044394 0.046027 0.108794 0.138521 0.193031 0.296442 0.971305 3.041383 5.525950 7.690701

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.646166 2.402485 1.929274 1.614042 1.497805 1.385949 1.229631 1.011834 0.870214 0.827721 0.699556 0.528175 0.423877 0.357786 0.296006 0.163094 0.032825 0.000635

1.321190 0.504219 0.041255 0.043345 0.047469 0.061556 0.044358 0.079283 0.146243 0.127776 0.161957 0.317223 0.511505 0.693811 0.925782 1.920036 5.122156 8.763804

222

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

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.631823 2.499929 2.369672 2.299741 2.154269 2.014168 1.920891 1.723300 1.562806 1.508123 1.419608 1.333158 1.222123 0.990589 0.718837 0.533342 0.456868 0.403304 0.351901 0.322195 0.244588 0.162032 0.128961 0.086856 0.049976 0.035601 0.014944 0.001441

1.175764 0.983423 0.059801 0.045698 0.041908 0.040449 0.040638 0.041974 0.044334 0.046831 0.053179 0.068082 0.044246 0.089437 0.157047 0.308671 0.437330 0.562787 0.712041 0.815992 1.201405 1.812798 2.193886 2.936724 3.930978 4.554562 6.133468 8.447354

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.575132 2.378846 2.089669 1.758479 1.282418 0.913297 0.786897 0.500456 0.229201 0.107009 0.034499 0.006394

1.215870 0.277608 0.041404 0.041682 0.048565 0.139003 0.130180 0.418826 1.301903 2.611818 4.691524 7.357926

Table A7 Discretised scattering phase function values for diffusely reflecting and transparent spheres for smallest weighting ordinate of S10 set. Interval no

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

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere


1.0
0.997235
0.947235
0.863395
0.701129
0.600000
0.464409
0.417174
0.369938
0.307672
0.254907
0.207671
0.171067 0.0 0.171067


0.997235
0.947235
0.863395
0.701129
0.600000
0.464409
0.417174
0.369938
0.307672
0.254907
0.207671
0.171067 0.0 0.171067 0.207671

2.665297 2.596816 2.438183 2.169735 1.903683 1.682524 1.520487 1.440248 1.349860 1.257832 1.180097 1.116499 0.967048 0.738897 0.611783

2.384899 1.038943 0.663439 0.042640 0.040666 0.042448 0.046228 0.050993 0.061647 0.043853 0.044696 0.045508 0.097831 0.146299 0.215101

223

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225 Table A7 (continued) Interval no

16 17 18 19 20 21 22 23 24 25 26

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere

0.207671 0.254907 0.307672 0.369938 0.417174 0.464409 0.600000 0.701129 0.863395 0.947235 0.997235

0.254907 0.307672 0.369938 0.417174 0.464409 0.600000 0.701129 0.863395 0.947235 0.997235 1.0

0.563330 0.507641 0.446491 0.390835 0.345116 0.263320 0.168844 0.083725 0.024020 0.004099 0.000047

0.265768 0.343648 0.459898 0.594675 0.735040 1.094248 1.756654 3.070123 5.330475 7.725988 9.185792

Table A8 Discretised scattering phase function values for diffusely reflecting and transparent spheres for 2nd smallest weighting ordinate of S10 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 29 30 31 32

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere


1.0
0.984969
0.954908
0.892641
0.845406
0.774339
0.721574
0.684969
0.648364
0.595600
0.524532
0.487928
0.451323
0.356852
0.210124
0.162888 0.0 0.162888 0.210124 0.356852 0.451323 0.487928 0.524532 0.595600 0.648364 0.684969 0.721574 0.774339 0.845406 0.892641 0.954908 0.984969


0.984969
0.954908
0.892641
0.845406
0.774339
0.721574
0.684969
0.648364
0.595600
0.524532
0.487928
0.451323
0.356852
0.210124
0.162888 0.0 0.162888 0.210124 0.356852 0.451323 0.487928 0.524532 0.595600 0.648364 0.684969 0.721574 0.774339 0.845406 0.892641 0.954908 0.984969 1.0

2.647130 2.590814 2.480792 2.355734 2.226952 2.097349 2.006852 1.934597 1.848500 1.732787 1.635122 1.570486 1.458248 1.26194 1.112385 0.961206 0.743964 0.615041 0.505985 0.380681 0.318183 0.285229 0.239274 0.189865 0.156842 0.131470 0.102770 0.067304 0.038350 0.017271 0.004311 0.000590

1.337422 0.996201 1.000999 0.05544 0.042861 0.041436 0.040465 0.040606 0.040779 0.041920 0.042832 0.044089 0.050114 0.050628 0.045579 0.100428 0.144019 0.212407 0.354315 0.626852 0.831709 0.970732 1.217608 1.558943 1.861269 2.160954 2.598785 3.392844 4.430369 5.839324 7.584977 8.779081

Table A9 Discretised scattering phase function values for diffusely reflecting and transparent spheres for 2nd largest weighting ordinate of S10 set. Interval no

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

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere


1.0
0.981698
0.945093
0.892328
0.845093
0.72730
0.619938
0.569938
0.488241
0.458179
0.421575
0.374339
0.171067


0.981698
0.945093
0.892328
0.845093
0.72730
0.619938
0.569938
0.488241
0.458179
0.421575
0.374339
0.171067
0.123832

2.642847 2.574811 2.46885 2.355109 2.177286 1.948663 1.797397 1.676385 1.576837 1.518914 1.447623 1.245703 1.054788

1.273917 1.030478 0.974097 0.055481 0.042314 0.040587 0.040932 0.042603 0.043856 0.046287 0.050491 0.050015 0.046439 (continued on next page)

224

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225

Table A9 (continued) Interval no

14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere


0.123832
0.071067
0.065538 0.0 0.065538 0.071067 0.123832 0.171067 0.374339 0.421575 0.458179 0.488241 0.569938 0.619938 0.727297 0.845093 0.892328 0.945093 0.981698


0.071067
0.065538 0.0 0.065538 0.071067 0.123832 0.171067 0.374339 0.421575 0.458179 0.488241 0.569938 0.619938 0.727297 0.845093 0.892328 0.945093 0.981698 1.0

0.982925 0.942843 0.893120 0.805646 0.758774 0.722993 0.661513 0.518516 0.386503 0.345899 0.314745 0.265542 0.210907 0.152342 0.080794 0.038492 0.018878 0.005794 0.000792

0.047641 0.049177 0.179363 0.127608 0.132098 0.144318 0.176485 0.342325 0.606493 0.731846 0.844578 1.070754 1.401932 1.929764 3.086146 4.423585 5.674822 7.279574 8.700867

Table A10 Discretised scattering phase function values for diffusely reflecting and transparent spheres for largest weighting ordinate of S10 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 29 30 31 32 33 34 35 36 37 38

Boundaries (cos h)

Scattering phase function values

Lower

Upper

Diffusely reflecting sphere

Transparent sphere


1.0
0.976382
0.950000
0.913902
0.875254
0.809716
0.755315
0.740285
0.674747
0.609716
0.547236
0.426382
0.408080
0.393049
0.325253
0.310223
0.158587
0.137520
0.018302 0.0 0.018302 0.137520 0.158587 0.310223 0.325253 0.393049 0.408080 0.426382 0.547236 0.609716 0.674747 0.740285 0.755315 0.809716 0.875254 0.913902 0.950000 0.976382


0.976382
0.950000
0.913902
0.875254
0.809716
0.755315
0.740285
0.674747
0.609716
0.547236
0.426382
0.408080
0.393049
0.325253
0.310223
0.158587
0.137520
0.018302 0.0 0.018302 0.137520 0.158587 0.310223 0.325253 0.393049 0.408080 0.426382 0.547236 0.609716 0.674747 0.740285 0.755315 0.809716 0.875254 0.913902 0.950000 0.976382 1.0

2.636188 2.574014 2.499430 2.413103 2.297300 2.168918 2.096651 2.015445 1.887330 1.766747 1.601209 1.480072 1.451992 1.383204 1.315769 1.185632 1.055537 0.955804 0.861409 0.836339 0.748012 0.660788 0.560551 0.468335 0.425515 0.383786 0.367660 0.303080 0.224153 0.174706 0.128749 0.102602 0.082237 0.050711 0.027653 0.014378 0.005752 0.001162

1.206835 1.019203 1.555368 0.051117 0.046365 0.042709 0.040365 0.040465 0.040698 0.041049 0.044095 0.048388 0.050043 0.061113 0.051962 0.044689 0.046487 0.106176 0.138988 0.128877 0.139021 0.176326 0.276775 0.414283 0.507274 0.613975 0.662481 0.902385 1.311905 1.691354 2.203863 2.589277 3.007042 3.925241 5.009902 6.072011 7.252854 8.548373

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.

T.H. Roos et al. / International Journal of Heat and Mass Transfer 101 (2016) 205–225 [5] B. Hunter, Z. Guo, Comparison of quadrature schemes in DOM for anisotropic scattering radiative transfer analysis, Numer. Heat Transfer, Part B 63 (2013) 485–507. [6] 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, Numerical Heat Transfer, Part B, 2012, pp. 203–222. [7] 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. [8] C. Caliot, Numerical Methods in Radiative Transfer: Introduction to DOM, FVM and MCM, Promes – CNRS, 11 June 2010. [Online]. Available: . [Accessed 30 January 2012] [9] T.H. Roos, T.M. Harms, A new radiative transfer scattering phase function discretisation approach with inherent energy conservation, Int. J. Heat Mass Transfer 73 (2014) 789–803. [10] 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.

225

[11] T.-K. Kim, H. Lee, Effect of anisotropic scattering on radiative heat transfer in two-dimensional rectangular enclosures, Int. J. Heat and Mass Transfer 31 (1988) 1711–1721. [12] W.J. Wiscombe, On initialisation, error and flux conservation in the doubling method, J. Quant. Spectrosc. Radiat. Transfer 18 (1976) 637–658. [13] 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. [14] 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. [15] B. Variot, T. Menigault, G. Flamant, Modelling and optimization of a two-slab selective volumetric solar receiver, Sol. Energy 53 (4) (1994) 359–368. [16] J.R. Howell, R. Siegel, M.P. Mengüç, Therm. Radiat. Heat Transfer, fifth ed., CRC Press, Boca Raton, 2011. [17] D. Griffith, Private Communication, CSIR, Pretoria, 2013. [18] G. Colomer, Numerical Methods for Radiative Heat Transfer, Ph.D. Thesis Universitat Politècnica de Catalunya, 2006. [19] D. Griffith, Private Communication, CSIR, Pretoria, 2014.