3D neutron diffusion computational code based on GFEM with unstructured tetrahedron elements: A comparative study for linear and quadratic approximations

Progress in Nuclear Energy 92 (2016) 119e132

Contents lists available at ScienceDirect

Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene

3D neutron diffusion computational code based on GFEM with unstructured tetrahedron elements: A comparative study for linear and quadratic approximations Seyed Abolfazl Hosseini Department of Energy Engineering, Sharif University of Technology, Tehran, 8639-11365, Iran

a r t i c l e i n f o

a b s t r a c t

Article history: Received 18 January 2016 Received in revised form 12 June 2016 Accepted 8 July 2016 Available online 28 July 2016

In the present study, the comparison between the results obtained from the linear and quadratic approximations of the Galerkin Finite Element Method (GFEM) for neutronic reactor core calculation was reported. The sensitivity analysis of the calculated neutron multiplication factor, neutron flux and power distributions in the reactor core vs. the number of the unstructured tetrahedron elements and order of the considered shape function was performed. The cost of the performed calculation using linear and quadratic approximation was compared through the calculation of the FOM. The neutronic core calculation was performed for both rectangular and hexagonal geometries. Both the criticality and fixed source calculations were done using the developed GFEM-3D computational code. An acceptable accuracy with low computational cost is the main advantage of applying the unstructured tetrahedron elements. The generated unstructured tetrahedron elements with Gambit software were used in the GFEM3D computational code via a developed interface. The criticality calculation was benchmarked against the valid data for IAEA-3D and VVER-1000 benchmark problems. Also, the neutron fixed source calculation was validated through the comparison with the similar computational code. The results show that the accuracy of the calculation for the both linear and quadratic approximations improves vs. the number of elements. Quadratic approximation gives acceptable results for almost all considered number of the elements, while the results obtained from the linear approximation have good accuracy for only high number of the elements. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Galerkin Finite Element Method Unstructured tetrahedron elements Neutron flux Power distribution Linear Quadratic

