Radiative transfer in two-dimensional anisotropic scattering media with collimated incidence

J. Quant. Spectrosc. Radiat. Transfer Vol. 42, No. 3, pp. 225-238, 1989 Printed in Great Britain

RADIATIVE

TRANSFER

ANISOTROPIC

IN TWO-DIMENSIONAL

SCATTERING

COLLIMATED

0022-4073/89 $3.00+ 0.00 Maxwell Pergamon Macmillanplc

MEDIA

WITH

INCIDENCE

TAE-KUK KIM a n d HAEOK SKARDA LEE Department of Mechanical Engineering, University of Minnesota, Minneapolis, MN 55455, U.S.A. (Received 20 September 1988)

A b s t r a c t - - A n accurate analysis is performed by using the S - N discrete ordinates method for

radiative transfer in two-dimensional rectangular enclosures exposed to an arbitrarily inclined, collimated incident beam. The enclosed gray medium absorbs and anisotropically scatters radiative energy. General anisotropic Mie-scattering phase functions are treated by Legendre polynomial expansions. A highly accurate S-14 approximation is used for the analysis. Benchmark results for the average incident radiations, the radiative fluxes and the source functions, obtained from the radiative intensity fields, are presented for purely scattering and absorbing-scattering media. Normal incidence and oblique incidence problems are analyzed. Anisotropy of the phase functions is found to have a strong effect on the radiative transfer. The reflected energy may differ by an order of magnitude when the phase function is changed from highly forward-scattering (phase function F I ) to backward-scattering (phase function B2). The transmitted energy changes by a factor of two for the two phase functions F I and B2. The effect of phase-function anisotropy on radiative transfer is found to be more significant for collimated incidence than for diffuse incidence.

INTRODUCTION There are many engineering applications in which collimated incident radiation plays an important role. In solar energy applications, the primary energy input is produced by collimated radiation. Collimated laser beams are often used as the radiation source in measurement systems. The analysis of radiative transfer for an incident collimated beam is also important in astrophysics and atmospheric science. Many astrophysics and atmospheric science studies deal with radiative transfer in one-dimensional planar media. Chandrasekhar ~presents a complete formulation and analysis of one-dimensional, collimated radiative transfer in anisotropic scattering media. General anisotropic Mie-scattering effects and azimuthal dependence have been included by Garcia and Siewert 2 and by Flicke. 3 Two-dimensional collimated incidence problems have been extensively studied by Crosbie and his co-workers. Crosbie and Schrenker4 treat a collimated incidence problem in a rectangular, isotropic-scattering medium. Crosbie and Koewing5 treat a collimated incidence problem in a medium of finite thickness, using Ambarzumian's method for the scaled isotropic equation of transfer. Radiative transfer in an anisotropically scattering, cylindrical medium has been studied numerically and experimentally by Nelson et al, 6 who also used a scaled isotropic equation of transfer for their analytical study. Crosbie and Dougherty7 treat the radiative transfer problem in a two-dimensional, linearly anisotropic-scattering, cylindrical medium exposed to a laser beam. Most of the specified studies contain results only at the boundaries of the enclosure. In this study, radiative transfer in a two-dimensional rectangular enclosure, with a collimated beam incident through a transparent top wall at an arbitrary incident angle, is modeled by using the S - N discrete ordinates method. Except for the transparent top wall, all other walls are considered to be opaque. Figure 1 shows a schematic diagram of the enclosure considered. The participating medium in the enclosure is gray, and it absorbs and anisotropically scatters radiative energy. The S - N computer code, written for a previous study,s has been modified for the collimated 225

TAE-KUK KIM and HAEOK SKARDA LEE

226

IC

wall 4](transparent wall) P

/,

'['yHor H

r

~

wall 2

Y

, L

wall 3

/

p

xx or x

z

X,.Lor L Fig. 1. The system geometry.

incidence problem. The radiative intensity field, the average incident radiation, the source function, and the radiative flux are obtained using the S-14 approximation. Basic formulations are presented for the analysis of two-dimensional, anisotropic-scattering problems with collimated incidence. A computational algorithm, developed for an effective treatment of general anisotropic scattering phase functions in symmetric two-dimensional problems, is discussed. The effects of anisotropy of the scattering phase functions, incidence angle, aspect ratio, optical thickness, scattering albedo, and boundary reflectivity on radiative transfer are examined. Both the pure scattering and the absorbing-scattering cases are investigated. RADIATIVE T R A N S F E R EQUATION When a collimated beam is incident upon a scattering medium, the angular distributions of the radiative intensities in the medium usually have very strong peaks. The extreme angular variation in the intensities is due to the unidirectional incident radiation, while the medium scatters radiation in all directions. The actual radiative intensity/.d can be separated into two components to facilitate the numerical analysis. One is the component which is scattered within the medium. The other component is the transmitted component, which travels along the direction of incidence and experiences an exponential decay. The radiative transfer equation, which describes the scattered component of the intensity I s, in a rectangular, absorbing, and anisotropically scattering medium is written by referring to Ref. 9, viz. bt

. + ~ ~z>.

