# Resolution of the Generalized Eigenvalue Problem in the Neutron Diffusion Equation Discretized by the Finite Volume Method

Authors: Álvaro Bernal; Rafael Miró; Damián Ginestar; Gumersindo Verdú
Journal: Abstract and Applied Analysis
 (2014)

---

## 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.


---

## Body

## 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 Section2.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. Section2 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]:
(1)
1
v
1
d
ϕ
1
(
r
→
,
t
)
d
t
=
-
∇
→
(
-
D
1
∇
ϕ
1
(
r
→
,
t
)
)
-
Σ
a
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
-
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
(
υ
Σ
f
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
υ
Σ
f
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
)
,
1
v
2
d
ϕ
2
(
r
→
,
t
)
d
t
=
-
∇
→
(
-
D
2
∇
ϕ
2
(
r
→
,
t
)
)
-
Σ
a
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
+
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
.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): D
1, D
2, Σ
a
,
1, Σ
a
,
2, υ
Σ
f
,
1, υ
Σ
f
,
2, and Σ
s
,
1
→
2.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(
k
1
) represents a measure of the steady state condition, numbering k in descending order:
(2)
∇
→
(
-
D
1
∇
ϕ
1
(
r
→
,
t
)
)
+
Σ
a
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
-
1
k
(
υ
Σ
f
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
υ
Σ
f
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
)
=
0
,
∇
→
(
-
D
2
∇
ϕ
2
(
r
→
,
t
)
)
+
Σ
a
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
-
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
=
0
.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(
V
i
) and the Divergence Theorem is applied, so that the divergence term could be avoided:
(3)
∫
(
-
D
1
∇
ϕ
1
(
r
→
,
t
)
)
d
S
V
i
+
∫
Σ
a
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
d
V
i
+
∫
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
d
V
i
=
1
k
∫
(
υ
Σ
f
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
υ
Σ
f
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
)
d
V
i
,
∫
(
-
D
2
∇
ϕ
2
(
r
→
,
t
)
)
d
S
V
i
+
∫
Σ
a
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
d
V
i
-
∫
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
d
V
i
=
0
.Finally, considering surface and volume averaged values and dividing by the volume of the element(
V
i
), the following equations will be obtained:
(4)
∑
j
(
S
j
V
i
u
j
(
-
D
1
V
i
∇
ϕ
1
-
j
)
)
+
Σ
a
,
1
V
i
ϕ
1
,
i
+
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
1
k
(
υ
Σ
f
,
1
V
i
ϕ
1
,
i
+
υ
Σ
f
,
2
V
i
ϕ
2
,
i
)
,
∑
j
(
S
j
V
i
u
j
(
-
D
2
V
i
∇
ϕ
2
-
j
)
)
+
Σ
a
,
2
V
i
ϕ
2
,
i
-
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
0
,
where j represents each face surrounding the element i, V
i is the volume of the element i, S
j is the area of the face j, u
j could be −1 or 1 depending on the direction of the face j with respect to the direction of the neutron flux gradient at this face, D
1
V
i and D
2
V
i are the first and second energy group diffusion coefficients volume averaged values for the element i, Σ
a
,
1
V
i and Σ
a
,
2
V
i are the first and second energy group absorption macroscopic cross-section volume averaged values for the element i, υ
Σ
f
,
1
V
i and υ
Σ
f
,
2
V
i are the first and second energy group nufission macroscopic cross-section volume averaged values for the element i, Σ
s
,
1
→
2
V
i is the scattering macroscopic cross-section volume averaged values for the element i, ϕ
1
,
i and ϕ
2
,
i are the first and second energy group volume averaged values of the neutron flux for the element i, and ∇
ϕ
1
-
j and ∇
ϕ
2
-
j are the first and second energy group surface averaged values of the neutron flux gradient for the face j.With regard to∇
ϕ
1
-
j and ∇
ϕ
2
-
j, although they are not known, moving Least-Squares reproducing kernel methods can express them as a weighted sum of ϕ
1
,
n and ϕ
2
,
n, respectively, by considering the elements n surrounding face j, as discussed elsewhere [6, 7]:
(5)
∇
ϕ
g
-
j
=
∑
n
k
n
,
j
grad
ϕ
g
,
n
,
g
=
1,2
,
where n is each element surrounding the face j and k
n
,
j
grad is the weighting factor of the element n with respect to the face j, which is called kernel by Harvie [6].Therefore, the final equations for each element will be(6)
∑
j
(
S
j
V
i
u
j
(
-
D
1
V
i
∑
n
k
n
,
j
grad
ϕ
1
,
n
)
)
+
Σ
a
,
1
V
i
ϕ
1
,
i
+
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
1
k
(
υ
Σ
f
,
1
V
i
ϕ
1
,
i
+
υ
Σ
f
,
2
V
i
ϕ
2
,
i
)
,
∑
j
(
S
j
V
i
u
j
(
-
D
2
V
i
∑
n
k
n
,
j
grad
ϕ
2
,
n
)
)
+
Σ
a
,
2
V
i
ϕ
2
,
i
-
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
0
.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,
(7)
ϕ
g
BC
=
0
,
g
=
1,2
,
(8)
∇
ϕ
g
BC
=
0
,
g
=
1,2
.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:(9)
ϕ
g
BC
=
∑
n
k
n
,
j
ϕ
g
,
n
=
0
,
g
=
1,2
,
(10)
∇
ϕ
g
BC
=
∑
n
k
n
,
j
grad
ϕ
g
,
n
=
0
,
g
=
1,2
,
where n represents each element surrounding the face j 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:
(11)
(
L
11
0
L
21
L
22
L
31
L
32
)
(
ϕ
1
ϕ
2
)
=
1
k
(
M
11
M
12
0
0
0
0
)
(
ϕ
1
ϕ
2
)
.L
11 entries correspond to the first equation left side of (6), M
11 and M
12 entries correspond to the first equation right side of (6), while L
21 and L
22 entries correspond to the second equation of (6), and L
31 and L
32 entries correspond to (9) and/or (10).On the other hand,ϕ
1 and ϕ
2 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 Figure1, an example of the virtual fluxes is shown. In this example the eigenvectors are
(12)
ϕ
g
=
(
ϕ
g
elements
,
ϕ
g
virtual
)
T
,
g
=
1,2
,
ϕ
g
elements
=
(
ϕ
g
1
,
ϕ
g
2
,
ϕ
g
3
,
ϕ
g
4
,
ϕ
g
5
,
ϕ
g
6
,
ϕ
g
7
,
ϕ
g
8
,
ϕ
g
9
)
T
,
g
=
1,2
,
ϕ
g
virtual
=
(
ϕ
g
10
,
ϕ
g
11
,
ϕ
g
12
,
ϕ
g
13
,
ϕ
g
14
,
ϕ
g
15
,
ϕ
g
16
,
ϕ
g
17
,
ϕ
g
18
,
ϕ
g
19
,
ϕ
g
20
,
ϕ
g
21
)
T
,
g
=
1,2
.Figure 1
Virtual fluxes example.In conclusion, the following formula represents the generalized eigenvalue problem:(13)
(
L
11
0
L
21
L
22
L
31
L
32
)
(
[
ϕ
1
elements
ϕ
1
virtual
]
[
ϕ
2
elements
ϕ
2
virtual
]
)
=
1
k
(
M
11
M
12
0
0
0
0
)
(
[
ϕ
1
elements
ϕ
1
virtual
]
[
ϕ
2
elements
ϕ
2
virtual
]
)
.