1. Introduction The finite element method has always been a fundamental numerical technique in reactor core calculations. It has been continuously improved over decades, starting from primal implementation in neutron diffusion equation, up to modern implementations with Raviart- Thomas, Hybrid, h-adaptivity, Response
bert, 2008), (Cavdar and Ozgener, 2004) and Matrix bases (He

Abbreviations: FEM, Finite Element Method; GFEM, Galerkin Finite Element Method; FA, Fuel assembly; fg(r), Forward flux in the energy group g; keff, Forward neutron multiplication factor; cg, Neutron spectrum in energy group g; Dg, Diffusion constant in the energy group g; St,g, Macroscopic total cross section in the energy group g; Sf,g, Macroscopic fission cross section in the energy group g; Ss;g’ /g , Macroscopic scattering cross section from energy group g’ to g; n, Fission neutron yield; V, The nabla operator. E-mail address: [email protected]. http://dx.doi.org/10.1016/j.pnucene.2016.07.006 0149-1970/© 2016 Elsevier Ltd. All rights reserved.

(Wang et al., 2009). Other numerical methods like Nodal (Zimin and Baturin, 2002) and Finite Difference (Ivanov et al., 1994) and may also be used to solve the neutron diffusion equation. In general, the Finite Element Methods (FEMs) is preferred in most applications to its principal alternative, the Finite Difference Method (FDM), due to its flexibility in the treatment of curved or irregular geometries and the high rates of convergence attainable by the use of high order elements. Several researchers have tried to develop convenient methods for solving multi-group neutron diffusion equations in the three dimensional geometries. Some selected researches are summarized as follows: 1. Ivanov and his collogues (Ivanov et al., 1994) developed an efficient algorithm for a more accurate solution of the multigroup coarse-mesh diffusion equations in hexagonal-z

120

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

geometry. The finite difference type nodal approach to solve the diffusion equations applied in neutronics design code HEXAB3D. 2. Cho and Lee (Cho and Lee, 2006) proposed the Analytic Function Expansion Nodal (AFEN) method to solve the multi-group neutron diffusion equation in hexagonal-z geometry. The COREDAX code implementing the AFEN method was verified testing on the VVER-440 benchmark problem. 3. Wang and his collogues (Wang et al., 2009) presented the three dimensional h-adaptivity for the multi-group neutron diffusion equations. Adaptive mesh refinement (AMR) showed to allow solving partial differential equations to significantly higher accuracy at reduced numerical cost. A novel algorithm for assembling the terms coupling shape functions from different meshes was proposed. In the present study, the Galerkin FEM (GFEM) (Hosseini and Vosoughi, 2013), a weighted residual method, is used to solve the multi-group neutron diffusion equation for rectangular and hexagonal three dimensional reactor cores. Unstructured tetrahedron elements generated by Gambit software (Feng et al., 2005) are used to discrete the equations. Indeed, a key advantage of the unstructured tetrahedron elements is their superiority in mapping the curved boundaries or material interfaces in three dimensional geometries. Also, the running time of computational code and the accuracy of the calculation may be optimized using the proper unstructured tetrahedron elements. For several reasons, such as desired precision yet being simple, the Galerkin method has been widely used in development of computational codes for solving the diffusion or transport equation in different geometries (Hosseini and Vosoughi, 2013) and (Maiani and Montagnini, 2004). The main advantage of GFEM is that the definition of boundary conditions in this method is easier than the other methods (Zhu et al., 2005). The mentioned reasons convinced us to use GFEM for solving multi-group forward/adjoint neutron diffusion equation in three dimensional geometries. An outline of the remainder of this contribution is as follows: In Section 2, we briefly introduce the mathematical formulation used to solve the3D neutron diffusion equation using GFEM by considering the linear and quadratic approximation of shape function. Section 3 presents the main specification of the IAEA-3D (Center, 1977) and VVER-1000 (Schulz, 1996) benchmark problems. Numerical results and a sensitivity analysis of the calculations to the number of elements are presented in Section 4. Also, in this section, we discuss on the results and advantages of applying the unstructured tetrahedron elements. Section 5 gives a summary and concludes the paper.

2.1.1. Linear approximation for shape function Eq. (1) is a linear partial differential equation which may be solved by different numerical methods. All of these methods transform the differential equation into a system of algebraic equations. Here, GFEM, a weighted residual method, is used to discrete the neutron diffusion equation. To start the discretization, the whole solution volume is divided into the unstructured tetrahedron elements as shown in Fig. 1. These elements have been generated using Gambit mesh generator. In the linear approximation of shape function, the neutron flux in each element might be considered as Eq. (2) (Zhu et al., 2005):

fðeÞ ðx; y; zÞ ¼ L1 ðx; y; zÞf1 þ L2 ðx; y; zÞf2 þ L3 ðx; y; zÞf3 þ L4 ðx; y; zÞf4 ;

(2)

where, Li,i ¼ 1,2,3,4. are the components of the shape function in Eq. (3): ðeÞ

N ðx; y; zÞ ¼ ½ L1 ðx; y; zÞ

L2 ðx; y; zÞ



L3 ðx; y; zÞ

L4 ðx; y; zÞ 

(3)

The shape function components are defined as Eq. (4):

Li ðx; y; zÞ ¼

ai þ bi x þ ci y þ di z ; i ¼ 1; 2; 3; 4: 6V

(4)

with

2

1 61 6 6V ¼ det4 1 1

x1 x2 x3 x4

y1 y2 y3 y4

3 z1 z2 7 7 z3 5 z4

(5)

in which, incidentally, the value V represents the volume of the tetrahedron. By expanding the other relevant determinants into their cofactors we have:

2

3 x2 y2 z2 a1 ¼ det4 x3 y3 z3 5 2x4 y4 z4 3 x2 1 z 2 c1 ¼ det4 x3 1 z3 5 x4 1 z 4

2

3 1 y2 z2 b1 ¼ det4 1 y3 z3 5 2 1 y4 z4 3 x2 y 2 1 d1 ¼ det4 x3 y3 1 5 x4 y 4 1

(6)

with the other constants defined by cyclic interchange of the subscripts in the order 1, 2, 3, 4. The components of the shape function satisfy the criterion given in Eq. (7) at all points of the domain:

2. Mathematical formulation 2.1. Fixed source calculation In the fixed source problem, the distribution of the neutron flux in the environment is obtained for the considered fixed neutron source. To this end, the neutron diffusion equation for fixed neutron source problem should be solved. The general matrix form of the fixed neutron source equation is as Eq. (1):

Ag fg ¼ Sg

(1)

P P where, A ¼ DgV2 þ St,g and Sg ¼ Sext;g þ G s;g’/g 4g’ g’¼1 PG g ðrÞ þ cg g’¼1 ng’ Sf ;g’ 4g’ ðrÞ are the loss operator and neutron source in energy group g, respectively. All quantities in Eq. (1) are defined as usual.

Fig. 1. The unstructured tetrahedron element by considering the linear approximation of the shape functions in each element (4 points in the vertices).

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

L1 ðx; y; zÞ þ L2 ðx; y; zÞ þ L3 ðx; y; zÞ þ L4 ðx; y; zÞ ¼ 1

(7)

The GFEM is a weighted residual method in which the purpose is to minimize the residual integral. In the weighted residual methods, the weighting function is considered as WðrÞ ¼ W T N ðrÞ. There are (at least) four sub-methods (Collocation, Sub-domain, Least Squares, Galerkin) for different functions of WT. Since WT is unit in the Galerkin method, the weighting function is considered as Eq. (8):

WðrÞ ¼ N ðrÞ

(8)

where, N ðrÞ is the global shape function. Multiplying Eq. (1) by the weighting function and integrating the results over the solution space, Eq. (9) is obtained:

Z

 dvWðrÞ  Dg V2 4g ðrÞ þ Sr;g 4g ðrÞ  Sext;g ðrÞ V X  Ss;g’ /g 4g’ ðrÞÞ ¼ 0

(9)

g ’ sg

In the above equation, the differential part may be transformed by applying the Divergence’s theorem:

Z V

  Z dvWðrÞ  Dg V2 fg ðrÞ ¼

When element matrices have to be evaluated it will follow that we are faced with integration of quantities defined in terms of volume coordinates over the tetrahedron region. The exact solution of the integral is given in Appendix A, see equation (A.1). The appendix contains mathematical details of the GFEM method. We have encountered with three types of integrals in solving the Eq. (14). The solution of the appeared integrals in the calculation of 3D neutron diffusion equation using linear approximation of the Galerkin finite element methods is given in Eqs. (A.2)-(A.3).

2.1.2. Quadratic approximation for shape function In the next step of the calculation, the quadratic approximation for shape function in each element is considered. As shown in Figs. 2 and 10 nodes including corners and mid-edge of tetrahedron are considered in the quadratic approximation. The shape function components for corner nodes are as Eq. (15):

Na ¼ ð2La  1ÞLa

V

a ¼ 1; 2; 3; 4

Z V

dvV$ðWðrÞVfg ðrÞ ¼

where,

vfg ðrÞ ¼ Vfg ðrÞ$n: vn

(11)

Here, n is the normal unit vector on the volumeV. Two types of boundary conditions (B.C.) are considered. The first B.C. is no incoming neutrons at vacuum boundaries (Marshak B.C.) which is expressed as Eq. (12):

vfg ðrÞ fg ðrÞ : ¼ vn 2Dg

(15)

where, La has been defined in Eq. (4). For mid-edge nodes, the shape function components are as Eq. (16):

Z dvVWðrÞ$Vfg ðrÞ 

121

Z V

dvVWðrÞ$Vfg ðrÞ 

N5 ¼ 4L1 L2 N6 ¼ 4L1 L3 N7 ¼ 4L1 L4 N8 ¼ 4L2 L3 N9 ¼ 4L3 L4 N10 ¼ 4L2 L4

dsWðrÞ A

vfg ðrÞ vn

(10)

(16)

Beginning from Eq. (9) and performing the same steps applied for linear approximation, one obtains the solution of integrals as Eq.

(12)

The second B.C. is zero net current or perfect reflective boundary condition which is described by Eq. (13):

vfg ðrÞ ¼ 0: vn

(13)

Substituting Eqs. (8), (10), (12) and (13), and converting the integration on the reactor domain to sum of the integration on finite elements, the final form of Eq. (9) is obtained as Eq. (14):

"Z E X e¼1

dvDg VN 

¼

ðrÞVN

TðeÞ



ðrÞfg

V

Z þ

ðeÞ

dsN A E h X e¼1

ðeÞ



ðeÞ

ðrÞN

TðeÞ



Z

Sext;g

dvN 

V

ðeÞ

ðeÞ

ðrÞ

ðeÞ

ðeÞ

Z

þ Sr;g

dvN 

ðeÞ

ðrÞN

TðeÞ



ðeÞ

ðrÞfg

V

3

fg 7 5 2

ðrÞ þ

g1 X g’ ¼1

Z Sg’ /g

dvN 

ðeÞ

ðrÞN 

TðeÞ

ðeÞ

ðrÞfg’

V

(14)

Fig. 2. The unstructured tetrahedron element by considering the quadratic approximation of the shape functions in each element (4 points in the vertices and 6 points in the middle of each edge).

122

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

(A.5) and (A.6) in the quadratic approximation for shape function. Finally, the solution of local appeared integral due to presence of volumetric external source for linear and quadratic approximations of shape functions are as Eqs. (A.7) and (A.8), respectively. The system of equations which is inverse problem may be obtained from assembling the calculated local matrixes for each element. The solution of the system of equations gives the neutron flux distribution.

2.2. Criticality calculation 2.2.1. Linear approximation for shape function In the absence of external neutron source, the multi-group neutron diffusion equation is as Eq. (17) (Duderstadt and Hamilton, 1976) and (Lamarsh, 1966). The lambda-mode is one of the possible calculations if the system has no fixed source and is at steady state.

Dg V2 fg ðrÞ þ St;g fg ðrÞ ¼

G cg X

keff þ

g ’ ¼1

nSf ;g’ fg’ ðrÞ

G X X

fg’ ðrÞ

; g ¼ 1; 2; :::; G

g’¼1 s;g’/g

(17) In operator presentation (Agfg ¼ Sg), Ag ¼ DgV þ St,g and PG P c PG Sg ¼ k g s;g’/g fg ’ ðrÞ are the loss and g’ ¼1 nSf ;g’ fg’ ðrÞ þ g’¼1 eff neutron source terms, respectively. In the criticality problem, three types of integral are appeared. The solutions of the appeared integrals in each element are as Eqs. (A.2)e(A.6). The system of equations which is an eigen-value problem is obtained from the assembling local matrixes into the global matrix. Here, the problem is solved using the power iteration method described as follow (Booth, 2006): 2

(1) Initialize the multi-group fluxes and the eigen-value with ð0Þ

fg ,

ð0Þ

keff ; compute the initial fission source PG ð0Þ ð0Þ Sf ðrÞ ¼ g¼1 nSf ;g ðrÞfg ðrÞ and set the iteration index to

n ¼ 1. ðnÞ (2) Solve for the multi-group fluxes fg using: ðnÞ

ðnÞ

Dg V2 fg ðrÞ þ Sr;g fg ðrÞ ¼

X

cg ðn1Þ ðn1Þ S þ Ss;g’ /g fg’ ðrÞ ðn1Þ f keff ’ g sg

3. Main specification of the benchmark problems 3.1. Neutron fixed source benchmark problem The GFEM-3D computational code may be applied to solve the neutron fixed source equation in anythree dimensional geometries created by Gambit. To validate the neutron fixed source calculation, benchmark problem, the results obtained from the GFEM-3D and CITATION (Fowler et al., 1971) are compared. Fig. 3 shows the multiregion cube with the mentioned materials and unit volumetric neutron source located in regions1. The dimension of considered cube is20 cm  20 cm  20 cm. The material cross section of each assembly for the considered cube is given in Table 1. The boundary condition in the all surfaces is the no incoming neutron current. 3.2. Criticality benchmark problems 3.2.1. IAEA-3D PWR The IAEA-3D PWR problem has been very important standard benchmark problem to measure the performance of neutronics calculation methods (Center, 1977). 177 fuel assemblies including 9 fully rodded fuel assemblies and 4 partially rodded fuel assemblies compose the core. 64 reflector assemblies surround the core. The fuel assembly pitch is 20 cm and the active height of a fuel assembly is 340 cm. The thickness of axial reflector is 20 cm. Fig. 4 displays the 1/8 of the IAEA-3D PWR. The boundary conditions of the reactor core comprise of no incoming current for the external boundaries and perfect reflective for the symmetry line boundaries. Table 2 represents the material cross section of each assembly for IAEA3D reactor core. 3.2.2. VVER-1000 The VVER-1000 is the second benchmark problem in this study (Schulz, 1996). The radial fuel assembly lattice pitch is 24.1 cm. This corresponds to the prototype VVER-1000 and is slightly different from the actual FA pitch of 23.6 cm in VVER-1000/V320, however it is acceptable for a mathematical benchmark. The core height is 355 cm, covered with axial and radial reflectors. The total height is 426 cm including 35.5 cm thick axial reflectors. Fig. 5 shows the 1/ 12 of the benchmark core configuration. The boundary conditions of the reactor core comprise of no incoming current for the external boundaries and perfect reflective for the symmetry line boundaries. The material cross section of each assembly for VVER-1000 is given in Table 3.

g

¼ 1; 2; :::; G: (18) P ðnÞ ðnÞ (3) Update the fission integral as: Sf ðrÞ ¼ G g¼1 nSf ;g ðrÞfg ðrÞ. (4) Update the eigen-value as: R ðnÞ ðn1Þ R ðnÞ ðn1Þ dU S ðrÞ= dU S ðrÞ. keff ¼ keff U U f f ðnÞ

ðn1Þ

(5) Compare keff with keff

ðnÞ

and fg

ðn1Þ

with fg

for all energy

groups. If the changes are greater than a prescribed tolerance, then set n ¼ n þ 1 and repeat the iteration stating at step-(2); otherwise the iteration is over. Applying the above mentioned algorithm to solve the system of equations results to the calculation of the neutron flux distribution in each energy group and the neutron multiplication factor.

Fig. 3. The considered cube for validation of neutron fixed source problem.

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132 Table 1 The material cross section of each assembly for considered cube in neutron fixed source problem. Region

D1(cm)

D2(cm)

Sa,1(cm1)

Sa,2(cm1)

Ss,1/2(cm1)

1 2 3 4 5 6 7 8 9

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

0.026 0.027 0.028 0.029 0.030 0.031 0.032 0.033 0.034

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.017245 0.017245 0.017245 0.017245 0.017245 0.017245 0.017245 0.017245 0.017245

123

comparisons between the calculated average neutron flux in each cube using GFEM-3D and CITATION computational codes vs. the number of elements are given in Tables 4e9. Tables 4e6 show the comparisons of the fast neutron flux distribution vs. the number of elements in each cube. Also, the comparisons of the thermal neutron flux distribution vs. the number of elements are given in Tables 7e9. As shown, the quadratic approximation gives the more accurate results in comparison with the linear approximation. The maximum error for both linear and quadratic approximations has been displayed as red color. Figs. 6e9 show the variation of the Relative Percent Error (RPE) of the calculated fast and thermal neutron flux distributions for both linear and quadratic approxi-

Fig. 4. The 1/8 of the IAEA-3D PWR reactor core.

Table 2 The material cross section of each section of the IAEA-3D reactor core. (1)

Fuel 1 D1(cm) D2(cm) nSf,1(cm1) nSf,2(cm1) Sa,1(cm1) Sa,2(cm1) S(cm1)s,1 / 2

1.500 0.400 0.000 0.135 0.010 0.080 0.020

(2)

(3)

Fuel 1 þ rod

Fuel 2

1.500 0.400 0.000 0.135 0.010 0.085 0.020

1.500 0.400 0.000 0.135 0.010 0.130 0.020

(4)

Reflector 2.000 0.300 0.000 0.000 0.000 0.010 0.040

(5)

mations. As shown from the mentioned tables and figures, the differences between the results obtained from the GFEM-3D and CITATION decreases vs. the number of elements. The presented quantity “RPE” in the tables is defined as Eq. (19):

Reflector þ rod 2.000 0.300 0.000 0.000 0.000 0.055 0.040

4. Numerical results and discussion First, we present the results obtained from the fixed neutron source calculation. To validate the performed calculation for neutron fixed source problem, the results of GFEM-3D and CITATION (Fowler et al., 1971) computational codes have been compared. Since the calculation was performed using small meshes in CITATION computational code, the obtained neutron flux distribution from CITATION may be considered as reference data. The

RPEð%Þ ¼

calculated value  refernce value  100: refernce value

(19)

To compare the cost of the performed calculation using CITATION, GFEM-3D (Linear) and GFEM-3D (Quadratic), the FOM has been calculated using Eq. (20):

FOM ¼ W1 T þ W2 M

(20)

where T is calculation time, M is memory footprint and W1 and W2 are weights. For example, W1 and W2 may be cost of one processor and cost of one memory, respectively. The considered values for W1andW2 in the present study are 2000 and 0.002, respectively. Fig. 10 shows the comparison between the calculated FOM for calculation performed using CITATION, GFEM-3D (Linear) and GFEM-3D (Quadratic) vs. accuracy of the calculation. As shown, FOM in quadratic approximation is less than two other ones. The lower FOM is considered to be better in the calculation.

124

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

Fig. 5. The 1/12 of the VVER-1000 reactor core (R: Reflector; 7/4/7: top and bottom section:material No. 7, middle section consists of material No. 4).

Table 3 The material cross section of each section of the VVER-1000 reactor core.

D1(cm) D2(cm) nSf,1(cm1) nSf,2(cm1) Sa,1(cm1) Sa,2(cm1) S(cm1)s,1 / 2

(1)

(2)

(3)

(4)

(5)

(6)

(7)

1.375 0.383 0.005 0.084 0.008 0.066 0.016

1.409 0.388 0.005 0.084 0.010 0.075 0.014

1.371 0.380 0.006 0.115 0.009 0.080 0.015

1.394 0.385 0.006 0.126 0.010 0.095 0.014

1.369 0.379 0.006 0.130 0.009 0.088 0.015

1.000 0.333 0.000 0.000 0.016 0.053 0.025

1.369 0.379 0.006 0.126 0.009 0.086 0.015

Here, we present the results obtained from the criticality calculation. Fig. 11 shows the IAEA-3D with unstructured tetrahedron elements. Table 10 shows the calculated neutron multiplication factors vs. number of the unstructured tetrahedron elements in the linear and quadratic approximation for shape function. In Table 11, power distribution in the reactor core has been compared with the reference data (Center, 1977). The data used as reference in the present study were obtained by extrapolation of results of the VENTURE computational code. Extrapolation of results is done on the basis of error dependence on the square of the mesh spacing

Table 4 Comparison between the calculated fast neutron flux distribution using CITATION and GFEM-3D computational codes for 5778 tetrahedron elements in the fixed source benchmark problem. Region

CITATION (#/cm2 s)

GFEM-3D (#/cm2 s) Linear

RPE (%)

Quadratic

RPE (%)

1 2 3 4 5 6 7 8 9

9.016Eþ00 1.053Eþ00 6.800E-03 1.048Eþ00 2.042E-01 3.600E-03 8.200E-03 3.900E-03 3.000E-04

8.586Eþ00 1.131Eþ00 6.100E-03 1.116Eþ00 2.027E-01 3.000E-03 7.200E-03 3.300E-03 2.000E-04

4.776Eþ00 7.426Eþ00 1.029Eþ01 6.528Eþ00 7.346E-01 1.667Eþ01 1.220Eþ01 1.538Eþ01

8.991Eþ00 1.058Eþ00 6.800E-03 1.052Eþ00 2.037E-01 3.600E-03 8.200E-03 4.000E-03 3.000E-04

2.751E-01 4.558E-01 0.000Eþ00 4.008E-01 2.449E-01 0.000Eþ00 0.000Eþ00

-3.333E+01

2.564E+00

0.000Eþ00

Table 5 Comparison between the calculated fast neutron flux distribution using CITATION and GFEM-3D computational codes for 9074 tetrahedron elements in the fixed source benchmark problem. Region

1 2 3 4 5 6 7 8 9

CITATION (#/cm2 s)

9.016Eþ00 1.053Eþ00 6.800E-03 1.048Eþ00 2.042E-01 3.600E-03 8.200E-03 3.900E-03 3.000E-04

GFEM-3D (#/cm2 s) Linear

RPE (%)

Quadratic

RPE (%)

8.701Eþ00 1.102Eþ00 6.300E-03 1.099Eþ00 2.053E-01 3.200E-03 7.700E-03 3.600E-03 3.000E-04

3.494Eþ00 4.653Eþ00 7.353Eþ00 4.839Eþ00 5.387E-01 1.111Eþ01 6.098Eþ00

8.997Eþ00 1.056Eþ00 6.800E-03 1.051Eþ00 2.040E-01 3.600E-03 8.200E-03 3.900E-03 3.000E-04

2.141E-01

-7.692E+00

0.000Eþ00

3.039E-01

0.000Eþ00 2.768E-01 9.794E-02 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

125

Table 6 Comparison between the calculated fast neutron flux distribution using CITATION and GFEM-3D computational codes for 69,246 tetrahedron elements in the fixed source benchmark problem. Region

1 2 3 4 5 6 7 8 9

CITATION (#/cm2 s)

9.016Eþ00 1.053Eþ00 6.800E-03 1.048Eþ00 2.042E-01 3.600E-03 8.200E-03 3.900E-03 3.000E-04

GFEM-3D (#/cm2 s) Linear

RPE (%)

Quadratic

RPE (%)

8.929Eþ00 1.067Eþ00 6.600E-03 1.061Eþ00 2.044E-01 3.500E-03 8.000E-03 3.900E-03 3.000E-04

9.660E-01 1.291Eþ00

9.003Eþ00 1.055Eþ00 6.700E-03 1.049Eþ00 2.044E-01 3.600E-03 8.100E-03 3.900E-03 3.000E-04

1.475E-01 1.519E-01

-2.941E+00

1.231Eþ00 9.794E-02 2.778Eþ00 2.439Eþ00 0.000Eþ00 0.000Eþ00

-1.471E+00

1.432E-01 9.794E-02 0.000Eþ00 1.220Eþ00 0.000Eþ00 0.000Eþ00

Table 7 Comparison between the calculated thermal neutron flux distribution using CITATION and GFEM-3D computational codes for 5778 tetrahedron elements in the fixed source benchmark problem. Region

CITATION (#/cm2 s)

1 2 3 4 5 6 7 8 9

1.154Eþ00 1.014E-01 4.000E-04 5.760E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

GFEM-3D (#/cm2 s) Linear

RPE (%)

Quadratic

RPE (%)

1.055Eþ00 1.090E-01 4.000E-04 6.260E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

8.645Eþ00 7.495Eþ00 0.000Eþ00

1.148Eþ00 1.018E-01 4.000E-04 5.800E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

5.371E-01 3.945E-01 0.000Eþ00

8.681E+00

0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00

6.944E-01

0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00

Table 8 Comparison between the calculated thermal neutron flux distribution using CITATION and GFEM-3D computational codes for 9074 tetrahedron elements in the fixed source benchmark problem. Region

CITATION (#/cm2 s)

1 2 3 4 5 6 7 8 9

1.154Eþ00 1.014E-01 4.000E-04 5.760E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

GFEM-3D (#/cm2 s) Linear

RPE (%)

Quadratic

RPE (%)

1.081Eþ00 1.059E-01 4.000E-04 6.140E-02 8.800E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

6.367Eþ00 4.438Eþ00 0.000Eþ00

1.150Eþ00 1.017E-01 4.000E-04 5.790E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

4.158E-01 2.959E-01 0.000Eþ00

6.597E+00

1.149Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00

5.208E-01

0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00

Table 9 Comparison between the calculated thermal neutron flux distribution using CITATION and GFEM-3D computational codes for 69246 tetrahedron elements in the fixed source benchmark problem. Region

1 2 3 4 5 6 7 8 9

CITATION (#/cm2 s)

1.154Eþ00 1.014E-01 4.000E-04 5.760E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

GFEM-3D (#/cm2 s) Linear

RPE (%)

1.134Eþ00 1.026E-01 4.000E-04 5.860E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

1.183Eþ00 0.000Eþ00 1.736Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00

(Center, 1977). Table 12 shows the calculated neutron multiplication factor from both linear and quadratic approximations. The comparisons between the averaged relative power calculated by GFEM-3D

-1.741E+00

Quadratic

RPE (%)

1.152Eþ00 1.015E-01 4.000E-04 5.770E-02 8.700E-03 1.000E-04 2.000E-04 1.200E-04 1.000E-05

9.862E-02 0.000Eþ00 1.736E-01 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00

-2.512E-01

computational code and reference value (Schulz, 1996) are given in Figs. 12 and 13. The reference values for the Schulz benchmark were calculated using the CRONOS computational code (Schulz, 1996) and with the extrapolated finite element solution of second

126

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

Fig. 6. The variation of the RPE (%) of calculated fast neutron flux in linear approximation of the shape function.

Fig. 7. The variation of the RPE (%) of calculated fast neutron flux in quadratic approximation of the shape function.

order with Lagrange polynomials on triangular-z meshes. CRONOS is a reactor code of CEA, which uses finite Q14 elements and nodal methods for homogenized diffusion and transport calculations; the code also has 3D kinetics and pin by-pin diffusion modules. The approximation in the hexagonal plane uses the Gauss-Legendre numerical quadrature corresponding to super convergent finite elements. The discretization was performed by considering 54 triangles per hexagon and Nz ¼ 24 meshes in the axial direction. The accuracy of the obtained results improves with increasing the

number of elements. The linear approximation gives the results with the accuracy of maximum RPEs (%) equal to 1.93420E-01 and 2.143Eþ00 for the neutron multiplication factor and power distribution, respectively. The same maximum RPEs (%) in quadratic approximation are 3.81000E-03 and 8.116E-01, respectively. As shown, the difference between the calculated neutron multiplication factor from linear and quadratic approximation, and reference value decreases as the number of elements is increased. For the quadratic approximation, the variation of the calculated

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

127

Fig. 8. The variation of the RPE (%) of calculated thermal neutron flux in linear approximation of the shape function.

Fig. 9. The variation of the RPE (%) of calculated thermal neutron flux in quadratic approximation of the shape function.

neutron multiplication factor vs. the number of the elements is slow. The results obtained from the linear approximation in the present study are in the range of other same reported results


lez-Pintor et al., 2009). (Maiani and Montagnini, 2004) and (Gonza Also as expected, the results obtained from the quadratic approximation are more accurate in comparison with the linear

128

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

Fig. 10. A view of IAEA-3D reactor core that has been gridded with unstructured tetrahedron elements.

Fig. 11. FOM for the performed calculation using CITATION, GFEM-3D (Linear) and GFEM-3D (Quadratic).

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

129

Table 10 The calculated neutron multiplication factor for IAEA-3D reactor core. Number of elements

Linear

RPE (%)

Quadratic

RPE (%)

36,629 41,691 43,262

1.03093Eþ00 1.03061Eþ00 1.03047Eþ00

1.84640E-01 1.53540E-01 1.39940E-01

1.02912Eþ00 1.02911Eþ00 1.02911Eþ00

8.75000E-03 7.77000E-03 7.77000E-03

*The reference effective multiplication factor is keff ¼ 1.02903Eþ00 (Center, 1977).

Table 11 The comparison of the calculated relative power distribution in the IAEA-3D reactor core by considering the 51,017 elements. FA

Ref. (Center, 1977)

Linear

RPE (%) in linear

Quadratic

RPE (%) in quadratic

1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 18 19 20 21 22 23 25 26 27 28 32 33 35

7.290E-01 1.281Eþ00 1.422Eþ00 1.193Eþ00 6.100E-01 9.530E-01 9.590E-01 7.770E-01 1.397Eþ00 1.432Eþ00 1.291Eþ00 1.072Eþ00 1.055Eþ00 9.760E-01 7.570E-01 1.368Eþ00 1.311Eþ00 1.181Eþ00 1.089Eþ00 1.000Eþ00 7.110E-01 1.178Eþ00 9.720E-01 9.230E-01 8.660E-01 7.000E-01 6.110E-01 5.970E-01

7.170E-01 1.239Eþ00 1.383Eþ00 1.157Eþ00 6.250E-01 9.330E-01 9.770E-01 8.210E-01 1.352Eþ00 1.386Eþ00 1.253Eþ00 1.038Eþ00 1.030Eþ00 9.920E-01 7.930E-01 1.332Eþ00 1.274Eþ00 1.159Eþ00 1.118Eþ00 1.025Eþ00 7.530E-01 1.157Eþ00 9.540E-01 9.360E-01 8.960E-01 7.160E-01 6.450E-01 6.300E-01

1.646Eþ00 3.279Eþ00 2.743Eþ00 3.018Eþ00 2.459Eþ00 2.099Eþ00 1.877Eþ00 5.663Eþ00 3.221Eþ00 3.212Eþ00 2.943Eþ00 3.172Eþ00 2.370Eþ00 1.639Eþ00 4.756Eþ00 2.632Eþ00 2.822Eþ00 1.863Eþ00 2.663Eþ00 2.500Eþ00 5.907Eþ00 1.783Eþ00 1.852Eþ00 1.408Eþ00 3.464Eþ00 2.286Eþ00 5.565Eþ00 5.528Eþ00

7.240E-01 1.274Eþ00 1.411Eþ00 1.185Eþ00 6.010E-01 9.490E-01 9.520E-01 7.730E-01 1.391Eþ00 1.424Eþ00 1.288Eþ00 1.066Eþ00 1.053Eþ00 9.800E-01 7.530E-01 1.363Eþ00 1.305Eþ00 1.178Eþ00 1.083Eþ00 1.006Eþ00 7.130E-01 1.172Eþ00 9.690E-01 9.270E-01 8.600E-01 7.020E-01 6.120E-01 5.940E-01

6.859E-01 5.464E-01 7.736E-01 6.706E-01 1.475Eþ00 4.197E-01 7.299E-01 5.148E-01 4.295E-01 5.587E-01 2.324E-01 5.597E-01 1.896E-01 4.098E-01 5.284E-01 3.655E-01 4.577E-01 2.540E-01 5.510E-01 6.000E-01 2.813E-01 5.093E-01 3.086E-01 4.334E-01 6.928E-01 2.857E-01 1.637E-01 5.025E-01

Table 12 The calculated neutron multiplication factor for VVER-1000 reactor core. Number of elements

Linear

RPE (%)

Quadratic

RPE (%)

30,148 35,412 47,612

1.05156Eþ00 1.05101Eþ00 1.05078Eþ00

1.93420E-01 1.41020E-01 1.19100E-01

1.04957Eþ00 1.04956Eþ00 1.04955Eþ00

3.81000E-03 2.86000E-03 7.77000E-03

*The reference effective multiplication factor is keff ¼ 1.04953Eþ00 (Schulz, 1996).

approximation. In the aforementioned calculation, the unstructured tetrahedron elements generated by Gambit software have been used. The advantage of these elements is their superiority in mapping the any three dimensional geometry created by Gambit. Also, as discussed by the author in the previously published work (Hosseini and Vosoughi, 2013), the results with RPE (%) less than 5% may be achieved using the unstructured elements. In the regions with high flux gradient, the number of the elements may be considered higher than the other regions. 5. Conclusion In the present study, the development of the GFEM-3D computational code for solving the multi-group neutron diffusion equation based on the Galerkin Finite Element Method was reported. The comparison between the considered linear and quadratic approximation for the shape function was performed.

Both the criticality and neutron fixed source calculation was done using the developed computational code. The calculations were performed using unstructured tetrahedron elements for both hexagonal and rectangular three dimensional geometries. The neutron multiplication factor, neutron flux and power distribution are the outputs of the performed calculation by GFEM-3D computational code. To validate the calculations, an analysis of three dimensional benchmark reactor cores was carried out and the obtained results have been compared with the reference data. The accuracy of the calculation for the linear approximation improves vs. the number of elements. The results obtained from the quadratic approximation have the accuracy less than 8.75000E-03 (neutron multiplication factor) and 1.475Eþ00 (power distribution) for all the number of the elements, while the results obtained from the linear approximation have the accuracy of the order 1.93420E-01 for the neutron multiplication factor and 5.907Eþ00 for the power distribution. The neutron fixed source

130

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

Fig. 12. The comparisons between the averaged relative power in each fuel assembly calculated by GFEM-3D (Linear) and reference value (Schulz, 1996).

Fig. 13. The comparisons between the averaged relative power in each fuel assembly calculated by GFEM-3D (Quadratic) and reference value (Schulz, 1996).

problem was also validated through comparison of the calculated neutron flux distributions by GFEM-3D and CITATION computational codes. A reader can conclude that the developed computational code may be considered as one of the suitable tools for the static calculations of the three dimensional reactor cores. The results with RPE (%) less than 5% may be obtained by considering the quadratic approximation for the shape function in the unstructured tetrahedron elements. The developed GFEM-3D computational code is applicable to full-core fuel management and design applications studies of nuclear reactor cores.

Appendix A The general form for exact solution for integrals that may be obtained in the calculation of 3D neutron diffusion equation is as Eq. (A.1):

Z V

La1 Lb2 Lc3 Ld4 dxdydz ¼

a!b!c!d! 6V ða þ b þ c þ d þ 3Þ!

 Linear approximation for shape function

(A.1)

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

131

The solution of the first integral in left side of Eq. (14) is as Eq. (A.2):

2 Z V

ðeÞ

TðeÞ





dvV N ðrÞV N ¼

ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b1 b1 þ c1 c1 þ d1 d1

ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b1 b2 þ c1 c2 þ d1 d2

ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b1 b3 þ c1 c3 þ d1 d3

ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b1 b4 þ c1 c4 þ d1 d4

3

6 ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ 7 7 1 6 6 b2 b1 þ c2 c1 þ d2 d1 b2 b2 þ c2 c2 þ d2 d2 b2 b3 þ c2 c3 þ d2 d3 b2 b4 þ c2 c4 þ d2 d4 7 ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ 7 36V ðeÞ 6 4 b3 b1 þ c3 c1 þ d3 d1 b3 b2 þ c3 c2 þ d3 d2 b3 b3 þ c3 c3 þ d3 d3 b3 b4 þ c3 c4 þ d3 d4 5 ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b4 b1 þ c4 c1 þ d4 d1

ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b4 b2 þ c4 c2 þ d4 d2

ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b4 b3 þ c4 c3 þ d4 d3

ðeÞ ðeÞ

ðeÞ ðeÞ

ðeÞ ðeÞ

b4 b4 þ c4 c4 þ d4 d4

(A.2)

ðeÞ

ðeÞ

ðeÞ

ðeÞ

where, parameters ai ; bi ; ci ; di was defined in Eq. (6). The solution of the second integral in left side of Eq. (14) may be calculated as Eq. (A.3):

2

Z V

1 6 60 6 6 6 1 6 ðeÞ TðeÞ 6 ðeÞ 6 120 dv N ðrÞ N ¼ 6V 6   6 1 6 6 120 6 4 1 120

1 120 1 60 1 120 1 120

1 120 1 120 1 60 1 120

3 1 120 7 7 7 1 7 7 120 7 7 1 7 7 7 120 7 7 1 5 60

(A.3)

Z V

1 6 420 6 6 1 6 6 6 630 6 6 1 6 6 6 2520 6 6 6 1 6 630 6 6 6 1 6 ðeÞ TðeÞ 6 420 ðeÞ 6 dv N ðrÞ N ¼ 6V  6   6 1 6 6 2520 6 6 6 1 6 6 630 6 6 6 1 6 420 6 6 6 1 6 6 420 6 4 1 2520

1  630 8 630 1  630 2 315 2 315 1  420 2 315 2 315 1 315 1  420

1 2520 1  630 1 420 1  420 1  630 1 2520 1  420 1  630 1  420 1 2520

1 6 12 6 6 6 1 Z ðeÞ TðeÞ ðeÞ 6 ds N ðrÞ N ¼ 2D 6 6 24   A 6 1 6 6 4 24 0

1 24 1 12 1 24 0

1 24 1 24 1 12 0

3 07 7 7 7 07 7 7 7 7 07 5 0

(A.4)

where, D(e) is the surface area of the boundary triangle element. Eq. (A.4) is the local boundary condition matrix for each element in no incoming current situation. Eq. (A.4) is zero matrix for the reflective boundary condition.  Quadratic approximation for shape function

The solution of the second integral in left side of Eq. (14) may be calculated as Eq. (A.3):

2

2

1 630 2 315 1  420 8 630 2 315 1  630 2 315 1 315 2 315 1  420 

1  420 2 315 1  630 2 315 8 630 1  630 1 315 2 315 2 315 1  420

The solution of the second integral in left side of Eq. (14) in quadratic approximation of the shape function may be calculated as Eq. (A.5):

1 2520 1  420 1 2520 1  630 1  630 1 420 1  420 1  420 1  630 1 2520

1  630 2 315 1  420 2 315 1 315 1  420 8 630 2 315 2 315 1  630

1  420 2 315 1  630 1 315 2 315 1  420 2 315 8 630 2 315 1  630

1 420 1 315 1  420 2 315 2 315 1  630 2 315 2 315 8 630 1  630 

3 1 2520 7 7 7 1 7 7  420 7 7 1 7 7 7 2520 7 7 1 7 7  420 7 7 7 1 7 7  420 7 7 1 7 7 7 2520 7 7 1 7 7  7 630 7 7 1 7 7  630 7 7 7 1 7 7  630 7 7 1 5 420

(A.5)

132

S.A. Hosseini / Progress in Nuclear Energy 92 (2016) 119e132

Z V

Z A

ðeÞ

TðeÞ





ds N ðrÞ N ¼ 0:5 3 1 1 1 1   0:0  0:0 7 6 30 180 180 45 7 6 7 6 1 1 1 1 7 6  0:0 0:0  7 6 6 180 30 180 45 7 7 6 7 6 1 1 1 1 6   0:0 0:0 7 7 6 180 30 45 7 ðeÞ 6 180 D 6 7 7 6 1 32 16 16 7 6 0:0 0:0  7 6 45 180 180 180 7 6 7 6 6 1 16 32 16 7 7 6 0:0 0:0 6 45 180 180 180 7 7 6 4 1 16 16 32 5 0:0 0:0  45 180 180 180 (A.6) 2

where, D(e) has been described earlier. The solution of the appeared integrals in the fixed neutron source calculation is as Eq. (A.7) and (A.8):

3 1 6 24 7 7 6 7 6 6 1 7 7 6 ðeÞ 6 24 7 7 dv N ðrÞ ¼ 6V ðeÞ 6 6 1 7  7 6 7 6 6 24 7 7 6 4 1 5 2

Z V

24

(A.7)

3 1  6 120 7 7 6 7 6 7 6 1 7 6 7 6 30 7 6 7 6 1 7 6 7 6 6 120 7 7 6 7 6 1 7 6 7 6 30 7 6 7 6 7 6 1 7 6 7 6 30 ðeÞ 7 dv SðrÞN ðrÞ ¼ 6 6 1 7  7 6 7 6 6 120 7 7 6 7 6 1 7 6 7 6 7 6 30 7 6 7 6 1 7 6 7 6 30 7 6 7 6 7 6 1 7 6 7 6 30 7 6 4 1 5  120 2

The solution of the first integral in left side of Eq. (14) for quadratic approximation of the shape function is so long and complicated and it is not given in the present paper. The no incoming boundary condition is inserted in the calculation via the Eq. (A.6):

(A.8)

References Booth, T.E., 2006. Power iteration method for the several largest eigenvalues and eigenfunctions. Nucl. Sci. Eng. 154, 48e62. Cavdar, S., Ozgener, H., 2004. A finite element/boundary element hybrid method for 2-D neutron diffusion calculations. Ann. Nucl. Energy 31, 1555e1582. Center, A.C., 1977. Benchmark Problem Book. Report ANL-7416 (Suppl. 2). Argonne National Laboratory, Argonne, IL. Cho, N.Z., Lee, J., 2006. Analytic function expansion nodal (AFEN) method in hexagonal-Z three-dimensional geometry for neutron diffusion calculation. J. Nucl. Sci. Technol. 43, 1320e1326. Duderstadt, J.J., Hamilton, L.J., 1976. Nuclear Reactor Analysis. Feng, C., Liu, C., Zhang, M., 2005. Application of Gambit and Fluent commercial software in chemical engineering. Comput. Appl. Chem. 3, 015. Fowler, T., Vondy, D., Cunningham, G., Jul. 1971. Nuclear Reactor Core Analysis Code; CITATION. Oak Ridge National Laboratory, Oak Ridge, Tenn.. ORNL-TM-2496, Rev 2. Gonz
alez-Pintor, S., Ginestar, D., Verdú, G., 2009. High order finite element method for the lambda modes problem on hexagonal geometry. Ann. Nucl. Energy 36, 1450e1462.
bert, A., 2008. A RaviarteThomaseSchneider solution of the diffusion equation in He hexagonal geometry. Ann. Nucl. Energy 35, 363e376. Hosseini, S.A., Vosoughi, N., 2013. Development of two-dimensional, multigroup neutron diffusion computer code based on GFEM with unstructured triangle elements. Ann. Nucl. Energy 51, 213e226. Ivanov, K., Manolova, M., Apostolov, T., 1994. An effective solution scheme of a three-dimensional reactor core model in hexagonal geometry. Comput. Phys. Commun. 82, 1e16. Lamarsh, J.R., 1966. Introduction to Nuclear Reactor Theory. Addison Wesley Publishing Company. Maiani, M., Montagnini, B., 2004. A Galerkin approach to the boundary elementresponse matrix method for the multigroup neutron diffusion equations. Ann. Nucl. Energy 31, 1447e1475. Schulz, G., 1996. Solutions of a 3D VVER-1000 benchmark. In: Proc. 6-th Symposium of AER on VVER Reactor Physics and Safety, Kirkkonummi, Finland. Wang, Y., Bangerth, W., Ragusa, J., 2009. Three-dimensional h-adaptivity for the multigroup neutron diffusion equations. Prog. Nucl. Energy 51, 543e555. Zhu, J., Taylor, Z., Zienkiewicz, O., 2005. The Finite Element Method: its Basis and Fundamentals. Butterworth-Heinemann. Zimin, V.G., Baturin, D.M., 2002. Polynomial nodal method for solving neutron diffusion equations in hexagonal-z geometry. Ann. Nucl. Energy 29, 1105e1117.