Finite element spherical harmonics (PN) solutions of the three-dimensional Takeda benchmark problems

A set of multi-group eigenvalue (K_eff) benchmark problems in three-dimensional homogenised reactor core configurations have been solved using the deterministic finite element transport theory code EVENT and the Monte Carlo code MCNP4C. The principal aim of this work is to qualify numerical methods and algorithms implemented in EVENT. The benchmark problems were compiled and published by the Nuclear Data Agency (OECD/NEACRP) and represent three-dimensional realistic reactor cores which provide a framework in which computer codes employing different numerical methods can be tested. This is an important step that ought to be taken (in our view) before any code system can be confidently applied to sensitive problems in nuclear criticality and reactor core calculations. This paper presents EVENT diffusion theory (P₁) approximation to the neutron transport equation and spherical harmonics transport theory solutions (P₃–P₉) to three benchmark problems with comparison against the widely used and accepted Monte Carlo code MCNP4C. In most cases, discrete ordinates transport theory (S_N) solutions which are already available and published have also been presented. The effective multiplication factors (K_eff) obtained from transport theory EVENT calculations using an adequate spatial mesh and spherical harmonics approximation to represent the angular flux for all benchmark problems have been estimated within 0.1% (100 pcm) of the MCNP4C predictions. All EVENT predictions were within the three standard deviation uncertainty of the MCNP4C predictions. Regionwise and pointwise multi-group neutron scalar fluxes have also been calculated using the EVENT code and compared against MCNP4C predictions with satisfactory agreements. As a result of this study, it is shown that multi-group reactor core/criticality problems can be accurately solved using the three-dimensional deterministic finite element spherical harmonics code EVENT.     

An exponential spectral nodal method for one-speed X,Y-geometry deep penetration discrete ordinates problems

Presented here is an exponential spectral nodal method applied to deep penetration X,Y-geometry heterogeneous neutron transport problems in the discrete ordinates (S_N) formulation. This numerical method uses the spectral Green's function (SGF) scheme for solving the one-dimensional transverse-integrated S_N exponential nodal equations with no spatial truncation error. Based on the physics of deep penetration problems, we approximate the transverse leakage terms by exponential functions. We show in two numerical experiments that the SGF-exponential nodal method (SGF-ExpN) generates very accurate results when compared to the conventional transport nodal methods for coarse-mesh deep penetration S_N problems, specially in highly absorbing media.     

Evaluation of dominant time-eigenvalues of neutron transport equation by Meyer’s sub-space iterations

The aim of this paper is to explore the use of Meyer’s sub-space iteration (SSI) method for the evaluation of dominant prompt time-eigenvalues of the neutron transport equation. The integro-differential form of the transport equation is considered. The SSI method is known to be an efficient technique to find the dominant eigenvalues of a non-symmetric matrix. It has been earlier used for eigenvalue problems in neutron diffusion theory. However, it does not seem to be tried in the transport theory case. Here, the use of SSI has been tested in transport theory for some 1-D mono-energetic homogeneous and heterogeneous benchmark problems. The space variable is discretised by finite differencing while neutron directions are discretised by discrete ordinates (S_n-) method. The SSI method needs frequent multiplication of the relevant matrix operator with vectors. As known from earlier works in this area, this can be achieved in terms of external source calculations for which a 1-D programme was developed and used. With the availability of more versatile S_n-method codes, it may perhaps be possible to extend use of SSI to more realistic cases.     

Axial SPN and Radial MOC Coupled Whole Core Transport Calculation