### 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 Figure2 is considered, the neutron current condition at the face separating these cells will be
(14)
J
g
+
=
J
g
-
,
g
=
1,2
.Figure 2
Adjacent cells.According to Fick’s Law,(15)
D
g
+
∇
ϕ
g
+
=
D
g
-
∇
ϕ
g
-
,
g
=
1,2
,
where D
g
+ is the g energy group diffusion coefficient for the left cell, D
g
- is the g energy group diffusion coefficient for the right cell, ∇
ϕ
g
+ is the g energy group partial neutron flux gradient for the left cell, and ∇
ϕ
g
- is the g energy group partial neutron flux gradient for the right cell.If both cells are composed of same material,D
g
+ and D
g
- 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 Figure3. As a result, a global diffusion coefficient (
D
g
j
) is needed in order to accomplish the neutron current condition by using the global gradient:
(16)
D
g
j
∇
ϕ
g
=
D
g
+
∇
ϕ
g
+
=
D
g
-
∇
ϕ
g
-
.Figure 3
Partial and global gradients.In this paper, 4 approaches of this global diffusion coefficient have been considered.(i)
Case
1 (CELL+): take the diffusion coefficient of the left cell as the global diffusion coefficient:
(17)
D
g
j
=
D
g
+
.
(ii)
Case
2 (CELL−): take the diffusion coefficient of the right cell as the global diffusion coefficient:
(18)
D
g
j
=
D
g
-
.
(iii)
Case
3 (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:
(19)
D
i
j
=
D
i
+
·
|
k
+
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
+
D
i
-
·
|
k
-
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
.
(iv)
Case
4 (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:
(20)
D
i
j
=
(
D
i
+
·
D
i
-
)
×
(
D
i
-
·
|
k
+
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
+
D
i
+
·
|
k
-
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
)
-
1
.k
+
grad 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.k
-
grad 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):
(21)
∑
j
(
S
j
V
i
u
j
(
-
D
1
j
∑
n
k
n
,
j
grad
ϕ
1
,
n
)
)
+
Σ
a
,
1
V
i
ϕ
1
,
i
+
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
1
k
(
υ
Σ
f
,
1
V
i
ϕ
1
,
i
+
υ
Σ
f
,
2
V
i
ϕ
2
,
i
)
,
∑
j
(
S
j
V
i
u
j
(
-
D
2
j
∑
n
k
n
,
j
grad
ϕ
2
,
n
)
)
+
Σ
a
,
2
V
i
ϕ
2
,
i
-
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
0
,
where D
1
j and D
2
j are the first and second energy group global diffusion coefficients for the face j.


### 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.


## 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]:
(1)
1
v
1
d
ϕ
1
(
r
→
,
t
)
d
t
=
-
∇
→
(
-
D
1
∇
ϕ
1
(
r
→
,
t
)
)
-
Σ
a
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
-
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
(
υ
Σ
f
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
υ
Σ
f
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
)
,
1
v
2
d
ϕ
2
(
r
→
,
t
)
d
t
=
-
∇
→
(
-
D
2
∇
ϕ
2
(
r
→
,
t
)
)
-
Σ
a
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
+
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
.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): D
1, D
2, Σ
a
,
1, Σ
a
,
2, υ
Σ
f
,
1, υ
Σ
f
,
2, and Σ
s
,
1
→
2.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(
k
1
) represents a measure of the steady state condition, numbering k in descending order:
(2)
∇
→
(
-
D
1
∇
ϕ
1
(
r
→
,
t
)
)
+
Σ
a
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
-
1
k
(
υ
Σ
f
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
υ
Σ
f
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
)
=
0
,
∇
→
(
-
D
2
∇
ϕ
2
(
r
→
,
t
)
)
+
Σ
a
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
-
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
=
0
.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(
V
i
) and the Divergence Theorem is applied, so that the divergence term could be avoided:
(3)
∫
(
-
D
1
∇
ϕ
1
(
r
→
,
t
)
)
d
S
V
i
+
∫
Σ
a
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
d
V
i
+
∫
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
d
V
i
=
1
k
∫
(
υ
Σ
f
,
1
(
r
→
)
ϕ
1
(
r
→
,
t
)
+
υ
Σ
f
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
)
d
V
i
,
∫
(
-
D
2
∇
ϕ
2
(
r
→
,
t
)
)
d
S
V
i
+
∫
Σ
a
,
2
(
r
→
)
ϕ
2
(
r
→
,
t
)
d
V
i
-
∫
Σ
s
,
1
→
2
(
r
→
)
ϕ
1
(
r
→
,
t
)
d
V
i
=
0
.Finally, considering surface and volume averaged values and dividing by the volume of the element(
V
i
), the following equations will be obtained:
(4)
∑
j
(
S
j
V
i
u
j
(
-
D
1
V
i
∇
ϕ
1
-
j
)
)
+
Σ
a
,
1
V
i
ϕ
1
,
i
+
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
1
k
(
υ
Σ
f
,
1
V
i
ϕ
1
,
i
+
υ
Σ
f
,
2
V
i
ϕ
2
,
i
)
,
∑
j
(
S
j
V
i
u
j
(
-
D
2
V
i
∇
ϕ
2
-
j
)
)
+
Σ
a
,
2
V
i
ϕ
2
,
i
-
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
0
,
where j represents each face surrounding the element i, V
i is the volume of the element i, S
j is the area of the face j, u
j could be −1 or 1 depending on the direction of the face j with respect to the direction of the neutron flux gradient at this face, D
1
V
i and D
2
V
i are the first and second energy group diffusion coefficients volume averaged values for the element i, Σ
a
,
1
V
i and Σ
a
,
2
V
i are the first and second energy group absorption macroscopic cross-section volume averaged values for the element i, υ
Σ
f
,
1
V
i and υ
Σ
f
,
2
V
i are the first and second energy group nufission macroscopic cross-section volume averaged values for the element i, Σ
s
,
1
→
2
V
i is the scattering macroscopic cross-section volume averaged values for the element i, ϕ
1
,
i and ϕ
2
,
i are the first and second energy group volume averaged values of the neutron flux for the element i, and ∇
ϕ
1
-
j and ∇
ϕ
2
-
j are the first and second energy group surface averaged values of the neutron flux gradient for the face j.With regard to∇
ϕ
1
-
j and ∇
ϕ
2
-
j, although they are not known, moving Least-Squares reproducing kernel methods can express them as a weighted sum of ϕ
1
,
n and ϕ
2
,
n, respectively, by considering the elements n surrounding face j, as discussed elsewhere [6, 7]:
(5)
∇
ϕ
g
-
j
=
∑
n
k
n
,
j
grad
ϕ
g
,
n
,
g
=
1,2
,
where n is each element surrounding the face j and k
n
,
j
grad is the weighting factor of the element n with respect to the face j, which is called kernel by Harvie [6].Therefore, the final equations for each element will be(6)
∑
j
(
S
j
V
i
u
j
(
-
D
1
V
i
∑
n
k
n
,
j
grad
ϕ
1
,
n
)
)
+
Σ
a
,
1
V
i
ϕ
1
,
i
+
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
1
k
(
υ
Σ
f
,
1
V
i
ϕ
1
,
i
+
υ
Σ
f
,
2
V
i
ϕ
2
,
i
)
,
∑
j
(
S
j
V
i
u
j
(
-
D
2
V
i
∑
n
k
n
,
j
grad
ϕ
2
,
n
)
)
+
Σ
a
,
2
V
i
ϕ
2
,
i
-
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
0
.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,
(7)
ϕ
g
BC
=
0
,
g
=
1,2
,
(8)
∇
ϕ
g
BC
=
0
,
g
=
1,2
.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:(9)
ϕ
g
BC
=
∑
n
k
n
,
j
ϕ
g
,
n
=
0
,
g
=
1,2
,
(10)
∇
ϕ
g
BC
=
∑
n
k
n
,
j
grad
ϕ
g
,
n
=
0
,
g
=
1,2
,
where n represents each element surrounding the face j 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:
(11)
(
L
11
0
L
21
L
22
L
31
L
32
)
(
ϕ
1
ϕ
2
)
=
1
k
(
M
11
M
12
0
0
0
0
)
(
ϕ
1
ϕ
2
)
.L
11 entries correspond to the first equation left side of (6), M
11 and M
12 entries correspond to the first equation right side of (6), while L
21 and L
22 entries correspond to the second equation of (6), and L
31 and L
32 entries correspond to (9) and/or (10).On the other hand,ϕ
1 and ϕ
2 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 Figure1, an example of the virtual fluxes is shown. In this example the eigenvectors are
(12)
ϕ
g
=
(
ϕ
g
elements
,
ϕ
g
virtual
)
T
,
g
=
1,2
,
ϕ
g
elements
=
(
ϕ
g
1
,
ϕ
g
2
,
ϕ
g
3
,
ϕ
g
4
,
ϕ
g
5
,
ϕ
g
6
,
ϕ
g
7
,
ϕ
g
8
,
ϕ
g
9
)
T
,
g
=
1,2
,
ϕ
g
virtual
=
(
ϕ
g
10
,
ϕ
g
11
,
ϕ
g
12
,
ϕ
g
13
,
ϕ
g
14
,
ϕ
g
15
,
ϕ
g
16
,
ϕ
g
17
,
ϕ
g
18
,
ϕ
g
19
,
ϕ
g
20
,
ϕ
g
21
)
T
,
g
=
1,2
.Figure 1
Virtual fluxes example.In conclusion, the following formula represents the generalized eigenvalue problem:(13)
(
L
11
0
L
21
L
22
L
31
L
32
)
(
[
ϕ
1
elements
ϕ
1
virtual
]
[
ϕ
2
elements
ϕ
2
virtual
]
)
=
1
k
(
M
11
M
12
0
0
0
0
)
(
[
ϕ
1
elements
ϕ
1
virtual
]
[
ϕ
2
elements
ϕ
2
virtual
]
)
.


