Radiative transfer in plane inhomogeneous media with exponentially varying albedo

Radiative transfer in plane inhomogeneous media with exponentially varying albedo

Ann. nuel.Energy,Vol. 12, No. 3, pp. 107-112, 1985 0306-4549/85 $3.00+0.00 Copyright © 1985Pergamon Press Ltd Printed in Great Britain. All rights r...

337KB Sizes 0 Downloads 32 Views

Ann. nuel.Energy,Vol. 12, No. 3, pp. 107-112, 1985

0306-4549/85 $3.00+0.00 Copyright © 1985Pergamon Press Ltd

Printed in Great Britain. All rights reserved

RADIATIVE T R A N S F E R IN PLANE I N H O M O G E N E O U S M E D I A WITH E X P O N E N T I A L L Y V A R Y I N G ALBEDO A. MAGNAVACCA,G. SPIGA and M. H. HAGGAG* Laboratorio di Ingegneria Nucleare dell'Universit~i di Bologna, Via dei Colli 16, 40136 Bologna, Italia

(Received 15 October 1984; accepted in revisedform 12 November 1984) A~traet--Accurate numerical results for the exit distributions and the global reflection and transmission coefficients relevant to radiative transfer in a stratified medium with exponentially varying albedo are obtained and compared to previous results. The semi-analytical solution of the linear transport equation is rigorously performed on the basis of a simple projectional method.

INTRODUCTION In this paper a rigorous solution of the benchmark problem of radiative transfer in finite atmospheres, recently investigated by Garcia and Siewert (1982) and Kelley (1983), is presented. The importance of exact benchmark solutions in transport theory, and their role in practical calculations for realistic configurations are well-known and well-recognized in the literature. The problem considered here is relevant to the analysis of radiative transfer with exponentially varying singlescattering albedo in a grey plane non-emitting medium with isotropic scattering and assigned radiation impinging on the boundaries. The angular radiation intensity I, as a function of the optical abscissa z ( - a, a) and the direction cosine It ~ [ - 1, 1],is governed by

