Numerical study of scattering legendre moments and effect of anisotropic scattering on S_N shielding calculation

Numerical study of scattering legendre moments and effect of anisotropic scattering on S N shielding calculation a first-author b a corresp Corresponding author, zhangbin@ncepu.edu.cn Corresponding author, zhangbin@ncepu.edu.cn zhangbin@ncepu.edu.cn c a Yi-Xue Chen a School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China; School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China; China Nuclear Power Engineering Co. Ltd., Beijing 100840, China; China Nuclear Power Engineering Co. Ltd., Beijing 100840, China; China Institute of Atomic Energy, Beijing 102413, China China Institute of Atomic Energy, Beijing 102413, China In neutron and photon transport problems, anisotropic scattering is of great importance for the particle flux, especially when the angular flux has a strong forward peak in shielding analyses. The conventional Legendre expansion is widely-used in discrete ordinates transport codes because of algebraic simplifications with spherical harmonics for the scattering source. However, negative cross-sections caused by the finitely truncated expansion may give rise to a negative source and flux. A simple method is adopted, based on integrating functions of scattering moments, to evaluate anisotropy and convergence of expanded functions. A series of problems were designed with angular fluxes of different anisotropy, and numerical simulations were performed using the ARES transport code to study different treatments and algorithms for scattering. Results show that the diagonal transport approximation is more stable, and obtains a similar accuracy with the extended approximation. A conservative fix-up for the negative source could ensure particle balance and improve computational accuracy significantly for photon transport. The effect of anisotropic scattering is problem-dependent, and no distinct differences among various methods are observed for volume source problems with a continuous energy source. For beam source problems, flux results are sensitive to negative scattering functions, and strictly non-negative cross-sections need to be implemented. Particle transport Shielding calculation Discrete ordinates method Anisotropic scattering Transport approximation Negative source fix-up 1. Introduction 1 The transport and slowing down process of particles in a nuclear system is described by the linear Boltzmann transport equation, and the discrete ordinates method is one of widely-used and highly-efficient deterministic methods for solving it. The anisotropic scattering effect is of importance both for reactor physics and the shielding calculation, as reactors’ heterogeneous effect cannot be neglected, and geometric structures become more complex. Appropriately handling scattering anisotropy has an important influence on accuracy and efficiency. Especially for shielding problems, highly-anisotropic angular flux distribution strengthens the influence of anisotropic scattering. Considering the fact that the influence on computational accuracy of scalar flux is problem-dependent, systematic analyses of numerical algorithms and cross-sections used in simulations are required. The conventional treatment for scattering in discrete ordinates transport codes is to expand the scattering cross-section with truncated Legendre polynomials [ 1 1 ], and the spherical harmonics expansion method is a usual option to calculate scattering source terms. The requirement for the Legendre expansion varies among diverse models, and significant experience shows that the P 5 order is sufficient to obtain accurate integral quantities [ 2 2 ], such as detector responses and fast neutron fluence rates for some problems of RPV analyses and fusion installations shielding. For the light nuclei’s elastic scattering reactions, merely increasing the expansion order cannot eliminate oscillations of approximated angular distributions. Therefore, Odom proposed the exact-kernel method [ 3 3 ] and modified the scattering source calculation method in deterministic codes, which generated much more accurate angular fluxes, but involved a larger storage requirement and longer execution time. The finite Legendre expansion generates negatives for differential scattering cross-sections. A constrained least-squares algorithm, which modified high-order scattering Legendre moments, and a maximum entropy method [ 4 4 ] to represent scattering cross-sections were developed by Dahl to correct this deficiency, and his methods could create a non-negative expansion. Another family of non-negative scattering cross-section generating methods, inspired by or based on the Monte Carlo theory, were applied in S N and multi-group MC codes, including the hybrid Monte Carlo-discrete elements method proposed by DelGrande [ 5 5 ] and the piecewise-average group cross section method proposed by Gerts [ 6 6 ]. Jong Woon Kim employed angular data of elastic scattering [ 7 7 ] from ENDF libraries as a complementary part for basic multi-group cross-sections, if needed. Moreover, transport approximations/corrections [ 8 8 , 9 9 ] and fix-up selections in the scattering source calculation [ 10 10 ] also play a significant role in flux results. Additionally, the Galerkin quadrature [ 11 11 , 12 12 ] proposed by Morel is a powerful technique for a highly anisotropic scattering source, and its application mainly concentrates on charged particle transport. Some numerical analyses on anisotropic scattering have been performed for reactor critical calculation [ 13 13 – 16 16 13–16 ] and analytical transport problems [ 17 17 , 18 18 ]. The purpose of this study is to quantitatively analyze the impact of scattering anisotropy on S N transport calculation using the conventional Legendre expansion method. The required Legendre expansion order depends on both materials’ physical properties and angular flux distributions. The analyses in this paper are beneficial for users to choose proper algorithms for realistic shielding problems. We produce multi-group cross-sections with the P 9 expansion in the LANL 30N12G and Vitamin-J 175N42G structures using NJOY2016 [ 19 19 ]. We studied the scattering cross section expansion convergence of actual materials. Effects of parameters on scattering for neutron and photon transport were investigated, including the expansion order, transport correction, scattering source fix-up, and other model features. The remainder of the paper is organized as follows. Section 2 briefly describes methods regarding the scattering source calculation for S N codes. The convergence analyses for water, iron, and other materials under truncated Legendre expansions are given in Section 3. We simulated a number of fix-source problems with simple geometry and homogeneous media, and the numerical results are summarized in Sect. 4. We draw some conclusions in Sect. 5. 2 Review of theory 1 ARES [ 20 20 – 22 22 20–22 ] is a multi-dimensional neutral particle transport code based on the discrete ordinates method, which is capable of solving critical and shielding problems. The general procedure and methods for calculating the scattering source in S N codes are briefly reviewed below. For simplicity, the time-independent transport equation without the fission source term in 1D slab geometry is considered as: <math display="block" id="M1"><mtable columnalign="left"><mtr><mtd><mi>μ</mi><mfrac><mrow><mo>∂</mo><mi>ψ</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>E</mi><mo>,</mo><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><msub><mi>Σ</mi><mtext>t</mtext></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo><mi>ψ</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>E</mi><mo>,</mo><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><mo>,</mo></mtd></mtr><mtr><mtd><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>∞</mi></msubsup><mrow><mtext>d</mtext><msup><mi>E</mi><mo>′</mo></msup></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><msup><mi>Ω</mi><mo>′</mo></msup></msub><mrow><msub><mi>Σ</mi><mtext>s</mtext></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><msup><mi>E</mi><mo>′</mo></msup><mo>→</mo><mi>E</mi><mo>,</mo><mover accent="true"><msup><mi>Ω</mi><mo>′</mo></msup><mo>→</mo></mover><mo>→</mo><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mi>ψ</mi></mrow></mrow></mstyle><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><msup><mi>E</mi><mo>′</mo></msup><mo>,</mo><mover accent="true"><msup><mi>Ω</mi><mo>′</mo></msup><mo>→</mo></mover><mo stretchy="false">)</mo><mi>d</mi><msup><mi>Ω</mi><mo>′</mo></msup><mo>+</mo><msub><mi>Q</mi><mrow><mtext>fixed</mtext></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>E</mi><mo>,</mo><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mtd></mtr></mtable></math> where <math id="M2"><mi>ψ</mi></math> is the angular flux; <math id="M3"><mi>μ</mi></math> is the direction cosine of <math id="M4"><mover accent="true"><mi>Ω</mi><mo>→</mo></mover></math> ; <math id="M5"><mrow><msub><mi>Q</mi><mrow><mtext>fixed</mtext></mrow></msub></mrow></math> is the fixed source term; <math id="M6"><mrow><msub><mi>Σ</mi><mtext>t</mtext></msub></mrow></math> and <math id="M7"><mrow><msub><mi>Σ</mi><mtext>s</mtext></msub></mrow></math> are macroscopic total and scattering cross-sections, respectively; x , E, and <math id="M8"><mover accent="true"><mi>Ω</mi><mo>→</mo></mover></math> are the spatial, energy, and direction variables, respectively. To solve the scattering integral after the energy and angular discretization, the scattering cross-section is a function of the scattering angle and can be expanded in the series of Legendre polynomials as Eq. (2). The angular flux can likewise be expanded in Legendre polynomials or the spherical harmonics functions. According to the addition theorem, the scattering source term for discrete direction m and group g is calculated as Eq. (3) using the L -order expansion. <math display="block" id="M9"><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><msub><mi>μ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mn>0</mn><mi>∞</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow><mrow/></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mi>P</mi><mi>l</mi></msub></mrow></mstyle><mo stretchy="false">(</mo><msub><mi>μ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>,</mo></mrow></math> <math display="block" id="M10"><mrow><msubsup><mi>Q</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi></mrow><mi>m</mi></msubsup><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>l</mi><mo>=</mo><mn>0</mn></mrow><mi>L</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>l</mi><mo>+</mo><mn>1</mn></mrow><mrow><mn>4</mn><mi>π</mi></mrow></mfrac></mrow></mstyle><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow><mrow/></msubsup><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mo>−</mo><mi>l</mi></mrow><mi>l</mi></munderover><mrow><msubsup><mi>ϕ</mi><mrow><mi>l</mi><mo>,</mo><mi>n</mi><mo>,</mo><mi>g</mi></mrow><mrow/></msubsup><msubsup><mi>Y</mi><mrow><mi>l</mi><mo>,</mo><mi>n</mi></mrow><mrow/></msubsup><mo stretchy="false">(</mo><msub><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mi>m</mi></msub><mo stretchy="false">)</mo></mrow></mstyle></mrow></mstyle><mo>,</mo></mrow></math> where <math id="M11"><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi></mrow></msub></mrow></math> is the l -th scattering Legendre moment; <math id="M12"><mrow><msub><mi>P</mi><mi>l</mi></msub></mrow></math> is the l -order Legendre polynomial; <math id="M13"><mrow><msubsup><mi>ϕ</mi><mrow><mi>l</mi><mo>,</mo><mi>n</mi></mrow><mrow/></msubsup></mrow></math> denotes the l -th, n -th flux moment; <math id="M14"><mrow><msubsup><mi>Y</mi><mrow><mi>l</mi><mo>,</mo><mi>n</mi></mrow><mrow/></msubsup></mrow></math> denotes the l -th, n -th spherical harmonics function, and G denotes the total number of energy groups. This spherical harmonics expansion method is easy to use because of its small data storage requirement. If the exact-kernel representation or the piecewise function representation for scattering cross sections is used, the formula for the scattering source is modified to Eq.(4). Hence, angular fluxes for each discrete angle, and transfer cross-sections for each possible scattering angle need to be stored. <math display="block" id="M15"><mrow><msubsup><mi>Q</mi><mrow><mi>s</mi><mo>,</mo><mi>g</mi></mrow><mi>m</mi></msubsup><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>m</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>w</mi><msup><mi>m</mi><mo>′</mo></msup></msub><msubsup><mi>Σ</mi><mrow><mi>s</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow><mrow><msup><mi>m</mi><mo>′</mo></msup><mo>→</mo><mi>m</mi></mrow></msubsup><msub><mi>ψ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>m</mi><mo>′</mo></msup></mrow></msub></mrow></mstyle></mrow></mstyle><mo>,</mo></mrow></math> where <math id="M16"><mrow><msubsup><mi>Σ</mi><mrow><mtext>s,</mtext><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow><mrow><msup><mi>m</mi><mo>′</mo></msup><mo>→</mo><mi>m</mi></mrow></msubsup></mrow></math> denotes the group-to-group and direction-to-direction scattering cross-section; <math id="M17"><mi>w</mi></math> denotes the discrete direction weight, and M denotes the total direction number. When the expansion order is not sufficiently high to adequately handle the anisotropy, the expanded cross-section may exhibit non-physical negative values among part of the angular range, and this may bring about a negative scattering source for certain directions. The negative source is one of reasons for negative angular flux besides spatial discretization and mesh size. In the iterations of transport sweep, treatments for the negative source include economic and conservative manners. The economic fix-up directly sets the negative total source to zero. The conservative method first sorts sources of all directions and then removes all negatives as well as some small positives until the 0-th source moment <math id="M18"><mrow><msub><mi>Q</mi><mn>0</mn></msub></mrow></math> in Eq.(5) is restored, which ensures particle conservation as far as possible. <math display="block" id="M19"><mrow><msubsup><mi>Q</mi><mn>0</mn><mrow/></msubsup><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mn>4</mn><mi>π</mi></mrow></msub><mrow><mi>Q</mi><mo stretchy="false">(</mo><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mtext>d</mtext><mi>Ω</mi><mo>≅</mo></mrow></mrow></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>w</mi><mi>m</mi></msub><msubsup><mi>Q</mi><mrow/><mi>m</mi></msubsup></mrow></mstyle></mrow></math> To avoid the negative scattering cross-section, increasing expansion orders, the Legendre moment adjustment, and piecewise constant function representation are all available methods. In addition, proper transport approximations [ 8 8 , 9 9 ] could improve accuracy with less cost for the Legendre expansion method. Transport cross-sections with isotropic scattering are usually used in core physics analyses as the treatment for scattering anisotropy [ 15 15 ]. Strictly speaking, multi-group total cross sections are angle-dependent to a certain extent in theory, as presented in the multi-group 1D P N equation as Eq. (6). <math display="block" id="M20"><mtable columnalign="left"><mtr><mtd><mi>μ</mi><mfrac><mrow><mo>∂</mo><msub><mi>ψ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>l</mi><mo>=</mo><mn>0</mn></mrow><mi>L</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msub><mi>P</mi><mi>l</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></mstyle><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>l</mi><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><mi>g</mi><mo>,</mo><mi>l</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo></mtd></mtr><mtr><mtd><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>l</mi><mo>=</mo><mn>0</mn></mrow><mi>L</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msub><mi>P</mi><mi>l</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup><mi>G</mi></munderover><mrow><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mi>N</mi></msub></mrow></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><mi>l</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo></mrow></mstyle><mo>+</mo><msub><mi>Q</mi><mrow><mtext>fix</mtext><mo>,</mo><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></mstyle></mtd></mtr></mtable></math> Eq. (6) can be transferred to the S N equation by moving the second term to the right side and then adding <math id="M21"><mrow><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>S</mi><mtext>N</mtext></msub></mrow></msubsup><msub><mi>ψ</mi><mi>g</mi></msub></mrow></math> term to both sides of the equation. However, angular-dependent total cross-sections or the high-order Legendre moment of total cross-sections are not used in the actual S N calculation, and this relationship is simplified by various transport approximations [ 19 19 ] by choosing total cross-sections and scattering moments as Eq. (7). <math display="block" id="M22"><mtable columnalign="left"><mtr><mtd><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow><mrow><msub><mi>S</mi><mtext>N</mtext></msub></mrow></msubsup><mo>=</mo><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup><mo>,</mo><mtext>     </mtext><msup><mi>g</mi><mo>′</mo></msup><mo>≠</mo><mi>g</mi></mtd></mtr><mtr><mtd><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi><mo>,</mo><mi>g</mi><mo>→</mo><mi>g</mi></mrow><mrow><msub><mi>S</mi><mtext>N</mtext></msub></mrow></msubsup><mo>=</mo><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>l</mi><mo>,</mo><mi>g</mi><mo>→</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup><mo>−</mo><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>l</mi><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup><mo>+</mo><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>S</mi><mtext>N</mtext></msub></mrow></msubsup></mtd></mtr><mtr><mtd><mtext>Consistent-</mtext><mi>P</mi><mo>:</mo><mtext> </mtext><msubsup><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>S</mi><mtext>N</mtext></msub></mrow></msubsup><mo>=</mo><msubsup><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup></mtd></mtr><mtr><mtd><mtext>        Diagonal</mtext><mo>:</mo><mtext> </mtext><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>S</mi><mtext>N</mtext></msub></mrow></msubsup><mo>=</mo><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>L</mi><mo>+</mo><mn>1</mn><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup><mo>−</mo><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>L</mi><mo>+</mo><mn>1</mn><mo>,</mo><mi>g</mi><mo>→</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup></mtd></mtr><mtr><mtd><mtext>            Extend</mtext><mo>:</mo><mtext> </mtext><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>S</mi><mtext>N</mtext></msub></mrow></msubsup><mo>=</mo><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>L</mi><mo>+</mo><mn>1</mn><mo>,</mo><mi>g</mi></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup><mo>−</mo><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>L</mi><mo>+</mo><mn>1</mn><mo>,</mo><mi>g</mi><mo>→</mo><msup><mi>g</mi><mo>′</mo></msup></mrow><mrow><msub><mi>P</mi><mtext>N</mtext></msub></mrow></msubsup></mrow></mstyle></mtd></mtr></mtable></math> The consistent-P approximation directly adopts the scalar flux-weighted total cross-sections and the L -order truncated Legendre moments of scattering cross-sections. Other approximations modify total cross-sections and within-group components of the scattering matrix using high-order information. The diagonal transport approximation considers within-group components in the ( L +1)-order scattering matrix and uses them to correct total cross sections. The extended transport approximation assumes the ( L +1)-order term of a scattering source for a given group equal to the scattering out of this group. These approximations reduce errors caused by disregarding high-order expansion terms. The cross-section data needed for transport corrections are generated by the GROUPR module of NJOY and stored in MATXS-format libraries for the processing code TRANSX [ 23 23 ]. 3 Convergence analyses of Legendre expansion 1 To date, the Legendre expansion has been the most widely-used method in deterministic codes because of its stability and usability for majority of neutron shielding problems. Here, a simple method is presented to compare the anisotropy of different materials. Results of convergence analyses could serve as auxiliaries and supplements for evaluating the relative strength of anisotropic scattering. For analyses of water and iron, multi-group cross-sections were generated based on the ENDF/B-VII.1 data by NJOY2016 with the P 9 expansion in LANL 30N12G and Vitamin-J 175N42G structures. Subsequently, TRANSX was used to produce materials’ macroscopic cross-sections, where the density of water and iron were set to 1.0 g·cm -3 and 7.87 g·cm -3 , respectively. Elastic neutron scattering of light nuclei for high energy and Compton scattering caused by photons exhibit a highly forward peak. We mainly focus on the anisotropy of macroscopic cross-sections of common shielding and structure materials. First, the normalized expansions for differential scattering cross-sections are adopted as Eq. (8) to investigate the relative angular distribution. <math display="block" id="M23"><mrow><msub><mi>F</mi><mi>L</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>l</mi><mo>=</mo><mn>0</mn></mrow><mi>L</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>s</mi><mo>,</mo><mi>l</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>Σ</mi><mrow><mi>s</mi><mo>,</mo><mn>0</mn></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo></mrow></mfrac><msub><mi>P</mi><mi>l</mi></msub></mrow></mstyle><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> Functions in Eq. (8) can be drawn directly to observe the distribution, and analyses were performed based on the differential scattering cross-section functions. Figure 1 Figure 1 shows some normalized differential neutron scattering cross-sections of water under the consistent P approximation in the LANL group structure. Results and experiences show that within-group scattering is generally forward-peaky, the peak place of group differential scattering cross-sections moves toward <math id="M24"><mrow><msub><mi>μ</mi><mn>0</mn></msub><mo>=</mo><mn>0</mn></mrow></math> as the energy loss increases, and scattering in the high-energy range is more anisotropic. Additionally, there is a relationship between self-shielding treatment and anisotropic scattering. The self-shielding effect affects not only total cross-sections in the resonance energy range, but also scattering cross-sections as well as the corresponding angular distribution. The Bondarenko method [ 23 23 ] was adopted as the treatment for the resonance self-shielding effect. Comparison of infinite-dilute and self-shielded cross-sections was performed based on the LANL-structure library for iron material, as shown in Figure 2 Figure 2 . According to the normalized within-group scattering functions of the 17 th group (2.48×10 -2 - 6.76×10 -2 MeV), anisotropy of within-group scattering is weakened with self-shielding treatment. Photon production cross-sections are also modified by the self-shielding effect, which is disregarded in this study because of its almost isotropic scattering property. Considering the investigation based on drawing scattering angular distributions is unrealistic for larger datasets, inspired by the scattering sampling functions of multi-group MC, we integrate normalized functions in the entire scattering angle interval and analyze anisotropy by the convergence of integrals. Factors <math id="M25"><mrow><msub><mi>f</mi><mi>L</mi></msub></mrow></math> are defined to describe the gap between L -order truncated and ( L -2)-order truncated expansions, as indicated by Eq.(9). <math display="block" id="M26"><mrow><msubsup><mi>f</mi><mi>L</mi><mrow/></msubsup><mo>=</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>1</mn></msubsup><mrow><mrow><mo>|</mo><mrow><msub><mi>F</mi><mi>L</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>F</mi><mrow><mi>L</mi><mo>−</mo><mn>2</mn></mrow></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow><mo>|</mo></mrow><mi>d</mi><mi>μ</mi></mrow></mrow></mstyle></mrow></math> Factors <math id="M27"><mrow><msub><mi>f</mi><mi>L</mi></msub></mrow></math> are integrals that measure the area of difference between normalized expansion functions, and the Legendre expansion is converged if this integral factor is small enough for L order. A more anisotropic scattering of materials leads to a larger integral <math id="M28"><mrow><msub><mi>f</mi><mi>L</mi></msub></mrow></math> for a given order. Considering that the true multi-group scattering angular distribution is unknown, these integrals evaluate the closeness of expanded functions of adjacent orders, and the intensity of anisotropy could be estimated. Integrals are computed including <math id="M29"><mrow><msubsup><mi>f</mi><mn>3</mn><mrow/></msubsup></mrow></math> , <math id="M30"><mrow><msubsup><mi>f</mi><mn>5</mn><mrow/></msubsup></mrow></math> , <math id="M31"><mrow><msubsup><mi>f</mi><mn>7</mn><mrow/></msubsup></mrow></math> , and <math id="M32"><mrow><msubsup><mi>f</mi><mn>9</mn><mrow/></msubsup></mrow></math> with the Gauss quadrature based on the library of LANL structure under the consistent P approximation, and factors for the within-group neutron and gamma scatters of water and iron are shown in Fig. 3 Fig. 3 . For neutron scattering of a given group, the factor value decreases as the order increases, which means that the Legendre expansion functions converge and negative ranges of cross-sections are gradually reduced. Within-group neutron scattering of water is more anisotropic than that of iron because of H elastic scattering. On the contrary to neutron cases, factors of gamma scattering increase slightly, which represents intensive oscillations and inevitable negatives in the differential cross-sections distribution. Negatives of expanded cross-sections could be neglected when factors have small values. Factors of scattering between different groups are difficult to converge, because the peak of differential cross-sections is close to the middle. The impact of group structures on the anisotropy was also investigated in this way. Factors of water’s neutron and gamma within-group scattering are listed in Table 1 Table 1 based on libraries with different group structures. As energy groups are refined, within-group scatters become more anisotropic, which is proved by the fact that larger <math id="M33"><mrow><msub><mi>f</mi><mi>L</mi></msub></mrow></math> factors of the same order are obtained for a finer group structure. 1 1   1 1 Group number 1 1 Upper energy(MeV) 1 1 f 3 1 1 f 5 1 1 f 7 1 1 f 9 1 1 Neutron 6 1 LANL 30N12G Group Structure   1 1 1 1 5 1 1 10.00 1 1 1.46 1 1 1.26 1 1 0.71 1 1 0.32 1 1 1 1 6 1 1 7.79 1 1 1.33 1 1 1.12 1 1 0.59 1 1 0.29 1 1 6 1 Vitamin-J 175N42G Group Structure   1 1 1 1 17 1 1 10.00 1 1 2.04 1 1 1.61 1 1 1.00 1 1 0.63 1 1 1 1 18 1 1 9.51 1 1 2.00 1 1 1.55 1 1 0.97 1 1 0.68 1 1 1 1 19 1 1 9.05 1 1 2.03 1 1 1.61 1 1 1.00 1 1 0.63 1 1 1 1 20 1 1 8.61 1 1 2.19 1 1 1.93 1 1 1.37 1 1 0.92 1 1 1 1 21 1 1 8.19 1 1 2.09 1 1 1.75 1 1 1.20 1 1 0.86 1 1 1 1 22 1 1 7.79 1 1 2.01 1 1 1.59 1 1 1.01 1 1 0.71 1 1 Gamma 6 1 LANL 30N12G Group Structure 1 1 1 1 3 1 1 8.0 1 1 2.74 1 1 3.46 1 1 3.93 1 1 4.24 1 1 1 1 4 1 1 7.0 1 1 2.73 1 1 3.41 1 1 3.84 1 1 4.07 1 1 6 1 Vitamin-J 175N42G Group Structure 1 1 1 1 7 1 1 8.0 1 1 2.76 1 1 3.53 1 1 4.09 1 1 4.52 1 1 1 1 8 1 1 7.5 1 1 2.76 1 1 3.52 1 1 4.07 1 1 4.48 1 1 1 1 9 1 1 7.0 1 1 2.76 1 1 3.51 1 1 4.04 1 1 4.44 1 1 1 1 10 1 1 6.5 1 1 2.75 1 1 3.49 1 1 4.01 1 1 4.38 The anisotropy of several materials was likewise investigated in this manner including water, polyethylene, concrete, graphite, iron, SS304, and lead based on the Fendl3.0MG libraries. For neutron scattering of materials, the scattering anisotropy is relatively strong for all materials from 13 MeV to 50 MeV, whose f 5 factors range from 1.0 to 2.6. Because of the hydrogen element and its elastic scattering, within-group scattering of water and polyethylene is remarkably more anisotropic than other materials’ scattering below 0.1 MeV, where their f 5 factors are greater than 1.0, whereas others are less than 0.1. Gamma scattering shows strong anisotropy similar to water and iron data above, and f 5 factors of these materials become 2.8–3.6 from 1.3 MeV to 50 MeV. Further analyses could be performed to predict accuracy and negatives of expanded scattering functions, and these comparisons may be helpful for users to pick key materials with strong anisotropy within a certain energy range. 4 Numerical results 1 The generation of multi-group cross-sections, algorithms of calculating the scattering source and model features determine the effects of scattering anisotropy and simulation accuracy together. To investigate these effects, a series of 2D fixed-source problems with simple geometry were simulated with the ARES code. Models of test problem 1 and 2 have similar structures shown in Fig. 4 Fig. 4 , where the two side lengths of the source region are 5 cm and 0.5 cm. Particles move more than 10 mean-free-paths from the source region to the system boundaries. The isotropic source only exists in the first group, because of the high-energy scattering’ strong anisotropy. Mesh sizes are set to 0.5 cm using weighted differencing spatial discretization, owing to the good resistance to negative fluxes. The Legendre Chebyshev quadrature sets of S 32 order are adopted to minimize the angular discretization error. The iteration convergence criterion for group flux is equal to 5E-5, and no iteration acceleration technique is employed. The behaviors of different self-shielding treatments are out of the scope of this paper. The scalar flux results of the P 9 expansion under consistent-P approximation with the negative source conservative fix-up serve as references. To evaluate the integral errors of group fluxes, dose rates are computed based on the ANSI/ANS-6.1.1-1977 conversion factors, and the output comprises root-mean-square errors in Eq.(10) and maximum errors in Eq.(11) for mesh dose rates. <math display="block" id="M34"><mrow><msub><mi>E</mi><mrow><mtext>RMS</mtext></mrow></msub><mo>=</mo><msqrt><mrow><mfrac><mn>1</mn><mi>N</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mi>j</mi><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>g</mi><mi>G</mi></munderover><mrow><msub><mi>C</mi><mi>g</mi></msub><msub><mi>ϕ</mi><mrow><mtext>calc</mtext><mo>,</mo><mi>i</mi><mi>j</mi><mi>k</mi><mo>,</mo><mi>g</mi></mrow></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>g</mi><mi>G</mi></munderover><mrow><msub><mi>C</mi><mi>g</mi></msub><msub><mi>ϕ</mi><mrow><mtext>ref</mtext><mo>,</mo><mi>i</mi><mi>j</mi><mi>k</mi><mo>,</mo><mi>g</mi></mrow></msub></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>g</mi><mi>G</mi></munderover><mrow><msub><mi>C</mi><mi>g</mi></msub><msub><mi>ϕ</mi><mrow><mtext>ref</mtext><mo>,</mo><mi>i</mi><mi>j</mi><mi>k</mi><mo>,</mo><mi>g</mi></mrow></msub></mrow></mstyle></mrow></mfrac><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></mstyle></mrow></msqrt></mrow></math> <math display="block" id="M35"><mrow><msub><mi>E</mi><mrow><mtext>MAX</mtext></mrow></msub><mo>=</mo><munder><mrow><mi>max</mi></mrow><mrow><mn>1</mn><mo>≤</mo><mi>i</mi><mi>j</mi><mi>k</mi><mo>≤</mo><mi>N</mi></mrow></munder><mrow><mo>|</mo><mrow><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>g</mi><mi>G</mi></munderover><mrow><msub><mi>C</mi><mi>g</mi></msub><msub><mi>ϕ</mi><mrow><mtext>calc</mtext><mo>,</mo><mi>i</mi><mi>j</mi><mi>k</mi><mo>,</mo><mi>g</mi></mrow></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>g</mi><mi>G</mi></munderover><mrow><msub><mi>C</mi><mi>g</mi></msub><msub><mi>ϕ</mi><mrow><mtext>ref</mtext><mo>,</mo><mi>i</mi><mi>j</mi><mi>k</mi><mo>,</mo><mi>g</mi></mrow></msub></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>g</mi><mi>G</mi></munderover><mrow><msub><mi>C</mi><mi>g</mi></msub><msub><mi>ϕ</mi><mrow><mtext>ref</mtext><mo>,</mo><mi>i</mi><mi>j</mi><mi>k</mi><mo>,</mo><mi>g</mi></mrow></msub></mrow></mstyle></mrow></mfrac></mrow><mo>|</mo></mrow></mrow></math> where <math id="M36"><mrow><msub><mi>ϕ</mi><mrow><mi>i</mi><mi>j</mi><mi>k</mi><mo>,</mo><mi>g</mi></mrow></msub></mrow></math> denotes the scalar flux of ijk -th mesh of g group, <math id="M37"><mrow><msub><mi>C</mi><mi>g</mi></msub></mrow></math> denotes the flux-dose conversion factor of g group, N denotes the total number of meshes. For neutron transport in water, the resulting errors obtained from model 1 and the computational time normalized by the shortest time are listed in Table 2 Table 2 . Generally, errors of mesh dose rates are quite small with P 3 or higher expansions. With expansions of the same order and economic fix-up, extended and diagonal approximations obtain better results compared to the consistent P method, because these two corrections implicitly contain high-order scattering moment information, which ameliorates the negatives of the within-group scattering angular distribution. Better results are acquired with no fix-up or conservative fix-up for the negative source compared to results obtained by the economic fix-up. The economic fix-up ignores the 0-th source moment conservation and damages particle balance. It exaggerates total scattering and leads to larger flux results. This effect may not be significant when negatives of scattering angular distribution are not severe, as in Table 2 Table 2 . Although good scalar quantities could be gained without the fix-up, the angular flux distribution would be disturbed at the permission of a negative source, which is discussed below. The situation for neutron transport in iron is similar, and the errors are smaller. From the error distribution of dose rates, relative errors of results with the economic fix-up retain positive signs, and error values rise gradually as neutron transport distance increases. No clear pattern is observed in the error distributions of the other two algorithms. The relative errors of dose rates are similar to those of high-energy group fluxes and slightly larger than those at low energy. The <math id="M38"><mrow><msub><mi>E</mi><mrow><mtext>RMS</mtext></mrow></msub></mrow></math> of dose rates could represent the average level of flux accuracy. *Negatives exist in 0-th scattering moments generated by TRANSX for within-group scattering with P 1 expansion. 1 1 Transport approximation 1 1 Expansion order 9 1 Negative source fix-up 1 1 1 1 3 1 No fix-up 3 1 Economic fix-up 3 1 Conservative fix-up 1 1 1 1 1 1 E RMS 1 1 E MAX 1 1 Time 1 1 E RMS 1 1 E MAX 1 1 Time 1 1 E RMS 1 1 E MAX 1 1 Time 1 1 Consistent-P 1 1 P 1 1 1 5.00% 1 1 7.92% 1 1 1.26 1 1 7.83% 1 1 10.52% 1 1 1.33 1 1 3.38% 1 1 5.24% 1 1 2.21 1 1 1 1 P 3 1 1 0.29% 1 1 0.68% 1 1 1.90 1 1 1.63% 1 1 2.13% 1 1 1.91 1 1 0.14% 1 1 0.42% 1 1 2.86 1 1 1 1 P 5 1 1 0.03% 1 1 0.09% 1 1 2.98 1 1 0.68% 1 1 0.92% 1 1 2.99 1 1 0.03% 1 1 0.09% 1 1 3.81 1 1 1 1 P 7 1 1 0.00% 1 1 0.01% 1 1 4.83 1 1 0.17% 1 1 0.25% 1 1 4.80 1 1 0.01% 1 1 0.02% 1 1 5.18 1 1 Extended 1 1 P 1 * 1 1 / 1 1 / 1 1   1 1 / 1 1 / 1 1   1 1 / 1 1 / 1 1   1 1 1 1 P 3 1 1 0.12% 1 1 0.41% 1 1 1.95 1 1 0.34% 1 1 0.57% 1 1 1.97 1 1 0.12% 1 1 0.48% 1 1 2.92 1 1 1 1 P 5 1 1 0.03% 1 1 0.17% 1 1 3.00 1 1 0.17% 1 1 0.31% 1 1 2.94 1 1 0.04% 1 1 0.16% 1 1 3.70 1 1 1 1 P 7 1 1 0.00% 1 1 0.01% 1 1 4.89 1 1 0.17% 1 1 0.24% 1 1 4.90 1 1 0.01% 1 1 0.02% 1 1 5.20 1 1 Diagonal 1 1 P 1 1 1 2.92% 1 1 4.52% 1 1 1.00 1 1 4.11% 1 1 5.61% 1 1 1.00 1 1 2.30% 1 1 3.56% 1 1 1.75 1 1 1 1 P 3 1 1 0.10% 1 1 0.35% 1 1 1.83 1 1 0.30% 1 1 0.47% 1 1 1.83 1 1 0.10% 1 1 0.40% 1 1 2.74 1 1 1 1 P 5 1 1 0.03% 1 1 0.17% 1 1 2.90 1 1 0.17% 1 1 0.31% 1 1 2.93 1 1 0.04% 1 1 0.17% 1 1 3.69 1 1 1 1 P 7 1 1 0.02% 1 1 0.09% 1 1 4.63 1 1 0.05% 1 1 0.11% 1 1 4.72 1 1 0.02% 1 1 0.08% 1 1 5.12 When low expansions are employed for neutron transport, relative errors of group fluxes and dose rates increase using the 175N42G library compared to the corresponding results using 30N12G library with the economic fix-up, while errors remain similar when using these two libraries with the conservative fix-up. Root-mean-square errors of flux for each group under the diagonal approximation are shown in Fig. 5 Fig. 5 , and the effect of group refinement is observed. The conservative fix-up method reduces differences of relative errors between fine and coarse libraries with the same expansion, which could compensate the effect of group refinement to a certain extent. Errors obtained for the gamma transport of iron are listed in Table 3 Table 3 . The extended approximation produces negative 0-th scattering moments for several groups at all expansion orders. The behavior of the economic fix-up becomes worse, as the negative situation becomes more severe for gamma scattering, which is indicated from the analyses in Sect.3. The diagonal approximation and conservative fix-up are more proper for gamma transport problems. Errors in the gamma transport of water are smaller than those in iron owing to the longer mean-free-path of photons in water, despite similar anisotropy. *One group of this case achieved the maximum number of inner iterations, and the iteration error of that group only reached 1e-4. 1 1 Transport approximation 1 1 Expansion order 9 1 Negative source fix-up 1 1 1 1 3 1 No fix-up 3 1 Economic fix-up 3 1 Conservative fix-up 1 1 1 1 1 1 E RMS 1 1 E MAX 1 1 Time 1 1 E RMS 1 1 E MAX 1 1 Time 1 1 E RMS 1 1 E MAX 1 1 Time 1 1  Consistent-P 1 1 P 1 1 1 4.40% 1 1 21.28% 1 1 1.43 1 1 4.32% 1 1 21.22% 1 1 1.44 1 1 6.59% 1 1 23.52% 1 1 4.13 1 1 1 1 P 3 1 1 0.32% 1 1 3.00% 1 1 2.10 1 1 4.06% 1 1 6.49% 1 1 2.07 1 1 0.17% 1 1 0.71% 1 1 3.22 1 1 1 1 P 5 1 1 0.09% 1 1 0.56% 1 1 3.22 1 1 2.43% 1 1 4.31% 1 1 3.21 1 1 0.11% 1 1 0.47% 1 1 4.33 1 1 1 1 P 7 1 1 0.05% 1 1 0.23% 1 1 7.69 * 1 1 1.37% 1 1 2.84% 1 1 5.17 1 1 0.04% 1 1 0.19% 1 1 5.81 1 1 Diagonal 1 1 P 1 1 1 3.00% 1 1 14.68% 1 1 1.00 1 1 2.96% 1 1 14.57% 1 1 1.01 1 1 4.49% 1 1 16.14% 1 1 1.84 1 1 1 1 P 3 1 1 0.22% 1 1 2.20% 1 1 1.70 1 1 2.59% 1 1 4.57% 1 1 1.73 1 1 0.14% 1 1 0.53% 1 1 2.64 1 1 1 1 P 5 1 1 0.10% 1 1 0.60% 1 1 2.86 1 1 1.38% 1 1 2.85% 1 1 2.81 1 1 0.11% 1 1 0.83% 1 1 3.78 1 1 1 1 P 7 1 1 0.09% 1 1 0.50% 1 1 3.78 1 1 0.72% 1 1 1.90% 1 1 4.66 1 1 0.08% 1 1 0.64% 1 1 5.28 The particle balance for the entire system is evaluated by comparing the sum of the total leakage and total collision, which are calculated based on fluxes of the i- th iteration, and total source, which is calculated based on fluxes of ( i -1)-th iteration without any fix-ups. The balance ratio is defined as Eq. (12), and particles are conserved if its value approaches zero. <math display="block" id="M39"><mrow><msub><mi>R</mi><mrow><mtext>balance</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mrow><mo>|</mo><mrow><mi>L</mi><mi>e</mi><mi>a</mi><mi>k</mi><mi>a</mi><mi>g</mi><mi>e</mi><mo>+</mo><mi>C</mi><mi>o</mi><mi>l</mi><mi>l</mi><mi>i</mi><mi>s</mi><mi>i</mi><mi>o</mi><mi>n</mi><mo>−</mo><mi>S</mi><mi>o</mi><mi>u</mi><mi>r</mi><mi>c</mi><mi>e</mi></mrow><mo>|</mo></mrow></mrow><mrow><mi>S</mi><mi>o</mi><mi>u</mi><mi>r</mi><mi>c</mi><mi>e</mi></mrow></mfrac></mrow></math> This ratio for simulations with no fix-up should be zero in theory, however its magnitude is approximately 1E-15 to 1E-14 because of round-off errors. For neutron transport problems, the imbalance caused by economic fix-up mainly happens in the first few groups and attenuates with increasing expansion order. Ratios of gamma transport in iron under the diagonal approximation are shown in Figure 6 Figure 6 , and ratios for 10-, 11- and 12-th groups are below 1E-14. The conservative method could improve the global balance by about four magnitudes compared to the economic fix-up. Besides monoenergic source problems, continuous energy source problems were calculated using the Watt fission spectrum of 235 U in model 1 involving neutron and photon coupled transport. The maximum E RMS of neutron dose rates is less than 0.1% among all cases for water and iron from P 3 to P 7 orders. The magnitude of error in the gamma dose rates is similar to those of neutron dose rates, considering that the accuracy of photon fluxes mainly depends on the local neutrons than on photons transporting from a distance. The differences among negative source fix-up methods can be disregarded, such that the economic fix-up is always recommended for shielding problems of RPV analyses [ 24 24 ]. For problems with a continuous energy source, the fixed source in each group weakens the effect of anisotropic scattering because of a larger proportion of uncollided flux in the scalar flux. More complicated scattering between groups may overshadow some effects observed in monoenergic source problems. The computational time of cases above was compared to evaluate the algorithm efficiency. With the P 3 or P 5 expansion, the diagonal approximation saves about 1%–15% of time in comparison to other approximations because of the smaller scattering ratio, which leads to fewer inner iterations to achieve convergence criteria. For negative source fix-ups, the conservative method would increase the computational time owing to its cost on sorting and source moment calculation. This influence on the computing time is related to the material property, expansion order, and group structure. When the LANL 30N12G or Vitamin-J 42G library is used, the conservative fix-up would increase the time by about 40% with the P 3 expansion and by about 20% with the P 5 expansion. For the Vitamin-J 175N library, this ratio decreases to about 10%, probably because the fix-up subroutine does not take effect in many low-energy groups. In general, the percentage of the conservative fix-up to the overall computational time would decrease with the expansion and group number increase. No fix-up may cause unacceptable disturbance for the iteration in some cases, preventing the iteration error to decrease stably and continuously. Model features likewise have an important role in the magnitude of anisotropic scattering effects. The size of the source region, angular dependence of source particles, geometric model, material composition, and other problem-related quantities determine the angular flux distribution. A more anisotropic angular flux leads to a stronger scattering influence. Monoenergic point source problems in model 2 were simulated. Compared to the model 1’s results, errors of model 2 are larger and the superiority of the conservative fix-up is more notable than the economic way. Error norms for photon transport in iron are shown in Table 4 Table 4 . To investigate the accuracy in the deep-penetration region, root-mean-square errors of mesh dose rates <math id="M40"><mrow><msubsup><mi>E</mi><mrow><mi>R</mi><mi>M</mi><mi>S</mi></mrow><mrow><mi>V</mi><mi>B</mi></mrow></msubsup></mrow></math> along the right vacuum boundary are likewise given. The angular flux is more forward-peaky for point source problems, and this makes the anisotropic scattering impact stronger. The small mean-free-path of photons encourages the accumulation of errors in spite of the small scattering ratio. 1 1 Transport approximation 1 1 Expansion order 9 1 Negative source fix-up 1 1 1 1 3 1 No fix-up 3 1 Economic fix-up 3 1 Conservative fix-up 1 1 1 1 1 1 E RMS 1 1 E MAX 1 1 <math id="M41"><mrow><msubsup><mi>E</mi><mrow><mtext>RMS</mtext></mrow><mrow><mtext>VB</mtext></mrow></msubsup></mrow></math> 1 1 E RMS 1 1 E MAX 1 1 <math id="M42"><mrow><msubsup><mi>E</mi><mrow><mtext>RMS</mtext></mrow><mrow><mtext>VB</mtext></mrow></msubsup></mrow></math> 1 1 E RMS 1 1 E MAX 1 1 <math id="M43"><mrow><msubsup><mi>E</mi><mrow><mtext>RMS</mtext></mrow><mrow><mtext>VB</mtext></mrow></msubsup></mrow></math> 1 1 Consistent-P 1 1 P 1 1 1 5.79% 1 1 20.05% 1 1 11.97% 1 1 5.40% 1 1 19.42% 1 1 10.95% 1 1 8.38% 1 1 24.70% 1 1 15.57% 1 1 1 1 P 3 1 1 0.64% 1 1 2.69% 1 1 1.04% 1 1 7.23% 1 1 13.93% 1 1 8.07% 1 1 0.68% 1 1 3.90% 1 1 1.10% 1 1 1 1 P 5 1 1 0.55% 1 1 2.41% 1 1 1.03% 1 1 5.95% 1 1 11.24% 1 1 6.98% 1 1 0.49% 1 1 2.27% 1 1 0.80% 1 1 1 1 P 7 1 1 0.42% 1 1 1.64% 1 1 0.70% 1 1 4.63% 1 1 8.01% 1 1 5.58% 1 1 0.21% 1 1 1.04% 1 1 0.35% 1 1 Diagonal 1 1 P 1 1 1 4.18% 1 1 20.90% 1 1 8.69% 1 1 3.88% 1 1 20.52% 1 1 7.90% 1 1 5.95% 1 1 24.58% 1 1 11.22% 1 1 1 1 P 3 1 1 1.24% 1 1 11.35% 1 1 2.44% 1 1 4.60% 1 1 10.66% 1 1 5.11% 1 1 1.09% 1 1 6.28% 1 1 1.90% 1 1 1 1 P 5 1 1 1.14% 1 1 6.66% 1 1 2.12% 1 1 3.39% 1 1 9.07% 1 1 4.00% 1 1 1.02% 1 1 6.56% 1 1 1.84% 1 1 1 1 P 7 1 1 1.01% 1 1 6.17% 1 1 1.88% 1 1 2.43% 1 1 7.53% 1 1 2.98% 1 1 0.91% 1 1 5.91% 1 1 1.66% The effect of source particles’ angular dependence is analyzed qualitatively by monoenergic surface incident source problems, whose model is shown in Fig. 7 Fig. 7 . A model with the isotropic incident ( <math id="M44"><mrow><msub><mi>ψ</mi><mrow><mtext>in</mtext></mrow></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mn>1.0</mn></mrow></math> ) and the other with exponential incident ( <math id="M45"><mrow><msub><mi>ψ</mi><mrow><mtext>in</mtext></mrow></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mrow><mi>exp</mi><mo stretchy="false">(</mo><mn>3.0</mn><mo>×</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow><mo>/</mo><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mi>exp</mi><mo stretchy="false">(</mo><mn>3.0</mn><mo>×</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>μ</mi></mrow></mrow></mstyle></mrow></mrow></mrow></math> ) about the direction cosine of discrete directions are adopted. The results of gamma transport in iron of model 3 are given in Table 5 Table 5 , where only the consistent P approximation is used, and results with the P 9 expansion for respective problems serve as references. The effect of anisotropic scattering is stronger for the exponential incident problem, because the angular flux of this problem is more anisotropic by our modification for the angular dependence of the incident source. 1 1 Source 1 1 Expansion order 6 1 Negative source fix-up 1 1 1 1 2 1 No fix-up 2 1 Economic fix-up 2 1 Conservative fix-up 1 1 1 1 1 1 E RMS 1 1 E MAX 1 1 E RMS 1 1 E MAX 1 1 E RMS 1 1 E MAX 1 1  Isotropic incident 1 1 P 3 1 1 0.34% 1 1 1.79% 1 1 1.83% 1 1 3.52% 1 1 0.40% 1 1 2.32% 1 1 1 1 P 5 1 1 0.05% 1 1 0.38% 1 1 0.48% 1 1 1.03% 1 1 0.09% 1 1 0.59% 1 1 1 1 P 7 1 1 0.03% 1 1 0.10% 1 1 0.23% 1 1 0.39% 1 1 0.03% 1 1 0.16% 1 1 Exponential incident 1 1 P 3 1 1 0.37% 1 1 2.45% 1 1 2.04% 1 1 3.20% 1 1 0.67% 1 1 3.46% 1 1 1 1 P 5 1 1 0.06% 1 1 0.51% 1 1 0.69% 1 1 1.45% 1 1 0.14% 1 1 0.75% 1 1 1 1 P 7 1 1 0.04% 1 1 0.12% 1 1 0.24% 1 1 0.49% 1 1 0.04% 1 1 0.23% To further evaluate the conservative fix-up and no fix-up methods, approximated slab problems with the isotropic incident from left side were solved, as shown in Fig. 8 Fig. 8 . A row of meshes is constructed with the side length equal to 1 cm. Although scalar flux results may converge quickly as the expansion order increases, it is extremely difficult to obtain accurate angular flux as references for each angle under the disturbance of spatial and angular discretization. Positivity is one of desirable qualities for the S N calculation, and the ratios of the number of directions with negative angular fluxes to the total number of discrete directions in the quadrature sets were computed for each mesh and each group. The incident source only exists in the first group, and the P 3 expansion is adopted under consistent P approximation. The differences between the two methods are not large for the neutron transport in water. However, the behaviors of these methods differ greatly for gamma transport in iron, and their results of negative angular flux ratios are shown in Fig. 9 Fig. 9 . The economic fix-up has the best ability to reduce occurrences of negative angular fluxes. This is not satisfactory for the economic fix-up for problems with highly anisotropic angular flux distributions, because it overestimates fluxes despite its elimination of the negative source. Almost all angular fluxes of backward directions are negative with no source fix-up, although good integral quantities could always be gained. Owing to discretization properties and improper mesh size, many spatial discretization schemes, such as diamond differencing and the linear finite element method, would produce non-physical negative angular flux despite of the positive incident and source term. There is no doubt that the negative source aggravates the appearance of negative angular fluxes, and it may cause iteration divergence for some spatial treatments based on the exponential assumption. Therefore, source fix-ups are recommended to obtain as few as possible disturbances on transport solving. Finally, monoenergic directional point source problems were simulated to analyze the most severe case of anisotropy. Unfortunately, no reliable results could be obtained with the existing algorithms in ARES. The differences between results with P 7 and P 9 expansions remain fairly large. The Galerkin quadrature and strictly non-negative scattering cross-section may be required. The conventional Legendre expansion method cannot deal with this high degree of anisotropic problems. The algorithms used to describe anisotropic scattering affects the flux results more significantly as the proportion of collided flux within the scalar flux increases. 5 Conclusion 1 The methodology employed to handle anisotropic scattering is reviewed for deterministic transport codes. Due to the insufficiency of the Legendre expansion, non-physical negatives of scattering functions may result in a negative source of transport solving for highly anisotropic scattering media. This affects computational accuracy and exacerbates negative angular fluxes, as well as the process of iteration calculation. Convergence analyses of scattering Legendre moments and a series of numerical tests were performed in this study. An integral factor method was proposed to evaluate the precision of expanded scattering angular distributions by integrating functions of scattering moments. By checking the convergence of integral factors, the severity of negatives could be predicted. Gamma scattering exhibits much stronger anisotropy than neutron scattering, considering normalized within-group scattering functions. Numerous 2D problems were simulated using the ARES transport code, and integral errors of dose rates were compared to analyze computational accuracy. The extended and diagonal transport approximations show better efficiency than the consistent P method for a given expansion order, which uses high-order information to reduce the negatives of scattering functions. However, the extended approximation may produce negative 0-th scattering moments in some cases. Different algorithms have stronger effects on the gamma transport than on the neutron transport. The economic fix-up overestimates total scattering and scalar fluxes, especially for monoenergic photon transport problems, and errors would accumulate and rise as the penetration distance increases. No fix-up and conservative fix-up show similar integral results, however, a serious situation exists with regard to negative angular fluxes if no fix-up is adopted. Few differences are detected among results using various algorithms for problems with the fission source, because the existence of a fixed source in every group increases the ratio of uncollided flux and weakens the scattering effect. The refinement of groups certainly makes within-group scattering more forward-peaky, however this effect may not be important for integral quantities of common problems. Simulations of beam point source problems are sensitive to the anisotropic scattering and require a non-negative cross-section to handle anisotropy more accurately. The proper treatments of scattering are important for computational accuracy, especially when angular flux is highly anisotropic or scattering is of strong anisotropy. This effect strengthens as the proportion of collided flux in the scalar flux increases. The analyses in this study assist to evaluate and determine options in S N codes. The Galerkin quadrature with the non-negative cross-section will be studied and implemented in the ARES code for beam source shielding problems in the future. Further comparisons will be performed for realistic engineering problems. References G.I. Bell , S. Glasstone , Nuclear reactor theory , US Atomic Energy Commission , Washington DC , 1970 . Nuclear reactor theory US Nuclear Regulatory Commission . Regulatory Guide 1.190: Calculational and Dosimetry Methods for Determining Pressure Vessel Neutron Fluence . US Nuclear Regulatory Commission, Office of Nuclear Regulatory Research , 2001 Regulatory Guide 1.190: Calculational and Dosimetry Methods for Determining Pressure Vessel Neutron Fluence J.P. Odom , J.K. Shultis , Anisotropic neutron transport without Legendre expansions , Nucl. Sci. Eng . 59 : 278 - 281 ( 1976 ). DOI: 10.13182/NSE76-5 http://doi.org/10.13182/NSE76-5 Anisotropic neutron transport without Legendre expansions J.A. Dahl , B.D. Ganapol , J.E. Morel , Positive Scattering Cross Sections Using Constrained Least Squares , in: Proceedings of ANS Mathematics and Computation Topical Meeting , 1999 . Positive Scattering Cross Sections Using Constrained Least Squares J.M. DelGrande , K.A. Mathews , Nonnegative anisotropic group cross sections: A hybrid Monte Carlo-discrete elements-discrete ordinates approach . Nucl. Sci. Eng . 139 : 33 - 46 ( 2001 ). DOI: 10.13182/NSE00-69 http://doi.org/10.13182/NSE00-69 Nonnegative anisotropic group cross sections: A hybrid Monte Carlo-discrete elements-discrete ordinates approach D.W. Gerts , K.A. Mathews , Non-negative Anisotropic Piecewise-Average Multigroup Cross Sections , in: Proceedings of International Conference on Mathematics, Computational Methods & Reactor Physics , Gatlinburg, Tennessee . 2003 . Non-negative Anisotropic Piecewise-Average Multigroup Cross Sections J.W. Kim , N.Z. Cho , An efficient deterministic method for generating non-negative scattering cross-sections , Ann. Nucl. Energy 34 : 967 - 976 ( 2007 ). DOI: 10.1016/j.anucene.2007.04.014 http://doi.org/10.1016/j.anucene.2007.04.014 An efficient deterministic method for generating non-negative scattering cross-sections G.I. Bell , G.E. Hansen , H.A. Sandmeier , Multitable treatments of anisotropic scattering in SN multigroup transport calculations , Nucl. Sci. Eng . 28 : 376 - 383 ( 1967 ). DOI: 10.13182/NSE67-2 http://doi.org/10.13182/NSE67-2 Multitable treatments of anisotropic scattering in SN multigroup transport calculations H. Brockmann , Treatment of anisotropic scattering in numerical neutron transport theory , Nucl. Sci. Eng , 77 : 377 - 414 ( 1981 ). DOI: 10.13182/NSE81-3 http://doi.org/10.13182/NSE81-3 Treatment of anisotropic scattering in numerical neutron transport theory M.B. Emmett , R.L. Childs , W.A. Rhoades , Repair for scattering expansion truncation errors in transport calculations , Oak Ridge Natl. Lab ., 1980 . No. CONF-791103-100 . Repair for scattering expansion truncation errors in transport calculations J.E. Morel , A hybrid collocation-Galerkin-SN method for solving the Boltzmann transport equation , Nucl. Sci. Eng . 101 : 72 - 87 ( 1989 ). DOI: 10.13182/NSE89-4 http://doi.org/10.13182/NSE89-4 A hybrid collocation-Galerkin-SN method for solving the Boltzmann transport equation R. Sanchez , J. Ragusa , On the construction of Galerkin angular quadratures , Nucl. Sci. Eng . 169 ( 2011 ): 133 - 154 . DOI: 10.13182/NSE10-31 http://doi.org/10.13182/NSE10-31 On the construction of Galerkin angular quadratures T. Ushio , T. Takeda , M. Mori , Neutron anisotropic scattering effect in heterogeneous cell calculations of light water reactors , J. Nucl. Sci. Technol . 40 : 464 - 480 ( 2003 ). DOI: 10.1080/18811248.2003.9715381 http://doi.org/10.1080/18811248.2003.9715381 Neutron anisotropic scattering effect in heterogeneous cell calculations of light water reactors T. Takeda , T. Okamoto , A. Inoue , et al ., Effect of anisotropic scattering in neutronics analysis of BWR assembly , Ann. Nucl. Energy 33 : 1315 - 1321 ( 2006 ). DOI: 10.1016/j.anucene.2006.08.008 http://doi.org/10.1016/j.anucene.2006.08.008 Effect of anisotropic scattering in neutronics analysis of BWR assembly A. Yamamoto , Y. Kitamura , Y. Yamane , Simplified treatments of anisotropic scattering in LWR core calculations , J. Nucl. Sci. Technol . 45 : 217 - 229 ( 2008 ). doi: 10.1080/18811248.2008.9711430 http://doi.org/10.1080/18811248.2008.9711430 Simplified treatments of anisotropic scattering in LWR core calculations S. Choi , K. Smith , H.C. Lee , et al ., Impact of inflow transport approximation on light water reactor analysis , J. Comput. Phys . 299 : 352 - 373 ( 2015 ). DOI: 10.1016/j.jcp.2015.07.005 http://doi.org/10.1016/j.jcp.2015.07.005 Impact of inflow transport approximation on light water reactor analysis E. Oblow , K. Kin , H Goldstein , et al ., Effects of highly anisotropic scattering on monoenergetic neutron transport at deep penetrations , Nucl. Sci. Eng . 54 : 72 - 84 ( 1974 ). Effects of highly anisotropic scattering on monoenergetic neutron transport at deep penetrations D.C. Sahni , R.G. Tureci , Discrete eigenvalues of Case spectrum with anisotropic scattering , Nucl. Sci. Eng . 191 : 1 - 15 ( 2018 ). DOI: 10.1080/00295639.2018.1463748 http://doi.org/10.1080/00295639.2018.1463748 Discrete eigenvalues of Case spectrum with anisotropic scattering R.E. Macfarlane , D.W. Muir , R.M. Boicourt , et al ., The NJOY Nuclear Data Processing System Version 2016 , Los Alamos Natl. Lab. , 2017 . No. LA-UR-17-20093 . The NJOY Nuclear Data Processing System Version 2016 Y. Chen , B. Zhang , L. Zhang et al ., ARES: a parallel discrete ordinates transport code for radiation shielding applications and reactor physics analysis , Science and Technology of Nuclear Installations , 2017 , Article ID 2596727 . doi: 10.1155/2017/2596727 http://doi.org/10.1155/2017/2596727 ARES: a parallel discrete ordinates transport code for radiation shielding applications and reactor physics analysis L. Zhang , B. Zhang , C. Liu , et al ., Evaluation of PWR pressure vessel fast neutron fluence benchmarks from NUREG/CR-6115 with ARES transport code , Nuclear Technology and Radiation Protection 32 : 204 - 210 ( 2017 ). Evaluation of PWR pressure vessel fast neutron fluence benchmarks from NUREG/CR-6115 with ARES transport code B. Zhang , L. Zhang , C. Liu , et al ., Goal-Oriented regional angular adaptive algorithm for the SN equations , Nucl. Sci. Eng . 189 : 120 - 134 ( 2018 ). DOI: 10.1080/00295639.2017.1394085 http://doi.org/10.1080/00295639.2017.1394085 Goal-Oriented regional angular adaptive algorithm for the SN equations R.E. MacFarlane , TRANSX 2: A code for interfacing MATXS cross-section libraries to nuclear transport codes , Los Alamos Natl. Lab ., 1992 . TRANSX 2: A code for interfacing MATXS cross-section libraries to nuclear transport codes B.G. Petrovic , A. Haghighat . Effects of SN method numerics on pressure vessel neutron fluence calculations . Nucl. Sci. Eng . 122 : 167 - 193 ( 1996 ). DOI: 10.13182/NSE96-3 http://doi.org/10.13182/NSE96-3 Effects of SN method numerics on pressure vessel neutron fluence calculations 10.1007/s41365-019-0695-y This work was supported by the National Natural Science Foundation of China (Nos. 11505059, 11575061) and the Fundamental Research Funds for Central Universities (No. 2017XS087). 2018-11-26 2019-05-25 2019-07-03 2019-10-21 2019-11-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2019 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-019-0695-y http://dx.doi.org/10.1007/s41365-019-0695-y . 2019 Nuclear Science and Techniques 11 30