Neutron transport problems in anisotropic media

Annals of Nuclear Energy 28 (2001) 1271±1286 www.elsevier.com/locate/anucene

Neutron transport problems in anisotropic media G. Lapenta a, P. Ravetto a,*, M.M. Rostagno a, R. Jacqmin b, F. Malvagi b, G. Rimpault b a

Politecnico di Torino, Dipartimento di Energetica, Corso Duca degli Abruzzi, 24- 10129 Torino, Italy b CEA, CEN-Cadarache-13108 Saint-Paul-lez-Durance, France Received 30 October 2000; accepted 14 November 2000

Abstract In nuclear reactor physics, neutron transport in media having anisotropic properties is of interest in modelling largely voided con®gurations, such as boiling-water and gas-cooled reactors. In order to treat such problems, neutronic methods need to be extended and adapted. The present work is devoted to the solution of transport problems for materials in which cross-sections are dependent on the direction of motion of the neutrons. A simple onedimensional plane con®guration is assumed, in order to be able to obtain analytical solutions and to gain physical insight into the transport phenomenon. Results can thus be considered as benchmarks for numerical algorithms employed in realistic cases. # 2001 Elsevier Science Ltd. All rights reserved.

1. Introduction The problem of anisotropic media is well established in many transport problems, especially in the radiative transfer ®eld. Recent developments can be found in the technical literature for application to plant canopies (Ganapol, 1991), and to the transport of radionuclides in fractured rocks (Cassell and Williams, 1992). The analytical solution of the transport equation for anisotropic media was proposed for in®nite con®gurations as an extension of the classical Case method (Williams, 1978). The phenomenon is not present in neutron transport, except in some applications considering very low energy neutrons. However, anisotropy in the total cross-section comes inevitably * Corresponding author. E-mail address: [email protected] (P. Ravetto). 0306-4549/01/$ - see front matter # 2001 Elsevier Science Ltd. All rights reserved. PII: S0306-4549(00)00126-2

1272

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

into play when constructing multigroup equations starting from the continuous energy version of the Boltzmann transport model and it is usually eliminated in some rather arti®cial manner (Bell and Glasstone, 1970). In fact, starting from the usual energy-dependent neutron transport equation, written using a conventional notation for an isotropic one-dimensional plane medium and for isotropic emissions: … @'…x; E; † 1  ‡  …x; E†'…x; E; † ˆ dE0 tr …x; E0 ! E†…x; E0 † ‡ S…x:E; †; …1† @x 2 the multigroup model is obtained by integration over a ®xed energy interval identi®ed by subscript g, ! „ @'g …x; † g dE …x; E†'…x; E; † ‡  'g …x; u† @x 'g …x; † ! „ …2† … X g0 dE0 tr …x; E0 ! E†…x; E0 † 1 dE ˆ g0 …x† ‡ Sg …x; †; 2 g g0 …x† g0 where: … 'g …x; † ˆ … Sg …x; † ˆ

… g

g

dE' …x; E; †;

…

g …x† ˆ d dE' …x; E; †; g

…3†

dES…x; E; †:

Group-average cross-sections are then de®ned accordingly, as: „ g dE …x; E†'…x; E; † ; g …x; † ˆ 'g …x; † „ „ 0 0 0 g dE g0 dE tr …x; E ! E†…x; E † ; g0 ! g …x† ˆ g0 …x†

…4†

and hence, obviously, anisotropy in the total cross-section comes into play for the set of multigroup integro-di€erential transport equations constituting a fully coupled system, while scattering properties turn out to be isotropic: 

@'g …x; † 1X ‡ g …x; †'g ˆ g0 ! g …x†g0 …x† ‡ Sg …x; †: @x 2 g0

…5†

Furthermore, recently, the problem of neutron transport in anisotropic media has received new attention in connection with neutronic safety studies of water reactors and with the design of gas-cooled systems. In situations presenting large voided regions, as the streaming in one direction is dominating with respect to other directions, the average properties of the homogenized material should physically account

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

1273

for such macroscopic anisotropy (Rimpault et al., 1999). For instance, in these cases it is suggested that cell calculations produce anisotropic average cross-sections, in the most general case also for scattering properties. It is therefore of interest to revisit previous work in the ®eld in order to make it suitable to current needs. Transport codes are being adapted to include anisotropic cross sections. An important aspect of code development is the validation of algorithms by analytical benchmarks. For that purpose, the present work is devoted to the development of reliable analytical solutions of the transport problem in slab geometry for the purpose of serving as benchmarks for numerical methods. For this reason spherical harmonics and discrete ordinate models are taken into consideration in this work. 2. Mathematical model and solution for the in®nite medium 2.1. General theory As a ®rst step, referring to Eq. (5) above, the following one-velocity transport problem is considered: 

