Application of Finite Element Method to Two-Dimensional Multi-Group Neutron Transport Equation in Cylindrical Geometry

The finite element method is applied to the spatial variables of multi-group neutron transport equation in the two-dimensional cylindrical (r, z) geometry. The equation is discretized using regular rectangular subregions in the (r, z) plane. The discontinuous method with bilinear or biquadratic Lagrange's interpolating polynomials as basis functions is incorporated into a computer code FEMRZ. Here, the angular fluxes are allowed to be discontinuous across the subregion boundaries.

Some numerical calculations have been performed and the results indicated that, in the case of biquadratic approximation, the solutions are sufficiently accurate and numerically stable even for coarse meshes. The results are also compared with those obtained by a diamond difference S _n code TWOTRAN-II. The merits of the discontinuous method are demonstrated through the numerical studies.     

A spatially adaptive grid-refinement approach for the finite element solution of the even-parity Boltzmann transport equation

A spatially adaptive grid-refinement approach has been investigated to solve the even-parity Boltzmann transport equation. A residual based a posteriori error estimation scheme has been utilized for checking the approximate solutions for various finite element grids. The local particle balance has been considered as an error assessment criterion. To implement the adaptive approach, a computer program ADAFENT (adaptive finite elements for neutron transport) has been developed to solve the second order even-parity Boltzmann transport equation using K⁺ variational principle for slab geometry. The program has a core K⁺ module which employs Lagrange polynomials as spatial basis functions for the finite element formulation and Legendre polynomials for the directional dependence of the solution. The core module is called in by the adaptive grid generator to determine local gradients and residuals to explore the possibility of grid refinements in appropriate regions of the problem. The a posteriori error estimation scheme has been implemented in the outer grid refining iteration module. Numerical experiments indicate that local errors are large in regions where the flux gradients are large. A comparison of the spatially adaptive grid-refinement approach with that of uniform meshing approach for various benchmark cases confirms its superiority in greatly enhancing the accuracy of the solution without increasing the number of unknown coefficients. A reduction in the local errors of the order of 10² has been achieved using the new approach in some 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.     

Applications of the multidimensional equations to complex fuel assembly problems

For the analysis of reactors with complex fuel assemblies or fine mesh applications as pin by pin neutron flux reconstruction, the usual approximation of the neutron transport equation by the multigroup diffusion equation does not provide good results. A classical approach to solve the neutron transport equation is to apply the spherical harmonics method obtaining a finite approximation known as the P_L equations. In this line, a nodal collocation method for the discretization of these equations on a rectangular mesh is used in this paper to analyse reactors with MOX fuel assemblies. Although the 3D P_L nodal collocation method becomes feasible due to the improvements in computer hardware, a complete treatment of the detailed structure of the fuel assemblies in actual three-dimensional geometry is still prohibitive, thus, an assembly homogenization method is necessary. A homogenization method compatible with our multidimensional P_L code is proposed and tested performing heterogeneous and homogenized calculations. In this work, we apply the method to 2D complex fuel assembly configurations.     

Coarse-Mesh Three-Dimensional Transport Calculation Method

A coarse-mesh 3-D (X-Y-Z, Hexagonal-Z) discrete ordinates transport calculation method has been developed. This method employs an weighted diamond difference approximation, the weight in which is a function of neutron direction and scalar flux, and can be easily incorpolated into conventional discrete ordinates transport codes. Results obtained in four-group S₄P₀ calculations on simple fast reactors indicate that, though the computer time of this method has increased by 30–50% compared with that of the conventional finite difference method using the same mesh, the errors of k_eff and the power distribution are reduced remarkably.     

A multigroup finite-element solution of the neutron transport equation—I: X-Y geometry

The result of extending a variational finite element method of solving the neutron transport equation, to energy dependence, is reported. Detailed results are given, in the form of tables and graphs, of P₁ and higher-order transport solutions to a number of benchmark problems in X-Y geometry. The accuracy and flexibility of the method are demonstrated. Some suggestions are made for the future development of the computer implementation of the method.     

Legendre expansion in multilayer slab geometry

The monoenergetic integral transport equation for a multilayer slab geometry has been solved by the Legendre expansion method. The method utilizes an expansion of the flux density in each layer in Legendre polynomials of the position co-ordinate. The use of these polynomials makes it possible to calculate most of the resulting matrix by means of recurrence formulae. These formulae have been obtained by a procedure which is an extension of Carlvik's method for a homogeneous slab. A code (MULREG) has been written for this purpose. Using MULREG a series of calculations have been performed for a homogeneous slab with vacuum boundary conditions. The slab has been divided into NR number of regions and in each region flux is expanded into Legendre polynomials of order NSI. For a particular value of NR, the NSI varies from 0 to 4. The effective multiplication factor k_eff of the slab is calculated. By comparing the computational time for all the cases, it is studied as how severe it is to consider the flat flux approximation (conventional collision-probability approach) as compared to a case when more terms in the flux expansion are considered per layer.     