## 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 Figure2 is considered, the neutron current condition at the face separating these cells will be
(14)
J
g
+
=
J
g
-
,
g
=
1,2
.Figure 2
Adjacent cells.According to Fick’s Law,(15)
D
g
+
∇
ϕ
g
+
=
D
g
-
∇
ϕ
g
-
,
g
=
1,2
,
where D
g
+ is the g energy group diffusion coefficient for the left cell, D
g
- is the g energy group diffusion coefficient for the right cell, ∇
ϕ
g
+ is the g energy group partial neutron flux gradient for the left cell, and ∇
ϕ
g
- is the g energy group partial neutron flux gradient for the right cell.If both cells are composed of same material,D
g
+ and D
g
- 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 Figure3. As a result, a global diffusion coefficient (
D
g
j
) is needed in order to accomplish the neutron current condition by using the global gradient:
(16)
D
g
j
∇
ϕ
g
=
D
g
+
∇
ϕ
g
+
=
D
g
-
∇
ϕ
g
-
.Figure 3
Partial and global gradients.In this paper, 4 approaches of this global diffusion coefficient have been considered.(i)
Case
1 (CELL+): take the diffusion coefficient of the left cell as the global diffusion coefficient:
(17)
D
g
j
=
D
g
+
.
(ii)
Case
2 (CELL−): take the diffusion coefficient of the right cell as the global diffusion coefficient:
(18)
D
g
j
=
D
g
-
.
(iii)
Case
3 (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:
(19)
D
i
j
=
D
i
+
·
|
k
+
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
+
D
i
-
·
|
k
-
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
.
(iv)
Case
4 (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:
(20)
D
i
j
=
(
D
i
+
·
D
i
-
)
×
(
D
i
-
·
|
k
+
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
+
D
i
+
·
|
k
-
grad
|
|
k
+
grad
|
+
|
k
-
grad
|
)
-
1
.k
+
grad 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.k
-
grad 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):
(21)
∑
j
(
S
j
V
i
u
j
(
-
D
1
j
∑
n
k
n
,
j
grad
ϕ
1
,
n
)
)
+
Σ
a
,
1
V
i
ϕ
1
,
i
+
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
1
k
(
υ
Σ
f
,
1
V
i
ϕ
1
,
i
+
υ
Σ
f
,
2
V
i
ϕ
2
,
i
)
,
∑
j
(
S
j
V
i
u
j
(
-
D
2
j
∑
n
k
n
,
j
grad
ϕ
2
,
n
)
)
+
Σ
a
,
2
V
i
ϕ
2
,
i
-
Σ
s
,
1
→
2
V
i
ϕ
1
,
i
=
0
,
where D
1
j and D
2
j are the first and second energy group global diffusion coefficients for the face j.


## 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:(22)
Power
error
(
%
)
=
|
Power
-
Power
reference
|
Power
reference
*
100
,
Eigenvalue
error
(
pcm
)
=
|
k
-
k
reference
|
k
reference
*
10
5
,
Power
=
∑
i
(
υ
Σ
f
,
1
V
i
ϕ
1
,
i
+
υ
Σ
f
,
2
V
i
ϕ
2
,
i
)
.In addition, power is normalized to attain that mean power equals the unity, calculated with the following formula:(23)
Mean
Power
=
∑
i
=
1
N
Power
i
V
i
∑
i
=
1
N
V
i
,
where V
i is volume or area of the element i, 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 Table1. In addition, it has a rectangular shape and its dimensions are 100
cm
×
60
cm. Furthermore, zero flux boundary conditions have been considered.Table 1
2D homogeneous reactor cross sections.
D
1 (cm)
D
2 (cm)
Σ
a
1 (cm−1)
Σ
a
2 (cm−1)
υ
Σ
f
1 (cm−1)
υ
Σ
f
2 (cm−1)
Σ
12 (cm−1)
1.28205128205
0.666667
0.01
0.1
0.01
0.109017634020268
0.075On 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. Table2 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.Table 2
2D homogeneous reactor results.
Structured mesh
Unstructured mesh
Mesh length (cm)
10
5
2
10
5
2
Number of cells
72
288
1620
234
666
4054
Computational time (s)
0.38
0.85
4.11
0.74
1.52
10.41
Eigenvalue-1 error (pcm)
139.903827
35.083823
5.733821
29.653822
12.923822
2.053821
Eigenvalue-2 error (pcm)
309.789649
77.398527
14.010938
117.773860
38.609075
6.124990
Eigenvalue-3 error (pcm)
1207.774879
246.016129
47.219525
437.256280
179.905306
18.743235
Eigenvalue-4 error (pcm)
1587.045813
462.009487
74.560259
221.291968
128.800577
20.948371
Eigenvalue-5 error (pcm)
1948.227429
492.735800
80.710670
436.132079
185.870476
29.279108Figure 4
2D structured mesh.Figure 5
2D unstructured mesh.Moreover, this problem has analytical solution and it will be the reference solution. The analytical eigenvaluesk 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 k is virtually 1.The results corresponding to computational time and eigenvalue errors are exposed in Table2. On the other hand, power errors are evaluated in 9 rectangular nodes of the same dimensions 33
.
33
cm
×
20
cm, 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.Table 3
Numbering of the 2D homogeneous reactor nodes.
7
8
9
4
5
6
1
2
3Table 4
2D homogeneous reactor power errors (%) corresponding to the first eigenvalue.
Node
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
1
0.009206
0.000081
0.000069
0.094629
0.131237
0.018789
2
0.005681
0.001217
0.000559
0.005927
0.002009
0.005058
3
0.002318
0.002675
0.001150
0.098431
0.136212
0.015750
4
0.010802
0.001601
0.000126
0.091223
0.125323
0.007585
5
0.014169
0.003062
0.000684
0.000774
0.005144
0.001788
6
0.017709
0.004348
0.001148
0.073893
0.111893
0.015444
7
0.009206
0.000081
0.000058
0.092636
0.096630
0.020416
8
0.005672
0.001208
0.000577
0.006823
0.000913
0.001457
9
0.002318
0.002686
0.001182
0.082769
0.141282
0.007689Table 5
2D homogeneous reactor power errors (%) corresponding to the second eigenvalue.
Node
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
1
0.010699
0.002026
0.000177
0.064440
0.131450
0.023185
2
—
—
—
—
—
—
3
0.011767
0.000900
0.000735
0.218325
0.056848
0.013567
4
0.013305
0.002670
0.000291
0.370735
0.142994
0.009263
5
—
—
—
—
—
—
6
0.012219
0.003781
0.001195
0.270132
0.000779
0.021504
7
0.010706
0.002026
0.000135
0.116602
0.091269
0.018453
8
—
—
—
—
—
—
9
0.011767
0.000900
0.000693
0.228939
0.057844
0.029325Figure 6
2D homogeneous reactor power corresponding to the first eigenvalue and the 2 cm structured mesh.Figure 7
2D homogeneous reactor power corresponding to the first eigenvalue and the 2 cm unstructured mesh.Figure 8
2D homogeneous reactor power corresponding to the second eigenvalue and the 2 cm structured mesh.Figure 9
2D homogeneous reactor power corresponding to the second eigenvalue and the 2 cm unstructured mesh.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 Table5. 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 Figure10 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.Table 6
Biblis reactor cross sections.
Material
D
1 (cm)
D
2 (cm)
Σ
a
1 (cm−1)
Σ
a
2 (cm−1)
Σ
12 (cm−1)
υ
Σ
f
1 (cm−1)
υ
Σ
f
2 (cm−1)
1
1.4360
0.3635
0.0095042
0.0750580
0.017754
0.0058708
0.0960670
2
1.4366
0.3636
0.0096785
0.0784360
0.017621
0.0061908
0.1035800
3
1.3200
0.2772
0.0026562
0.0715960
0.023106
0.0
0.0
4
1.4389
0.3638
0.0103630
0.0914080
0.017101
0.0074527
0.1323600
5
1.4381
0.3665
0.0100030
0.0848280
0.017290
0.0061908
0.1035800
6
1.4385
0.3665
0.0101320
0.0873140
0.017192
0.0064285
0.1091100
7
1.4389
0.3679
0.0101650
0.0880240
0.017125
0.0061908
0.1035800
8
1.4393
0.3680
0.0102940
0.0905100
0.017027
0.0064285
0.1091100Figure 10
Biblis reactor: geometry and material composition.Moreover, the same meshes as in Section3.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 4
*
4 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 (
k
1
) is 1.025110.Table 7
Biblis reactor meshes.
Structured mesh
Unstructured mesh
Mesh length (cm)
10
5
2
10
5
2
Number of cells
1028
2416
9252
1808
4542
23972Table 8
Biblis reactor computational time (minutes:seconds).
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
CELL+
0:3.17
0:6.78
0:40.31
0:4.02
0:12.79
3:33.1
CELL−
0:2.73
0:6.75
0:40.48
0:4.03
0:12.6
3:42.44
HOM
0:2.71
0:6.79
0:40.74
0:4.43
0:12.64
3:34.98
LIN
0:2.73
0:6.94
0:40.15
0:3.99
0:12.9
3:44.03Table 9
Biblis reactor eigenvalue-1 error (pcm).
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
CELL+
97.882178
42.210104
11.208553
49.516637
19.919813
4.536099
CELL−
40.941948
7.784530
3.521573
15.305675
0.097551
1.951010
HOM
67.641521
23.968160
3.463043
31.411263
9.452644
1.160851
LIN
69.182819
24.797339
3.726429
32.250198
9.862356
1.258402In 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 Figure10, 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.Table 10
Numbering of the Biblis reactor nodes.
49
50
51
52
53
54
55
56
—
41
42
43
44
45
46
47
48
—
33
34
35
36
37
38
39
40
—
25
26
27
28
29
30
31
32
—
18
19
20
21
22
23
24
—
—
11
12
13
14
15
16
17
—
5
6
7
8
9
10
—
—
1
2
3
4
—
—
—
—
—
—
—
—Only power errors corresponding to HOM are exposed in Table11, 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.Table 11
Biblis reactor power errors (%) corresponding to HOM.
Node
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
1
2.284672
0.001124
0.367491
0.879946
0.030794
0.180845
2
2.304411
0.135778
0.367554
0.824867
0.142893
0.191021
3
2.376446
0.134061
0.376089
0.864341
0.162517
0.180412
4
2.376446
0.134061
0.376089
0.864341
0.162517
0.180412
5
2.488581
0.003299
0.489467
0.429311
0.390778
0.228766
6
1.074203
0.220179
0.038275
0.588963
0.054508
0.060416
7
0.934102
0.176358
0.083588
0.479943
0.043084
0.086734
8
1.469022
0.483910
0.021256
0.658828
0.161010
0.032085
9
0.027694
0.558208
0.369578
0.276270
0.457187
0.183278
10
1.969739
0.070576
0.335412
0.446243
0.289564
0.156511
11
2.755899
0.126755
0.433999
0.793181
0.334266
0.202575
12
1.849153
0.858648
0.268538
1.049713
0.424649
0.110011
13
0.442460
0.372698
0.185346
0.282549
0.266245
0.089598
14
1.537896
0.814179
0.251776
0.831365
0.408647
0.095628
15
1.077245
0.592678
0.167885
0.500974
0.281762
0.078132
16
0.221870
0.218571
0.199870
0.032188
0.187027
0.063343
17
2.075160
0.300133
0.189122
0.755988
0.026896
0.045875
18
2.755899
0.126755
0.433972
0.840364
0.366588
0.198116
19
0.014891
0.540713
0.404170
0.298276
0.443251
0.188347
20
2.219547
0.768456
0.099073
1.019476
0.302488
0.036096
21
0.295151
0.629288
0.353496
0.274481
0.436990
0.164800
22
1.672244
0.729190
0.169991
0.865452
0.342074
0.059898
23
0.882870
0.491431
0.140994
0.374551
0.237709
0.074700
24
0.221876
0.218560
0.199864
0.043428
0.188407
0.065455
25
1.969739
0.070576
0.335412
0.473057
0.317933
0.151979
26
2.189996
0.490979
0.095543
0.948384
0.057490
0.051558
27
0.351107
0.582043
0.464038
0.056930
0.475770
0.206069
28
2.368009
0.709176
0.016244
1.057085
0.220311
0.005560
29
0.223202
0.671285
0.391755
0.271355
0.448080
0.172726
30
1.672244
0.729190
0.169991
0.855210
0.346483
0.059600
31
1.077245
0.592678
0.167868
0.535406
0.269561
0.079436
32
0.027702
0.558224
0.369601
0.231229
0.486828
0.179438
33
2.488581
0.003299
0.489478
0.525898
0.439935
0.231726
34
0.717648
0.496628
0.533502
0.053316
0.522542
0.251418
35
2.951741
0.772375
0.049890
1.260491
0.198974
0.034441
36
0.382809
0.717441
0.545718
0.150071
0.595033
0.251171
37
2.368009
0.709176
0.016217
1.041303
0.221315
0.004546
38
0.295151
0.629279
0.353451
0.301348
0.431097
0.163978
39
1.537890
0.814179
0.251782
0.790510
0.415010
0.094983
40
1.469022
0.483910
0.021268
0.718255
0.135214
0.031668
41
2.376439
0.134061
0.376089
1.033266
0.223136
0.182892
42
3.220990
0.751593
0.089823
1.413971
0.171813
0.054183
43
0.764614
0.595606
0.606721
0.003153
0.592573
0.276228
44
2.951741
0.772375
0.049916
1.244259
0.210189
0.034626
45
0.351107
0.582043
0.464013
0.098390
0.470016
0.204001
46
2.219547
0.768446
0.099092
0.984197
0.309352
0.036133
47
0.442460
0.372698
0.185365
0.357257
0.261884
0.086875
48
0.943520
0.185705
0.074237
0.575157
0.042333
0.074883
49
2.304417
0.135772
0.367517
1.003710
0.178850
0.186330
50
0.861206
0.553378
0.631204
0.085637
0.536703
0.290387
51
3.220990
0.751593
0.089823
1.405356
0.184740
0.055653
52
0.717648
0.496628
0.533502
0.092312
0.518655
0.248030
53
2.189996
0.490979
0.095543
0.910271
0.061734
0.046650
54
0.014881
0.540713
0.404161
0.357925
0.437896
0.182018
55
1.849153
0.858642
0.268519
0.983297
0.433686
0.112597
56
2.284672
0.001134
0.367482
1.015056
0.047411
0.185290Figure 11
Biblis reactor power corresponding to the first eigenvalue and the 2 cm structured mesh and HOM.Figure 12
Biblis reactor power corresponding to the first eigenvalue and the 2 cm unstructured mesh and HOM.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:100
cm
×
60
cm
×
180
cm. 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 Table12. On the other hand, power errors are evaluated in 54 parallelepiped nodes of the same dimensions 33.33
cm
×
20
cm
×
30
cm, 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.Table 12
3D homogeneous reactor results.
Structured mesh
Mesh length (cm)
10
5
2
Number of cells
1080
8640
135000
Computational time (h:min:s)
0:0:8.56
0:2:51.75
12:41:9
Eigenvalue-1 error (pcm)
148.492003
37.141902
5.927090
Eigenvalue-2 error (pcm)
174.642897
43.473709
6.912315
Eigenvalue-3 error (pcm)
276.887015
68.855621
10.939931
Eigenvalue-4 error (pcm)
390.546564
97.622804
15.555099
Eigenvalue-5 error (pcm)
416.479277
103.762637
16.498151Table 13
Numbering of 3D homogeneous reactor nodes.
1st axial plane
2nd axial plane
3rd axial plane
4th axial plane
5th axial plane
6th axial plane
7
8
9
16
17
18
25
26
27
34
35
36
43
44
45
52
53
54
4
5
6
13
14
15
22
23
24
31
32
33
40
41
42
49
50
51
1
2
3
10
11
12
19
20
21
28
29
30
37
38
39
46
47
48Table 14
3D homogeneous reactor power errors (%) corresponding to the first eigenvalue.
Node
Structured mesh
10 cm
5 cm
2 cm
1
0.982818
0.283377
0.044883
2
1.021300
0.279512
0.044440
3
0.982822
0.283386
0.044883
4
1.067796
0.290603
0.045257
5
0.991369
0.277265
0.044169
6
1.067802
0.290610
0.045211
7
0.982818
0.283364
0.044932
8
1.021300
0.279545
0.044217
9
0.982822
0.283373
0.044914
10
1.020863
0.283123
0.044656
11
1.011269
0.282063
0.044683
12
1.020863
0.283113
0.044617
13
1.042381
0.284688
0.044685
14
1.043681
0.285467
0.044857
15
1.042373
0.284664
0.044693
16
1.020863
0.283142
0.044694
17
1.011261
0.282119
0.044569
18
1.020863
0.283103
0.044712
19
1.015342
0.281914
0.044589
20
1.017327
0.283325
0.044599
21
1.015328
0.281928
0.044560
22
1.036297
0.283471
0.044532
23
1.050299
0.286768
0.044747
24
1.036291
0.283471
0.044526
25
1.015350
0.281900
0.044489
26
1.017315
0.283308
0.044682
27
1.015342
0.281893
0.044503
28
1.015342
0.281921
0.044582
29
1.017327
0.283337
0.044599
30
1.015328
0.281921
0.044567
31
1.036297
0.283453
0.044502
32
1.050299
0.286774
0.044753
33
1.036291
0.283471
0.044526
34
1.015350
0.281893
0.044483
35
1.017315
0.283290
0.044688
36
1.015350
0.281900
0.044510
37
1.020863
0.283142
0.044656
38
1.011269
0.282063
0.044674
39
1.020863
0.283113
0.044626
40
1.042381
0.284688
0.044677
41
1.043681
0.285467
0.044857
42
1.042381
0.284672
0.044685
43
1.020863
0.283132
0.044712
44
1.011261
0.282103
0.044569
45
1.020863
0.283093
0.044712
46
0.982822
0.283373
0.044888
47
1.021306
0.279519
0.044453
48
0.982809
0.283377
0.044883
49
1.067796
0.290610
0.045244
50
0.991362
0.277277
0.044169
51
1.067802
0.290616
0.045197
52
0.982822
0.283364
0.044914
53
1.021306
0.279545
0.044210
54
0.982822
0.283373
0.044905Table 15
3D homogeneous reactor power errors (%) corresponding to the second eigenvalue.
Node
Structured mesh
10 cm
5 cm
2 cm
1
0.996377
0.284899
0.045155
2
1.019201
0.279652
0.044237
3
0.996370
0.284906
0.045097
4
1.062091
0.289487
0.045069
5
1.008268
0.280020
0.044300
6
1.062091
0.289473
0.045105
7
0.996377
0.284871
0.045155
8
1.019187
0.279666
0.044165
9
0.996377
0.284856
0.045155
10
1.037403
0.284821
0.044523
11
1.008552
0.282182
0.044436
12
1.037403
0.284835
0.044473
13
1.034915
0.283264
0.044419
14
1.064620
0.288742
0.044959
15
1.034915
0.283235
0.044389
16
1.037410
0.284806
0.044658
17
1.008552
0.282277
0.044626
18
1.037403
0.284813
0.044729
19
1.031311
0.283454
0.044715
20
1.015537
0.283716
0.044796
21
1.031304
0.283440
0.044693
22
1.027935
0.281797
0.044396
23
1.072478
0.290229
0.045138
24
1.027935
0.281790
0.044459
25
1.031311
0.283412
0.044765
26
1.015558
0.283695
0.044655
27
1.031318
0.283412
0.044686
28
1.031311
0.283462
0.044715
29
1.015537
0.283716
0.044782
30
1.031304
0.283440
0.044686
31
1.027928
0.281804
0.044403
32
1.072478
0.290224
0.045144
33
1.027928
0.281797
0.044452
34
1.031311
0.283419
0.044772
35
1.015558
0.283695
0.044655
36
1.031311
0.283419
0.044679
37
1.037403
0.284813
0.044516
38
1.008552
0.282182
0.044442
39
1.037395
0.284828
0.044466
40
1.034915
0.283258
0.044419
41
1.064620
0.288751
0.044959
42
1.034915
0.283235
0.044389
43
1.037403
0.284806
0.044650
44
1.008552
0.282289
0.044626
45
1.037403
0.284799
0.044729
46
0.996377
0.284892
0.045155
47
1.019208
0.279659
0.044244
48
0.996384
0.284913
0.045097
49
1.062091
0.289473
0.045076
50
1.008250
0.280026
0.044300
51
1.062091
0.289458
0.045119
52
0.996377
0.284871
0.045162
53
1.019194
0.279659
0.044151
54
0.996377
0.284871
0.045155Figure 13
3D homogeneous reactor power corresponding to the first eigenvalue and the 2 cm structured mesh.Figure 14
3D homogeneous reactor power corresponding to the second eigenvalue and the 2 cm structured mesh.The analytical eigenvaluesk 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 k 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 Table16. 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.Table 16
Langenbuch reactor cross sections.
Material
D
1 (cm)
D
2 (cm)
Σ
a
1 (cm−1)
Σ
a
2 (cm−1)
Σ
12 (cm−1)
υ
Σ
f
1 (cm−1)
υ
Σ
f
2 (cm−1)
Comb.1
1.423913
0.356306
0.01040206
0.08766217
0.0175555
0.006477691
0.1127328
Comb.2
1.425611
0.350574
0.01099263
0.09925634
0.01717768
0.007503284
0.1378004
Absorbent
1.423913
0.356306
0.01095206
0.09146217
0.0175555
0.006477691
0.11273228
Reflector
1.634227
0.264002
0.002660573
0.04936351
0.02759693
0.0
0.0Figure 15
Axial plane of Langenbuch reactor.Figure 16
Frontal plane of Langenbuch reactor.In this case, only structured meshes have been used, in particular the 10 cm and 5 cm mesh length. Table17 shows the number of cells for each mesh.Table 17
Langenbuch reactor meshes.
Structured mesh
Mesh length (cm)
10
5
Number of cells
2800
18720Moreover, 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 (k
1 and k
2) are 0.994881227 and 0.948210698.Regarding the results, Table18 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.Table 18
Langenbuch reactor computational time (h:min:s).
Structured mesh
10 cm
5 cm
CELL+
0:0:25.69
0:10:54.5
CELL−
0:0:25.68
0:11:23.11
HOM
0:0:25.77
0:10:55.77
LIN
0:0:25.73
0:10:43.59Table 19
Langenbuch reactor eigenvalue-1 error (pcm).
Structured mesh
10 cm
5 cm
CELL+
83.3771
12.5047
CELL−
212.8518
63.0611
HOM
143.1450
22.6874
LIN
150.0886
25.2636Table 20
Langenbuch reactor eigenvalue-2 error (pcm).
Structured mesh
10 cm
5 cm
CELL+
226.1219
7.7211
CELL−
350.5162
80.2113
HOM
280.9989
39.9871
LIN
292.1948
44.2362Table 21
Numbering of the Langenbuch reactor nodes.
2nd axial plane
5th axial plane
—
—
—
—
—
—
—
—
—
—
21
22
23
24
—
—
45
46
47
48
—
—
16
17
18
19
20
—
40
41
42
43
44
—
11
12
13
14
15
—
35
36
37
38
39
—
6
7
8
9
10
—
30
31
32
33
34
—
1
2
3
4
5
—
25
26
27
28
29
—Table 22
Langenbuch reactor power errors (%) corresponding to the first eigenvalue and HOM.
Node
Structured mesh
10 cm
5 cm
1
0,88732759
0,083512937
2
0,413636006
0,169843064
3
0,285898499
0,26806574
4
1,650354128
0,463970268
5
7,308155289
2,443154619
6
0,413578042
0,169901321
7
0,218154684
0,211939733
8
0,89582742
0,279213792
9
2,440629291
0,498489677
10
8,245123598
2,368554084
11
0,285938524
0,268114276
12
0,895860783
0,279229236
13
1,485778349
0,264874283
14
3,250361067
0,503872312
15
8,254395078
1,996055827
16
1,650484996
0,46408882
17
2,440737883
0,498585195
18
3,250407371
0,5039109
19
6,540873894
1,269229237
20
10,49036132
3,000440824
21
7,308460522
2,44345405
22
8,245387401
2,368777266
23
8,25461267
1,99626084
24
10,49053455
3,000606939
25
5,466731896
0,931224354
26
4,140848872
0,716286587
27
3,657496267
0,698504374
28
2,237579532
0,474563288
29
3,189694202
1,268772356
30
4,140836529
0,716252343
31
2,96775873
0,633299552
32
2,476161123
0,618678432
33
0,910178355
0,364113213
34
4,548866137
1,200011218
35
3,657467753
0,698474984
36
2,476131091
0,618647924
37
2,967442105
0,56874817
38
0,292240082
0,412798339
39
4,267581802
0,724868524
40
2,237496931
0,474473118
41
0,910102518
0,364043283
42
0,292200002
0,412751132
43
2,930647399
0,281476818
44
6,41914469
1,657371938
45
3,189876024
1,268969435
46
4,549015577
1,2001458
47
4,267727516
0,724995924
48
6,419233002
1,657446439Table 23
Langenbuch reactor power errors (%) corresponding to the second eigenvalue and HOM.
Node
Structured mesh
10 cm
5 cm
1
2,748470704
0,472742419
2
1,064854909
0,505072706
3
0,379827902
0,520035451
4
0,846862035
0,527032189
5
6,400564495
2,241810405
6
1,064870054
0,505064685
7
0,500339771
0,504531289
8
1,13600946
0,496430539
9
2,54996541
0,525233749
10
8,338474277
2,124761105
11
0,379868741
0,519985959
12
1,135973249
0,496377465
13
1,719135726
0,413130915
14
3,303282255
0,426080951
15
8,233916013
1,633618147
16
0,846698163
0,526850063
17
2,549801498
0,525079406
18
3,303171711
0,426002352
19
6,461736033
0,930426868
20
10,39493074
2,481116214
21
6,400329869
2,24158495
22
8,338216136
2,124537285
23
8,233749447
1,633473241
24
10,39486316
2,481044325
25
8,712874235
0,839916928
26
6,327127
0,784501724
27
5,950753659
0,862124448
28
4,44225414
0,601601562
29
0,062727543
0,461900883
30
6,327225132
0,78459806
31
4,313496832
0,803439342
32
3,885076646
0,873376122
33
2,432393528
0,789928103
34
2,461696614
0,258876027
35
5,951087173
0,86245633
36
3,885366412
0,873674968
37
4,2237795
0,779967598
38
1,90610207
0,962348979
39
2,22742644
0,28534924
40
4,44301111
0,602394933
41
2,433102829
0,790642934
42
1,906524369
0,962753237
43
0,944598463
0,649099892
44
4,201989073
0,422463819
45
0,061633966
0,460821104
46
2,460613715
0,257845259
47
2,226705013
0,286064257
48
4,20167533
0,422178482Figure 17
Langenbuch reactor power corresponding to the first eigenvalue and the 5 cm structured mesh and HOM.Figure 18
Langenbuch reactor power corresponding to the second eigenvalue and the 5 cm structured mesh and HOM.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.


## 3.1. 2D Homogeneous Reactor

The reactor considered is composed only of one material, whose cross sections are shown in Table1. In addition, it has a rectangular shape and its dimensions are 100
cm
×
60
cm. Furthermore, zero flux boundary conditions have been considered.Table 1
2D homogeneous reactor cross sections.
D
1 (cm)
D
2 (cm)
Σ
a
1 (cm−1)
Σ
a
2 (cm−1)
υ
Σ
f
1 (cm−1)
υ
Σ
f
2 (cm−1)
Σ
12 (cm−1)
1.28205128205
0.666667
0.01
0.1
0.01
0.109017634020268
0.075On 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. Table2 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.Table 2
2D homogeneous reactor results.
Structured mesh
Unstructured mesh
Mesh length (cm)
10
5
2
10
5
2
Number of cells
72
288
1620
234
666
4054
Computational time (s)
0.38
0.85
4.11
0.74
1.52
10.41
Eigenvalue-1 error (pcm)
139.903827
35.083823
5.733821
29.653822
12.923822
2.053821
Eigenvalue-2 error (pcm)
309.789649
77.398527
14.010938
117.773860
38.609075
6.124990
Eigenvalue-3 error (pcm)
1207.774879
246.016129
47.219525
437.256280
179.905306
18.743235
Eigenvalue-4 error (pcm)
1587.045813
462.009487
74.560259
221.291968
128.800577
20.948371
Eigenvalue-5 error (pcm)
1948.227429
492.735800
80.710670
436.132079
185.870476
29.279108Figure 4
2D structured mesh.Figure 5
2D unstructured mesh.Moreover, this problem has analytical solution and it will be the reference solution. The analytical eigenvaluesk 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 k is virtually 1.The results corresponding to computational time and eigenvalue errors are exposed in Table2. On the other hand, power errors are evaluated in 9 rectangular nodes of the same dimensions 33
.
33
cm
×
20
cm, 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.Table 3
Numbering of the 2D homogeneous reactor nodes.
7
8
9
4
5
6
1
2
3Table 4
2D homogeneous reactor power errors (%) corresponding to the first eigenvalue.
Node
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
1
0.009206
0.000081
0.000069
0.094629
0.131237
0.018789
2
0.005681
0.001217
0.000559
0.005927
0.002009
0.005058
3
0.002318
0.002675
0.001150
0.098431
0.136212
0.015750
4
0.010802
0.001601
0.000126
0.091223
0.125323
0.007585
5
0.014169
0.003062
0.000684
0.000774
0.005144
0.001788
6
0.017709
0.004348
0.001148
0.073893
0.111893
0.015444
7
0.009206
0.000081
0.000058
0.092636
0.096630
0.020416
8
0.005672
0.001208
0.000577
0.006823
0.000913
0.001457
9
0.002318
0.002686
0.001182
0.082769
0.141282
0.007689Table 5
2D homogeneous reactor power errors (%) corresponding to the second eigenvalue.
Node
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
1
0.010699
0.002026
0.000177
0.064440
0.131450
0.023185
2
—
—
—
—
—
—
3
0.011767
0.000900
0.000735
0.218325
0.056848
0.013567
4
0.013305
0.002670
0.000291
0.370735
0.142994
0.009263
5
—
—
—
—
—
—
6
0.012219
0.003781
0.001195
0.270132
0.000779
0.021504
7
0.010706
0.002026
0.000135
0.116602
0.091269
0.018453
8
—
—
—
—
—
—
9
0.011767
0.000900
0.000693
0.228939
0.057844
0.029325Figure 6
2D homogeneous reactor power corresponding to the first eigenvalue and the 2 cm structured mesh.Figure 7
2D homogeneous reactor power corresponding to the first eigenvalue and the 2 cm unstructured mesh.Figure 8
2D homogeneous reactor power corresponding to the second eigenvalue and the 2 cm structured mesh.Figure 9
2D homogeneous reactor power corresponding to the second eigenvalue and the 2 cm unstructured mesh.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 Table5. 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 Figure10 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.Table 6
Biblis reactor cross sections.
Material
D
1 (cm)
D
2 (cm)
Σ
a
1 (cm−1)
Σ
a
2 (cm−1)
Σ
12 (cm−1)
υ
Σ
f
1 (cm−1)
υ
Σ
f
2 (cm−1)
1
1.4360
0.3635
0.0095042
0.0750580
0.017754
0.0058708
0.0960670
2
1.4366
0.3636
0.0096785
0.0784360
0.017621
0.0061908
0.1035800
3
1.3200
0.2772
0.0026562
0.0715960
0.023106
0.0
0.0
4
1.4389
0.3638
0.0103630
0.0914080
0.017101
0.0074527
0.1323600
5
1.4381
0.3665
0.0100030
0.0848280
0.017290
0.0061908
0.1035800
6
1.4385
0.3665
0.0101320
0.0873140
0.017192
0.0064285
0.1091100
7
1.4389
0.3679
0.0101650
0.0880240
0.017125
0.0061908
0.1035800
8
1.4393
0.3680
0.0102940
0.0905100
0.017027
0.0064285
0.1091100Figure 10
Biblis reactor: geometry and material composition.Moreover, the same meshes as in Section3.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 4
*
4 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 (
k
1
) is 1.025110.Table 7
Biblis reactor meshes.
Structured mesh
Unstructured mesh
Mesh length (cm)
10
5
2
10
5
2
Number of cells
1028
2416
9252
1808
4542
23972Table 8
Biblis reactor computational time (minutes:seconds).
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
CELL+
0:3.17
0:6.78
0:40.31
0:4.02
0:12.79
3:33.1
CELL−
0:2.73
0:6.75
0:40.48
0:4.03
0:12.6
3:42.44
HOM
0:2.71
0:6.79
0:40.74
0:4.43
0:12.64
3:34.98
LIN
0:2.73
0:6.94
0:40.15
0:3.99
0:12.9
3:44.03Table 9
Biblis reactor eigenvalue-1 error (pcm).
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
CELL+
97.882178
42.210104
11.208553
49.516637
19.919813
4.536099
CELL−
40.941948
7.784530
3.521573
15.305675
0.097551
1.951010
HOM
67.641521
23.968160
3.463043
31.411263
9.452644
1.160851
LIN
69.182819
24.797339
3.726429
32.250198
9.862356
1.258402In 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 Figure10, 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.Table 10
Numbering of the Biblis reactor nodes.
49
50
51
52
53
54
55
56
—
41
42
43
44
45
46
47
48
—
33
34
35
36
37
38
39
40
—
25
26
27
28
29
30
31
32
—
18
19
20
21
22
23
24
—
—
11
12
13
14
15
16
17
—
5
6
7
8
9
10
—
—
1
2
3
4
—
—
—
—
—
—
—
—Only power errors corresponding to HOM are exposed in Table11, 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.Table 11
Biblis reactor power errors (%) corresponding to HOM.
Node
Structured mesh
Unstructured mesh
10 cm
5 cm
2 cm
10 cm
5 cm
2 cm
1
2.284672
0.001124
0.367491
0.879946
0.030794
0.180845
2
2.304411
0.135778
0.367554
0.824867
0.142893
0.191021
3
2.376446
0.134061
0.376089
0.864341
0.162517
0.180412
4
2.376446
0.134061
0.376089
0.864341
0.162517
0.180412
5
2.488581
0.003299
0.489467
0.429311
0.390778
0.228766
6
1.074203
0.220179
0.038275
0.588963
0.054508
0.060416
7
0.934102
0.176358
0.083588
0.479943
0.043084
0.086734
8
1.469022
0.483910
0.021256
0.658828
0.161010
0.032085
9
0.027694
0.558208
0.369578
0.276270
0.457187
0.183278
10
1.969739
0.070576
0.335412
0.446243
0.289564
0.156511
11
2.755899
0.126755
0.433999
0.793181
0.334266
0.202575
12
1.849153
0.858648
0.268538
1.049713
0.424649
0.110011
13
0.442460
0.372698
0.185346
0.282549
0.266245
0.089598
14
1.537896
0.814179
0.251776
0.831365
0.408647
0.095628
15
1.077245
0.592678
0.167885
0.500974
0.281762
0.078132
16
0.221870
0.218571
0.199870
0.032188
0.187027
0.063343
17
2.075160
0.300133
0.189122
0.755988
0.026896
0.045875
18
2.755899
0.126755
0.433972
0.840364
0.366588
0.198116
19
0.014891
0.540713
0.404170
0.298276
0.443251
0.188347
20
2.219547
0.768456
0.099073
1.019476
0.302488
0.036096
21
0.295151
0.629288
0.353496
0.274481
0.436990
0.164800
22
1.672244
0.729190
0.169991
0.865452
0.342074
0.059898
23
0.882870
0.491431
0.140994
0.374551
0.237709
0.074700
24
0.221876
0.218560
0.199864
0.043428
0.188407
0.065455
25
1.969739
0.070576
0.335412
0.473057
0.317933
0.151979
26
2.189996
0.490979
0.095543
0.948384
0.057490
0.051558
27
0.351107
0.582043
0.464038
0.056930
0.475770
0.206069
28
2.368009
0.709176
0.016244
1.057085
0.220311
0.005560
29
0.223202
0.671285
0.391755
0.271355
0.448080
0.172726
30
1.672244
0.729190
0.169991
0.855210
0.346483
0.059600
31
1.077245
0.592678
0.167868
0.535406
0.269561
0.079436
32
0.027702
0.558224
0.369601
0.231229
0.486828
0.179438
33
2.488581
0.003299
0.489478
0.525898
0.439935
0.231726
34
0.717648
0.496628
0.533502
0.053316
0.522542
0.251418
35
2.951741
0.772375
0.049890
1.260491
0.198974
0.034441
36
0.382809
0.717441
0.545718
0.150071
0.595033
0.251171
37
2.368009
0.709176
0.016217
1.041303
0.221315
0.004546
38
0.295151
0.629279
0.353451
0.301348
0.431097
0.163978
39
1.537890
0.814179
0.251782
0.790510
0.415010
0.094983
40
1.469022
0.483910
0.021268
0.718255
0.135214
0.031668
41
2.376439
0.134061
0.376089
1.033266
0.223136
0.182892
42
3.220990
0.751593
0.089823
1.413971
0.171813
0.054183
43
0.764614
0.595606
0.606721
0.003153
0.592573
0.276228
44
2.951741
0.772375
0.049916
1.244259
0.210189
0.034626
45
0.351107
0.582043
0.464013
0.098390
0.470016
0.204001
46
2.219547
0.768446
0.099092
0.984197
0.309352
0.036133
47
0.442460
0.372698
0.185365
0.357257
0.261884
0.086875
48
0.943520
0.185705
0.074237
0.575157
0.042333
0.074883
49
2.304417
0.135772
0.367517
1.003710
0.178850
0.186330
50
0.861206
0.553378
0.631204
0.085637
0.536703
0.290387
51
3.220990
0.751593
0.089823
1.405356
0.184740
0.055653
52
0.717648
0.496628
0.533502
0.092312
0.518655
0.248030
53
2.189996
0.490979
0.095543
0.910271
0.061734
0.046650
54
0.014881
0.540713
0.404161
0.357925
0.437896
0.182018
55
1.849153
0.858642
0.268519
0.983297
0.433686
0.112597
56
2.284672
0.001134
0.367482
1.015056
0.047411
0.185290Figure 11
Biblis reactor power corresponding to the first eigenvalue and the 2 cm structured mesh and HOM.Figure 12
Biblis reactor power corresponding to the first eigenvalue and the 2 cm unstructured mesh and HOM.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:100
cm
×
60
cm
×
180
cm. 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 Table12. On the other hand, power errors are evaluated in 54 parallelepiped nodes of the same dimensions 33.33
cm
×
20
cm
×
30
cm, 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.Table 12
3D homogeneous reactor results.
Structured mesh
Mesh length (cm)
10
5
2
Number of cells
1080
8640
135000
Computational time (h:min:s)
0:0:8.56
0:2:51.75
12:41:9
Eigenvalue-1 error (pcm)
148.492003
37.141902
5.927090
Eigenvalue-2 error (pcm)
174.642897
43.473709
6.912315
Eigenvalue-3 error (pcm)
276.887015
68.855621
10.939931
Eigenvalue-4 error (pcm)
390.546564
97.622804
15.555099
Eigenvalue-5 error (pcm)
416.479277
103.762637
16.498151Table 13
Numbering of 3D homogeneous reactor nodes.
1st axial plane
2nd axial plane
3rd axial plane
4th axial plane
5th axial plane
6th axial plane
7
8
9
16
17
18
25
26
27
34
35
36
43
44
45
52
53
54
4
5
6
13
14
15
22
23
24
31
32
33
40
41
42
49
50
51
1
2
3
10
11
12
19
20
21
28
29
30
37
38
39
46
47
48Table 14
3D homogeneous reactor power errors (%) corresponding to the first eigenvalue.
Node
Structured mesh
10 cm
5 cm
2 cm
1
0.982818
0.283377
0.044883
2
1.021300
0.279512
0.044440
3
0.982822
0.283386
0.044883
4
1.067796
0.290603
0.045257
5
0.991369
0.277265
0.044169
6
1.067802
0.290610
0.045211
7
0.982818
0.283364
0.044932
8
1.021300
0.279545
0.044217
9
0.982822
0.283373
0.044914
10
1.020863
0.283123
0.044656
11
1.011269
0.282063
0.044683
12
1.020863
0.283113
0.044617
13
1.042381
0.284688
0.044685
14
1.043681
0.285467
0.044857
15
1.042373
0.284664
0.044693
16
1.020863
0.283142
0.044694
17
1.011261
0.282119
0.044569
18
1.020863
0.283103
0.044712
19
1.015342
0.281914
0.044589
20
1.017327
0.283325
0.044599
21
1.015328
0.281928
0.044560
22
1.036297
0.283471
0.044532
23
1.050299
0.286768
0.044747
24
1.036291
0.283471
0.044526
25
1.015350
0.281900
0.044489
26
1.017315
0.283308
0.044682
27
1.015342
0.281893
0.044503
28
1.015342
0.281921
0.044582
29
1.017327
0.283337
0.044599
30
1.015328
0.281921
0.044567
31
1.036297
0.283453
0.044502
32
1.050299
0.286774
0.044753
33
1.036291
0.283471
0.044526
34
1.015350
0.281893
0.044483
35
1.017315
0.283290
0.044688
36
1.015350
0.281900
0.044510
37
1.020863
0.283142
0.044656
38
1.011269
0.282063
0.044674
39
1.020863
0.283113
0.044626
40
1.042381
0.284688
0.044677
41
1.043681
0.285467
0.044857
42
1.042381
0.284672
0.044685
43
1.020863
0.283132
0.044712
44
1.011261
0.282103
0.044569
45
1.020863
0.283093
0.044712
46
0.982822
0.283373
0.044888
47
1.021306
0.279519
0.044453
48
0.982809
0.283377
0.044883
49
1.067796
0.290610
0.045244
50
0.991362
0.277277
0.044169
51
1.067802
0.290616
0.045197
52
0.982822
0.283364
0.044914
53
1.021306
0.279545
0.044210
54
0.982822
0.283373
0.044905Table 15
3D homogeneous reactor power errors (%) corresponding to the second eigenvalue.
Node
Structured mesh
10 cm
5 cm
2 cm
1
0.996377
0.284899
0.045155
2
1.019201
0.279652
0.044237
3
0.996370
0.284906
0.045097
4
1.062091
0.289487
0.045069
5
1.008268
0.280020
0.044300
6
1.062091
0.289473
0.045105
7
0.996377
0.284871
0.045155
8
1.019187
0.279666
0.044165
9
0.996377
0.284856
0.045155
10
1.037403
0.284821
0.044523
11
1.008552
0.282182
0.044436
12
1.037403
0.284835
0.044473
13
1.034915
0.283264
0.044419
14
1.064620
0.288742
0.044959
15
1.034915
0.283235
0.044389
16
1.037410
0.284806
0.044658
17
1.008552
0.282277
0.044626
18
1.037403
0.284813
0.044729
19
1.031311
0.283454
0.044715
20
1.015537
0.283716
0.044796
21
1.031304
0.283440
0.044693
22
1.027935
0.281797
0.044396
23
1.072478
0.290229
0.045138
24
1.027935
0.281790
0.044459
25
1.031311
0.283412
0.044765
26
1.015558
0.283695
0.044655
27
1.031318
0.283412
0.044686
28
1.031311
0.283462
0.044715
29
1.015537
0.283716
0.044782
30
1.031304
0.283440
0.044686
31
1.027928
0.281804
0.044403
32
1.072478
0.290224
0.045144
33
1.027928
0.281797
0.044452
34
1.031311
0.283419
0.044772
35
1.015558
0.283695
0.044655
36
1.031311
0.283419
0.044679
37
1.037403
0.284813
0.044516
38
1.008552
0.282182
0.044442
39
1.037395
0.284828
0.044466
40
1.034915
0.283258
0.044419
41
1.064620
0.288751
0.044959
42
1.034915
0.283235
0.044389
43
1.037403
0.284806
0.044650
44
1.008552
0.282289
0.044626
45
1.037403
0.284799
0.044729
46
0.996377
0.284892
0.045155
47
1.019208
0.279659
0.044244
48
0.996384
0.284913
0.045097
49
1.062091
0.289473
0.045076
50
1.008250
0.280026
0.044300
51
1.062091
0.289458
0.045119
52
0.996377
0.284871
0.045162
53
1.019194
0.279659
0.044151
54
0.996377
0.284871
0.045155Figure 13
3D homogeneous reactor power corresponding to the first eigenvalue and the 2 cm structured mesh.Figure 14
3D homogeneous reactor power corresponding to the second eigenvalue and the 2 cm structured mesh.The analytical eigenvaluesk 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 k 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 Table16. 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.Table 16
Langenbuch reactor cross sections.
Material
D
1 (cm)
D
2 (cm)
Σ
a
1 (cm−1)
Σ
a
2 (cm−1)
Σ
12 (cm−1)
υ
Σ
f
1 (cm−1)
υ
Σ
f
2 (cm−1)
Comb.1
1.423913
0.356306
0.01040206
0.08766217
0.0175555
0.006477691
0.1127328
Comb.2
1.425611
0.350574
0.01099263
0.09925634
0.01717768
0.007503284
0.1378004
Absorbent
1.423913
0.356306
0.01095206
0.09146217
0.0175555
0.006477691
0.11273228
Reflector
1.634227
0.264002
0.002660573
0.04936351
0.02759693
0.0
0.0Figure 15
Axial plane of Langenbuch reactor.Figure 16
Frontal plane of Langenbuch reactor.In this case, only structured meshes have been used, in particular the 10 cm and 5 cm mesh length. Table17 shows the number of cells for each mesh.Table 17
Langenbuch reactor meshes.
Structured mesh
Mesh length (cm)
10
5
Number of cells
2800
18720Moreover, 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 (k
1 and k
2) are 0.994881227 and 0.948210698.Regarding the results, Table18 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.Table 18
Langenbuch reactor computational time (h:min:s).
Structured mesh
10 cm
5 cm
CELL+
0:0:25.69
0:10:54.5
CELL−
0:0:25.68
0:11:23.11
HOM
0:0:25.77
0:10:55.77
LIN
0:0:25.73
0:10:43.59Table 19
Langenbuch reactor eigenvalue-1 error (pcm).
Structured mesh
10 cm
5 cm
CELL+
83.3771
12.5047
CELL−
212.8518
63.0611
HOM
143.1450
22.6874
LIN
150.0886
25.2636Table 20
Langenbuch reactor eigenvalue-2 error (pcm).
Structured mesh
10 cm
5 cm
CELL+
226.1219
7.7211
CELL−
350.5162
80.2113
HOM
280.9989
39.9871
LIN
292.1948
44.2362Table 21
Numbering of the Langenbuch reactor nodes.
2nd axial plane
5th axial plane
—
—
—
—
—
—
—
—
—
—
21
22
23
24
—
—
45
46
47
48
—
—
16
17
18
19
20
—
40
41
42
43
44
—
11
12
13
14
15
—
35
36
37
38
39
—
6
7
8
9
10
—
30
31
32
33
34
—
1
2
3
4
5
—
25
26
27
28
29
—Table 22
Langenbuch reactor power errors (%) corresponding to the first eigenvalue and HOM.
Node
Structured mesh
10 cm
5 cm
1
0,88732759
0,083512937
2
0,413636006
0,169843064
3
0,285898499
0,26806574
4
1,650354128
0,463970268
5
7,308155289
2,443154619
6
0,413578042
0,169901321
7
0,218154684
0,211939733
8
0,89582742
0,279213792
9
2,440629291
0,498489677
10
8,245123598
2,368554084
11
0,285938524
0,268114276
12
0,895860783
0,279229236
13
1,485778349
0,264874283
14
3,250361067
0,503872312
15
8,254395078
1,996055827
16
1,650484996
0,46408882
17
2,440737883
0,498585195
18
3,250407371
0,5039109
19
6,540873894
1,269229237
20
10,49036132
3,000440824
21
7,308460522
2,44345405
22
8,245387401
2,368777266
23
8,25461267
1,99626084
24
10,49053455
3,000606939
25
5,466731896
0,931224354
26
4,140848872
0,716286587
27
3,657496267
0,698504374
28
2,237579532
0,474563288
29
3,189694202
1,268772356
30
4,140836529
0,716252343
31
2,96775873
0,633299552
32
2,476161123
0,618678432
33
0,910178355
0,364113213
34
4,548866137
1,200011218
35
3,657467753
0,698474984
36
2,476131091
0,618647924
37
2,967442105
0,56874817
38
0,292240082
0,412798339
39
4,267581802
0,724868524
40
2,237496931
0,474473118
41
0,910102518
0,364043283
42
0,292200002
0,412751132
43
2,930647399
0,281476818
44
6,41914469
1,657371938
45
3,189876024
1,268969435
46
4,549015577
1,2001458
47
4,267727516
0,724995924
48
6,419233002
1,657446439Table 23
Langenbuch reactor power errors (%) corresponding to the second eigenvalue and HOM.
Node
Structured mesh
10 cm
5 cm
1
2,748470704
0,472742419
2
1,064854909
0,505072706
3
0,379827902
0,520035451
4
0,846862035
0,527032189
5
6,400564495
2,241810405
6
1,064870054
0,505064685
7
0,500339771
0,504531289
8
1,13600946
0,496430539
9
2,54996541
0,525233749
10
8,338474277
2,124761105
11
0,379868741
0,519985959
12
1,135973249
0,496377465
13
1,719135726
0,413130915
14
3,303282255
0,426080951
15
8,233916013
1,633618147
16
0,846698163
0,526850063
17
2,549801498
0,525079406
18
3,303171711
0,426002352
19
6,461736033
0,930426868
20
10,39493074
2,481116214
21
6,400329869
2,24158495
22
8,338216136
2,124537285
23
8,233749447
1,633473241
24
10,39486316
2,481044325
25
8,712874235
0,839916928
26
6,327127
0,784501724
27
5,950753659
0,862124448
28
4,44225414
0,601601562
29
0,062727543
0,461900883
30
6,327225132
0,78459806
31
4,313496832
0,803439342
32
3,885076646
0,873376122
33
2,432393528
0,789928103
34
2,461696614
0,258876027
35
5,951087173
0,86245633
36
3,885366412
0,873674968
37
4,2237795
0,779967598
38
1,90610207
0,962348979
39
2,22742644
0,28534924
40
4,44301111
0,602394933
41
2,433102829
0,790642934
42
1,906524369
0,962753237
43
0,944598463
0,649099892
44
4,201989073
0,422463819
45
0,061633966
0,460821104
46
2,460613715
0,257845259
47
2,226705013
0,286064257
48
4,20167533
0,422178482Figure 17
Langenbuch reactor power corresponding to the first eigenvalue and the 5 cm structured mesh and HOM.Figure 18
Langenbuch reactor power corresponding to the second eigenvalue and the 5 cm structured mesh and HOM.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.



---

*Source: 913043-2014-01-23.xml*