@'…x; † 1 ‡  …x; †'…x; † ˆ s …x†…x† ‡ S…x; †: @x 2

…6†

As from Eq. (4), the scattering cross-section is assumed to be isotropic, and hence anisotropy in the total cross-section is due to absorption only. This is also the case when in voided cell calculations anisotropic leakage is simulated by a capture crosssection modi®cation. However, anisotropy in the scattering proprieties can be also handled rather easily and it is done in the latest part of the work. The objective of the present paper is the determination of the solution of the transport equation and of the corresponding spherical harmonics and discrete-ordinate approximate models for a materially homogeneous system. A fully analytical procedure is to be employed, so that the results obtained can constitute a reliable benchmark for numerical techniques. The problem (6) can be handled rigorously in the in®nite medium by following the previous work by Williams (1978). Assuming an isotropic external source having strength S0 , and taking the Fourier transform of Eq. (6), one obtains for the transform of the total ¯ux ~ the following expression: S0 I…k† ~…k† ˆ ; 2 1  s I… k † 2

…7†

where I…k† is de®ned as: I…k† ˆ

…1 1

d

1 : ik ‡  …†

…8†

1274

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

The explicit solution can be obtained by taking the inverse Fourier transform of Eq. (7). Following the standard technique which resorts to a line integration in the complex plane performed by the residue theorem (Case et al., 1953; Case and Zweifel, 1976), the full solution can be written in a form which separates the contributions of the residues and of the branch points of ~…k†, as: …x† ˆ asy …x† ‡ tr …x†:

…9†

In the above formula the asymptotic term dominant at large distances from the source is explicitly:   asy …x† ˆ 2iRes ~…k†; ik0 e 0 jxj ; …10† where i0 is one of the two imaginary conjugate poles of ~…k†, i.e. the imaginary solution of the transcendental equation: 1

s I…k† ˆ 0: 2

…11†

The transient portion of the solution is due to the branch cut of ~…k†, and it can be written down as: …1 tr …x† ˆ dA…†e kjxj ; …12† t

where A … † ˆ

1 h~ dx …† 2i

i ~ sx …† ; 

…13†

~dx …† and ~sx …k† are the values obtained by the limit calculation for the function (7) as k ! ik from either side of the branch cut and the essential singularity is located at ikt ; kt being given by (Williams, 1978):    …† ; withjj41: …14† t ˆ min  Since t > 0 , it is evident that the transient part of the solution gives a contribution that can be relevant near the source but fades away by far more rapidly than the pole contribution as the distance from the source increases. 2.2. Speci®c applications For full calculations two explicit forms of the angular dependence of the total cross-section are assumed, either an elliptic form, as (Fig 1):

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

1275

Fig. 1. Graphs of  …† (a) and  …†= (b) using an elliptic form of the cross-section, Eq. (15), for di€erent values of the anisotropy coecient a.



2 1 2  … † ˆ 2 ‡ A T2



1=2

…15†

or a parabolic form as (Fig 2): T †2 ‡ T ;

 …† ˆ …A

…16†

to interpolate between a given longitudinal value A and a transverse value T of the total cross-section. Both expressions show an even parity with respect to the direction cosine . Such a symmetry is often supposed in the study of voided systems. In the following calculations it is assumed that T ˆ 1. Therefore, distances are measured in terms of the transverse mean free path (t.m.f.p.). An anisotropy coecient can be introduced, in the ®rst case as: q a ˆ 1 …T =A †2 ; …17† which can take imaginary values, and in the second case as: b ˆ A =T

1

…18†

Assuming the cross-section takes on the form (15), the above integral (8) can be calculated as: p! p! 1 1 2a2 ‡ i 1 a2 1 1 ‡ ik 1 a2 p ‡ p ; log log I…k† ˆ …19† 2i 2ik 1 2a2 i 1 a2 1 ik 1 a2

1276

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

Fig. 2. Logarithmic scale graphs of  …† (a) and  …†= (b) using a parabolic form of the cross-section, Eq. (16), for di€erent values of the anisotropy coecient b.

where ˆ becomes:

