Abstract
Numerical methods are usually required to solve the neutron diffusion equation applied to nuclear reactors due to its heterogeneous nature. The most popular numerical techniques are the Finite Difference Method (FDM), the Coarse Mesh Finite Difference Method (CFMD), the Nodal Expansion Method (NEM), and the Nodal Collocation Method (NCM), used virtually in all neutronic diffusion codes, which give accurate results in structured meshes. However, the application of these methods in unstructured meshes to deal with complex geometries is not straightforward and it may cause problems of stability and convergence of the solution. By contrast, the Finite Element Method (FEM) and the Finite Volume Method (FVM) are easily applied to unstructured meshes. On the one hand, the FEM can be accurate for smoothly varying functions. On the other hand, the FVM is typically used in the transport equations due to the conservation of the transported quantity within the volume. In this paper, the FVM algorithm implemented in the ARB Partial Differential Equations solver has been used to discretize the neutron diffusion equation to obtain the matrices of the generalized eigenvalue problem, which has been solved by means of the SLEPc library.
1. Introduction
The neutron diffusion equation is used to calculate the neutron flux distribution, which is one of the most important variables in a Nuclear Power Reactor (NPR). This equation is a simplification of the neutron transport equation using Fick’s Law, as discussed by Stacey [1]. Nevertheless, the use of the neutron diffusion equation is justified by the lower computational time and relatively low heterogeneity of commercial NPR.
In order to calculate the spatial distribution of the neutron flux, the steady state of the neutron diffusion equation is considered by transforming the neutron diffusion equation into a generalized eigenvalue problem, explained in Section 2.1.
The greatest eigenvalue is the most important one and it has a special interest for nuclear reactors safety. As a result, most methods used to calculate only this eigenvalue and utilizing iterative methods to avoid solving the generalized eigenvalue problem. Nevertheless, the calculation of several eigenvalues and eigenvectors is important for different applications as the modal analysis of nuclear reactors and BWR instabilities analysis, as discussed elsewhere [2, 3].
However, the resolution of this generalized eigenvalue problem could be a difficult task due to the large and sparse nature of the matrices. In this paper, the SLEPc library solves this problem. Actually, the emphasis of this library is on methods appropriate for problems in which the associated matrices are sparse, such as those arising after the discretization of partial differential equations, as discussed by Hernandez et al. [4].
On the other hand, numerical methods are usually required to solve the neutron diffusion equation applied to nuclear reactors due to its heterogeneous nature, discretizing the partial differential terms. The most popular numerical techniques are the Finite Difference Method (FDM), the Coarse Mesh Finite Difference Method (CFMD), the Nodal Expansion Method (NEM), and the Nodal Collocation Method (NCM), used virtually in all the neutronic diffusion codes, which give accurate results in structured meshes. However, the application of these methods in unstructured meshes to deal with complex geometries is not straightforward and it may cause problems of stability and convergence of the solution, as discussed by Hoffmann and Chiang [5]. In fact, the use of unstructured meshes is justified by the thermal hydraulic-neutronic coupled calculation, which sometimes uses the unstructured meshes. By contrast, the Finite Element Method (FEM) and the Finite Volume Method (FVM) are easily applied to unstructured meshes. On the one hand, the FEM can be accurate for smoothly varying functions. On the other hand, the FVM is typically used in the transport equations due to the conservation of the transported quantity within the volume.
In this paper, the FVM algorithm implemented in the ARB Partial Differential Equations solver has been used to discretize the neutron diffusion equation to obtain the matrices of the generalized eigenvalue problem. The strength of ARB is that a fully implicit numerical formulation is generated and formulated easily, for an arbitrary set of equations, which are input by the user using pseudomathematical expressions, as discussed by Harvie [6].
The outline of the paper is as follows. Section 2 presents the discretization of the equations and the methodology used. Section 3 describes the reactors used to validate the method and their results. Section 4 contains few comments and conclusions about the results.
2. Materials and Methods
2.1. Discretization of the Diffusion Equation Using the Finite Volume Method
Although several approaches in terms of energy could be applied to the neutron diffusion equation, the most widely used for Light Water Reactors (LWR) is the 2-energy group neutron diffusion approximation, as discussed by Stacey [1]:
In order to calculate the spatial distribution of the neutron flux, the steady state of the previous differential equations system is considered by putting the temporal derivation of (1) equal to 0. However, if the geometry of the problem is fixed, the steady state will be accomplished only for certain set of diffusion coefficients and cross sections, and vice versa. These nuclear parameters depend on the materials and are the following coefficients of (1): , , , , , , and .
Therefore, the problem is transformed into an eigenvalue problem to assure the steady state accomplishment, where the eigenvectors are the spatial distribution of the neutron flux and the inverse of the first eigenvalue represents a measure of the steady state condition, numbering in descending order:
Moreover, the geometry should be discretized in elements where the neutron diffusion equation will be applied, due to the reactor heterogeneity and the cross-section spatial dependence. Consequently, a set of equations will be obtained for each element. Then, these equations are integrated in each element volume and the Divergence Theorem is applied, so that the divergence term could be avoided:
Finally, considering surface and volume averaged values and dividing by the volume of the element , the following equations will be obtained: where represents each face surrounding the element , is the volume of the element , is the area of the face , could be −1 or 1 depending on the direction of the face with respect to the direction of the neutron flux gradient at this face, and are the first and second energy group diffusion coefficients volume averaged values for the element , and are the first and second energy group absorption macroscopic cross-section volume averaged values for the element , and are the first and second energy group nufission macroscopic cross-section volume averaged values for the element , is the scattering macroscopic cross-section volume averaged values for the element , and are the first and second energy group volume averaged values of the neutron flux for the element , and and are the first and second energy group surface averaged values of the neutron flux gradient for the face .
With regard to and , although they are not known, moving Least-Squares reproducing kernel methods can express them as a weighted sum of and , respectively, by considering the elements surrounding face , as discussed elsewhere [6, 7]: where is each element surrounding the face and is the weighting factor of the element with respect to the face , which is called kernel by Harvie [6].
Therefore, the final equations for each element will be
The boundary conditions typically used in nuclear reactors are face boundary conditions, and consequently these boundary conditions should be expressed in terms of the volume averaged values of the neutron flux, which are the unknown values. The boundary conditions most widely used are zero flux or reflective flux expressed mathematically by (7) and (8), respectively,
Therefore, the same procedure as the neutron flux gradient could be used to obtain these face averaged values in terms of the volume averaged values: where represents each element surrounding the face that is part of the boundary condition BC.
As a result, the previous equations will compose a generalized eigenvalue problem considering (6) for each element and (9) and/or (10) for each boundary:
entries correspond to the first equation left side of (6), and entries correspond to the first equation right side of (6), while and entries correspond to the second equation of (6), and and entries correspond to (9) and/or (10).
On the other hand, and are vectors containing the flux volume averaged values in each element for the first and second energy group, respectively. Nonetheless, the addition of boundary conditions equations implies an equations excess, and therefore more unknown values have to be taken into account. In order to solve this problem, some virtual fluxes are considered, in particular, the same number as the boundary conditions equations. These virtual fluxes have no physical meaning and (6) is not applied to them, but they are needed to evaluate the boundary conditions equations and obtain the same number of equations and unknowns. In addition, these virtual fluxes have the same centroid as the faces where they are defined, and their volume values are the area values of these faces.
In Figure 1, an example of the virtual fluxes is shown. In this example the eigenvectors are
In conclusion, the following formula represents the generalized eigenvalue problem:
2.2. Neutron Current Condition
In the neutron diffusion theory, the partial neutron current calculated at the face separating 2 cells must be the same in both cells. Therefore, if Figure 2 is considered, the neutron current condition at the face separating these cells will be
According to Fick’s Law, where is the energy group diffusion coefficient for the left cell, is the energy group diffusion coefficient for the right cell, is the energy group partial neutron flux gradient for the left cell, and is the energy group partial neutron flux gradient for the right cell.
If both cells are composed of same material, and will be the same and the partial neutron flux gradient will be the same; therefore, the neutron current condition will be accomplished. Nevertheless, if the materials are different, the accomplishment of this condition is not assured.
On the other hand, the method proposed in this paper does not include the calculation of the partial gradients, only the global gradient calculation, and it depends on the face but not on the cell. These gradients are exposed in Figure 3. As a result, a global diffusion coefficient is needed in order to accomplish the neutron current condition by using the global gradient:
In this paper, 4 approaches of this global diffusion coefficient have been considered.(i) (CELL+): take the diffusion coefficient of the left cell as the global diffusion coefficient: (ii) (CELL−): take the diffusion coefficient of the right cell as the global diffusion coefficient: (iii) (HOM): take a weighted sum of the diffusion coefficients of both cells as the global diffusion coefficient. In this case, the weighting factors are based on the kernels mentioned in Section 2.1: (iv) (LIN): by considering the neutron flux is a line function of the distance between the centroids of the cell and face. In addition, the distances have been replaced by the kernels mentioned in Section 2.1 due to the fact that kernels depend on distances:
is the kernel of the left cell to calculate the contribution to the neutron flux gradient at the interface. It may be negative, so the absolute value is considered to avoid a null dividing.
is the kernel of the right cell to calculate the contribution to the neutron flux gradient at the interface. It may be negative, so the absolute value is considered to avoid a null dividing.
Therefore, (6) is transformed into (21), and this last equation should be used instead of (6) to obtain the matrices entries of (13): where and are the first and second energy group global diffusion coefficients for the face .
2.3. Calculation Methodology
The following steps of the calculation have been done by the next codes:(i)geometry meshing by means of Gmsh developed by Geuzaine and Remacle [8];(ii)discretization of the neutron diffusion equations with the Finite Volume Method by means of Arb developed by Harvie [6];(iii)solution of the generalized eigenvalue problem by means of SLEPc developed by Hernández et al. [4, 9].
A code has been developed to do these steps automatically.
3. Results and Discussions
In this section, 4 different cases are exposed: homogeneous and heterogeneous 2D and 3D reactors, with the aim of showing the capabilities of the method for 2D and 3D reactors. In fact, the homogeneous reactors check the discretization of the equations without taking into account the use of the global diffusion coefficient and the heterogeneous reactors check the different approaches of the global diffusion coefficient developed in this study. Regarding the calculation, 5 eigenvalues have been calculated in each case. With respect to the results, the following magnitudes will be used so as to evaluate them:
In addition, power is normalized to attain that mean power equals the unity, calculated with the following formula: where is volume or area of the element , depending on 3D or 2D geometry, respectively, and considering only the elements with not null power.
3.1. 2D Homogeneous Reactor
The reactor considered is composed only of one material, whose cross sections are shown in Table 1. In addition, it has a rectangular shape and its dimensions are . Furthermore, zero flux boundary conditions have been considered.
On the other hand, this study includes the results for structured and unstructured meshes, each one with several sizes. In particular, 3 different mesh lengths have been used: 10 cm, 5 cm, and 2 cm. Table 2 contains the number of cells of structured and unstructured meshes. Furthermore, Figures 4 and 5 show the structured and unstructured mesh, respectively, for the 5 cm mesh length.
Moreover, this problem has analytical solution and it will be the reference solution. The analytical eigenvalues are 0.99999996179, 0.94332347153, 0.85966257142, 0.85451196366, and 0.81030009136. In this case, the reactor is in steady state since the greatest eigenvalue is virtually 1.
The results corresponding to computational time and eigenvalue errors are exposed in Table 2. On the other hand, power errors are evaluated in 9 rectangular nodes of the same dimensions , which are shown in Table 3, and they are presented in Tables 4 and 5, but only those corresponding to the first and second eigenvalues. In addition, the power corresponding to the first and second eigenvalues and the 2 cm structured and unstructured meshes is shown in Figures 6, 7, 8, and 9.
For the second eigenvalue, there is not error in nodes 2, 5, and 8, since the power in these nodes is 0, so the error in these nodes has been represented with “—” in Table 5. Actually, it can be seen in Figures 8 and 9 that power in nodes 2, 5, and 8 are virtually 0. Moreover, there are high errors in structured mesh of 10 cm in length for the third and fourth eigenvalue, and consequently this mesh is not acceptable. It can be noted that unstructured mesh errors are higher than structured mesh ones. Furthermore, the results show an error decrease as the mesh is finer. In any case, the maximum power error is 0.370735%, corresponding to the second eigenvalue, 10 cm unstructured mesh, and node 4.
3.2. Biblis Reactor
Biblis is a 2D heterogeneous reactor composed of 8 materials. Its geometry and cross sections are described in Figure 10 and Table 6, respectively. A quarter of the reactor has been simulated, and therefore reflective flux has been assumed at west and north boundaries, and zero flux at east and south boundaries.
Moreover, the same meshes as in Section 3.1 have been used.
On the other hand, the reference solution considered was the solution obtained by Müller and Weiss [10]. This solution was obtained by means of PANIC analytic nodal code using a mesh, but only the results for the first eigenvalue are available. Therefore, only the results for the first eigenvalue will be exposed in this section, although the computational time corresponds to the calculation of five eigenvalues. Number of elements, computational times and eigenvalue errors are shown in Tables 7, 8, and 9, respectively. The reference solution of the first eigenvalue is 1.025110.
In this case, the 4 approaches of the global diffusion coefficient have been used since this reactor is heterogeneous.
With respect to power errors, the results are evaluated at the nodes of Figure 10, but without taking into account the nodes of material 3, since the power in these nodes is 0. In particular, the nomenclature exposed in Table 10 will be used.
Only power errors corresponding to HOM are exposed in Table 11, because they are the most accurate results. The power corresponding to the first eigenvalue and the 2 cm structured and unstructured meshes is shown in Figures 11 and 12.
The finer the mesh, the lower the errors. In addition, unstructured mesh errors are lower than structured mesh ones. Regarding the different approaches of the global diffusion coefficient, the lowest power errors correspond to HOM and LIN. However, the lowest eigenvalue error corresponds to CELL− in coarse meshes, yet this error is virtually the same for all the approaches in the finest mesh.
3.3. 3D Homogeneous Reactor
The reactor considered is composed of the same material as the 2D homogeneous reactor, it is a parallelepiped of the next dimensions: . In addition, the boundary conditions applied are zero flux. Regarding the meshes, the same mesh lengths as in the previous sections have been used again, but only the structured mesh.
Moreover, the analytical solution exists, and therefore it will be the reference solution. On the one hand, eigenvalue errors and computational time are shown in Table 12. On the other hand, power errors are evaluated in 54 parallelepiped nodes of the same dimensions , which are exposed in Table 13, and they are shown in Tables 14 and 15, but only those corresponding to the first and second eigenvalues. In addition, the power corresponding to the first and second eigenvalues and the 2 cm structured mesh is shown in Figures 13 and 14.
The analytical eigenvalues are 0.99391916952, 0.97602952377, 0.94734259138, 0.93778667636, and 0.92148598185. In this case, the reactor is in steady state since the greatest eigenvalue is virtually 1.
In conclusion, the finer the mesh, the more accurate the results. However, the maximum error of the coarse mesh is about 1%, which is acceptable. On the other hand, the computational time of the finest mesh is not practical.
3.4. Langenbuch Reactor
Langenbuch is a 3D heterogeneous reactor composed of 4 materials. Its cross sections are exposed in Table 16. Its geometry is described in Figures 15 and 16. A quarter of the reactor has been simulated, and therefore reflective flux has been assumed at west and south boundaries, and zero flux at east, north, top, and bottom boundaries.
In this case, only structured meshes have been used, in particular the 10 cm and 5 cm mesh length. Table 17 shows the number of cells for each mesh.
Moreover, the reference solution has been obtained with VALKIN code developed by Miró et al. [3], which was called MODKIN in the cited reference. VALKIN is a nodal modal code that is able to calculate several eigenvalues and their respective eigenvectors, which are the neutronic fluxes. In this case, the reference solution was obtained for 5 eigenvalues at nodes of Figures 15 and 16. Furthermore, the 4 approaches of the global diffusion coefficient have been used owing to the heterogeneous nature of this reactor. The reference solution of the first and second eigenvalues ( and ) are 0.994881227 and 0.948210698.
Regarding the results, Table 18 exposes the computational time for each mesh. Moreover, eigenvalue errors are shown in Tables 19 and 20, but only those corresponding to the first and second eigenvalues. On the other hand, power errors have been calculated in the nodes of Figures 15 and 16, but only the nodes numbered in Table 21 are exposed in this paper due to the results extent. These power errors are shown in Tables 22 and 23, but only those corresponding to the first and second eigenvalues and HOM global diffusion coefficient approach. In addition, the power corresponding to the first and second eigenvalues and the 5 cm structured mesh is shown in Figures 17 and 18.
With respect to the meshes used, the 10 cm mesh length implies the highest errors and the 5 cm mesh length the lowest ones. Furthermore, the highest errors are located near the reflector. Regarding the global diffusion coefficient approaches, HOM and LIN give better results than CELL+ and CELL−. It is important to remark the computational time dependence on the mesh that implies 2 cm mesh is not practical owing to the fact that it has an order of magnitude of hours.
4. Conclusions
A method has been developed to solve the steady state of the 2-energy group neutron diffusion equation for LWR in any reactor configuration, using the Finite Volume Method and calculating several eigenvalues.
This method supplies accurate results for 2D reactors and low computational times, about seconds. However, 3D reactors computational times are higher that could be about hours for fine meshes. Moreover, 3D reactor results are less accurate than 2D reactor ones.
If the global diffusion coefficient had not been used, the power errors corresponding to the first eigenvalue would have been about 15–20%, although this study does not show these results. Consequently, the global diffusion coefficient has to be used to obtain acceptable results.
With reference to future work, the method will include the parallelization of both geometry processing and eigenvalue calculation to reduce the computational time. Furthermore, more global diffusion coefficients approaches will be developed, and another alternatives to evaluate the face averaged gradient flux are being considered as the implementation of high-order schemes. Regarding the nuclear applications, the transitory state calculation will be developed in order to evaluate any reactor condition. Finally, the advanced thermal-hydraulic coupling will be the final step to take into account the thermal hydraulic influence in the neutronic calculation.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work has been partially supported by the Spanish Ministerio de Ciencia e Innovación under Projects ENE2011-22823 and ENE2012-34585, the Generalitat Valenciana under Projects PROMETEO/2010/039 and ACOMP/2013/237, and the Universitat Politècnica de València under Project UPPTE/2012/118.