I l,(z.,.,Zy, f~) = S(z.~, "Cy,f~).

(I)

The source function is expressed as

S(zx, zy, f~) = -~n

./~(Zx,Zy, f~')~(t)'; f~) dtT + ~ Ic(zx, tyL, f~)~(D~; fl) exp{ -- (~yn -- ~y)/l¢cl}, (2)

where the collimated beam Ic is incident through the top boundary at an angle tic = (/~c, ¢c). The ' indicates the incident direction and the (/~, ~) are the direction cosines in the x and y coordinates, respectively, that describe the direction f~. All other variables are explained in the Nomenclature.

RT in 2-D anisotropic scattering media

227

I([~ i ) ~[~(['~B;~F )

Planc o f Symmetry

~'~..~.w'~k~_.........." ~

(xyplane,

. •

Fig. 2. Scattered intensity diagram for the front and back hemispheres. The actual intensity is related to the scattered component of the intensity through the relation, Ia(Zx, "rv, ta) = Is(zx, zy, ta) + [c('rx, "ryrl,tac) e x p { - (Tya -- ~y)/l~01}~(ta - ~ ) .

(3)

The scattering phase function ¢(ta'; ta) may be approximated by a finite series of Legendre polynomials as in a previous study, 8 i.e. K

~(ta'; ta) = ¢(cos ~k) = ~ C, Pt(cos tl'),

(4)

1=0

where K + 1 is the number of terms used to approximate the phase function and ~k is the scattering angle. The Ct are the expansion coefficients. The integrand of the source function [Eq. (2)] is rewritten to take advantage of the symmetry condition implied in the two-dimensional problem. Figure 1 shows the x - y plane as the plane of symmetry. The relation between the radiative intensities in the front and back hemispheres, with the x - y plane serving as the base plane for the hemispheres, is illustrated in Fig. 2. Figure 2 shows a bottom view of the enclosure outlined in Fig. 1, with the front hemisphere in the positive Z direction indicated by the subscript F and the back hemisphere in the negative Z direction indicated by the subscript B. The magnitudes of the intensities in the corresponding mirror-image directions must be equal at the plane of symmetry, i.e., I(ta') = I(ta~:) = I(ta~), where the position dependence (zx, Zy) is simply implied. This relation holds as long as the collimated incidence is parallel to the x - y plane. Figure 2 also illustrates the relation between the front- and the back-hemisphere scattered intensities. For a point on the symmetry plane (Xx, ~y), the intensities in two incident directions ta~ in the front hemisphere and its mirror image ta~ in the back hemisphere are shown. The scattered intensities into the directions taF and its mirror image are shown as I(ta~:)-¢(ta~; taF) and I(ta~,) • ~ ( a . ,t . riB). The total scattered intensity to the front hemisphere is written as (I • ¢)nF and is given by

( t - ~).~ = i(ta;~)~(n~; fir) + z(ta~)~(ta~,; tar) = I(ta')[C~(ta~; n r ) + ¢(ta~; nr)].

(Sa)

The total scattered intensity to the back hemisphere is written as ( I . ¢)nB and is given by

( l . ~)n. = I(f~D~(tar, taB) + l(taB)~(taB, taB) i

t

.

i

t

.

/. = I(ta / )[~(tar,t . taB) + ~(taB, taB)].

(Sb)

Since ~(tab; taF) = ¢(ta~; taB) and ¢(ta~:; taB) = ¢(ta~; tar) for azimuthally-symmetric phase functions, Eqs. (5a) and (5b) are equivalent and the following single expression remains:

(I. ~),~ = (I. ~)n, = I(ta~)~(ta~; tar) + I(ta~)~(ta~; tar) t

t ,

= I(ta )['/~(tar, tar) + ¢(ta~; tar)].

(6)

228

TAE-KUK KIM and HAEOK SKARDA LEE

This expression is the integrand of the source function [Eq. (2)]. If scattering is isotropic, Eq. (6) becomes ( I . ¢)nv = (I" ¢)n, = 2I(IT).

(7)

Using Eq. (6), the range of integration for Eq. (2) becomes 2n in only one of the hemispheres, either in front or in the back of the symmetry plane. The enclosure walls are all gray and diffusely reflecting with the same reflectivity, including the transparent top wall. The intensities from the walls are written as

In" l)'[Lw(Zx, zy, •') dfl' + P-P-In" fldI=(zx, Zyn, fl¢)exp{ - ( z y . - Zy)/l~cl}

Zy, l)) = p- ~ Jn .if<0

TC

for

n.~>0,

n.D'<0

and

n.lac<0,

(8)

where p is the diffuse wall reflectivity and n is the inward normal to the boundary wall. The only radiative source is applied at the transparent top wall from an external collimated source, viz.

l¢(zx, "CyH,]-Lc, ~c)

~-

I¢ = constant for tp = 0,

(9)

where the incident direction described by kt¢ and ~¢ is parallel to the x - y plane since ~p = 0. The ranges of #¢ and ~¢ are 0 ~< #c < 1 and - l < ~c < 0. The direction, #, = 0 and ~ = - l, indicates normal incidence to the top wall. Once the intensity field is obtained, the average incident radiation G(Tx, ~y) and the radiative fluxes Q~+, Q~-, Qy+, and Q7 can be obtained as

]

I,(z~, Ty, ta) d n + I< exp{ -(Ty. - *y)/l¢ol} ,

G(zx, Zy) = ~ Q ~ (T~, ~:,) =

#I,(~x, ~ , n ) d n + Io.~

exp{-

(zyH - ~y)/I ~d },

(10) (1 l a)

>0

Q;(Tx, ~Y) = f ,
(1 lb)

Q~-(~x. %.) = f

(12a)

~I~(~x.Ty. n) dfL >0

Qy(zx, zy) = f

<0

~l~(zx, Zy, [2) dD + I¢~cexp{-(Zyn -

z,)/l~0i}.

(12b)

For isotropic scattering media, the source function and the average incident radiation can be related by comparing Eqs. (2) and (10) as

S(z x, Zy) = coG(zx, Zy).

(13)

The net radiative fluxes are obtained by summing the positive and negative components of the Q as follows: Qx(~., ~,.) =

. + Q.~-(Tx, L.), Qx+ (z.~, T,.)

(14a)

Q:,(zx, Zy) = Q.,+.(zx, zy) + Qy(zx, Zy).

(14b)

The fluxes Q.; and Q,:- are always negative in view of the definition given in Eqs, (11) and (12).

NUMERICAL CALCULATIONS The S - N discrete ordinates method replaces the radiative transfer equation with a coupled set of equations for a finite number of M ordinate directions.~ The integral appearing in the source term [Eq. (2)] is replaced by a quadrature of order M with the appropriate angular weights WI, where Z~= ~ 2w,~ = 4~. For a specific ordinate direction m defined by fire = (#m, ~m), the source

RT in 2-D anisotropic scattering media

229

function is written by using the modified expression of the integrand given in Eq. (6), i.e., s( x,

n.) =

n.t) + m'=l (.0

+ ~ Id0.5{a'(~¢F; t2,,~) + ~(t2¢B; t2mF)}]exp{-- (~yx -- ~y)/l~¢l}-

(15)

The source term is assumed to be a constant in each control volume. Spatially averaging the collimated incidence term over a control volume results in better accuracy. For the two-dimensional rectangular geometry shown in Fig. 1, the average of the exponential term over Azy can be expressed as

exp[{ -(zyH - zy)/l~d}l,ve= ~

exp{-(zyH - ~,)/1~¢1} d-Cy y.s

= Azy [exp{ - - ( T y

H --

Zy.n)/l~d} -- e x p { - - (ZyH -- T,)/1¢¢1}],

(16)

where the subscripts n and s designate the north and the south faces of the control volume, respectively. The phase function appearing in the source function expression [Eq. (2) or Eq. (15)] is re-normalized to avoid normalization error, which is the major error source in the energy balance. The re-normalization is performed by using the S - N quadrature set and the following expression: ~(f~,,,; tim) = 1

J~(flm'; fire)

(17)

E Wm'~(Om';~'~ra)

4-~ m'= 1

The re-normalized ~ differs only slightly from the original • and leads to a better overall energy balance in the calculation. The finite difference discrete ordinates equation of transfer is solved over the rectangular enclosure, subdivided into small control volumes by a M X × M Y mesh. Using the weighted diamond difference scheme, the radiative intensity at the center of the control volume is related to the intensities at each face as in Ref. 8. The spatial weights in the weighted difference scheme are obtained from a positive scheme suggested by Lathrop) ° The required ordinates set/~,, and ~, and the corresponding angular weights w,, are taken from the TWOTRAN-II code. I~ The radiative intensity field is solved iteratively, as explained in Ref. 8, by using the spatial-angular iteration procedure suggested by Lathrop and Brinkley." The boundary conditions given in Eqs. (8) and (9) are used to compute the wall intensities and the collimated source term. In order to check the energy balance in each control volume and accelerate convergence of the numerical solution, the energy balance is computed during every iteration. An expression for the divergence of the radiative flux is derived, which includes the collimated incidence term, i.e., V.Q=

~

(1

~

o )[fo z,(z.,Zy, fl')d['I'+Icexp{-(zyH-zy)/,,¢,}l

=-4n(1-to)G(%,

zy).

(18)

For the (i,j) control volume, a re-balance f a c t o r f j is obtained from a balance of all the incoming and outgoing energies, which gives the V. Q expression given in Eq. (18). If the outgoing energies are multiplied byrd, and the incoming energies are multiplied by the corresponding f, the re-balance factor is obtained as f~,j = Azy(Q~f~+ ,,j + O +xwfi-Lj) + Azx(O~o~,j+ , + O+ f~,j- , ) , Azy(Q~%+ Qfw) + Axx(Qg + Q~) + SCA,,/

(19)

where

= 41tAV(1 - og)G(~x,-, ~yj).

(20)

230

TAE-KuKK[M and HAEOKSKARDALEE 1.0

i I I =' I -. . . . . . .

,i II

I '! ' f i,J+l -.

'° ntl; n Cl ÷

XW i

I

~ ~ SCA~,I.~ fn-'1,1~,~

........

i O~w--

....

.....

fl,,J ~

,i

o;j!*o;.

t u

f ~ ....

I,j-1

I i !

0.9

,%°

o/O

A.LX" "

I

Q+ !I Ixo

f

!

O/O/

o

.o /

0.8 @

Phase function : F2

@'@

O--O : 1-D

! I I I I

0.7

a=l

0--@

:

z~--Zx

: 0=2.5

&--A

: am5

n--n

: a=10

I

Z I

Fig. 3. Energy components and grids for the energy re-balance factors.

0.6 0.0

i

I

I

I

I

0.2

0.4

0.6

0.8

1.0

T* y Fig. 4. Effect of aspect ratios on the y-component of the centerline net flux (p = 0 and co = I).

The energy components and grid points for computing the re-balance factors are indicated in Fig. 3. The exponential is evaluated by using Eq. (16) for better accuracy. When the energy is balanced in a control volume, the re-balance factor for that control volume should be equal to one. The re-balance factors are used to accelerate the iteration convergence, by multiplying the Ip by f~j before a new source term is evaluated. If any of the re-balance factors have negative values or do not converge after a given iteration limit, a representative re-balance factor for the entire system is obtained by balancing the overall energy for the system. A similar, more detailed discussion of the re-balance technique can be found in Ref. 11. The convergence of the numerical solution is checked during the iteration process by using the cell average intensities (Ip) at current and previous iterations. It is assumed that convergence is obtained when the maximum change of the intensities is <0.0001%.

RESULTS A numerical study has been performed for some selected cases to study the effect of anisotropic scattering, Most of the cases consider purely scattering media (co = l) for which scattering is most significant. The collimated source described by Eq. (9) is applied at the top wall (z* --- zy/z m = 1.O), as is shown in Fig. 1. The scattering phase functions used for this study have been described previously in Ref. 8 (Table 1 and Fig. 2), except for the F3. The phase functions marked F1, F2, and F3 are the forward scattering phase functions and have 13 terms (the asymmetry factor defined 12 by C~/3 is 0.84534), nine terms (C~/3 = 0.66972), and three terms (C1/3 = 0.40000) in the Legendre polynomial series expansion, respectively. The B1 and B2 are the backward scattering phase functions. The B 1 has six terms (C~/3 = -0.18841) in the expansion, while B2 has three terms (Ct/3 = -0.40000). The F1 and F2 are the phase functions for spheres with a refractive index of 1.33 and size parameters of 5 and 2, respectively. The B 1 is cited from Ozi~ik ~ and B2 is from Siegel and Howell? 4 The phase function F3 is the forward scattering equivalent of B2 with a positive expansion coefficient C~. The isotropic scattering phase function used for comparison has an asymmetry factor of zero. The dimensionless quantities G*=4rcG/l~cI~l, S*=4ttS/l~oIcl, Q*=Qx/(~J~), Q~+*--Q+~ /1~cI¢[, Q2* = Q2/(~clc), Q* = Qy/(~clc), Q-~* = Q-~ /1~I~[, and Q f * = Q f /(~¢I¢) are presented by using the normalized coordinates r* = Zx/Z~Land '~y0 = *Cy/'~y H • By normalizing with respect to either (¢~I~) or I~¢I¢1, most of the dimensionless quantities are positive. Results for the normal incidence problems ( # c = 0 , ¢ ¢ = - 1 ) are presented in the x - y graphs of Figs. 4-12 and in

R T in 2 - D a n i s o t r o p i c s c a t t e r i n g m e d i a , , , ,

0.5

0.4

: ISOTROPIC : ISOTROPIC [ 4 ] : F1

l--I

: F'2

vmv

: F3

A--A r~mD

: B1 :B2

_..o.~o...o--o--o--o--o--o-.o~.o...o ~

0.3

Io~

o~o e~e A--a,

231

0.2

I-('

A s A~"

.A"

wO ~

:-.-

A'~A.. A

o _ o - O - O ~ ~

",_~ "4

, o -°

- ~ 1 ~ "~

0,1 rV.v..V~V •

•

•

•

•

•

•

•

•

•

0.0 'P"~"I'f-xSx-z-x-xTz---'--=~ffi . . . . 0.0

0.2

0.4.

0.6

V~V.~$,

V~

"~-- "~;~"~ 0.8

1.0

q-* X

Fig. 5. Effect o f p h a s e f u n c t i o n s o n the reflected flux c o m p o n e n t (p = 0, co = l, a n d ZxL = ~yH = !).

1.0

/

0.8

A"A

~'A

o--o e--e &--& I1--11

: : : : :

".,

ISOTROPIC IS0TROPIC [ 4 ] F1 F2

v--v F3 ~ _ _ s , : B1 r~D : B2

0.2

0.0

,

0.0

I

i

0.2

I

I

0.4.

I

I

0.6

I

,

0.8

1.0

7.* X

Fig. 6. Effect o f p h a s e f u n c t i o n s o n the t r a n s m i t t e d flux c o m p o n e n t (p = 0, co = 1, a n d TxL = zy. = l). 0.30

,

0.25

,

,

o~o e~e AmI1~11

: ISOTROPIC : ISOTROPIC [4.] : F1 :F2 :F3 : B1 : B2

v~v ~'~ C~

0.20

I

0.15

A~A Omn

,

O. I 0

0.05

0.00 0.0

0.2

0.4

0.6

0.8

1.0

Fig. 7. Effect o f p h a s e f u n c t i o n s o n t h e side wall e n e r g y loss (p = O, to = 1, a n d Z~L = zy. = 1).

232

HAEOK SKAKDALEE

T A E - K u K KIM and 1.0

'

0.8 - , ~.4-~"

'

,,,A- A

_,...-"

_m_m....m.~'m

m...m-'m .,~i.~"mi

0--O"

m~.m ..-.,

-

~.~O ~

O.0-00~ A.A -A

0.6 A*,.*,..,~r'l~ DO.re.C]-rl

0.2

0,0

i

0.0

: ISOTROPIC : F1 :F2 :F3 : B1 : 82

o--o I--I z~--z~ l--I a--a m--El

0.4

J

0.2

I

i

t

0.4

0.6

0.8

i

1.0

"T* Y Fig. 8. Effect o f phase functions on the y - c o m p o n e n t centerlin¢ net flux (p = 0, co = l , and ~ L = ~yH = 1)' 1.0

.

.

.

.

.

.

~

~0

~00"v

o--O __ - ' 'u-- O-O~'~1~ ,.,a

~0

i/I

0t0 "la ~-~

0.6

-

•O ~

0.4

~ ~'a'~

o.O ~O~

~'~

~.z~ ~

0.2

I

O.0

•

~/

Phase 0--0 I--I a--Z~

0.0

Z~"

_/@

,

I

0.2

,

function : F2 : t ~ = 1.0 : ~ : 0.5 : ~ = 0.I

i

0.4,

,

,,I

0.6

0.8

1.0

Y Fig. 9. Effect o f scattering albedos on the y - c o m p o n e n t centerlin¢ net f l u x (p = 0 and ~xL = ~yH = l). 1.0

.

~

.

,

.

.

,

~0..O - ' O ~ O - ' O ' ' 0 0.8

@" •

• O'°'O"O

~0.0.0.00

~..

0.6

',~>'

0.4,

z¢~.~'

0.2

&.A"

'

,

~&z~~

/,.A

_.Q.O"

&.4"

.a.a a.a -a"

.,,i

,&"

Phase function : F2

&.A °

ua~aaaa a a.

O ~ O

:I"

l--I

: "rxL = ' r ,),nu =

t~,, a~a

0.0

'

0.0

~

0.2

' 0.4

x-.m'r

,,=10

yN

2,5

"rxL - "ryH- 5.0 -rxL = "ryH= I0.0

'

0.6

0.8

1.0

T* Y Fig. 10. Effect o f optical thickness on the y - c o m p o n e n t centerlin¢ net flux (p = 0, co = 1, and a = 1).

R T in 2 - D anisotropic scattering media

233

1.0 ....~....0..~.0~0....0__0~0._

0..~0 ' ' 0 0.8

i..l.011

)o.O.O O 0

/I/l~l~rl

0/0/0/0 A/A

~Q/

/ A "/

0.0.O •

n/&/A

A~

0.4

Z~A ~

"~AI~" A A

O~0 ~

0.6

.

0 0-0-0

A~

Phase function : F2

0 - - O :p=0.0

~

"~'n~

O--o

: p = 0.5

~--~

: p = 0.9

0.2

0.0

I 0.2

I 0.0

I 0.4

I

I 0.6

I 0.8

1.0

I"*

Y Fig. 11. Effect of surface reflectivity on the y-component centerline net flux (to = 0 a n d ZxL = ~y8 = l).

Tables 1-3. Oblique incidence results (0
0.8

/

A/

Collimated

:

0--o

: ISOTROPIC

Q--D

:B2

~z~ \ "

0,6

": IE~ :~

0.0.0_.0~0---0--0~0~0---0--0~0--0--0._0_0..~.^

0.4

& A & & & & A A & & & & & & & & & & ,t, & ji~& &

I~.-O.O

, 0.0

- Diffuse . . . . . L8] I 0.2

: o--e &--&

~ 0.4

II--II ,

: ISOTROPIC : F1

:B2 I 0.6

-~n~

0.8

1.0

I"* X

Fig. 12. Comparison of transmitted fuxes for the collimated and the diffuse incidences (p = 0, to = 1, and TxL = TyH = 1). ~RT 42/~E

234

TAE-KUKKIM and HAEOK SKARDALEE Table 1. Source function comparisons for isotropic scattering (p=O, c o = l , andzxL=ZyH=I), T w

S*(0,x;)

S*(0.5,T*)

Y

Y

current Ref. 0.00000 0.00611 0.02037 0.04251 0.07216 0.10884 0.15194 0.20076 0.25449 0.31225 0.37309 0.43602 0.50000 0.56398 0.62691 0.68775 0.74551 0.79924 0.84806 0.89116 0.92784 0.95749 0.97963 0.99389 1.00000

0.5377 0.5446 0.5583 0.5798 0.6066 0.6391 0.6748 0.7142 0.7559 0.8008 0.8473 0.8961 0.9449 0.9943 1.0415 1.0875 1.1287 1.1661 1.1951 1.2180 1.2305 1.2392 1.2402 1.2415 1.2388

[4]

(0.5362) (0.5440) (0.5591) (0.5803) (0.6066) (0.6378) (0.6731) (0.7122) (0.7547) (0.7999) (0.8472) (0.8959) (0.9452) (0.9940) (1.0413) (1.0860) (1.1269) (1.1628) (1.1930) (1.2164) (1.2329) (1.2424) (1.2455) (1.2437) (1.2396)

current 0.6478 0.6593 0.6842 0.7214 0.7650 0.8140 0.8664 0.9228 0.9821 1.0438 1.1065 1.1693 1.2303 1.2883 1.3412 1.3881 1.4266 1.4568 1.4761 1.4853 1.4805 1.4660 1.4424 1.4249 1.4145

Ref.

[4]

(0.6472) (0.6615) (0.6875) (0.7224) (0.7644) (0.8122) (0.8648) (0.9214) (0.9809) (1.0425) (1.1051) (1.1676) (1.2286) (1.2865) (1.3398) (1.3868) (1.4257) (1.4552) (1.4743) (1.4827) (1.4807) (1.4697) (1.4520) (1.4319) (1.4161)

comparison for the source function. The centerline and edge comparisons for a pure scattering square enclosure show excellent agreement. The source function is related to the average incident radiation by Eq. (13). Two-dimensional anisotropic Mie-scattering results are verified by checking the energy balance in each control volume. The maximum control volume energy balance error, which is equal to the maximum of (1 - f 0 ) , is of the order 10 -6 with the re-normalized phase functions. The overall energy balance error is of the order l0 -7 for all the cases considered in this study.

Table 2. The net y-direction fluxes for different phase functions (p = 0, ~o = 1, and ZxL= 5'H = 1). F1

T*

y

0.00000 0.00282 0.01222 0.02905 0.05301 0.08375 0.12080 0.16360 0.21149 0.26377 0.31962 0.37820 0.43864 0.50000 0.56136 0.62180 0.68038 0.73623 0.78851 0.83640 0.87920 0.91625 0.94699 0.97095 0.98778 0.99718 1.00000

x* = 0 x

0.6376 0.6387 0.6408 0.6448 0.6506 0.6582 0.6676 0.6788 0.6918 0.7066 0.7231 0,7411 0.7605 0.7811 0.8026 0.8248 0.8474 0.8699 0.8917 0.9124 0.9315 0.9486 0.9631 0.9746 0.9826 0.9869 0.9881

F2

0.5 0.9029 0.9034 0.9048 0.9073 0.9107 0.9150 0.9200 0.9257 0.9317 0.9378 0.9438 0.9497 0.9552 0.9603 0.9650 0.9690 0.9725 0.9752 0.9774 0.9791 0.9803 0.9810 0.9814 0.9815 0.9813 0.9808 0.9803

x* = 0

B1 0.5

x

0.5835 0.5845 0.5870 0.5918 0.5987 0.6080 0.6195 0.6332 0,6491 0.6671 0.6869 0.7084 0.7313 0.7553 0.7801 0.8054 0.8308 0.8557 0.8797 0.9022 0.9227 0.9407 0.9558 0.9676 0.9757 0.9800 0.9811

0.7543 0.7550 0.7576 0.7623 0.7693 0.7785 0.7900 0.8034 0.8183 0.8345 0.8515 0.8687 0.8855 0.9016 0.9163 0.9293 0.9405 0.9496 0.9570 0.9626 0.9668 0.9697 0.9715 0.9723 0.9724 0.9721 0,9717

Zx = 0 0.4380 0.4386 0.4407 0.4448 0.4508 0.4589 0.4692 0.4817 0.4965 0.5136 0.5329 0.5543 0.5777 0.6028 0.6294 0.6570 0.6853 0.7136 0.7413 0.7677 0.7920 0.8135 0.8316 0.8460

B2 0.5 0.4933 0.4934 0.4942 0.4960 0.4990 0.5035 0.5097 0.5175 0.5270 0.5384 0.5515 0.5663 0.5826 0.5999 0.6180 0.6362 0.6540 0.6710 0.6868 0.7010 0.7133 0.7236 0.7315 0.7368

0.8562 i.7402 0.8617 0.8636

.7420 .7424

x'=0

0.5

x

0.4198 0.4204 0.4223 0.4261 0.4318 0.4394 0.4491 0.4609 0.4750 0.4913 0.5098 0.5305 0.5531 0,5775 0,6035 0.6307 0,6586 0.6867 0.7143 0.7408 0.7654 0.7873 0.8059 0.8208 0.8314 0.8373 0.8394

0.4625 0.4625 0.4631 0.4644 0.4667 0.4703 0.4753 0.4818

~.4899 .4997 0.5112 0.5245 0.5392 ! 0.5552 0.5719 0.5891 0.6061 0.6224 0 6378 0.6517 0.6640 0.6743 0.6822 0.6878 0.6914 0.6934 0.6939

RT in 2-D anisotropic scattering media

235

Table 3. The average incidence radiations for different phase functions (p = 0 , o~ = 1, and ZxL ZyH I). =

F1 T*

y

0.00000 0.00282 0.01222 0.02905 0.05301 0.08375 0.12080 0.16360 0.21149 0.26377 0.31962 0.37820 0.43864 0.50000 0.56136 0.62180 0.68038 0.73623 0.78851 0.83640 0.87920 0.91625 0.94699 0.97095 0.98778 0.99718 1.00000

T* x

= 0

0.5

0.6943 0.6959 0.6990 0.7045 0.7123 0.7221 0.7336 0.7466 0.7609 0.7765 0.7933 0.8112 0.8300 0.8495 0.8695 0.8898 0.9099 0.9297 0.9486 0.9661 0.9821 0.9961 1.0078 1.0168 1.0230 1.0262 1.0270

=

F2

0.9937 0.9949 0.9985 1.0045 1.0122 1.02061 1.0294 1.0381 1.0466 1.0545 1.0617 1.0678 1.0727 1.0762 1.0783 1.0789 1.0780 1.0758 1.0724 1.0682 1.0634 1.0584 1.0533 1.0486 1.0448 1.0422 1.0410

T"

X

= 0

B1 0.5

T* X

0.7056 0.9372 0.7080 0.9402 0.7143 0.9495 0.7250 0.9650 0.7394 0.9846 0.7566 1.0064 0.7758 1.0296 0.7964 1.0537 0.8182 1.0783 0.8408 1.1029 0.8640 1.1268 0.8875 1.1490 0.9109 1.1686 0.9337 1.1847 0.9556 1.1964 0.9760 1.2030 0.9946 1.2042 1.0111 1.2000 1.0250 1.1908 1.0360 1.1773 1.0437 1.1605 1.04831.1411 1.0504 1.1200 1.0508ii.0991 1.0500 1.0820 1.0488 1.0714 1.0481 1.0676

=

B2

0.5

0

0.5030 0.5061 0.5155 0.5317 0.5542 0.5824 0.6151 0.6514 0.6910 0~7336 0.7788 0.8263 0.8755 0.9256 0.9760 1.0256 1.0736 1.1186 1.1596 1.1948 1.2232 1.2441 1.2584 1.2673 1.2723 1.2741 1.2739

0.5935 0.5989 0.6161 0.6460 0.6856 0.7319 0.7829 0.8383 0.8976 0.9602 1.0251 1.0911 1.1572 1.2220 1.2841 1.3420 1.3945 1.4401 1.4778 1.5065 1.5252 1.5327 1.5285 1.5155 1.5006 1.4907 1.4869

~*X

= 0

0.5

0.4710 0.5418 0.4741 0.5474 0.4835 0.5649 0.499810.5956 0.5226 0.6367 0.5512 0.6853 0.5847 0.7396 0.6222 0.7990 0.6632 0.8630 0.7075 0.9308 0.7547 1.0013 0.8045 1.0733 0.8562 1.1456 0.9093 1.2169 0.9628 1.2856 1.0158 1.3504 1.0675 1.4098 1.1166 1.4624 1.1617 1.5070 1.2014 1.5425 1.2344 1.5677 1.2600 1.5812 1.2788 1.5825 1.2919 1.5744 1.3003 1.5632 1.3044 1.5556 1.3049 1.5525

Our solution is further verified when the effects of aspect ratio are examined with a medium with the F2 phase function. Figure 4 shows the y-component of the net radiative heat flux along the y-direction centerline shown in Fig. 1. For an aspect ratio of unity, the energy losses through the side walls significantly affect the radiative flux distribution. The flux is large near the top wall a)

Incident angle (XxL-- ~yH= I) _

: 0 deg.

b)

Incident angle ( ~ = ZyH= I)

: 30 deg.

i c)

Incident angle (%L = ZyH= 5)

: 30

deg.

d)

Incident

w.

angle

: 60 deg.

Q

% Fig. 13. The net flux maps for different incident angles (p = 0, m = 1, a = 1, and F2).

236

TAE-KUK KIM and HAEOK SKARDALEE

(~* = 1), where the collimated energy is incident, and it is small near the bottom wall (T* = 0). As the aspect ratio is increased, the effect of the side wall energy loss is rapidly reduced at the center region of the enclosure. When the aspect ratio is increased to 10, the two-dimensional centerline analysis closely predicts the one-dimensional result. The two-dimensional result is a nearly constant value of about 0.9085, whereas the corresponding one-dimensional flux, obtained from our one-dimensional version of the discrete ordinates code, is 0.90688. The effect of the different scattering phase functions is the focus of our study. In Figs. 5 and 6, the reflected [Q~'*(¢*, 1)] and transmitted [Q~-*(¢*, 0)] components of the radiative flux along the top and bottom walls (walls 4 and 3) are presented for the purely scattering media. The isotropic scattering results are again compared with the results presented in Ref. 4. These figures show that the forward-scattering media transmit more radiative energy and reflect less energy back, than the backward-scattering media. The relative augmentation of the reflected or transmitted energy is a strong function of the phase function anisotropy. The reflected energy can differ by an order of magnitude when the phase function is changed from the F1 to B2. The transmitted energy changes by a factor of two when the two phase functions F1 and B2 are considered. The two oppositely shaped phase functions F3 and B2 reduce and augment the reflected component of the radiative fluxes by nearly the same amount relative to the isotropic-scattering medium. For the transmitted component, the F3 phase function results in a slightly larger amount of augmentation than the reduction caused by B2. Figure 7 shows the side wall energy losses [Q.~+*(1, T*) or - Q f*(0, ~*)] for the different phase functions. The side wall energy loss is also a strong function of the phase function anisotropy. The largest side wall loss is found for an isotropic scattering medium, since the sideways scattering is equally as strong as the scattering to the forward or backward direction (total loss through one side wall is equal to 0.1570). The scattering in either the forward or backward direction dominates as the absolute magnitude of the asymmetry factor gets larger, and the side wall loss is reduced to 0.0673 for F1 and 0.1392 for B2. The side wall loss for the F3 (total loss of 0.1547) is slightly larger than that for the B2. The location of the maximum energy loss gradually moves from the bottom-side (~* = 0) to the top-side (~* -- 1), as the asymmetry factor is changed from + I (forward scattering) to - 1 (backward scattering). The purely scattering results presented in Figs. 5-7 can be used to check for flux conservation. The sum of all the energies leaving the four boundary walls has to equal the incident energy for this non-absorbing case. The estimated errors in the radiative flux conservation here are of the order l 0 -7 (7.753 x 10 -7, 4.706 x 10 -7, 7.598 x 10 -s, 2.586 x 10 -7, 3.380 x 10 -7, and 4.296 x 10-7 for F1, F2, F3, ISOTROPIC, B1, and B2, respectively). Figure 8 and Table 2 show the y-component of the net radiative flux along the y-direction centerline of a square enclosure containing a purely scattering medium of unit optical thickness. The flux values along the enclosure edge are also given in Table 2. These results also show that energy transmission is a strong function of the phase function asymmetry factor. The phase functions F3 and B2, which are opposite in shape, produce about the same amount of augmentation and reduction in the fluxes compared to the isotropic flux, In Table 3, the centerline and edge distributions of the average incident radiation G* are presented for the different phase functions. The scattering albedo is also an important parameter to consider. In Fig. 9, the y-direction centerline net radiative flux is shown for a square enclosure of unit optical thickness with the phase function F2. Scattering albedos of 1.0, 0.5 and 0.1 result in significant changes in the net radiative flux distributions. The flux is conserved for the purely scattering case (co = 1), but the side-wall losses result in a nonuniform flux distribution along the centerline. For an absorbing-scattering medium (co < 1), some of the radiative energy is absorbed by the medium. As the result of absorption and side-wall losses, the net radiative flux distributions for absorbing-scattering media are considerably lower than for purely scattering media. Near the top wall (~* --- 1.0), however, the net flux for the smaller co is higher than for the larger co, since the back-scattered radiation is smaller. The optical thickness is another important parameter in radiative transfer. It is difficult to obtain fast convergence and good accuracy for large optical thicknesses. The number of grid points is increased with the optical thickness. For ZxL= TyX= 1, 26 X 26 grids are used in the x- and y-directions. For ¢xL= Zy8 = 2.5, 5 and 10, 46 point Lobatto quadrature points are used.

RT in 2-D anisotropic scattering media

237

In Fig. 10, some of the net flux results for the various optical thicknesses are presented for a purely scattering square enclosure with the phase function F2. The centerline distribution of the y component of the net flux shows a strong dependence on optical thickness. At an optical thickness of 1, more than 75% of the collimated energy entering the enclosure at ~* = 1.0 is transmitted. As the optical thickness is increased from 1 to 10, the percentage of the transmitted energy rapidly reduces to 15%. The reflected energy, although not shown, is low for small optical thicknesses ( < 5% for an optical thickness of 1) and is high for large optical thicknesses (about 60% for an optical thickness of 10). Different optical depths affect the amount of scattering, which then affects the reflected energy and the side-wall energy losses. The effect of the diffuse wall reflectivity on radiative transfer is examined for a purely scattering square enclosure with the phase function F2. As shown in Fig. 11, the net radiative flux is significantly affected by different surface reflectivities. For p = 0, no radiative energy is reflected back into the enclosure, and the net flux remains large. As p is increased from 0 to 0.9, the reflected energy is increased, and the net flux is significantly reduced. The angular distribution of the intensity becomes nearly isotropic with high p. In Fig. 12, a transmitted energy comparison of the collimated and diffuse incidence problems is made. For this comparison, equivalent radiative flux is incident on the same medium. The phase functions considered are F1, B2 and isotropic. The collimated incidence case is shown to transmit more radiative energy through the medium than the diffuse incidence. The effect of phase function anisotropy appears to be more significant for collimated incidence. The F1 transmits about 77% more energy than the isotropic medium for collimated incidence, while it transmits about 65% more energy than the isotropic medium for diffuse incidence. In Figs. 13(a~l), the three-dimensional wire frame graphs of Q* over z* and z* are shown for different incident angles of the collimated radiation beam. These graphs show the net radiative fluxes in square enclosures containing purely scattering media with the phase function F2. Comparing Figs. 13(a) (0 ° incidence), 13(b) (30 ° incidence), and 13(d) (60 ° incidence) for the same optical depth clearly shows regions of small flux for oblique incidence. The opaque wall at z* = 0 (wall 1) shades much of the two-dimensional region when the incidence angle increases to 60 °. The radiative flux in the shaded region is due to the scattered radiation from the illuminated region, and from the diffuse reflection at the bottom or the side wall. Comparing Figs. 13(b) (z = 1) and 13(c) (z = 5) shows the effect of optical thickness. Although the overall level of the flux is lower, the scattering through a larger optical depth smoothes the transmitted flux at "Cy* = 0. That is, the shaded region of the enclosure is less distinct. The negative Q*, shown near ~* = 0 and ~y* = 1 in Fig. 13(c), are due to the large back-scattered flux from the thick medium. The dip in the Qy* distribution at the top wall is also due to the significant reflected flux in that region. Acknowledgements--This work was supported, in part, by National Science Foundation Grant. No. NSF/CBT-8451076. A grant from the Minnesota Supercomputer Institute is also gratefullyacknowledged. REFERENCES S. Chandrasekhar, Radiative Transfer, pp. 149-150, Dover, New York, NY (1960). R. D. M. Garcia and C. E. Siewert, Transp. Theory Statist. Phys. 14, 437 (1985). C. Flicke, JQSRT 20, 429 (1978). A. L. Crosbie and R. G. Schrenker, JQSRT 33, 101 (1985). A. L. Crosbie and J. W. Koewing, JQSRT 21, 573 (1979). H. F. Nelson, D. C. Look, Jr., and A. L. Crosbie, J. Heat Transfer 108, 619 (1986). A. L. Crosbie and R. L. Dougherty, JQSRT 33, 487 (1985). T.-K. Kim and H. Lee, Int. J. Heat Mass Transfer 31, 1711 (1988). M. P. Mengiic and R. Viskanta, JQSRT 33, 533 (1985). K. D. Lathrop, J. Comput. Phys. 4, 475 (1968). K. D. Lathrop and F. W. Brinkley, "TWOTRAN-II: An Interfaced, Exportable Version of the TWOTRAN Code for Two-Dimensional Transport," Los Alamos Scientific Laboratory Report No. LA-4848-MS, Los Alamos, NM (1973). 12. W. M. Irvine, Bull. Astr. Inst. Neth. 17, 176 (1963). 13. M. N. 0zi~ik, Radiative Transfer, p. 83, Wiley-Interscience, New York, NY (1973). 14. R. Siegel and J. R. Howell, Thermal Radiation Heat Transfer, 2nd edn, p. 591, McGraw-Hill, New York, NY (1981). 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.

238

TAE-KuK K1M and HAEOK SI~RDA LEE

NOMENCLATURE a = Aspect ratio = L/H or Z~L/Zya CI = E x p a n s i o n coefficients o f p h a s e function = Energy re-balance factor G = Average incidence r a d i a t i o n o b t a i n e d from la H = Height o f the rectangular enclosure = Scattered c o m p o n e n t o f the m o n o c h r o m a t i c intensity /= = Actual m o n o c h r o m a t i c intensity, Eq. (3) K = O r d e r of phase function e x p a n s i o n L = L e n g t h o f the rectangular enclosure M = Total n u m b e r of o r d i n a t e directions n = I n w a r d n o r m a l vector to enclosure walls N = O r d e r o f S - N discrete o r d i n a t e s a p p r o x i m a t i o n p, --- Legendre p o l y n o m i a l s o f order l Q = Radiative flux vector o b t a i n e d f r o m la ' Q.~, Q.; = Positive a n d negative c o m p o n e n t s o f Q in x direction Q.~, a,. = Positive a n d negative c o m p o n e n t s o f Q in y direction Qx~ Or = N e t radiative fluxes in x a n d y directions S = Source function, Eq. (2) AV --- V o l u m e o f the control volume Wm = A n g u l a r weights Extinction coefficient, fl = x + as D i r a c delta function Scattering p h a s e function A b s o r p t i o n coefficient Direction cosine in x direction = cos 0 Direction cosine in y direction = sin 0 cos tp p = Reflectivity of the enclosure wall Optical coordinates, z~ = fix a n d ~y = fly "CxL, TyH Overall optical thicknesses, ZxL = flL a n d Zyn = flH Scattering angle Scattering coefficient O= Polar angle
Superscripts *= + = = '= -

Dimensionless variable Positive direction Negative direction Incident direction

Subscripts c F, B i, j rn n, s, e, a n d w p w

= = = = = = =

Collimated c o m p o n e n t F r o n t a n d b a c k hemispheres divided by the Spatial indices A n g u l a r index C o m p a s s directions C o n t r o l volume center p o i n t Wall

x - y s y m m e t r y plane