Transport calculations of neutron energy spectra in a graphite cylinder with a D-T source

Neutron energy spectra resulting from the transport of 14.7 MeV neutrons from a collimated D-T source through a graphite cylinder, have been calculated with the discrete-ordinates 2-D transport dot 4.2 code, with multigroup cross-sections generated using the njoy code from the ENDF/B (IV & V) libraries. The results confirm the conclusion of Goldfeld et al. (1985), that energy spectra at mesh points close to the axis of the system, in front of the collimated beam, consist mainly of one-collision contributions of elastically or inelastically scattered neutrons. Investigation of the dependence of the calculated spectra on the order of truncation of the Legendre polynomials expansion of the flux and of the cross-sections (i.e. the order of scattering) leads to the following observations:

• 1.(a) the P₆ or P₇ approximations seem to be adequate enough for flux calculations, with less than 3% error, in spite of the high degree of the source and the cross-section's anisotropy;
• 2.(b) the calculation error is reduced significantly by increasing the order of scattering from P₄ to P₇, mainly in mesh points close to the axis and of those energies in which the anisotropy of the elastic and discrete level inelastic scattering processes is most pronounced.

Finally, the dot 4.2 calculations are compared with Monte Carlo mcnp calculations; both calculated spectra are in a good agreement.     

Spatially adaptive eigenvalue estimation for the SN equations on unstructured triangular meshes

Adaptive mesh refinement is a powerful method for efficiently solving physical problems described by partial differential equations at reduced computational cost. In this paper we present a new adaptive algorithm for estimating the effective multiplication factor, k_eff, for the neutron transport equation. The method is based on a dual weighted residual approach where an appropriate adjoint problem is solved to obtain the importance of residual errors to the multiplication factor. The forward residuals and the importance are then combined to give accurate error measures which are used to design economical finite element meshes. We illustrate the effectiveness of our methods by applying them to two 2-dimensional reactor problems by comparing the quality of the error estimator which is the basis for adaptivity and the overall efficiency which is judged by the number of elements required for a given accuracy.     

An Improvement of the Transverse Leakage Treatment for the Nodal SN Transport Calculation Method in Hexagonal-Z Geometry

To overcome the divergent behavior of the NSHEX code, a nodal S_N code for hexagonal geometry, for some transport calculations, an improvement has been made in the calculation of the axial leakage. The axial leakage, previously calculated by using the quadratic transverse leakage approximation (QLA), is calculated by a new method of analytically treating the spatial distribution within a node, based on the axial homogeneity of the ordinary core. The verification tests were performed for the KNK-II model geometry of the NEACRP 3-D Neutron Transport Benchmarks and the large assembly-size KNK-II model. It is found that k_ett values obtained by introducing the new method agree with the reference Monte Carlo calculation results within 0.1% Δk/k for the KNK-II model, although the QLA method did not converge for two cases. Furthemore the new method succeeded in calculations for the large assembly-size model, in which the QLA method failed for all cases. Thus the new method has been found accurate and convergence achieved for the cases in which the QLA method failed.     

Finite element analysis of crack propagation in three-dimensional solids under cyclic loading

This paper is concerned with an efficient computational procedure for analyzing crack propagation in solids. The method is general; however, its application to semi-elliptical surface cracks in thick plates is discussed in particular. The strain energy release rate G for a crack in mode I is a function of the crack geometry, the direction of crack propagation and the state of loading. When G is known, the stress intensity factor K_I can easily be obtained. In this paper the strain energy of the plate is computed numerically for a wide range of crack geometries using the finite element method. A 20-node isoparametric solid element is employed in modelling the structure. Certain special techniques for increasing the computational efficiency of the method, such as multilevel subdivision of the structure (substructuring) and condensation of degrees of freedom that are not needed in the crack propagation analysis, are emphasized. In fact, analysis of a large number of crack geometries requires only insignificantly more computational efforts than treating a single crack. Certain other aspects of the finite element modelling are also discussed.Two methods for replacing the computed discrete values of strain energy by continuous functions are presented. These functions are expressed in terms of the two half-axes defining the geometry of the elliptical crack and they are determined using a least square technique. G and K_I are easily deduced from these functions. As an example, a semi-elliptical, part-through, surface crack in a thick nickel steel plate is analyzed. The crack is subjected to a combination of axial and bending loading, applied cyclically. From the finite element calculations of the strain energy and the stress intensity factors which are computed accordingly, crack propagation along the two half-axes of the ellipse is calculated by utilization of a formula suggested by Paris. The results are checked against laboratory fatigue tests. The method has proved to be very efficient and accurate, and due to its generality it can also be applied to complicated geometries and complex states of loading.     