The Simplified P_N (SP_N) method is applied to the axial solution of the two-dimensional (2-D) method of characteristics (MOC) solution based whole core transport calculation. A sub-plane scheme and the nodal expansion method (NEM) are employed for the solution of the one-dimensional (1-D) SP_N equations involving a radial transverse leakage. The SP_N solver replaces the axial diffusion solver of the DeCART direct whole core transport code to provide more accurate, transport theory based axial solutions. In the sub-plane scheme, the radial equivalent homogenization parameters generated by the local MOC for a thick plane are assigned to the multiple finer planes in the subsequent global three-dimensional (3-D) coarse mesh finite difference (CMFD) calculation in which the NEM is employed for the axial solution. The sub-plane scheme induces a much less nodal error while having little impact on the axial leakage representation of the radial MOC calculation. The performance of the sub-plane scheme and SP_N nodal transport solver is examined by solving a set of demonstrative problems and the C5G7MOX 3-D extension benchmark problems. It is shown in the demonstrative problems that the nodal error reaching upto 1,400 pcm in a rodded case is reduced to 10pcm by introducing 10 sub-planes per MOC plane and the transport error is reduced from about 150pcm to 10pcm by using SP₃. Also it is observed, in the C5G7MOX rodded configuration B problem, that the eigenvalues and pin power errors of 180 pcm and 2.2% of the 10 sub-planes diffusion case are reduced to 40 pcm and 1.4%, respectively, for SP₃ with only about a 15% increase in the computing time. It is shown that the SP₅ case gives very similar results to the SP₃ case.     

A new spectral nodal method based on Larsen’s extended diamond approach

In this article, we describe a new spectral nodal method for solving discrete ordinates (S_N) neutron transport problems with anisotropic scattering for arbitrary order N of angular quadrature. The key to our new spectral nodal method is a consistent derivation of nonstandard auxiliary equations that relate angular neutron fluxes only in the upwind directions. These nonstandard equations are angularly coupled extensions of very basic auxiliary equations proposed by Edward W. Larsen in his extended diamond scheme of solving S₂ problems in the presence of scattering and free from spatial truncation error. The resulting method here is also free from spatial truncation error and, in contrast to previously developed spectral nodal methods, it is compatible with an efficient use of iteration on the scattering source and is free from the storage of cell-edge angular fluxes.     

Cellular neural network to the spherical harmonics approximation of neutron transport equation in x–y geometry. Part I: Modeling and verification for time-independent solution

This paper describes a novel method based on using cellular neural networks (CNN) coupled with spherical harmonics method (P_N) to solve the time-independent neutron transport equation in x–y geometry. To achieve this, an equivalent electrical circuit based on second-order form of neutron transport equation and relevant boundary conditions is obtained using CNN method. We use the CNN model to simulate spatial response of scalar flux distribution in the steady state condition for different order of spherical harmonics approximations. The accuracy, stability, and capabilities of CNN model are examined in 2D Cartesian geometry for fixed source and criticality problems.     

Solution of the multiplying binary stochastic media based on L-P equation in 1D spherical geometry

In order to obtain the ensemble average k_eff for the binary stochastic media system, a statistical transport equation with eigenvalue is derived based on the L-P equation. Combined method of statistical and deterministic is proposed to deal with this problem. In the first step, Monte Carlo approach is employed to calculate the mean chord length of the background scattering material, deterministic transport method using diamond difference and source iteration is applied to solve the eigenvalue L-P equation. Test problems of different scattering ratios, different chord-path ratios and different number of random fissile lumps and simple UO₂-H₂O problem are calculated and compared with references. Results show that the eigenvalue L-P equation can provide the mean value of k_eff for the given stochastic systems in most cases.     

Solution of some 2D transport problems by a high order AN–SP2N−1 method

A collection of classical 2D transport problems (the escape probability from prisms of various shapes, the current-to-flux ratio of a wedge-shaped reflector, the transport and asymptotic flux as well as the extrapolation length near a corner) are solved by means of the boundary element version of a high order A_N method, an equivalent form of the odd order simplified spherical harmonics (SP_2N−1) method. The use of a high order approximation is motivated by the fact that all the above problems can be made to fulfil the condition of constant total mean free path, which makes A_N–SP_2N−1 to be equivalent, in turn, to the classical odd order spherical harmonics (P_2N−1) method, so that for these problems A_N–SP_2N−1 shares with the latter method the property that, by increasing the order 2N − 1, the error can be made as small as we want. A second purpose of the paper is to show that the boundary element approach can handle such highly singular boundary integrals as those implied by the partial derivatives of the asymptotic flux at the boundary.     

On the integration scheme along a trajectory for the characteristics method