p 4a2 ‡ k2 , while, assuming the cross-section takes on the form (16), it

1 1 I…k† ˆ p log 2 i 4b ‡ k 1

p! b ‡ i 4b ‡ k2 p : b i 4b ‡ k2

In both cases it can be veri®ed that in the limit of isotropy …A ˆ T †:   1 1 ‡ ik log ; for a; b ! 0; I…k† ! ik 1 ik

…20†

…21†

and hence the classical result is obtained. Furthermore, from Eq. (14) the branch point is located on the imaginary axis: ( p p 1= 1 a2 at  ˆ 1  if a41= 2  p p t ˆ …22† 2a at  ˆ 1= a 2 ; if 1= 2 < a < 1 using Eq. (15) while,  1p ‡b at  ˆ 1 p if 14b < 1 t ˆ 2 b at  ˆ 1= b if b > 1 using Eq. (16). Finally, in the ®rst case it is possible to write:

…23†

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

  (h S0 1 1 ‡ A … † ˆ 1 2 2 2

 ) 1 i2   1 s 1 2 s ‡ I … † ‡ ; 2 2 2 2

1277

…24†

and in the second case: ( h S0 1 A…† ˆ p 1 2 2 4b

2 ) 1 i2  s s I…† ‡ p : 2 2 2 4b

…25†

These formulae [Eqs. (22)±(25)] constitute a special case of the results obtained by Williams (1978) from which they can be directly worked out; from them tr …x† can be determined making use of Eq. (12). 3. Analytical solution for the approximate models In applications to real systems the transport problem is treated by approximate models, such as the expansion in terms of spherical harmonics or the angular discretisation through the discrete ordinate technique. It is then worthwhile to develop analytical solutions for such models. In the following some considerations on the spherical harmonics method and the di€usion approximation and full analytical solution for the discrete ordinate model are presented. 3.1. Spherical harmonics model For plane geometry, the unknown angular ¯ux is expressed using an expansion in terms of the Legendre polynomials, namely: '…x; † ˆ

N X 2n ‡ 1

2

nˆ0

Pn …†n …x†;

…26†

where the moments n …x† of the expansion constitute the new unknowns of the problem. The equations for the unknown moments are obtained substituting Eq. (26) into Eq. (6) and using the orthogonality properties of the Legendre polynomials to obtain: N X n‡1 0 n n‡1 …x† ‡ 0n 1 …x† ‡ I nm m …x† ˆ 0n ‰S0 …x† ‡ s 0 …x†Š; 2n ‡ 1 2n ‡ 1 nˆ0

…27†

where I nm ˆ

2n ‡ 1 2

…1 1

dPn …†Pm …† …†:

…28†

1278

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

Due to the even parity of the cross-section, I nm ˆ 0 for all odd values of …n ‡ m† However, unlike usual PN equations where the coupling introduced by the streaming term involves only consecutive moments, a full coupling appears in Eq. (27). A consistent di€usion model can be obtained from the P1 approximation: 8 d > < 1 ‡ I 00 0 ˆ S0 …x† ‡ s 0 …x† dx …29† > : 2 d 0 ‡ I 11 1 ˆ 0: 3 dx The equation for the ®rst moment provides a P1 de®nition for the di€usion coecient, assuming the form of the cross-section (15) it becomes: D … a† ˆ

2 a3 p ; 9 arcsin…a† a 1 a2

…30†

and with the form (16), it is: D … b† ˆ

5 : 3…5 ‡ 3b†

…31†

In the limit of an angular independent cross-section …a; b ! 0†; the above di€usion coecients reduce to the usual di€usion coecient for isotropic media and isotropic emissions. 3.2. Discrete ordinate model Eq. (6) is written for a set of discrete directions: n

N X d 1 1 'n …x† ‡ n 'n …x† ˆ s …x† wn 'n …x† ‡ S0 …x†; dx 2 2 nˆ1

…32†

where 'n …x† ˆ '…x; n †; n ˆ  …n † and wn are the weights of the quadrature set. The formulation of the SN method is conveniently expressed in a matrix form:  d  ' ˆ A ' ‡ jSi; dx

…33†

 
T 
T 1 where vectors ' ˆ '1 …x†; . . . ; 'N …x† ; jSi ˆ S0 …x† 1=1 . . .; 1=N are intro2 duced. The linear system of di€erential equations above is solved using the eigenstate method which was previously employed in analytical discrete ordinate calculations (Coppa and Ravetto, 1980). Previous applications to anisotropic media were proposed by Lapenta et al. (2000). The solution is found as a superposition of eigen-

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

1279

vectors juk i of A and its corresponding eigenvalues k determine the spatial behaviour of the angular ¯uxes. In the present case, the matrix is not symmetric and the basis of eigenvectors must be supplemented with the adjoint eigenvectors huk j. It is worth mentioning that the fundamental eigenvalue is connected with the spatially most persisting component of the solution, and it is usually referred to as the inverse of the di€usion length. The other eigenvalues come into play in simulating the continuous portion of the spectrum of the transport operator. The unknown vector is expressed as a superposition of eigenvectors: N  X ' ˆ ck …x†juk i;

…34†

kˆ1

where the coecients ck …x† are the new unknowns. Substituting the expansion in system (33) and projecting upon the basis of adjoint eigenvectors huh j, the equations for the unknown coecients follow directly: d ch …x† ˆ h ch …x† ‡ sh …x†; dx

…35†

where sh …x† ˆ huh jSi: These equations are uncoupled and can be solved easily. The general solution reads: …x 0 h x ch …x† ˆ ch0 e ‡ dx0 sh …x0 †eh …x x † : …36† 0

The coecients ch0 are determined solving a complete algebraic system determined by the boundary conditions imposed on the angular ¯uxes given by Eq. (34). A full analytical solution can be also obtained for the most general case when a known angular dependence of the scattering cross-section is assumed, and thus s …x† is substituted by s …x; †: In this situation Eq. (32) is simply substituted by: n

N d 1X 1 'n …x† ‡ n 'n …x† ˆ s;n …x†wn 'n …x† ‡ S0 …x†; dx 2 nˆ1 2

…37†

where s;n …x† ˆ s …x; n †: The structure of the solution procedure is exactly the same as the one presented above. Obviously, the characteristic matrix A of the problem is changed and, consequently, also eigenvalues and eigenvectors are changed. 4. Results In this section, some results are presented and discussed. Firstly, Tables 1 and 2 report the values of the eigenvalues for di€erent orders of the discrete ordinate approximation and for two values of the scattering cross-section. Fig. 3 represents pictorially the localization of the poles for the same material con®guration considered in

1280

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

Table 1 Eigenvalue spectrum fh g for di€erent orders of the discrete-ordinate model and for s ˆ 0:5 T =A

N

1.5

2 4 8 16

0.92595 0.65700 0.62517 0.62306

2.24089 0.83920 0.69694

1.47030 0.78769

4.83742 0.95642

1.25916

1.84133

3.21846

9.96540

2 4 8 16

1.22474 0.98876 0.95949 0.95753

2.44271 1.16090 1.02879

1.73824 1.11432

4.95218 1.27094

1.54867

2.08432

3.38497

10.02535

2 4 8 16

1.42415 1.34476 1.35860 1.35834

2.55610 1.52315 1.50428

1.91974 1.51107

5.00762 1.57238

1.76650

2.23104

3.46952

10.05296

1

0.5

Table 2 Eigenvalue spectrum for s ˆ 0:9 T =A

N

1

2 4 8 16

0.54772 0.52558 0.52548 0.52552

2.05519 1.10894 1.02069

1.61565 1.09350

4.55497 1.23671

1.49808

2.00833

3.25380

9.62602

2 4 8 16

0.83482 0.91742 0.91617 0.91621

2.21826 1.51883 1.50401

1.83262 1.51072

4.61923 1.56270

1.73543

2.16737

3.34480

9.65561

2 4 8 16

1.07943 1.46832 1.38901 1.38759

2.51726 2.11114 2.01456

3.29612 2.18914

4.68982 2.33929

2.95444

3.43081

5.00651

9.68067

0.5

0.1

the previous tables. From the results presented, it can be seen how higher order poles simulate the continuous portion of the transport spectrum. The higher order poles or the continuum determine the transport-like portion of the solution, while the fundamental pole is associated with the spatially most persistent portion of the solution, which is often referred to as asymptotic or di€usive contribution. Note also that for higher values of the scattering cross-section, the separation between the fundamental pole and the higher order ones (or the continuum, in exact transport) increases. The physical consequence is that the contribution of the continuum is less relevant with respect to the fundamental pole, as well known in conventional transport theory.

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

1281

Fig. 3. Spectra for di€erent discrete ordinate approximations and for exact transport. Crosses (+) indicate the fundamental pole, while circles (*) indicate higher order poles. Di€erent orders of the discrete ordinate approximation and exact transport (TR) are shown. Anisotropy amounts to T =A ˆ 0:7.

Fig. 4. Fundamental pole for some values of the anisotropy coecient a ˆ values of the scattering cross-sections.

q 1 …T A †2 and di€erent

1282

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

Fig. 5. Ratio of asymptotic to exact solution as a function of the distance from the source for s ˆ 0:7, using an elliptic form of the total cross-section.

Fig. 6. Ratio of asymptotic to exact solution as a function of the distance from the source for s ˆ 0:7, using a parabolic form of the total cross-section.

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

1283

Fig. 7. Exact solution of the discrete ordinate model (S16), using an elliptic form for the total cross-section.

Fig. 8. Exact solution of the discrete ordinate model (S16), using a parabolic form for the total cross-section.

1284

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

Fig. 9. Exact solution of the discrete ordinate model (S16), using an elliptic form for the total and the scattering cross-sections, having the ratio s = ˆ 0:7.

Fig. 10. Exact solution of the discrete ordinate model (S16), using an elliptic form for the total and the scattering cross-sections, having the ratio s = ˆ 0:9.

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

1285

Fig. 4 shows the dependence of the exact transport fundamental pole  0 upon the anisotropy coecient a for the elliptic form of the total cross section and for di€erent values of s . Figs. 5 and 6 report the ratios of the asymptotic to the full solution for an in®nite system injected by a localised plane source at x ˆ 0 for the two forms of the total crosssection. The contribution of the asymptotic portion of the solution is larger for smaller values of T =A . This is due to the fact that neutrons travelling around the transverse direction are less absorbed and thus have a higher weight in the neutron balance. As a consequence, such neutrons undergo scattering and become di€usive closer to the source. Figs. 7 and 8 illustrate the behaviour of the total neutron ¯ux for di€erent anisotropies in ®nite slab injected by a symmetrically localised source. The following Figs. 9 and 10 concern the case where the same type of material anisotropy as for the total cross-section is assumed to a€ect also the scattering cross-section. In particular, the two ®gures report results for an elliptic form of the anisotropy, for both the total and the scattering cross-sections, for two di€erent values of the scattering ratio s =: 5. Conclusions In the present work, an analytical solution is derived for transport problems with anisotropic media. The solution is obtained for in®nite systems both in exact transport and in the discrete ordinate approximation. Using the spherical harmonics technique, the PN model is obtained and the consistent di€usion equation is derived, including the proper de®nition of the di€usion coecient. Numerical calculations are performed, using the analytical solution in order to obtain results that can be used as benchmarks for numerical method validation. The interpretation of the results allows one to gain a physical insight into the phenomena involved in neutron transport in anisotropic media. Acknowledgements This work was performed in the framework of a collaboration between CEACadarache and Politecnico di Torino and it was funded by CEA.

References Bell, G.I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand Reinhold, New York. Case, K.M., Zweifel, P.F., 1967. Linear Transport Theory. Addison-Wesley, Reading. Case, K.M., De Ho€mann, F., Placzek, G., 1953. Introduction to the Theory of Neutron Di€usion. Los Alamos Scienti®c Laboratory, Los Alamos, New Mexico. Cassell, J.S., Williams, M.M.R., 1992. Discrete eigenvalues of the one-speed transport equation for asymmetric scattering. Annals of Nuclear Energy 19, 403. Coppa, G., Ravetto, P., 1980. Green function for the discrete ordinate approximation to the Boltzmann transport equation. Transport Theory and Statistical Physics 9, 185.

1286

G. Lapenta et al. / Annals of Nuclear Energy 28 (2001) 1271±1286

Ganapol, B., 1991. Radiative Transfer in One-Dimensional Plant Canopies. Conference on Transport Theory, Albuquerque. Lapenta, G., Ravetto, P., Rostagno, M.M., 2000. Transactions of the American Nuclear Society 82, 146. Rimpault, G., Smith, P.J., Newton, T.D., 1999. Advanced Method for Treating Heterogeneity and Streaming E€ects in Gas Cooled Fast Reactors, Proceedings of the International Conference on Mathematics and Computation, Reactor Physics and Environmental Analysis in Nuclear Applications, 568, Madrid. Williams, M.M.R., 1978. Transport theory in anisotropic media. Math. Proc. Camb. Phil. Soc. 84, 549.