Application of ‘natural coordinate system’ in the finite element solution of multigroup neutron diffusion equation

A simple method of generating stiffness matrices for the solution of multigroup diffusion equation by ‘natural coordinate system’ has been presented. A comparative study has been made using triangular elements with linear model, triangular elements with quadratic model and rectangular elements with bilinear model to demonstrate their relative efficiencies. The quadratic interpolation model has been shown to be superior to linear and bilinear models with respect to computing time, computer storage and relative error in K_eff for a two group diffusion example. The flexibility of the finite element treatment has been demonstrated by the calculation of the reactivity of a partially inserted control rod. Good agreement has been obtained with a perturbation calculation.     

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.     

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.     

Analysis of Novo-Voronezh NPP unit 3 radiation lifetime of the WWER-440 reactor pressure vessel given samplings

The analysis of the residual radiation lifetime of the Novo-Voronezh NPP Unit 3 reactor pressure vessel which had spherical samplings after annealing was performed for the spectrum of the ‘worst’ modes of the emergency situation category. For the residual radiation lifetime estimation within the given study, two approaches to determine stress intensity factors, K_I have been used simultaneously. The first approach included a direct numeric modelling of postulated cracks in the cut-out zone with the use of the 3D finite element method. The second approach included K₁ calculation using 3D weight functions calculated with the use of the boundary element method. For K_I, calculation flaws have been postulated as surface longitudinal semielliptical flaws located in the deepest point of a cut-out. The results of K_I calculations obtained using different methods were practically the same. The allowable critical brittleness temperature was determined as 175°C that permitted the extension of the radiation lifetime by up to 6 years after annealing.     

Transport Calculations of MOX and UO2 Pin Cells by the Method of Characteristics

A new transport theory code for two-dimensional calculations of both square and hexagonal fuel lattices by the method of characteristics has been developed. The ray tracing procedure is based on the macroband method, which permits more accurate spatial integration in comparison to the equidistant method of tracing. The neutron source within each region is approximated by a linear function and linearly anisotropic scattering can be optionally accounted for. Efficient new techniques for both azimuthal and polar integration are presented. The spatial discretization problem in case of P ₁-scattering has been studied. Detailed analyses show that the P ¹-scattering in case of regular infinite array of fuel cells is significant, especially for MOX fuel, while the transport correction is inadequate in case of real geometry multi-group calculations. Finally, the complicated nature of the angular flux in MOX and UO₂ fuel cells is demonstrated.     

Recent progress in the application of the finite element method to the neutron transport equation

The finite element technique is applied to the even-parity form of the neutron transport equation by means of a variational functional. A computer code, FELICIT, has been developed which solves criticality and source problems in R–Z and X–Y geometries assuming isotropic scattering in the lab-system. The code is compared with a number of exact analytic solutions in the P₁ approximation and in exact transport theory thus enabling an estimate of its accuracy to be made. Many tables and graphs are given which illustrate the usefulness of the method for a variety of source and criticality problems in complicated geometries.     

Fast Neutron Spectra from Natural Uranium Target Bombarded with High Energy Electrons

The numerical solution of the transport equation has the errors caused by the approximations used in the computational method. In the past estimations of these errors have been performed experimentally. In the present study, formulas to estimate the errors have been derived on the basis of the perturbation theory. This method enables us to deterministically estimate the numerical errors due to the iteration, spatial discretization and Legendre polynomial expansion of scattering transfer cross sections.

Using the error estimation method developed in the present study, two examples of error analyses were carried out to confirm its validity and applicability to error estimation for a practical purpose. The errors of the calculated tritium breeding ratio for ⁷Li in a infinite slab geometry were estimated, and they agreed well with the values predicted from direct calculation. As the second example, error analysis was carried out for one-dimensional nuclear calculations on two types of commercial fusion reactor blankets. In this analysis the tritium breeding ratio and the fast neutron leakage flux from the inboard shield were investigated, and the errors from different causes were quantitatively compared.