The issue of the integration scheme along a trajectory which appears for all tracking-based transport methods is discussed from the point of view of the method of characteristics. The analogy with the discrete ordinates method in slab geometry is highlighted along with the practical limitation in transposing high-order S_N schemes to a trajectory-based method. We derived an example of such a transposition starting from the linear characteristic scheme. This new scheme is compared with the standard flat-source approximation of the step characteristic scheme and with the diamond differencing scheme. The numerical study covers a 1D analytical case, 2D one-group critical and fixed-source benchmarks and finally a realistic multigroup calculation on a BWR-MOX assembly.     

The criticality problems with the FN method for the FBIS model

In one-speed, time-independent, neutron transport theory, the F_N method is used for the FBIS (forward-backward-isotropic scattering) model to reinvestigate the behaviour of the critical size in plane and spherical geometries. For the FIS (forward-isotropic scattering) model the numerical results are compared with previously obtained variational results and it is shown that they are in agreement. For the BIS (backward-isotropic scattering) model exact results are obtained and compared with the first-order approximate results obtained using the method of elementary solutions.     

Coolant void reactivity adjustments in advanced CANDU lattices using adjoint sensitivity technique

Coolant void reactivity (CVR) is an important factor in reactor accident analysis. Here we study the adjustments of CVR at beginning of burnup cycle (BOC) and k_eff at end of burnup cycle (EOC) for a 2D Advanced CANDU Reactor (ACR) lattice using the optimization and adjoint sensitivity techniques. The sensitivity coefficients are evaluated using the perturbation theory based on the integral neutron transport equations. The neutron and flux importance transport solutions are obtained by the method of cyclic characteristics (MOCC). Three sets of parameters for CVR-BOC and k_eff-EOC adjustments are studied: (1) Dysprosium density in the central pin with Uranium enrichment in the outer fuel rings, (2) Dysprosium density and Uranium enrichment both in the central pin, and (3) the same parameters as in the first case but the objective is to obtain a negative checkerboard CVR-BOC (CBCVR-BOC). To approximate the EOC sensitivity coefficient, we perform constant-power burnup/depletion calculations using a slightly perturbed nuclear library and the unperturbed neutron fluxes to estimate the variation of nuclide densities at EOC. Our aim is to achieve a desired negative CVR-BOC of −2 mk and k_eff-EOC of 0.900 for the first two cases, and a CBCVR-BOC of −2 mk and k_eff-EOC of 0.900 for the last case. Sensitivity analyses of CVR and eigenvalue are also included in our study.     

Linear and quadratic octahedral wavelets on the sphere for angular discretisations of the Boltzmann transport equation

In this paper, two new wavelet bases are developed for discretising the angular term of the first-order Boltzmann transport equation. The wavelets proposed are based on Sweldens second generation wavelets [Sweldens, W., 1993. The lifting scheme: a construction of second generation wavelets. SIAM J. Math. 1, 54], which are constructed through the lifting procedure [Sweldens, W., 1995. The lifting scheme: a new philosophy in biorthogonal wavelet construction. Wavelet Applications in Signal and Image Processing III]. In this paper, the wavelets are built on an octahedral domain, Fig. 2, and the angular flux approximation takes the form of finite element linear and quadratic representations. Full details of the meshing over the octahedron and derivation of the wavelet functions are given. The wavelets discussed are similar to the wavelets developed in Buchan [Buchan, A., 2003 c. Angular discretisation of the first order Boltzmann transport equation. Part 2: linear spherical wavelets. Technical Report, Imperial College, London, Dep. Earth Sci. Eng.] and [Buchan, A., 2003b. Angular discretisation of the first order Boltzmann transport equation. Part 3: quadratic spherical wavelets. Technical Report, Imperial College, London, Dep. Earth Sci. Eng.], in this paper the bases use a new fundamental amendment for mitigating the inaccuracies observed with the earlier bases. The performance of the new angular discretisation techniques are demonstrated using 2 one-dimensional and 4 two-dimensional test problems. These problems demonstrate the accuracy and susceptibility to ray effects of the proposed methods. Comparisons of all calculations are made with the conventional S_N and P_N approximations. Benchmark solutions are provided by the established code EVENT.     