Ol o9 It ~ + l(r, #) = ~ exp ( -- z/s)lo(z ) CO

(" 1

=sexp(-T/s)J_

Equation (1) may be formally solved upon the boundary conditions (2) to yield, for # > 0,

ll(z,It,)dff (1)

t(a,--it)=G(#),

xlo(z')dz'+exp(-~f-)F(it)

(3a)

and

It>0,

(2)

where F and G are the assigned incident distributions, co represents the albedo in the midplane z = 0 and s is a real non-zero parameter. The constant albedo case corresponds to the limit s ~ oo. In the present rigorous approach a Fredholm integral equation of the second kind is derived for I o, and its solution is performed

*CNR fellow; permanent address : Department of Physics, University of Mansoura, Mansoura, Egypt. 12:3-A

METHOD

l(z,it)=~ f~ exp(-z'/s)exp(-LS~-)

with boundary conditions

l(--a, it)=F(it),

analytically by a projection procedure with polynomial basis functions, which is particularly suitable for optically not too thick layers--as is the case for most layers of practical interest in technology (e.g. nuclear engineering). The only task to be accomplished numerically is the solution of a linear set of algebraic equations. Explicit results are reported and comparison is made with the other semi-analytical work of Garcia and Siewert (1982), where the problem has been treated in terms of the elementary solutions of Mullikin and Siewert (1980) and solved by the F N method.

I(z,--it)=~ffexp(--z'/s)exp(--5~5) xlo(z')dz'+exp(--a~--f)G(it).

For the only remaining unknown, Io, an integral equation is easily derived by a suitable integration of equations (3a, b) with respect to It. Setting K(~, z') =

107

(3b)

1 2

coexp(-z'/s)Ex(lz-z'D

(4a)

A. MAGNAVACCAet

108

al.

and

with matrix elements

S(r)=2rcf/exp(-~-)F(lx)dlx

C'"=[(2m+lX2n+l)]mf~2a a P'(~)

+ 2n f / e x p ( - ~ ) G ( # )

d/x, (4b)

where E. denotes the nth exponential integral function, one ends up with the Fredholm equation in canonical form :

co exp (a/lsl) = ~Oo ~< 1

(6)

(as it occurs in any physical situation), so that the spectral radius of the integral operator generated by K is less than unity, and we are guaranteed that a unique L 2 solution to equation (5) exists and, in addition, it is the limit in the Lz norm of the sequence of approximate solutions obtained by projectional methods. Once Io has been found, the angular intensity follows from equations (3a, b); the quantities of outstanding physical meaning are the exit distributions I ( - a , - # ) and l(a,/t), and the corresponding global parameters for the slab :

(10a)

and

4a

aran - -

#o(z)=f~oK(z,z')Io(z')dz'+S(z),(5) where it is easy to check that S eL2(-a,a) and K belongs to L2[( - a, a) x ( - a, a)]. To be on the safe side we will further assume that

dz

ka/

-a Pm

a

t

adz ,[.,

.[ - - ,It

×f_
(10b)

and inhomogeneous terms

(2m+

1~ 1/2

8m=\--~a)

fa S(z)pm(~)exp(z/s)dz. .I-o

(11)

The symmetric matrix C,.. can be obtained by the recursion formula / 2 n + 1 ~1/2

(m+l)t

)

["2n+ 1"~1/2

C.+,,.+mtT-U:-7) C._,..

( 2 m + l ~ '/2

=(n+l)\-~-~-.]

['2m+l'~ '/2 C,...+, + nt-~--~ ) Cm..-, (12a)

starting from the first column only, Cmo,which in turn is simply given by (Luke, 1969)

R = {f:,r,(,)+ G(,)1d,}-' f:,,(-<.. -,)d,

C,.o = ~1~ t ~

)

s'12I,.+1/2

,

(12b)

(7a) and

T={;:#[F(#x)+G(#x)]dla}-t;:lxl(a, lx)dv.

(7b)

where I v is a modified Bessel function of the first kind. Another diamond recurrence relation can be derived for Amn,namely

{2m+ 1'~ 1/2

2m+ 1"~ 1/2 Choosing as the expansion basis the orthogonal set P. of the Legendre polynomials, we write the Nth-order approximate solution as

L {2n+l'k 112 N {z'~

)

(2n+l~ '/z

=\~/

Am.n+l

{2n+lk 1/z

Am+l,.-t~)

A._,,.,

(8)

(13)

as a projection basis it proves convenient to take, according to the Petrov-Galerkin method (Krasnoselskii, 1972), the same set of polynomials multiplied by exp (T/s). In this way the linear algebraic system for the calculation of the N + I expansion coefficients ~u reads as

which allows the computation of all elements starting from the first two columns and the first row. Anyway, using a result of Kschwendt (1970), the general Am, can be cast in the form a

n-ra+ln-ml

A,.. = ~ [(2rn + 1) (2n + 1)] 1/2(_ 1)

z

N

E n=0

(Cmn - -

~°A,,,.)¢u. = B,., m = O, 1..... N,

(9)

x

~ v=O

" r ( - 1),,,+"D$-.t-D~-], Z Y.

(14)

Radiative transfer in finite atmospheres [--(a + z)//~0], the inhomogeneous terms read as

where

D~ = f~ x"El(ax ) exp ( _+ax/s) dx

(15)

--

a,

co exp(--a//~) ~ (_1).(2n+1)1/2 4x/-~P .= o

-- #)

N

s#

×~n('ff-~)

1/2

B,, = 2• 3/2 exp (-a/#o)(-- 1)m(2m+ 1) 1/2

×( SIlO ~l/2lm+l/2(aS--#°~,

and fl~" is the Kschwendt coefficient, as defined by Cupini et al. (1969). For each approximation of order N the approximate values of the sought quantities are analytically expressed by

Is(

s+#

I,+xlz(aT)+G(I.z)exp(--2allO

k s - #o /

\

M{2m+ 1"~1/2 B. = 2~a exp(--a/s)(--1) t~-a ) H+.

Is(a, #) = - -co 4.4~-

= _1 o9

1).(2n+ 1~ 1/2 ~ffH~- (21) \~,/

~,

and

exp(--a/#) ~ (2n+l)X/2

s / st* \~/2

/~

lo9

.=o

N (2n+1~1/2

TN = ~ ~ - a exp (-- a/s).~o \

[ s-#'~

x~. t~_~ ) ,.+l/2ta~)+FO, Oexp(-2al#x) (16b) and, for the global quantities, o9 Izlt~(-T-a,-T-p)d# = ~aexp( +_a/s) ~N

(20)

Once the vector {~} has been computed for N increasing up to the desired accuracy, the exit distributions are provided by equations (16a,b), with the appropriate specializations, and the slab reflectivity and transmissivity are given, respectively, by /.2 (-RN ct ~ a exp (a/s).=o

and

(19)

Spo /

whilst in the case ofisotropicincident radiation we have S(z) = 27rE2(a + ~) and thus

(16a)

fO

109

(-T-l)"

~ H + + fl; (22)

with a = t% and fl = exp ( - 2a/Po) in the first case, and = 1/2 and fl = 2E3(2a ) in the second. All quantities needed for the numerical procedure are now in closed analytical form, except D.~ and H.~. However, on integrating by parts equation (15) yields, after some manipulation,

/s'"+tF/

L'/"+''

°:

n=O

2a )

2a)

]

"~1,2 fo/t exp (-- 2a/#) F(p) G(p)d#' ~ H~ +

(23)

x(2n+ 1 \ 2a ]

where (17) Q~= E,(~-2a)-E,(2a'+ln

with

H.~ = f~ P.(1--x)E2(ax)exp( +ax/s)dx.

(18)

& (-T-l)"

(-Vl)RESULTS

The above theory has been processed numerically on the VAX 11/780 computer of the University of Bologna for the case where G(#) = 0, so that the parameters R and T actually represent the reflection and transmission coefficients, respectively, and for two different options for F(#), namely F(#) = 3(#--#o) where 6 is the Dirac delta function and 0
1

,,=12" m, ~ =,.=.+1

(

( ~ --~) sT-1

"~

Ttm'----~ 2a) 1

[

sTl~'~

m! ~ T t m ' ~ - z a )

(24a)

and y denotes an incomplete F function. When s = + 1 (s = - 1) we must use for Q+ (for Q~-) Q.~ = - El(Za)- c - In (2a) ~, (-T-1)m (2a)" -Y- 7 1 mf m =_+

~ m=n+l

(T- 1)m (2a)= --, m! m

s = +1,

(24b)

where c is the Euler constant. Of the two alternative

110

A. MAGNAVACCA

expressions given in equations (24a,b), the infinite series has to be preferred to the closed form (involving a finite sum) in all cases when the latter gives rise to nonnegligible cancellations of significant figures by difference. As regards H,~, we have initially u. • =

(-1) k

Some of our numerical results are given in Tables 1-7 for different values of the input parameters a, s, coo and #o (including s = ~ , which has been run separately, with greatly simplified expressions) and the variable #. Negative values of s correspond to an exponentially increasing albedo in the direction of propagation ofthe impinging radiation. The convergence rate is quite satisfactory; as expected, it improves substantially when Isl increases or a decreases, and is excellent for optically very thin layers (2a << 1). Within approximation orders of about 20, six stable significant figures are achieved for the functions I(+_ a, +_t~) in all the cases considered with 2a < 10; if 2a increases beyond this limit, the number of converged digits typically

(n+k)!

k=O k ~ 2k(n-k)!k! U~

(25a)

and then by recursion

U~ = j~2 xkE2(ax) exp ( +

=

ax/s)dx

s-T-1 ± (1 s y+lT(k+l,~s__2a)_aDk+l.

Table 1. l(-a, p

s = ~

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.129401E+ 0.134259E+ 0.136758E + 0.137898E+ 0.138141E+ 0.137746E + 0.136880E + 0.135656E + 0.134158E + 0.132452E +

et al.

(25b)

-,u) for negative values of s, when e) o = s = 1 1 1 1 1 1 1 1 1 1

10 3

0.126465E + 0.130996E + 0.133248E + 0.134198E+ 0.134296E+ 0.133793E + 0.132849E+ 0.131571E+ 0.130041E + 0.128320E +

s = 1 1 1 1 1 1 1 1 1 1

1, 2a = 5 and F(#) = 6 ~ - 0 . 9 )

10 2

0.105912E+ 0.108313E+ 0.109000E + 0.108773E + 0.107991E + 0.I06846E + 0.105454E + 0.103891E+ 0.102208E+ 0.100444E +

1 1 1 1 1 1 1 1 1 1

s = -- 10

s = -- 1

0.412082 0.399804 0.385648 0.371383 0.357612 0.344549 0.332246 0.320689 0.309835 0.299633

0.334714E- 2 0.331464E- 2 0.328337E- 2 0.325377E- 2 0.322633E- 2 0.320126E--2 0.317831E--2 0.315682E--2 0.313597E--2 0.311501E--2

Table 2. l ( a , # ) for positive values of s, when 090 = 1, 2a = 5 and F(/0 = 6 ( # - 0 . 9 ) ,u 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

s=

1

s=

0.212259E 4 0.292138E-4 0.470821E 4 0.102175E-3 0.268549E- 3 0.650347E - 3 0.132760E-2 0.232491E - 2 0.361612E - 2 0.514440E- 2

10

0.144975E0.173068E0.207231E0.250539E0.305785E0.374222E 0.454844E0.544866E0.640732E0.738937E-

1 1 1 1 1 1 1 1 1 1

s = 102

s = 103

0.147725 0.172372 0.196559 0.220967 0,245791 0.270889 0.295861 0,320191 0,343389 0,365066

0.226076 0.262914 0.297875 0.331902 0.365250 0.397794 0.429157 0.458866 0.486492 0.511727

Table 3. l(a,,u) for negative values of s, when e~o = 1, 2a = 5 and F(#) = 6(#--0.9) /t

S = O0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.238416 0.277156 0.313779 0.349269 0.383896 0.417544 0.449841 0.480324 0.508578 0.534308

S=

10 3

0.228063 0.265067 0.300036 0.333929 0.367018 0.399203 0.430133 0.459364 0.486492 0.511227

S = --10 2

S = --10

S = --1

0.161083 0.186868 0.211165 0.234740 0.257874 0.280568 0.302601 0.323648 0.343389 0.361575

0.315909E-- I 0.359690E- 1 0.399477E- 1 0.437909E - 1 0.476653E- 1 0.516580E- 1 0.557702E-- 1 0.599391E- 1 0.640732E- 1 0.680812E- 1

0.412604E- 2 0.411132E-2 0.404153E- 2

0.395935E- 2 0.387867E - 2 0.380393E- 2 0.373585E- 2 0.367367E- 2 0.361612E-2 0.356198E- 2

111

Radiative transfer in finite atmospheres

Table 4.

l(-a, -#)

for positive values of s, when to o = I, 2a = 5 and F(/0 = 1

#

s = 1

s = 10

s = 102

s = 103

s =00

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.531121 0.443281 0.380307 0.332965 0.296091 0.266564 0.242390 0.222235 0.205174 0.190547

0.733977 0.686320 0.644181 0.606474 0.572568 0.541972 0.514267 0.489096 0.466153 0.445170

0.850283 0.825231 0.801810 0.779464 0.757995 0.737304 0.717335 0.698051 0.679431 0.661459

0.883186 0.864099 0.846057 0.828567 0.811431 0.794552 0.777879 0.761394 0.745105 0.729035

0.887836 0.869581 0.852~99 0.835503 0.818996 0.802676 0.786493 0.770429 0.754496 0.738721

Table 5. The reflection coefficient R for selected values of 090, s and 2a, when F(#) = 6(/*- 0.9) tOo

s

2a = 10 -2

2a = 10 -1

2a = 0.5

2a = 1

2a = 5

0.7

1 10 102 103 oo -103 - 102 -10 - 1

0.381553E-2 0.383288E-2 0.383462E-2 0.383480E--2 0.383482E--2 0.383480E--2 0.383462E--2 0.383284E--2 0.381514E--2

0.330938E-1 0.346827E- 1 0.348478E- 1 0.348644E-1 0.348662E--1 0.348641E--1 0.348453E- 1 0.346478E-1 0.328571E- 1

0.966545E-1 0.121451 0.124583 0.124904 0.124940 0.124886 0.124402 0.119699 0.838771E- 1

0.115005 0.168417 0.177022 0.177937 0.178039 0.177841 0.176069 0.159655 0.703596E- 1

0.118853 0.196979 0.216745 0.219185 0.219464 0.217494 0.200917 0.102994 0.123364E-2

1.0

1 10 102 103 oo - 103 -102 -10 1

0.549702E-2 0.552224E-2 0.552477E-2 0.552503E-2 0.552506E-2 0.552503E-2 0.552477E-2 0.552219E-2 0.549646E - 2

0.498765E-1 0.524177E-1 0.526825E-1 0.527091E-1 0.527121E--1 0.527087E - 1 0.526787E-1 0.523797E-1 0.495168E - 1

0.160287 0.212226 0.219254 0.219982 0.220063 0.219949 0.218920 0.209020 0.138379

0.197101 0.334012 0.361677 0.364743 0.365087 0.364523 0.359508 0.314834 0.117145

0.205174 0.466153 0.679431 0.745105 0.754496 0.732695 0.584276 0.187838 0.177247E-2

Table 6. The transmission coefficient T for selected values of too, s and 2a, when F(#) = 3(/~- 0.9) tOo

s 1 10 102 /103

0.7 | |

- 103 - 102

L -10

1 10 102 103 1.0 ]--103 ]--102 [--_1~

2a = 10 -2

2a = 10 -1

2a = 0.5

2a = 1

2a = 5

0.992765 0.992783 0.992785 0.992785 0.992785 0.992785 0.992785 0.992783 0.992765

0.927586 0.929322 0.929503 0.929521 0.929523 0.929521 0.929504 0.929334 0.927697

0.654983 0.685146 0.689071 0.689475 0.689520 0.689481 0.689125 0.685668 0.658773

0.394446 0.457548 0.468911 0.470135 0.470272 0.470155 0.469108 0.459352 0.403411

0.444660E - 2 0.101850E- 1 0.169413E- 1 0.182531E- 1 0.184121E- 1 0.182769E- 1 0.171560E- 1 0A11466E-- 1 0.504880E -- 2

0.994446 0.994472 0.994475 0.994475 0.994475 0.994475 0.994475 0.994472 0.994447

0.944201 0.946967 0.947256 0.947285 0.947288 0.947285 0.947257 0.946985 0.944370

0.709104 0.770418 0.778951 0.779838 0.779937 0.779848 0.779052 0.771380 0.715644

0.442193 0.596084 0.630603 0.634477 0.634913 0.634524 0.631068 0.600175 0.458570

0.486873E-- 2 0.291545E--1 0.163099 0.234432 0.245504 0.234955 0.166668 0.341662E--1 0.596036E -- 2

112

A. MAGNAVACCA et al. Table 7. The reflection coefficient R for selected values of tOo, s and 2a, when F(#) = 1 tOo

s 1 10 102 103

0.7 1-103 [-102 1--10 -

1 10 102 /103 1.0 ]-103 | - 102 L-10

2a = 10 -2

2a = 10 -1

2a = 0.5

2a = 1.0

2a = 5.0

0.671647E--2 0.674682E- 2 0.674987E - 2 0.675017E-2 0.675020E-2 0.675017E-2 0.674985E-2 0.674671E- 2 0.671534E-2

0.530284E-1 0.554877E-- 1 0.557431E- 1 0.557687E-1 0.557716E-I 0.557681E-1 0.557372E-1 0.554289E- 1 0.524694E-1

0.133802 0.164641 0.168507 0.168904 0.168948 0.168872 0.168188 0.161556 0.111241

0.152071 0.211600 0.220966 0.221959 0.222070 0.221815 0.219535 0.198425 0.841249E-1

0.155240 0.235414 0.254017 0.256264 0.256519 0.254299 0.235572 0.122769 0.137761E-2

0.967639E- 2 0.972050E- 2 0.972493E- 2 0.972537E- 2 0.972542E-2 0.972537E-2 0.972491E - 2 0.972034E- 2 0.967477E- 2

0.799031E- 1 0.838411E- 1 0.842514E- 1 0.842925E- 1 0.842971E- 1 0.842916E- 1 0.842423E - 1 0.837514E- 1 0.790539E- 1

0.221113 0.286097 0.294827 0.295730 0.295831 0.295672 0.294241 0.280479 0.182584

0.258891 0.412506 0.442863 0.446217 0.446594 0.445893 0.439658 0.384099 0.138235

0.265892 0.531182 0.725972 0.784073 0.792343 0.771288 0.626532 0.216804 0.197547E- 2

d e c r e a s e s , d e p e n d i n g also o n t h e v a l u e s o f s, COoand/~o. T h u s , we believe t h a t t h e r e s u l t s s h o w n in T a b l e s 1 - 7 a r e a c c u r a t e to _ 1 in t h e last figure. O u r r e s u l t s a r e also in perfect a g r e e m e n t w i t h t h o s e o f G a r c i a a n d Siewert (1982), w h e r e v e r c o m p a r i s o n is possible. T h e i n f l u e n c e o f s a n d o f different d i s t r i b u t i o n s o f t h e i m p i n g i n g r a d i a t i o n o n t h e exit d i s t r i b u t i o n c a n be s e e n f r o m T a b l e s 1-4. A q u a n t i t a t i v e d e s c r i p t i o n o f t h e v a r i a t i o n s o f t h e reflection a n d t r a n s m i s s i o n coefficients for v a r y i n g s, coo a n d a a r e g i v e n in T a b l e s 5-7. I n a n a l y s i n g different v a l u e s o f a for a g i v e n coo a n d for a fixed n e g a t i v e s, o n e m u s t b e a r in m i n d t h a t t h e a l b e d o o f t h e i l l u m i n a t e d s u r f a c e (at t h e inlet) is e q u a l to coo exp(2a/s), w h i c h is n o t c o n s t a n t a n d d e c r e a s e s r a p i d l y w h e n a increases. T h i s e x p l a i n s , for i n s t a n c e , t h e s u d d e n d e c r e a s e o f R in T a b l e s 5 a n d 7, for s = - 10 a n d s = - 1, w h e n 2 a i n c r e a s e s f r o m 1 to 5 m.f.p. T h i s effect is a b s e n t for p o s i t i v e s s i n c e coo itself t h e n r e p r e s e n t s t h e a l b e d o at t h e inlet.

Acknowledgements--We are grateful to Dr P. Vestrucci for several helpful suggestions and enlightening discussions in relation to this paper. M.H.H. would like to thank the Nuclear Engineering Laboratory of the University of Bologna for the creative atmosphere and research opportunities provided, and for the kind hospitality extended to him while this paper was being written. This work was supported by the CNR. REFERENCES

Cupini E., De Matteis A., P r e m u d a F. and Trombetti T. (1969) C N E N Report RT/FI-(69)-45. Garcia R. D. M. and Siewert C. E. (1982) J. quant. Spectrosc. Radiat. Transfer 27, 141. Kelley C. T. (1983) Trans. theor, statist. Phys. 12, 183. Krasnoselskii M. A. (1972) Approximate Solution of Operator Equations. Noordhoff, Groningen, The Netherlands. Kschwendt K. (1970) Atomkernenergie 15, 122. Luke Y. L. (1969) The Special Functions and Their Approximation. Academic Press, New York. Mullikin T. W. and Siewert C. E. (1980) Ann. nucl. Energy 7, 205.