Abstract and Applied Analysis

Volume 2014 (2014), Article ID 913043, 15 pages

http://dx.doi.org/10.1155/2014/913043

Research Article

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

Álvaro Bernal,¹ Rafael Miró,¹ Damián Ginestar,² and Gumersindo Verdú¹

¹Institute for Industrial, Radiophysical and Environmental Safety (ISIRYM), Universitat Politècnica de València, 46022 Valencia, Spain

²Institute for Multidisciplinary Mathematics, Universitat Politècnica de València, 46022 Valencia, Spain

Received 31 October 2013; Accepted 22 November 2013; Published 23 January 2014

Academic Editor: Benito Chen-Charpentier

Copyright © 2014 Álvaro Bernal et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

Figure 1: Virtual fluxes example.

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

Figure 2: Adjacent cells.

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:

Figure 3: Partial and global gradients.

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.

Table 1: 2D homogeneous reactor cross sections.

(cm)
(cm)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)

1.28205128205
0.666667
0.01
0.1
0.01
0.109017634020268
0.075

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.

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

Figure 4: 2D structured mesh.

Figure 5: 2D unstructured mesh.

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.

Table 3: Numbering of the 2D homogeneous reactor nodes.

7
8
9

4
5
6

1
2
3

Table 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.007689

Table 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.029325

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

Table 6: Biblis reactor cross sections.

Material
(cm)
(cm)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)

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

Figure 10: Biblis reactor: geometry and material composition.

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.

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
23972

Table 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.03

Table 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.258402

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.

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

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

Figure 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: . 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.

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

Table 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
48

Table 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.044905

Table 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.045155

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

Table 16: Langenbuch reactor cross sections.

Material
(cm)
(cm)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)
(cm⁻¹)

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

Figure 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. Table 17 shows the number of cells for each mesh.

Table 17: Langenbuch reactor meshes.

Structured mesh

Mesh length (cm)
10
5

Number of cells
2800
18720

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.

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

Table 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.2636

Table 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.2362

Table 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,657446439

Table 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,422178482

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

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.

References

W. M. Stacey, Nuclear Reactor Physics, John Wiley & Sons, New York, NY, USA, 2001.

G. Verdú, D. Ginestar, V. Vidal, and J. L. Muñoz-Cobo, “3D λ-modes of the neutron-diffusion equation,” Annals of Nuclear Energy, vol. 21, no. 7, pp. 405–421, 1994.

R. Miró, D. Ginestar, G. Verdú, and D. Hennig, “A nodal modal method for the neutron diffusion equation. Application to BWR instabilities analysis,” Annals of Nuclear Energy, vol. 29, no. 10, pp. 1171–1194, 2002.

V. Hernandez, J. E. Roman, and V. Vidal, “SLEPc: a scalable and flexible toolkit for the solution of eigenvalue problems,” ACM Transactions on Mathematical Software, vol. 31, no. 3, pp. 351–362, 2005.

K. A. Hoffmann and S. T. Chiang, Computational Fluid Dynamics, vol. 2, Engineering Education System, Wichita, Kan, USA, 4th edition, 2000.

D. J. E. Harvie, “An implicit finite volume method for arbitrary transport equations,” ANZIAM Journal, vol. 52, pp. C1126–C1145, 2012.

L. Cueto-Felgueroso, I. Colominas, X. Nogueira, F. Navarrina, and M. Casteleiro, “Finite volume solvers and moving least-squares approximations for the compressible Navier-Stokes equations on unstructured grids,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 45–48, pp. 4712–4736, 2007.

C. Geuzaine and J.-F. Remacle, “Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009.

V. Hernández, J. E. Román, and V. Vidal, “SLEPc: scalable library for eigenvalue problem computations,” in High Performance Computing for Computational Science—VECPAR 2002, vol. 2565 of Lecture Notes in Computer Science, pp. 377–391, Springer, Berlin, Germany, 2003.

E. Z. Müller and Z. J. Weiss, “Benchmarking with the multigroup diffusion high-order response matrix method,” Annals of Nuclear Energy, vol. 18, no. 9, pp. 535–544, 1991.