The finite-element method for multigroup neutron transport: Anisotropic scattering in 1-D slab geometry

Proof-tests on 1-D multigroup neutron transport problems are reported for strong anisotropic scattering. These tests have been undertaken as part of the validation of the 3-D multigroup finite-element transport code fel tran for ansisotropic scattering media. To illustrate the treatment of within-group and intergroup anisotropic scattering in the finite-element method the relevant theory is outlined. Ingroup scattering is checked using the backward-forward-isotropic (BFI) scattering law for source and eigenvalue problems. With this law anisotropic scattering problems can be transformed into equivalent isotropic scattering problems. In this way the well-validated isotropic scattering version of fel tran is used to validate the anisotropic version. Intergroup scattering effects are checked by solving few-group source problems for P₁ and P₃ scattering and the BFI scattering law. For P₁ and P₃ scattering checks are made with the discrete-ordinate finite-difference code anisn and the spherical harmonics finite-difference code marc/pn. For the BFI scattering law comparison is made with two-group exact solutions of Williams (1985) for 1-D systems.     

Solution of one-speed neutron transport equation for strongly anisotropic scattering by TN approximation: Slab criticality problem

In this study, a recently proposed version of Chebyshev polynomial approximation which was used in spectrum and criticality calculations by one-speed neutron transport equation for slabs with isotropic scattering is further developed to slab criticality problems for strongly anisotropic scattering. Backward–forward-isotropic model is employed for the scattering kernel which is a combination of linearly anisotropic and strongly backward–forward kernels. Further to that, the common approaches of using the same functional form for scattering and fission kernels or embedding fission kernel into the scattering kernel even in strongly anisotropic scattering is questioned for T_N approximation via taking an isotropic fission kernel in the transport equation. As a starting point, eigenvalue spectrum of one-speed neutron transport equation for a multiplying slab with different degrees of anisotropy in scattering and for different cross-section parameters is obtained using Chebyshev method. Later on, the spectra obtained for different degree of anisotropies and cross-section parameters are made use of in criticality problem of bare homogeneous slab with strongly anisotropic scattering. Calculated critical thicknesses by Chebysev method are almost in complete agreement with literature data except for some limiting cases. More importantly, it is observed that using a different kernel (isotropic) for fission rather than assuming it equal to the scattering kernel which is a more realistic physical approach yields in deviations in critical sizes in comparison with the values presented in literature. This separate kernel approach also eliminates the slow convergency and/or non-convergent behavior of high-order approximations arising from unphysical eigenspectrum calculations.     

The Davidson method as an alternative to power iterations for criticality calculations

The Davidson method is implemented within the neutron transport core solver parafish to solve k-eigenvalue criticality transport problems. The parafish solver is based on domain decomposition. It uses spherical harmonics (P_N method) for angular discretization, and non-conforming finite elements for spatial discretization. The Davidson method is compared to the traditional power iteration method in this context.     

An overview of the Boltzmann transport equation solution for neutrons, photons and electrons in Cartesian geometry

Questions regarding accuracy and efficiency of deterministic transport methods are still on our mind today, even with modern supercomputers. The most versatile and widely used deterministic methods are the P_N approximation, the S_N method (discrete ordinates method) and their variants. In the discrete ordinates (S_N) formulations of the transport equation, it is assumed that the linearised Boltzmann equation only holds for a set of distinct numerical values of the direction-of-motion variables. In this paper, looking forward to confirm the capabilities of deterministic methods in obtaining accurate results, we describe the recent advances in the class of deterministic methods applied to one and two dimensional transport problems for photons and electrons in Cartesian Geometry. First, we describe the Laplace transform technique applied to S_N two dimensional transport equation in a rectangular domain considering Compton scattering. Next, we solved the Fokker-Planck (FP) equation, an alternative approach for the Boltzmann transport equation, assuming a mono-energetic electron beam in a rectangular domain. The main idea relies on applying the P_N approximation, a recent advance in the class of deterministic methods, in the angular variable, to the two dimensional Fokker-Planck equation and then applying the Laplace Transform in the spatial x variable. Numerical results are given to illustrate the accuracy of deterministic methods presented.     

