Abstract
A neutronic module for the solution of two-dimensional steady-state multigroup diffusion problems in nuclear reactor cores is developed. The module can produce both direct fluxes as well as adjoints, that is, neutron importances. Different numerical schemes are employed. A standard finite-difference approach is firstly implemented, mainly to serve as a reference for less computationally challenging schemes, such as nodal methods and boundary element methods, which are considered in the second part of the work. The validation of the methods proposed is carried out by comparisons of results for reference structures. In particular a critical problem for a homogeneous reactor for which an analytical solution exists is considered as a benchmark. The computational module is then applied to a fast spectrum system, having physical characteristics similar to the proposed lead-cooled ELSY project. The results show the effectiveness of the numerical techniques presented. The flexibility and the possibility to obtain neutron importances allow the use of the module for parametric studies, design assessments, and integral parameter evaluations as well as for future sensitivity and perturbation analyses and as a shape solver for time-dependent procedures.
1. Introduction and Objective of the Work
The European ELSY Project focuses on the assessment of the fast spectrum lead-cooled reactor technology [1]. Various sophisticated numerical tools are today available for the detailed neutronic design of such a system. However, a simple multidimensional computational tool may in any case be rather useful for its flexibility and the possibility to retrieve information on the physical features of the system and to carry out extensive parametric studies. Furthermore, it can be appropriate for an easy and efficient coupling with a thermal-hydraulic module to allow investigations on nonlinear effects. For the same reasons, while it may be very complicate to couple high-performance computational codes to a kinetic module within a quasistatic framework, the tool developed by the authors in this work can be easily used for the evaluation of neutron shapes for system dynamic computations.
The neutronic module is based on a coarse-mesh approach for the calculation of the neutron distribution of the full core. While several nodal methods are available for neutronic calculations, it is interesting to compare and benchmark their performance and to assess their convergence trends. The present work aims at the evaluation of the performance of different computational schemes, starting from a standard finite volume approach in two-dimensional multigroup diffusion on a fine-mesh structure and then proceeding to a nodal formulation of the same problem on a coarse-mesh discretization [2], with the inclusion of a computational option based on the boundary element method (BEM) [3, 4] on the same coarse mesh. The interest of such methods is due to the possibility to develop a computational algorithm based on a response matrix formulation, thus allowing an efficient treatment of a full core geometry including the possibility to easily improve the physical model described (i.e., passing from diffusion to transport), retaining the same structure of the numerical code. The comparison of the performances of the two algorithms (nodal and BEM) for the evaluation of the flux distribution in a full core, both in terms of accuracy of prediction and efficiency of calculation, constitutes an important step in the assessment of coarse-mesh methods for the application to fast reactors. The introduction of the BEM approach in the code developed can also be of interest for further developments that may lead to the possibility of solving the problem with a higher-order neutronic model with respect to diffusion, such as a second-order form of the transport equation [5, 6].
As pointed out, the availability of a steady-state flux solver is rather important also for the development of a time-dependent module which is needed for safety evaluations and assessments, within a quasi-static treatment [7]. For that purpose an adjoint solution is also needed. The adjoint can also be used for sensitivity analyses which are important for the preliminary design of a new nuclear system.
In a first stage a fine-mesh finite-difference module in 2D cartesian geometry is developed, as a reference validation tool for coarse mesh and boundary element algorithms. On its turn, such module is validated against analytical results for homogeneous configurations.
In a second stage the implementation of a nodal method is carried out. A polynomial representation of the neutron flux within a node is used up to second and third order. Eigenvalue and flux results are reported for different number of meshes. The results show the good performance of the coarse-mesh technique associated to a reduced computational effort.
The introduction of a boundary element formulation is then performed and results are presented and compared with the previous techniques. The flexibility of the scheme could be of great importance for future applications to different configurations of lead-cooled fast reactors.
2. The Physical Model
For the neutronic design of a nuclear reactor core the solution of the system of multigroup diffusion equations is usually sufficient. The physical characteristics of a lead-cooled reactor call for a multidimensional treatment. In this paper the diffusion model is adopted and the equations are solved for a two-dimensional configuration. The neutron steady-state multigroup diffusion equations without and with an external independent source, respectively, are the following:
where is the diffusion coefficient of group ; is the removal cross-section of group ; is the scattering cross-section from group to group ; is the mean number of neutrons per fission times the fission cross-section of group is the neutron spectrum of group ; is the neutron scalar flux for group ; is the neutron source for group ; is the effective multiplication constant.
For the homogeneous equation (1) the standard effective multiplication constant is introduced as an eigenvalue to guarantee solvability.
3. Finite-Difference and Nodal Approach
The multigroup diffusion equations without independent source in two-dimensional Cartesian coordinate systems are considered. The solution of the source-driven problem is derived as one iteration step in the power method used to obtain the eigenvalue. The equations can be written explicitly in the following form:
A numerical solution can be readily obtained by the standard finite-difference method, once the domain with -dimension and -dimension is subdivided into meshes that may be chosen as and . The finite-difference solution is used here as a reference for the assessment of the coarse-mesh schemes. The multiplication eigenvalue is determined by the classical power method [8].
Nodal methods have proved to be very effective for the accurate solution of reactor core physics problems in nuclear engineering applications. The general principle of nodal diffusion theory methods require to subdivide the reactor core into a relatively small number of regions, or nodes, which present dimensions larger than the physical characteristic length associated to the diffusion process. The detailed flux distribution within each node is approximated by the superposition of a suitable set of spatial functions. The global flux distribution is then determined by the application of a weighted residual technique associated to proper conservation principles and coupling conditions for the nodes. In the following the derivation of the nodal scheme is presented for a general 3D system. The node is defined by its sides , , and . For each spatial coordinate a one-dimensional equation is obtained by integration over the transverse coordinates. For instance, by integration over and , the following equation is obtained for node , centered at point (, and group :
where Following a well-assessed nodal procedure [9, 10] the unknown transverse-averaged flux is expressed as a superposition of suitable polynomial functions as
If the polynomials herewith appearing are chosen according to the following formulae:
the coefficient of the zeroth-order terms turns out to play the role of node-averaged flux. Furthermore, since , the values of the fluxes at the node boundary can be related to the expansion coefficients by
Outgoing surface currents can also be related to the expansion coefficients and to the incoming currents:
The weighted residual method is now applied for the determination of a relationship involving the unknown coefficients . This can be achieved by spatial integration after the application of the weighting function for any group, with being any positive integer. The following expression is readily obtained:
It is worth taking the limit for , yielding
The physical meaning of this process has been discussed in the previous reference works [9, 10], where the present nodal approach is described in full detail.
Combining (11) with (8)-(9), a relationship connecting incoming and outgoing partial currents with the average flux is obtained
Alternatively the above formulae can be given the following form:
with appropriate definition of the coefficients. The formula above clearly shows the connection between boundary partial currents and node average flux. The algebraic problem is then solved by standard techniques and outer iterations are carried out for the determination of the multiplication eigenvalue, as discussed before.
4. Boundary Elements Approach
The Boundary Elements Method (BEM) [11] constitutes an innovative approach for the effective solution of problems in applied sciences. In the previous works the method has been exploited for applications to the solution of multigroup diffusion problems [3, 4]. Also applications to transport problems in space second-order forms have been recently proposed. Results have proven the effectiveness of the technique and its excellent properties [5, 6].
In the following, the basics of a formulation of the boundary element method for the diffusion problem are outlined. The neutron diffusion equation in a homogeneous region with a closed smooth boundary surface is written as
The fission production and the group-to-group transfer in the multigroup model may be included in the source term. The fundamental solution, that is, the Green function, is found by solving the problem:
Using the properties of the Green function and integrating over the volume equation (15), the following integral formulation is obtained:
where
and is equal to , , or , depending on whether is outside , on the boundary of or inside , respectively. Equation (17) allows to obtain and at any point inside , provided that the flux and the normal current at the points of the boundary surface are known. If we take , (17) becomes an integral equation
in terms of the boundary values of and . Introducing the partial currents
into (19), one has
4.1. Multigroup Boundary Integral Equations
The multigroup equations can be cast into matrix form as
where and are -dimensional vectors. In the case of a subcritical system driven by an external source the following definition holds:
where the fission operator takes the form
On the other hand, for the homogeneous problem,
The eigenequation
is assumed to have real distinct eigenvalues . is defined as the matrix which has the normalized eigenvectors as its columns. Thus, the following equality is verified:
where .
Now setting: after left multiplication by , (22) reads By redefining , the following set of equations is obtained for each component of vector : The corresponding Green fundamental solutions of the following equations: are easily found out as together with their derivatives
where and are the modified Bessel functions of orders and , respectively.
Boundary integral equations can be generated with the same procedure leading to (17), obtaining
By left multiplying the set of (34) written in matrix form by , one gets back to relations involving , namely
where
Let be taken on the boundary, . Then, finally, the multigroup boundary integrated equations can be explicitly written down:
where
4.2. Reduction of the Domain Integrals of the Source Term
The domain integrals involving the independent and flux dependent fission sources in (37) can be reduced into a boundary integrated form as
The source distributions are now assumed to be either spatially constant or first-order polynomials. In the latter case for 2D geometry the following expression is assumed:
where, obviously, is the average source on the domain and . The following relationship holds:
As in the presence of fission emissions is evaluated by the flux of the previous generation on the boundary, the coefficients of expansion (40) can be evaluated by using a least-square approach. The problem can be nicely written in matrix form
having set the following definitions:
where are points on the boundary at which the source takes the value . By left multiplication of (42) by , the normal equations are obtained for the problem or, explicitly
The above system of equations yields the value of the two-dimensional vector .
4.3. Discretization of the Boundary Integral Equations
In this section the numerical solution of (37) is outlined. The boundary is partitioned into elements of length . After denoting the two extremes and the middle point by , , and , respectively, a local dimensionless coordinate is introduced to denote each point along as
Thus, it is assumed that
Then, assuming a linear behaviour on each boundary mesh of the integrand functions, the boundary integrals are approximated by finite sums as
By indicating
one has
The integrals herewith appearing are obtained by application of a numerical integration formula, namely
where and are the quadrature weights and abscissae, respectively. In the present paper a standard Gauss-Legendre formula is used.
At last the problem is cast into the following algebraic response-matrix form:
where
The problem for the whole system where nodes are coupled is then solved by numerical iteration schemes.
5. Verification of Numerical Consistency and Convergence of Methods
It is worth going through a verification of the performance of various numerical approaches by comparison of results with highly accurate benchmark results that may be obtained directly by an analytical formula. Therefore a homogeneous rectangular system is now considered, whose dimensions are 160 cm and 140 cm. The material data for a three-group model are reported in Table 1. A critical calculation is carried out and errors with respect to the analytical value of the effective multiplication constant are reported in Table 2. The various models are identified as follows: FD indicates the finite-difference scheme, N is nodal, BEM denotes a distribution of order for the source (0, constant, or 1, linear, in the present work), the number of points used for the discretization of each edge, and the order of the quadrature formula for integration. The number reported in the node column indicates the number of nodes assumed for each coordinate; hence the number of nodes covering the 2D domain increases as the square of the indicated number. The performance of the nodal method is quite remarkable, while good results are obtained by BEM only with a first-order approximation of the source distribution. Also, it is important to call the attention on the important role associated to the number of points used for the discretization of the boundary.
Furthermore, it is interesting to observe the convergence pattern for the methods studied, see Figure 1. While FD shows a monotonic behaviour, Nodal and some BEM methods are characterized by nonmonotonic trends and reach a saturation. Such properties for the convergence of nodal methods are well-known to exist [12] and would deserve deeper investigation. They are connected to the assumption of the spatial distribution within the node, which, in general, may not be consistent with the physical model as described by the differential equations to be solved.
); BEM 1 10 24 (*); BEM 1 10 240 (•); BEM 1 20 240 ().Some final considerations concerning the computational effort are in order. Considering Figure 2, it is clear that to obtain relative errors below 10 pcm nodal and BEM are by far more effective.
A nonhomogeneous system is also considered, to test the performance of the schemes in the presence of a spatial interface. The symmetrical core is 180 cm wide and characterized by the data reported in Table 1. The data for the 90 cm thick reflector are specified in Table 3. The system is reflected only in the -direction. The number of nodes herewith indicated refers to the discretization in the -direction only, since along only two nodes are used with zero-current on the boundary condition, to simulate an infinite medium. Results concerning values are reported in Table 4. A very fine finite-difference solution (320 nodes), to be assumed as reference, yields the value . The enhanced performance of nodal and boundary element methods is clearly visible.
6. Results for a Fast-Spectrum Configuration
A realistic fast-spectrum system is considered, with physical characteristics of the same type as the ones characterizing the ELSY project. The core is constituted by three regions, with increasing multiplicativities passing from the inner to the outermost regions. Lead is assumed for the coolant and the reflector. Cross-sections are provided by spectrum and homogenization calculations carried out by the ENEA group participating to the ELSY project. The sketch of the system and the nuclear data for the core regions and for the reflector are reported in Figure 3 and Table 5. Although only three energy groups are used, the test is useful to verify the performance of the computational schemes in a highly heterogeneous system where spectral effects are important.
As already anticipated, adjoints may well serve in many interesting and useful applications in nuclear core analysis. Once the physical problem is cast into an algebraic form, the adjoint problem is easily constructed by taking the mathematical adjoint of the matrix operators appearing in the model. It is then possible to compute neutron importances that may be used for sensitivity and perturbation analyses and for the construction of integral parameters for time-dependent calculations. The results of some such calculations are now presented in Figures 4, 5, 6, 7, 8, and 9. For brevity, only results for the nodal method are presented, having previously proved its consistency with BEM.
7. Conclusions and Future Developments
A neutronic module for multidimensional multigroup diffusion calculations is developed. Different discretization algorithms are implemented and tested. Some benchmark calculations allow to assess the consistency and applicability of various techniques. Nodal methods are shown to perform adequately for the simulation of fast spectrum systems, such as lead-cooled reactors. Calculations are carried out for a system characterized by properties similar to the system adopted as the basis of the ELSY project.
As a future development, it is foreseen to couple the neutronic module developed in the frame of this work with a thermohydraulic code. This would allow to have a flexible and efficient computational tool for the assessment of the ELSY project and for parametric calculations aimed at the neutronic-thermal-hydraulic optimization. Furthermore, the treatment of full three-dimensional Cartesian geometries is envisaged.
A further future development will consider the extension of the computational tool to allow time-dependent evaluations within a quasi-static approach. To that end, the module can perform already adjoint calculations. The extension is to be done following the same philosophy as the one already employed for other recent applications [13].
Acknowledgments
The work is performed with the support of the ELSY project (VI European Framework Program). The authors are grateful to ENEA-Bologna for providing the cross-section data utilized for the simulations and to A. Cencio for the help in carrying out part of the numerical work. One of the authors (S. Han) acknowledges the support received by the Zhong Guo program and by the Asia Link Project “An EU-China Campus for Energy and Environment” for making available the opportunity of a PhD program on Energy Engineering at Politecnico di Torino.