An overview of the Boltzmann transport equation solution for neutrons,photons and electrons in Cartesian geometry

Questions regarding accuracy and efficiency of deterministic transport methods are still on our mind today, even with modern supercomputers. The most versatile and widely used deterministic methods are the P_N approximation, the S_N method (discrete ordinates method) and their variants. In the discrete ordinates (S_N) formulations of the transport equation, it is assumed that the linearised Boltzmann equation only holds for a set of distinct numerical values of the direction-of-motion variables. In this paper, looking forward to confirm the capabilities of deterministic methods in obtaining accurate results, we describe the recent advances in the class of deterministic methods applied to one and two dimensional transport problems for photons and electrons in Cartesian Geometry. First, we describe the Laplace transform technique applied to S_N two dimensional transport equation in a rectangular domain considering Compton scattering. Next, we solved the Fokker–Planck (FP) equation, an alternative approach for the Boltzmann transport equation, assuming a mono-energetic electron beam in a rectangular domain. The main idea relies on applying the P_N approximation, a recent advance in the class of deterministic methods, in the angular variable, to the two dimensional Fokker–Planck equation and then applying the Laplace Transform in the spatial x variable. Numerical results are given to illustrate the accuracy of deterministic methods presented.     

A nodal SN transport method for three-dimensional triangular-z geometry

A new nodal S_N transport method has been developed to perform accurate transport calculation in three-dimensional triangular-z geometry, where arbitrary triangles are transformed into regular triangles via a coordinate transformation. The transverse integration procedure is applied to treat the neutron transport equation in the regular triangle. The neutron angular distributions of intra-node fluxes are represented using the S_N quadrature set, and the spatial distributions of neutron fluxes and sources are approximated by a quadratic polynomial. The nodal-equivalent finite difference algorithm for 3D triangular geometry is applied to establish a stable and efficient iterative scheme. The present method was tested on four 3D Takeda benchmark problems published by the nuclear data agency (NEACRP), in which the first three problems are in XYZ geometry and the last one is in hexagonal-z geometry. The results of the present method agree well with those of the reference Monte-Carlo calculation method, the difference in k_eff being less than 0.1%. This shows that multi-group reactor core/criticality problems can be accurately and effectively solved using the present method.     

Two-dimensional nodal transport method for triangular geometry

The advanced nodal method for solving the multi-group neutron transport equation in two-dimensional triangular geometry is developed. To apply the transverse integration procedure, an arbitrary triangular node is transformed into a regular triangular node using coordinate transformation. The angular distributions of intra-node neutron fluxes and its transverse-leakage are represented by the S_N quadrature set. The spatial distributions of neutron flux and source in the regular triangle are given approximately by an orthogonal quadratic polynomial, and the spatial expansion of transverse-leakage is approximated by a second-order polynomial. To establish a stable and efficient iterative scheme, the improved nodal-equivalent finite difference algorithm is used. The results for several benchmark problems demonstrate the higher capability of the method to yield the accurate results in significantly smaller computing times than those required by the standard finite difference method and the finite element spherical-harmonics method.     

A scheme for the evaluation of dominant time-eigenvalues of a nuclear reactor

This paper presents a scheme to obtain the fundamental and few dominant solutions of the prompt time eigenvalue problem (referred to as α-eigenvalue problem) for a nuclear reactor using multi-group neutron diffusion theory. The scheme is based on the use of an algorithm called Orthomin(1). This algorithm was originally proposed by Suetomi and Sekimoto [Suetomi, E., Sekimoto, H., 1991. Conjugate gradient like methods and their application to eigenvalue problems for neutron diffusion equations. Ann. Nucl. Energy 18 (4), 205–227] to obtain the fundamental K-eigenvalue (K-effective) of nuclear reactors. Recently, it has been shown that the algorithm can be used to obtain the further dominant K-modes also. Since α-eigenvalue problem is usually more difficult to solve than the K-eigenvalue problem, an attempt has been made here to use Orthomin(1) for its solution. Numerical results are given for realistic 3-D test case.