- Email: [email protected]
1D-MOC on pin-resolved neutron transport calculations of PWR
Annals of Nuclear Energy 138 (2020) 107227
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Comparison of the performance of heterogeneous variational nodal method and 2D/1D-MOC on pin-resolved neutron transport calculations of PWR Yongping Wang ⇑, Chen Zhao, Hongchun Wu, Liangzhi Cao Xi’an Jiaotong University, Xi’an, Shaanxi, China
a r t i c l e
i n f o
Article history: Received 27 August 2019 Received in revised form 24 September 2019 Accepted 16 November 2019
Keywords: Variational nodal method Method of characteristic Pin-resolved PWR
a b s t r a c t Heterogeneous Variational Nodal Method (HVNM) and the 2D/1D Method of Characteristics (2D/1D-MOC) are two promising methods for high-fidelity neutronics calculation. HVNM starts from even-parity form of neutron transport equation and employs iso-parameteric finite elements, orthogonal polynomials and piece-wise constants for spatial discretization. Additionally, spherical harmonics are applied for the angular expansion. In terms of 2D/1D-MOC, it starts from first-order transport equation and splits the three-dimensional transport equation into one 2D and one 1D equations which are coupled by the radial and axial leakage terms. Then space and angle are discretized by characteristic lines in radial plane and 1D Sn calculation is employed in axial direction. Both HVNM and 2D/1D-MOC have the capability of dealing with pin-resolved geometry and strong transport effect in PWR core problem. Thus, we developed two codes PANX and NECP-X which based on the methods of HVNM and MOC respectively. In this paper, to compare their performance on pin-resolved problems in PWR, the KAIST problem, the NuScale problem and the Beavrs problem are calculated. Both computational accuracy and efficiency are evaluated. The numerical results shown that both PANX and 2D/1D-MOC can obtain accurate eigenvalues with the error less than 100 pcm. However, PANX obtains more accurate pin power distribution and achieves higher computational efficiency than 2D/1D-MOC. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction The two-step framework (Smith, 1986) has been widely applied to neutronics analysis of reactor core since 1980s, including lattice calculation for homogenized few-group cross sections and wholecore diffusion calculation for eigenvalue and power distribution. However, in recent years, high-fidelity neutronics calculation has become the research focus since the computational power has made a significant progress. Among the deterministic methods, the Method of Characteristics (MOC) (Boyd et al., 2014; Talamo, 2013; Marin-Lafleche et al., 2013) is one of the promising methods for high-fidelity calculation of reactor core. The basic idea of MOC is to generate a set of parallel rays for every discretized angle and solve the one-dimensional neutron transport equation along the rays. Thus, MOC has strong capability of treating complicated geometry. Obviously, the generation of the rays and computing the coordinates of the intersection points are two important processes in MOC. Generally, the characteristics lines can be divided into two categories: long characteristics lines (Chen et al., 2008) and modular characteristics lines (Tang and Zhang, 2009). Long characteristics lines go through the whole problem, while modular characteristics lines are generated in one of the repeated geometric modules and get connected when repetitive geometric modules constitutes the
whole problem domain. To guarantee the computational accuracy, the distance between the parallel characteristics lines should be adjusted to ensure every material region has been passed through by several lines. In two-dimensional geometry (like 2D lattice calculation), MOC has been widely used as it achieves high accuracy and efficiency, such as AutoMOC (Chen et al., 2008). If we directly apply MOC to 3D problem (generate the rays in 3D geometry), we call it 3D-MOC (Liu, 2013). 3D-MOC also obtains high accuracy, however, the computational cost becomes unacceptable because the sharp increase of the number of rays in 3D geometry. Therefore, a compromise approach which employs the MOC in the lateral plane combined with other diffusion or transport methods in the axial direction is carried out, referred as 2D/1D-MOC. It carries out an iteration between MOC calculation in the radial plane and diffusion (or transport) calculation in the axial direction, coupled by the leakage along the radial and axial surfaces respectively. This scheme has been successfully applied to 3D calculations and many popular codes have been developed, such as MPACT (Kochunas et al., 2013); nTRACER (Ryu et al., 2015) and DECART (Joo et al., 2004). Although 2D/1D-MOC reduces the computational cost of 3D-MOC to make it practical, in the meantime it introduces some errors in the coupling leakage term which sometimes even cause diverge of the 2D-1D iteration (Kelley and Larsen, 2013; Yuk and Cho, 2015).
2
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Variational Nodal Method (VNM) was first proposed by Elmer Lewis of Northwestern University in 1980s. The first VNM code VARIANT (Palmiotti et al., 1995) was then developed by Argonne National Laboratory (ANL) in 1990 s and it was embedded as a transport solver in REBUS (Toppel, 1983) code system of ANL and ERANOS (Doriath et al., 1994) code system of Atomic Energy Commission. In 2007, ANL released a new version of VARIANT which named as NODAL. Then in 2011, Idaho National Laboratory (INL) also developed the code INSTANT (Wang et al., 2011) based on VNM. In 2015, Xi’an Jiaotong University developed VIOLET-HEX and VIOLET (Li et al., 2015) for fast reactor and PWR core calculations respectively. However, the VNM codes mentioned above are all homogeneous nodal methods which requires homogenized cross section within each spatial mesh (nodes). From 1997, some research has made efforts to extend the VNM to treat heterogeneous material within the nodes, trying to eliminate the homogenization procedure. The first code (Fanning and Palmiotti, 1997) based on heterogeneous VNM was developed by University of Wisconsin and ANL. The basic idea is to divide the heterogeneous node into several homogeneous sub-domains and calculate the response matrices within each sub-domains. Then the nodal response matrices is a combination of all the sub-domains. In this method, orthogonal polynomials are applied to expand the spatial variable which becomes the main error source. Then in 2003, ANL has extended VARIANT (Smith et al., 2003) to treat heterogeneous nodes in which finite element is employed to describe complicated material boundaries within the nodes. However, the two methods mentioned above are only studied and evaluated in 2D geometry, but have not been applied to 3D core-analysis. In 2018, we proposed a 3D Heterogeneous Variational Nodal Method (HVNM) for pin-resolved problem. It treats each pin cell as one node and uses iso-parameteric finite element to explicitly describe the pin-resolve geometry. Spherical harmonics is employed for angular expansion, and polynomials and piece-wise constants are applied for radial and axial leakage expansion. Some tests have been done in our previous publication (Zhang et al., 2017). As both HVNM and 2D/1D-MOC tend to perform high-fidelity neutronics calculation without homogenization, it is necessary to
Vacuum
Fig. 2. The finite elements (black lines) and zones (yellow lines) in a pin cell. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Table 1 Eigenvalue results of KAIST problem. Methods
HVNM
2D/1D-MOC
MC
keff Error /pcm
1.13652 57
1.13660 49
1.13709 ± 0.00003 /
make a comparison between them, including accuracy and efficiency. The results can be a reference for choosing the transport solver in high-fidelity analysis. Two codes, PANX (Zhang et al., 2017) and NECP-X (Chen et al., 2018) were developed by Nuclear Engineering Computational Physics laboratory of Xi’an Jiaotong University based on HVNM and 2D/1D-MOC respectively. The paper is arranged as follows: in section II, the basic formulations of the HVNM and 2D/1D-MOC is presented; in section III, results are presented for the KAIST problem (Cho), NuScale problem (Haugh and Mohamed, 2012) and for a large three-dimensional model of the Beavrs (2013) reactor core. In section IV, the results are discussed, conclusions presented, and implications for future work outlined. 2. Theory 2.1. Formulations of HVNM
Moderator
UOX-1
Heterogeneous variational nodal methods starts from the evenparity form of neutron transport equation with isotropic scattering (the energy group and spatial index is omitted):
UOX-1
MOX-1
UOX-2 (CR)
UOX-2 (CR)
MOX-1 (BA8)
þ þ X rR1 t X rw þ Rt w ¼ Rs / þ q
Vacuum
Reflective
Baffle
UOX-1
UOX-2 (CR)
UOX-1
where the even- and odd-parity components of angular flux are defined by:
w ðXÞ ¼
UOX-2 (CR)
MOX-1
UOX-1
Reflective Fig. 1. Quarter core of KAIST problem.
1 ½wðXÞ wðXÞ 2
ð2Þ
and the scalar flux is the integral of angular flux:
Z
/¼ UOX-2 (BA16)
ð1Þ
wðXÞdX
ð3Þ
wðXÞ, Rt and Rs are the angular flux, total cross section and within-group scattering cross section respectively. In addition, q is the source term including the fission source and scattering source from other energy groups. Then the functional of the whole domain can be written as a summation of nodal functional:
3
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Fig. 3. Pin power error distribution of HVNM of KAIST problem /%.
Fig. 4. Pin power error distribution of 2D/1D-MOC of KAIST problem /%.
F½wþ ; w ¼
X m
F m ½wþ ; w
ð4Þ
Where the nodal functional is: Z Z h þ 2 i þ 2 2 F m ½wþ ; w ¼ dV dX R1 X r w þ R w R / 2/q t s t v Z Z þ 2 dA dX nC Xwþ w C
ð5Þ
The first term is an integral within the volume and the second term presents the integral on the nodal surfaces, where C and nC stands for the surface and the normal vector attached to the surface, respectively. Then, we expand the even-parity flux as:
wþ ðr; XÞ ¼ f ðzÞ g T ðx; yÞwðXÞ T
and expand the odd-parity component on the radial surface as:
w ðr; XÞ ¼ f ðzÞyT ðXÞ f ðnÞvc ; n ¼ x; y T
Pin power error/%
CPU hour
Max. power Max. error AVG RMS MRE 100
T
ð7Þ
The expansion on the axial surface is different from radial surface:
Table 2 Results of pin power and CPU time of KAIST problem. Methods
ð6Þ
HVNM
2D/1D-MOC
MC
0.50 3.34 0.50 0.01 0.59 120
1.10 6.44 1.18 0.02 1.17 160
±0.08 ±0.14 ±0.09 ±0.01 ±0.05
w ðr; XÞ ¼ yT ðXÞ h ðx; yÞvz T
ð8Þ
where f ðuÞ; u ¼ x; y; z is the vector of one dimensional orthogonal polynomials; g ðx; yÞ is a vector of the trial functions of iso-parameteric finite elements; hðx; yÞ is a vector of functions which equals to a constant in some region and equals to zero in
4
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227 þ
/ ¼ Vq Cðj j Þ qg ¼
X
Rs;g0 g /g0 þ
g 0 –g
BANK 2
þ
BANK 2
BANK 2
BANK 2
BANK 1
BANK 2
BANK 2
BANK 1
BANK 1
BANK 2
1 X v v Rf ;g0 /g0 k g g0
ð11Þ
where j and j stand for the vectors of the expansion coefficients of outgoing and incoming currents along the nodal surfaces. B, R, V and C are the nodal response matrices defined by geometry and material. Typical power iteration method is then applied to solve these three equations. The detailed derivation of this method can be found in our previous research paper (Zhang et al., 2017). HVNM solved the 3D transport equation directly without any coupling procedure. Thus, the main error source of HVNM is the spatial and angular discretization which has been analyzed in detail in reference (Zhang et al., 2017). The calculation parameters used in this paper is based on previous analysis.
BANK 2
BANK 1
ð10Þ
BANK 2
BANK 2
BANK 2
2.2. Formulations of 2D/1D-MOC BANK 2
The 2D/1D Method of Characteristic starts from the first-order neutron transport equation (the index of energy group is omitted):
X wðr; XÞ þ Rt ðr Þwðr; XÞ ¼ Q ðrÞ
ð12Þ
For a specific angle Xm , the equation can be written as: Fig. 5. Layout of NuScale problem.
nm other regions on the axial surface; yðXÞ is the vector of spherical harmonics for angular expansion. Then the expansion form of scalar flux and group source can be derived based on the expansion equations of even- and odd-parity flux. Submitting all the expansions into the nodal functional and employing the variational process finally results in the following set of three equations: þ
j ¼ Bq þ R j
ð9Þ
@wðr; Xm Þ @wðr; Xm Þ @wðr; Xm Þ þ gm þ lm þ Rt ðr Þwðr; Xm Þ @x @y @z
¼ Q ðr Þ
ð13Þ
where nm , gm and lm are the projection of Xm in the x, y and z direction respectively; Q is the isotropic source term including fission source and scattering source. Integrating the equation in z direction from zk-1/2 to zk+1/2 (the lower and upper surface of an axial mesh k) results in:
nm
@wm;k ðx;yÞ @x
þ gm
@wm;k ðx;yÞ @y
h i þ Dlzm wm;kþ1=2 ðx; yÞ wm;k1=2 ðx; yÞ k
þRt;k ðx; yÞwm;k ðx; yÞ ¼ Q k ðx; yÞ
ð14Þ
where
Table 3 Eigenvalue results of NuScale problem.
Dzk ¼ zk
Methods
HVNM
2D/1D-MOC
MC
keff Error /pcm
1.01337 90
1.01320 73
1.01247 ± 0.00003 /
þ 1=2
wm;k ðx; yÞ ¼
zk1=2
1 Dz k
Z
zkþ1=2
ð15Þ wðr; Xm Þdz
zk1=2
Fig. 6. Pin power error distribution of HVNM of NuScale problem /%.
ð16Þ
5
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Fig. 7. Pin power error distribution of 2D/1D-MOC of NuScale problem /%.
wm;k1=2 ðx; yÞ ¼ w x; y; zk1=2 ; Xm Q k ðx; yÞ ¼
1 Dz k
Z
zkþ1=2
ð17Þ
Q ðrÞdz
ð18Þ
zk1=2
TLAxial m;k ðx; yÞ ¼
lm h
Dzk
wm;kþ1=2 ðx; yÞ wm;k1=2 ðx; yÞ
i
ð19Þ
we can obtain:
@wm;k ðx; yÞ @wm;k ðx; yÞ þ gm þ Rt;k ðx; yÞwm;k ðx; yÞ nm @x @y
By defining the leakage term as:
¼ Q k ðx; yÞ TLAxial m;k ðx; yÞ
ð20Þ
0.31 -0.32 0.10 -0.12
0.11 -0.44
-0.29 -0.31
0.12 0.05
0.30 0.09
0.00
0.10 -0.50
0.34 0.13
0.09 -0.19
0.17 -0.35
0.26 -0.05
-0.06 -0.63
0.26 0.16
0.21 0.07
0.00
0.27 -0.36
0.27 -0.18
0.00
0.34 -0.02
0.18 -0.38
-0.06 -0.20
0.36 0.04
0.07 -0.23
0.36 0.17
0.20 -0.25
0.00 -0.49
-0.09 -0.31
-0.21 -0.51
-0.16 -0.36
0.13 -0.62
0.07 -0.07
0.07 -0.34
0.16 -0.23
-0.07 -0.63
-0.06 -0.43
-0.11 -0.73
0.24 -0.26
0.00
0.05 -0.32
0.10 -0.45
0.00
0.08 -0.56
0.12 -0.46
C B A
HVNM 2D/1D-MOC
0.00 Fig. 8. The locations of A, B and C assmblies.
Fig. 9. Pin power error distribution in assembly A /%.
6
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
This equation can be solved by 2D MOC by converting to this form:
dwðs; Xm Þ þ Rt ðsÞwðs; Xm Þ ¼ Q 0 ðsÞ ds
wm;y1=2 ðzÞ ¼
ð21Þ
wm;k ðzÞ ¼
where
Q 0 ðsÞ ¼ Q ðsÞ TLAxial m;k ðsÞ
ð22Þ
Q k ðzÞ ¼
Similarly, integrating Eq. (13) in x-y direction in the area xk1=2 to xkþ1=2 and yk1=2 to ykþ1=2 (the four surfaces of a radial mesh k) results in:
i g h i nm h wm;xþ1=2 ðzÞ wm;x1=2 ðzÞ þ m wm;yþ1=2 ðzÞ wm;y1=2 ðzÞ Dx Dy dwm;k ðzÞ þ Rt;k ðzÞwm;k ðzÞ þ lm dz ¼ Q k ðzÞ Z
ykþ1=2
w xk1=2 ; y; z; Xm dy
1 DxDy
TLRadial m;k ðzÞ ¼
Z
w x; yk1=2 ; z; Xm dx
xkþ1=2
ð25Þ
xk1=2
Z
Z
ykþ1=2
yk1=2
ykþ1=2 yk1=2
xkþ1=2
wðr; Xm Þdxdy
ð26Þ
xk1=2
Z
xkþ1=2
Q ðrÞdxdy
ð27Þ
xk1=2
i nm h wm;xþ1=2 ðzÞ wm;x1=2 ðzÞ Dx i g h þ m wm;yþ1=2 ðzÞ wm;y1=2 ðzÞ Dy
ð28Þ
the equation can be written as:
ð23Þ
lm 1 Dy
1 D xD y
Z
If we define the leakage term as:
where
wm;x1=2 ðzÞ ¼
1 Dx
dwm;k ðzÞ þ Rt;k ðzÞwm;k ðzÞ ¼ Q k ðzÞ TLRadial m;k ðzÞ dz
This is a one-dimensional transport equation which can be easily solved by any diffusion or transport methods. In our code, we choose Sn method as the 1D solver.
ð24Þ
yk1=2
0.80 0.71
0.45 0.20
0.64 0.19
0.07 -0.31
0.27 0.31
0.15 -0.04
0.00
0.37 -0.07
0.67 0.15
0.18 -0.19
0.12 0.01
0.51 0.12
0.58 0.37
0.35 0.28
0.18 -0.16
0.00
0.62 0.41
0.70 0.37
0.00
0.17 0.10
0.15 -0.20
0.14 0.23
0.23 -0.21
0.20 -0.28
0.44 0.15
0.43 -0.25
0.02 -0.47
0.30 0.17
0.24 0.26
0.18 -0.06
0.11 -0.02
0.01 0.07
0.36 0.07
0.47 0.18
0.24 0.18
0.26 -0.16
0.00
0.07 -0.60
0.37 0.17
0.00
0.23 -0.15
0.32 -0.23
0.00
-0.07 -0.42
0.28 -0.24
0.23 0.05
0.63 0.09
0.06 -0.32
-0.21 -0.23
0.16 -0.44
0.07 -0.28
0.14 0.11
0.42 -0.32
-0.21 -0.59
-0.02 0.08
0.09 0.15
0.01 -0.32
0.22 -0.54
0.00 0.08
0.02 -0.19
0.09 -0.28
0.16 0.03
0.21 0.29
0.37 -0.08
-0.05 -0.55
0.02 0.14
0.00
0.19 -0.06
-0.02 -0.44
0.00
0.31 -0.21
0.37 0.15
0.00
0.37 -0.30
0.48 0.28
0.00
0.23 -0.45
0.09 -0.23
0.01 -0.41
0.17 -0.31
-0.14 -0.15
-0.02 -0.17
0.37 -0.20
0.33 0.00
0.34 0.11
0.65 0.17
0.17 -0.01
0.30 0.11
0.21 -0.02
0.15 -0.37
0.24 -0.09
0.00
-0.02 -0.67
0.27 -0.36
0.12 -0.09
0.26 0.02
0.38 -0.29
0.47 0.04
0.18 -0.25
0.47 0.60
0.30 -0.34
0.00
0.40 -0.15
0.26 0.10
0.05 -0.28
-0.08 -0.36
0.15 -0.55
0.22 -0.04
0.00
0.23 -0.10
0.19 -0.33
0.00
0.27 -0.37
0.56 0.19
0.00
0.08 -0.47
0.23 -0.03
0.10 -0.34
0.20 0.04
0.29 0.21
0.10 -0.08
0.25 -0.06
0.35 -0.28
0.13 -0.02
0.25 -0.14
-0.04 -0.21
-0.06 -0.10
-0.01 -0.60
-0.02 -0.30
0.20 0.21
0.21 -0.01
0.09 0.06
0.23 -0.33
0.01 -0.32
0.35 0.10
0.41 0.19
0.28 0.10
0.19 0.01
0.13 -0.46
0.02 -0.19
0.07 -0.34
-0.06 -0.57
0.01 -0.50
0.35 0.15
0.17 -0.07
0.19 0.24
0.18 0.03
0.43 -0.05
0.44 0.25
0.20 -0.19
0.41 0.06
0.14 -0.49
HVNM 2D/1D-MOC
0.08 -0.36
ð29Þ
Fig. 10. Pin power error distribution in assembly B /%.
7
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
2D/1D-MOC respectively. Three pin-resolved problem are calculated, including the KASIT problem, the NuScale problem and the Beavrs problem. Both accuracy and efficiency are evaluated. The reference is generated by Monte Carlo (MC) calculation. All the calculations are run on the same machine and same convergence criteria are applied for fair comparisons.
We can find that the 2D and 1D equations are coupled by the leakage term. Therefore, an iteration between the 2D and 1D calculation is carried out to solve the three-dimensional equation, and we call this method as 2D/1D MOC. The derivation of MOC and the treatment of leakage terms can be detailed in Ref. (Chen et al., 2018). From above derivations, we can find that the 2D/1D MOC equations is obtained by integrating the three-dimensional neutron transport equation along the axial and radial direction, coupling by the leakage term. The derivation process does not introduce any approximation. Thus, if the integral region is small enough, the strict transport equation can be solved and very accurate results can be obtained. However, in practical calculations, considering the computational efficiency and memory requirements, approximations are usually applied which introduce some errors. In this paper, the source of error of the 2D/1D MOC is mainly from two aspects: first, the homogenization procedure in radial plane before doing the axial 1D calculation; second, the flat approximation of the leakage term on the surfaces of the coupling units (usually the pin-cell size in radial and several centimeters in axial).
3.1. Kaist problem This problem is based on the KAIST (Cho) benchmark problem. The quarter core configuration is shown in Fig. 1. The detailed description can be found in our previous research paper (Wang et al., 2019). In PANX, we treat each pin cell as one node (1.26 cm 1.26 cm) and divide the whole problem into 20 layer in axial evenly, with the height of 10 cm for each layer. Thus, the total number of nodes is 85 85 20. Within the pin-cell, the radius of the fuel rod (or guide tube) is 0.475 cm and 80 isoparametric finite elements are applied to explicitly describe the geometry of fuel rod and surroundered moderator. On the axial surface, those finite elements are divided into 8 zones for currents which are coupled between axial layers: 4 for the fuel and 4 for the moderator, as shown in Fig. 2. The expansion order for orthgonal ploynomials on the radial surfaces is 4th. In addition, the angular expansion order is P23
3. Results To compare the performance of HVNM and 2D/1D-MOC, two codes, PANX and NECP-X, were developed based on HVNM and
0.84 2.97
0.94 3.16
0.68 3.03
-0.06 1.56
-0.11 1.83
0.21 2.32
0.00
0.72 2.04
0.62 2.46
0.13 1.91
0.06 0.99
0.52 1.42
0.30 1.58
0.95 2.81
0.21 2.04
0.00
0.38 1.08
0.22 1.08
0.00
0.89 2.45
0.69 2.38
-0.01 0.60
0.23 0.68
0.08 0.95
0.04 0.97
1.25 2.25
0.95 2.56
0.62 2.49
0.09 0.41
-0.08 0.42
-0.14 0.10
-0.04 0.75
0.16 1.13
0.45 1.40
0.90 2.58
0.95 2.58
0.00
0.49 0.66
-0.29 -0.07
0.00
0.59 0.87
0.19 0.87
0.00
0.39 1.65
0.35 1.93
-0.09 0.15
-0.25 -0.19
-0.12 0.20
-0.01 0.29
0.56 0.85
0.17 0.86
0.22 1.13
0.17 1.06
-0.04 1.43
0.11 1.75
0.19 0.36
0.14 0.35
0.17 0.15
0.00 0.23
-0.08 0.22
0.31 0.44
0.12 0.64
0.57 1.46
0.04 0.81
0.31 1.74
0.42 1.97
0.00
0.42 0.28
0.42 0.28
0.00
0.44 0.51
-0.03 0.04
0.00
-0.28 0.03
0.36 0.77
0.00
0.19 1.25
0.48 1.91
0.22 0.33
0.52 0.36
0.28 0.37
0.29 0.40
0.40 0.39
0.19 0.37
0.01 0.25
-0.02 0.16
-0.34 0.18
0.10 0.44
-0.04 0.72
0.17 1.45
0.20 1.52
0.00
0.42 0.21
0.25 -0.01
0.18 0.33
0.28 0.48
0.45 0.37
0.19 0.42
0.32 0.59
0.16 0.15
-0.08 0.15
0.00
0.55 1.10
0.35 1.32
0.22 1.39
-0.16 -0.12
-0.09 -0.36
-0.01 -0.20
0.00
0.22 -0.01
0.43 0.16
0.00
0.40 0.25
0.27 0.25
0.00
0.31 0.45
0.13 0.39
-0.05 0.78
-0.16 0.88
0.18 1.22
0.45 0.37
0.14 0.05
0.44 0.37
0.01 0.03
0.14 -0.06
0.28 0.26
0.27 0.34
0.29 0.24
-0.05 0.04
0.16 0.34
0.12 0.16
-0.16 0.17
-0.06 0.35
-0.18 0.55
-0.13 0.76
0.02 0.79
0.07 -0.13
0.11 0.10
0.20 0.10
0.18 0.13
0.05 -0.02
0.32 0.36
0.23 0.32
0.30 0.39
0.02 0.09
0.13 0.26
0.31 0.58
0.11 0.41
0.14 0.56
-0.19 0.36
-0.08 0.61
0.01 0.78
HVNM 2D/1D-MOC
0.23 0.01
Fig. 11. Pin power error distribution in assembly C /%.
8
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
(Zhang et al., 2017) and P3 on the radial and axial surfaces, respectively. For NECP-X calculation, the problem is also divided into 20 layers in axial. The distance between the characteristic lines is 0.03 cm and the numbers of discreted polar and prone angle in each octant are 8 and 3 respectively. For Monte Carlo calculation, the partical number is 8 million in each cycle and 600 active cycles are carried out with first 300 cycles abandoned. Table 1 indicates both 2D/1D-MOC and HVNM obtain very accurate keff, the errors are around 50 pcm. Figs. 3. and 4 have shown the error of pin power distribution. We can find from Figs. 3 and 4 that: 1. the max. nagative error is 3.0% for HVNM and 6.0% for 2D/1D-MOC, which locates at the burnable rod in assembly UOX-2(BA16); 2. the max. positive error is 1.8% for HVNM and 6.0% for 2D/1D-MOC which locates at the interface of active core and reflector; 3. both HVNM and 2D/1D-MOC underestimate the pin power of MOX assemblies, while the assembly power error is about 1.0% and 2.0% for HVNM and 2D/1D-MOC respectively; 4. for UOX-2(CR) assemblies with control rods inserted, the max. error of is 2.0% and 1.0% for HVNM and 2D/1D-MOC at the center of the assemblies. Thus, we can conclude that the large errors of HVNM and 2D/1D-MOC usually occurs at strong absorbable rods, burnable rods and the outer surface of active core with sharp flux gradient; the error is less than 1.0% for other positions. To show the comparison more clearly, we summarize the results in Table 2. The error of the pin cell with max. power is 0.5% and the max. pin power error is 3.34% for HVNM, while the corresponding errors for 2D/1D-MOX is 1.1% and 6.44%. In addition, the AVG, RMS and MRE of HVNM is 0.50%, 0.01% and 0.59% respectively; the corresponding errors of 2D/1D-MOC is 1.18%, 0.02% and 1.17%. Therefore, we can find that HVNM achieved more accurate pin power distribution than 2D/1D-MOC. Additionally, the CPU time of HVNM is also 20% less than 2D/1DMOC.
1.6% UO2
2.4% UO2
3.1% UO2
Fig. 12. Layout of Beavrs problem.
Fuel pin cell Guide tube pin cell
3.2. NuScale core problem
Fission chamber pin cell
This problem is based on NuScale (Haugh and Mohamed, 2012) core design as shown in Fig. 5. The detailed description of the geometry and materila can be found in Ref. (Wang et al., 2019). The whole core is modeled in PANX. Each pin cell is treated as one node and the core is divided into 12 axial segaments. Thus, the total number of nodes is 239292. 96 isoparameteric finite elements are utilized to model the pin-resolved geometry. All other spatial and angular discretion is same as the KAIST problem. For NECP-X calculation, it employed same axial mesh size as PANX, the characteristic line distance and angular discretion is same as that in the KAIST problem. For MC calculation, same number of particals and active cycles is employed as before. Table 3 shows the results of eigenvalue. HVNM and 2D/1D-MOC acheves similar accuracy with the errors less than 100 pcm. The error of pin power distribution is shown in Figs. 6 and 7. We can find that: 1. most of the errors of HVNM are within 1.0% to 1.0%, a small number of errors close to the reflector are larger than 1.0% but still less than 1.5%; however, in 2D/1D-MOC results, a Table 4 Results of pin power and CPU time of NuScale problem. Methods Pin power error/%
CPU hour
Max. power Max. error AVG RMS MRE 120
Reflector
HVNM
2D/1D-MOC
MC
0.22 1.25 0.30 0.01 0.21 150
0.36 3.76 0.48 0.01 0.35 635.8
±0.08 ±0.17 ±0.07 ±0.01 ±0.06
Fig. 13. Configuration of fuel assembly of Beavrs problem.
larger number of errors adjacent to the reflector exceed 1.0% with 3.5% as the max. value. We compared the explicit pin power of three assemblies A, B, and C whose locations are marked in Fig. 8. The comparison is shown in Figs. 9–11. The upper and lower values in the squares are the normalized pin power obtained by HVNM and 2D/1DMOC, respectively. We can find the errors of assemblies A and B are all less than 1.0%. In assembly C which is adjacent to the reflector, there is only one pin cell has the error larger than 1.0% for HVNM. However, 41 pin power errors exceed 1.0% for 2D/1D-MOC calculation with 3.16% as the max. error. Obviously, HVNM obatin better pin power results than 2D/1D-MOC.
Table 5 Eigenvalue results of Beavrs problem. Methods
HVNM
2D/1D-MOC
MC
keff Error /pcm
1.26770 45
1.26766 41
1.26725 ± 0.00002 /
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
9
Fig. 14. Pin power distribution of Beavrs problem.
Table 4 summarizes the pin power comparison. The max. error, error of max. pin power, AVG, RMS and MRE of HVNM are 0.22%, 1.25%, 0.22%, 0.30% and 0.21%, and they are all less than those of 2D/1D-MOC. In terms of CPU time, the HVNM is 20% more efficient than 2D/1D-MOC. 3.3. Beavrs problem This problem is based on BEAVRS (Horelik et al., 2013) benchmark problem. It has 193 UO2 fuel assemblies which three enrichment: 1.6%, 2.4% and 3.1%. The active core is surrounded by the
reflector with the thick of one assembly. The configuration of the core is shown in Fig. 12. The layout of the assembly is typical 17 17 design as Fig. 13. The size of each pin cell is 1.26 cm and the radius of fuel rod (or guide tube) is 0.54 cm. The height of the core is 380 cm with 40 cm reflector on the top and at the bottom. The 7-group cross sections of the materials are listed in the appendix. Upper half core is calculated by PANX, NECP-X and MC code and the number of axial meshes is 12 for PANX and NECP-X. The total number of nodes for PANX is 891276. All other spatial and angular discretions are same as those of NuScale problem. As for MC
Fig. 15. Pin power error distribution of HVNM of Beavrs problem /%.
10
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Fig. 16. Pin power error distribution of 2D/1D-MOC of Beavrs problem /%.
calculation, the number of partical is 12,000,000 and 700 active cycles are carried out with first 300 cycles abandoned. Table 5 shows HVNM and 2D/1D-MOC obain accurate eigenvalue with the errors less than 50 pcm. The reference of pin power distribution calculated by MC code is shown in Fig. 14. It can be told that pin power in the inner core in relatively flat, while sharp power gradient occurs at the interface of the assemblies with high enrichment (3.1%) and low enrichment. The power peak is also locates at assemblies with high enrichment at the outer core. The error of pin power distribution of HVNM is shown in Fig. 15, compared with MC calculation. It can be found that the majority of errors are within 1.0%. The pin cells on the interface of high enrichment and low enrichment assemblies has relatively larger errors. The max. error occurs at the boundary of the active core which can be 2.0% ~ 2.5%. Overall, HVNM obtain accurate pin power, the average error is 0.45%. Fig. 16 shows the pin power error distribution of 2D/1D-MOC. Same as HVNM, most errors are less than 1.0%. However, there are some errors exceed 1.0% in the assemblies at the innermost core. Similarly, large negative errors occurs within the high enrichment assemblies. Table 6 lists the comparison of power distribution. It indicates that the error of the max. pin power calculated by 2D/1D-MOC is much larger than HVNM. They are 1.65% and 0.63% respectively. The max. error is similar for HVNM and 2D/1D-MOC, in the range of 2.0%–2.5%, both of them locate at the outermost
Table 6 Results of pin power and CPU time of Beavrs problem. Methods Pin power error /%
CPU hour
Max. power Max. error AVG RMS MRE 180
layer of the active core. The AVG and MRE of HVNM is lightly smaller than those of 2D/1D-MOC, which means HVNM acheves better accuracy. In terms of CPU time, the HVNM cost only half of the time of 2D/1D-MOC, but they are far less than MC. 4. Conclusions This paper compared the performance of HVNM and 2D/1DMOC in pin-resolved Neutron Transport Calculations of PWR. HVNM starts from even-parity form of neutron transport equation and employs iso-parameteric finite elements and spherical harmonics for spatial and angular discretization respectively. However, in 2D/1D-MOC, spatial and angular variables are discretized by characteristic lines in radial plane and 1D Sn calculation is employed in axial direction. PANX and NECP-X are developed based on the two methods. In numerical results section, the KASIT problem, the NuScale problem and the Beavrs problem are calculated for comparison. To summarize the numerical results, both HVNM and 2D/1DMOC can obtain accurate eigenvalue for pin-resolved problem with the error less than 100 pcm. However, HVNM obtains more accurate pin powers than 2D/1D-MOC. For the small reactor core KAIST problem, the error of the pin with maximum power, the maximum error and average error of 2D/1D-MOC are doubled compared with HVNM. Additionally, HVNM achieves relative higher computational efficiency than 2D/1D-MOC, especially for large reactor core like BEAVRS problem, the CPU time of HVNM is less than half of the time of 2D/1D-MOC. Acknowledgement
HVNM
2D/1D-MOC
MC
0.63 2.4 0.45 0.01 0.55 400
1.65 2.1 0.63 0.01 0.67 1250
±0.10 ±0.22 ±0.12 ±0.03 ±0.10
This work is financially supported by the National Natural Science Foundation of China (11775169 and 11735011). Appendix Tables A1–A6
11
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227 Table A1 Cross section of 1.6%UO2 in Beavrs problem (cm1). Group
1
2
3
4
5
6
7
Rgt Rgf
1.68490E1 1.39392E2
3.68396E1 8.50310E4
5.59531E1 9.78541E3
4.17033E1 1.74050E2
5.02485E1 7.66481E2
5.76455E1 1.28384E1
7.68239E1 2.58837E1
3.34539E2
2.04075E3
2.34850E2
4.17721E2
1.83956E1
3.08122E1
6.21208E1
5.92520E1 1.00870E1
4.07140E1 5.44760E2
3.31930E4 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
3.60500E1
2.36670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.86770E1
2.38670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.68140E1
1.56420E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.92050E3
3.54520E1
3.29690E2
7.63850E4
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.97670E2
3.20770E1
3.57300E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.44270E3
7.54980E2
3.33210E1
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
Table A2 Cross section of 2.4% UO2 in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt Rgf
1.68630E1 1.40792E2
3.68404E1 8.58856E4
5.59629E1 9.88375E3
4.17207E1 1.75800E2
5.03255E1 7.74185E2
5.77745E1 1.29674E1
7.70841E1 2.61439E1
3.37901E2
2.06126E3
2.37210E2
4.21919E2
1.85804E1
3.11218E1
6.27452E1
5.92520E1 1.00870E1
4.07140E1 5.44760E2
3.31930E4 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
3.60500E1
2.36670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.86770E1
2.38670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.68140E1
1.56420E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.92050E3
3.54520E1
3.29690E2
7.63850E4
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.97670E2
3.20770E1
3.57300E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.44270E3
7.54980E2
3.33210E1
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
Table A3 Cross section of 3.1% UO2 in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt Rgf
1.68581E1 1.42282E2
3.67807E1 1.30784E3
5.68419E1 1.56254E2
4.32233E1 2.85611E2
5.58721E1 1.24018E1
6.73176E1 2.11689E1
9.59349E1 4.24566E1
3.41476E2
3.13881E3
3.75008E2
6.85466E2
2.97644E1
5.08052E1
1.01896E+0
5.92520E1 1.00780E1
4.07140E1 5.43980E2
3.31930E4 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
3.59480E1
2.33750E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.87360E1
2.30110E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.69670E1
1.53950E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.50720E3
3.58350E1
3.01000E2
6.92940E4
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.05850E2
3.22200E1
3.51460E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.47810E3
7.66610E2
3.33390E1
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
Table A4 Cross section of fission chamber in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt Rgf
1.26032E1 4.79002E9
2.93160E1 5.82564E9
2.84250E1 4.63719E7
2.81020E1 5.24406E6
3.34460E1 1.45390E7
5.65640E1 7.14972E7
1.17214E+0 2.08041E6
1.32340E8
1.43450E8
1.12860E6
1.27630E5
3.53850E7
1.74010E6
5.06330E6
5.87910E1 6.61659E2
4.11760E1 5.90700E2
3.39060E4 2.83340E4
1.17610E7 1.46220E6
0.00000E+0 2.06420E8
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
2.40377E1
5.24350E2
2.49900E4
1.92390E5
2.98750E6
4.21400E7
0.00000E+0
0.00000E+0
1.83425E1
9.22880E2
6.93650E3
1.07900E3
2.05430E4
0.00000E+0
0.00000E+0
0.00000E+0
7.90769E2
1.69990E1
2.58600E2
4.92560E3
0.00000E+0
0.00000E+0
0.00000E+0
3.73400E5
9.97570E2
2.06790E1
2.44780E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
9.17420E4
3.16774E1
2.38760E1
0.00000E+0
0.00000E+0
0.00000E + 0
0.00000E+0
0.00000E+0
4.97930E2
1.09910E+0
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
12
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Table A5 Cross section of guide tube in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt R1g s R2g s R3g s R4g s R5g s R6g s R7g s
1.26032E1 6.61659E2
2.93160E1 5.90700E2
2.84240E1 2.83340E4
2.80960E1 1.46220E6
3.34440E1 2.06420E8
5.65640E1 0.00000E+0
1.17215E+0 0.00000E+0
0.00000E+0
2.40377E1
5.24350E2
2.49900E4
1.92390E5
2.98750E6
4.21400E7
0.00000E+0
0.00000E+0
1.83297E1
9.23970E2
6.94460E3
1.08030E3
2.05670E4
0.00000E+0
0.00000E+0
0.00000E+0
7.88511E2
1.70140E1
2.58810E2
4.92970E3
0.00000E+0
0.00000E+0
0.00000E+0
3.73330E5
9.97372E2
2.06790E1
2.44780E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
9.17260E4
3.16765E1
2.38770E1
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.97920E2
1.09912E+0
Table A6 Cross section of moderator and reflector in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt R1g s R2g s R3g s R4g s R5g s R6g s R7g s
1.59206E1 4.44777E2
4.12970E1 1.13400E1
5.90310E1 7.23470E4
5.84350E1 3.74990E6
7.18000E1 5.31840E8
1.25445E+0 0.00000E+0
2.65038E+0 0.00000E+0
0.00000E+0
2.82334E1
1.29940E1
6.23400E4
4.80020E5
7.44860E6
1.04550E6
0.00000E+0
0.00000E+0
3.45256E1
2.24570E1
1.69990E2
2.64430E3
5.03440E4
0.00000E+0
0.00000E+0
0.00000E+0
9.10284E2
4.15510E1
6.37320E2
1.21390E2
0.00000E+0
0.00000E+0
0.00000E+0
7.14370E5
1.39138E1
5.11820E1
6.12290E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.21570E3
6.99913E1
5.37320E1
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
1.32440E1
2.48070E+0
References Cho NZ. Benchmark problems in reactor and particle transport physics[R/OL]. [2000-06-23]. http://nurapt.kaist.ac.kr/benchmark. Boyd, W., Shaner, S., Li, L., et al., 2014. The OpenMOC method of characteristics neutral partical transport code. Ann. Nucl. Energy 68, 43–52. Chen, J., Liu, Z.Y., Chen, Z., et al., 2018. A new high-fidelity neutronics code NECP-X. Ann. Nucl. Energy 116, 417–428. Chen, Q., Wu, H., Cao, L., 2008. Auto MOC—A 2D neutron transport code for arbitrary geometry based on the method of characteristics and customization of AutoCAD. Nucl. Eng. Design 238 (10), 2828–2833. Doriath, J.Y., Ruggieri, J.M., Buzzi, G., et al., 1994. Reactor analysis using the variational nodal method implemented in the ERANOS system. Physor, Knoxville, TN USA. Apirl 11-15 III:464-471. Fanning, T.H., Palmiotti, G., 1997. Variational nodal transport methods with heterogeneous nodes. Nucl. Sci. Eng. 127, 154–168. Haugh B, Mohamed A. 2012. Comparison of a NuScale SMR Conceptual Core Design Using CASMO5/SIMULATE5 and MCNP5[C]. Physor, Knoxville, Tennessee, USA, April 15–20. Horelik N, Herman B, Forget B, et al. 2013. Benchmark for evaluation and validation of reactor simulations (BEAVRS)[C]. International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Sun Valley, Idaho. Joo, H.G., Cho, Y., Kim, K.S., et al., 2004. Methods and performance of a threedimensional whole-core transport code DeCART. American Nuclear Society, Physor, USA. Kelley, B.W., Larsen, E.W., 2013. 2D/1D approximations to the 3D neutron transport equation. I: Theory May 5-9. In: Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering. American Nuclear Society, USA, pp. 1346–1357. Kochunas, B., Collins, B., Jabaay, D., et al., 2013. Overview of development and design of MPACT: Michigan parallel characteristics transport code. Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering. American Nuclear Society, USA, pp. 42–53. Li, Y., Wang, Y., Liang, B., et al., 2015. Partitioned-Matrix acceleration to the FissionSource iteration of the Variational Nodal Method. Progr. Nucl. Energy 85, 640– 647.
Liu Z. 2013. Research on Three-dimensional modular method of characteristics of neutron transport equation. Xi’an Jiaotong University, China. Marin-Lafleche A, Smith MA, and Lee CH. 2013. PROTEUS-MOC: A 3D Deterministic Solver Incorporating 2D Method of Characteristics. In: Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Sun Valley, ID, May 5–9. Palmiotti, G., Lewis, E.E., Carrico, C.B., 1995. VARIANT: VARIational Anisotropic Nodal Transport for multidimensional Cartesian and hexagonal geometry calculation. Argonne National Laboratory, Argonne, IL USA. ANL-95/40. Ryu, M., Jung, Y.S., Cho, H.H., et al., 2015. Solution of the BEAVRS benchmark using the nTRACER direct whole core calculation code. J. Nucl. Sci. Technol. 52, 961–969. Smith, K.S., 1986. Assembly homogenization techniques for light water reactor analysis. Progr. Nucl. Energy 17 (3), 303–335. Smith, M.A., Tsoulfanidis, N., Lewis, E.E., et al., 2003. A finite subelement generalization of the variational nodal method. Nucl. Sci. Eng. 144, 36–46. Talamo, A., 2013. Numerical solution of the time dependent neutron transport equation by the method of the characteristics. J. Comput. Phys. 240, 248–267. Tang, C.T., Zhang, S.H., 2009. Development and verification of an MOC code employing assembly modular ray tracing and efficient acceleration techniques. Ann. Nucl. Energy 36 (8), 1013–1020. Toppel, B.J., 1983. The Fuel Cycle Analysis Capability REBUS-3. Argonne National Laboratory, Argonne, IL USA. ANL-83-2. Wang YQ, Rabiti C, Palmiotti G. Krylov 2011. Solvers preconditioned with the loworder red-black algorithm for the PN hybrid FEM for the INSTANT code. In: International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Rio de Janeiro, RJ, Brazil,May 8–12. Wang, Y., Zhang, T., Lewis, E.E., et al., 2019. Three-dimensional variational nodal method parallelization for pin resolved neutron transport calculations. Progr. Nucl. Energy 117. Yuk, S., Cho, N.Z., 2015. Whole-core transport solutions with 2-D/1-D fusion kernel via p-CMFD acceleration and p-CMFD embedding of nonoverlapping local/global iterations. Nucl. Sci. Eng. 181 (1), 1–6. Zhang, T., Wang, Y., Lewis, E.E., et al., 2017. A three-dimensional variational nodal method for pin resolved neutron transport analysis of pressurized water reactors. Nucl. Sci. Eng. 188, 160. Zhang, T., Lewis, E.E., Smith, M.A., et al., 2017. A variational nodal approach to 2D/1D pin-resolved neutron transport for pressurized water reactors. Nucl. Sci. Eng. 186.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Comparison of the performance of heterogeneous variational nodal method and 2D/1D-MOC on pin-resolved neutron transport calculations of PWR Yongping Wang ⇑, Chen Zhao, Hongchun Wu, Liangzhi Cao Xi’an Jiaotong University, Xi’an, Shaanxi, China
a r t i c l e
i n f o
Article history: Received 27 August 2019 Received in revised form 24 September 2019 Accepted 16 November 2019
Keywords: Variational nodal method Method of characteristic Pin-resolved PWR
a b s t r a c t Heterogeneous Variational Nodal Method (HVNM) and the 2D/1D Method of Characteristics (2D/1D-MOC) are two promising methods for high-fidelity neutronics calculation. HVNM starts from even-parity form of neutron transport equation and employs iso-parameteric finite elements, orthogonal polynomials and piece-wise constants for spatial discretization. Additionally, spherical harmonics are applied for the angular expansion. In terms of 2D/1D-MOC, it starts from first-order transport equation and splits the three-dimensional transport equation into one 2D and one 1D equations which are coupled by the radial and axial leakage terms. Then space and angle are discretized by characteristic lines in radial plane and 1D Sn calculation is employed in axial direction. Both HVNM and 2D/1D-MOC have the capability of dealing with pin-resolved geometry and strong transport effect in PWR core problem. Thus, we developed two codes PANX and NECP-X which based on the methods of HVNM and MOC respectively. In this paper, to compare their performance on pin-resolved problems in PWR, the KAIST problem, the NuScale problem and the Beavrs problem are calculated. Both computational accuracy and efficiency are evaluated. The numerical results shown that both PANX and 2D/1D-MOC can obtain accurate eigenvalues with the error less than 100 pcm. However, PANX obtains more accurate pin power distribution and achieves higher computational efficiency than 2D/1D-MOC. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction The two-step framework (Smith, 1986) has been widely applied to neutronics analysis of reactor core since 1980s, including lattice calculation for homogenized few-group cross sections and wholecore diffusion calculation for eigenvalue and power distribution. However, in recent years, high-fidelity neutronics calculation has become the research focus since the computational power has made a significant progress. Among the deterministic methods, the Method of Characteristics (MOC) (Boyd et al., 2014; Talamo, 2013; Marin-Lafleche et al., 2013) is one of the promising methods for high-fidelity calculation of reactor core. The basic idea of MOC is to generate a set of parallel rays for every discretized angle and solve the one-dimensional neutron transport equation along the rays. Thus, MOC has strong capability of treating complicated geometry. Obviously, the generation of the rays and computing the coordinates of the intersection points are two important processes in MOC. Generally, the characteristics lines can be divided into two categories: long characteristics lines (Chen et al., 2008) and modular characteristics lines (Tang and Zhang, 2009). Long characteristics lines go through the whole problem, while modular characteristics lines are generated in one of the repeated geometric modules and get connected when repetitive geometric modules constitutes the
whole problem domain. To guarantee the computational accuracy, the distance between the parallel characteristics lines should be adjusted to ensure every material region has been passed through by several lines. In two-dimensional geometry (like 2D lattice calculation), MOC has been widely used as it achieves high accuracy and efficiency, such as AutoMOC (Chen et al., 2008). If we directly apply MOC to 3D problem (generate the rays in 3D geometry), we call it 3D-MOC (Liu, 2013). 3D-MOC also obtains high accuracy, however, the computational cost becomes unacceptable because the sharp increase of the number of rays in 3D geometry. Therefore, a compromise approach which employs the MOC in the lateral plane combined with other diffusion or transport methods in the axial direction is carried out, referred as 2D/1D-MOC. It carries out an iteration between MOC calculation in the radial plane and diffusion (or transport) calculation in the axial direction, coupled by the leakage along the radial and axial surfaces respectively. This scheme has been successfully applied to 3D calculations and many popular codes have been developed, such as MPACT (Kochunas et al., 2013); nTRACER (Ryu et al., 2015) and DECART (Joo et al., 2004). Although 2D/1D-MOC reduces the computational cost of 3D-MOC to make it practical, in the meantime it introduces some errors in the coupling leakage term which sometimes even cause diverge of the 2D-1D iteration (Kelley and Larsen, 2013; Yuk and Cho, 2015).
2
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Variational Nodal Method (VNM) was first proposed by Elmer Lewis of Northwestern University in 1980s. The first VNM code VARIANT (Palmiotti et al., 1995) was then developed by Argonne National Laboratory (ANL) in 1990 s and it was embedded as a transport solver in REBUS (Toppel, 1983) code system of ANL and ERANOS (Doriath et al., 1994) code system of Atomic Energy Commission. In 2007, ANL released a new version of VARIANT which named as NODAL. Then in 2011, Idaho National Laboratory (INL) also developed the code INSTANT (Wang et al., 2011) based on VNM. In 2015, Xi’an Jiaotong University developed VIOLET-HEX and VIOLET (Li et al., 2015) for fast reactor and PWR core calculations respectively. However, the VNM codes mentioned above are all homogeneous nodal methods which requires homogenized cross section within each spatial mesh (nodes). From 1997, some research has made efforts to extend the VNM to treat heterogeneous material within the nodes, trying to eliminate the homogenization procedure. The first code (Fanning and Palmiotti, 1997) based on heterogeneous VNM was developed by University of Wisconsin and ANL. The basic idea is to divide the heterogeneous node into several homogeneous sub-domains and calculate the response matrices within each sub-domains. Then the nodal response matrices is a combination of all the sub-domains. In this method, orthogonal polynomials are applied to expand the spatial variable which becomes the main error source. Then in 2003, ANL has extended VARIANT (Smith et al., 2003) to treat heterogeneous nodes in which finite element is employed to describe complicated material boundaries within the nodes. However, the two methods mentioned above are only studied and evaluated in 2D geometry, but have not been applied to 3D core-analysis. In 2018, we proposed a 3D Heterogeneous Variational Nodal Method (HVNM) for pin-resolved problem. It treats each pin cell as one node and uses iso-parameteric finite element to explicitly describe the pin-resolve geometry. Spherical harmonics is employed for angular expansion, and polynomials and piece-wise constants are applied for radial and axial leakage expansion. Some tests have been done in our previous publication (Zhang et al., 2017). As both HVNM and 2D/1D-MOC tend to perform high-fidelity neutronics calculation without homogenization, it is necessary to
Vacuum
Fig. 2. The finite elements (black lines) and zones (yellow lines) in a pin cell. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Table 1 Eigenvalue results of KAIST problem. Methods
HVNM
2D/1D-MOC
MC
keff Error /pcm
1.13652 57
1.13660 49
1.13709 ± 0.00003 /
make a comparison between them, including accuracy and efficiency. The results can be a reference for choosing the transport solver in high-fidelity analysis. Two codes, PANX (Zhang et al., 2017) and NECP-X (Chen et al., 2018) were developed by Nuclear Engineering Computational Physics laboratory of Xi’an Jiaotong University based on HVNM and 2D/1D-MOC respectively. The paper is arranged as follows: in section II, the basic formulations of the HVNM and 2D/1D-MOC is presented; in section III, results are presented for the KAIST problem (Cho), NuScale problem (Haugh and Mohamed, 2012) and for a large three-dimensional model of the Beavrs (2013) reactor core. In section IV, the results are discussed, conclusions presented, and implications for future work outlined. 2. Theory 2.1. Formulations of HVNM
Moderator
UOX-1
Heterogeneous variational nodal methods starts from the evenparity form of neutron transport equation with isotropic scattering (the energy group and spatial index is omitted):
UOX-1
MOX-1
UOX-2 (CR)
UOX-2 (CR)
MOX-1 (BA8)
þ þ X rR1 t X rw þ Rt w ¼ Rs / þ q
Vacuum
Reflective
Baffle
UOX-1
UOX-2 (CR)
UOX-1
where the even- and odd-parity components of angular flux are defined by:
w ðXÞ ¼
UOX-2 (CR)
MOX-1
UOX-1
Reflective Fig. 1. Quarter core of KAIST problem.
1 ½wðXÞ wðXÞ 2
ð2Þ
and the scalar flux is the integral of angular flux:
Z
/¼ UOX-2 (BA16)
ð1Þ
wðXÞdX
ð3Þ
wðXÞ, Rt and Rs are the angular flux, total cross section and within-group scattering cross section respectively. In addition, q is the source term including the fission source and scattering source from other energy groups. Then the functional of the whole domain can be written as a summation of nodal functional:
3
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Fig. 3. Pin power error distribution of HVNM of KAIST problem /%.
Fig. 4. Pin power error distribution of 2D/1D-MOC of KAIST problem /%.
F½wþ ; w ¼
X m
F m ½wþ ; w
ð4Þ
Where the nodal functional is: Z Z h þ 2 i þ 2 2 F m ½wþ ; w ¼ dV dX R1 X r w þ R w R / 2/q t s t v Z Z þ 2 dA dX nC Xwþ w C
ð5Þ
The first term is an integral within the volume and the second term presents the integral on the nodal surfaces, where C and nC stands for the surface and the normal vector attached to the surface, respectively. Then, we expand the even-parity flux as:
wþ ðr; XÞ ¼ f ðzÞ g T ðx; yÞwðXÞ T
and expand the odd-parity component on the radial surface as:
w ðr; XÞ ¼ f ðzÞyT ðXÞ f ðnÞvc ; n ¼ x; y T
Pin power error/%
CPU hour
Max. power Max. error AVG RMS MRE 100
T
ð7Þ
The expansion on the axial surface is different from radial surface:
Table 2 Results of pin power and CPU time of KAIST problem. Methods
ð6Þ
HVNM
2D/1D-MOC
MC
0.50 3.34 0.50 0.01 0.59 120
1.10 6.44 1.18 0.02 1.17 160
±0.08 ±0.14 ±0.09 ±0.01 ±0.05
w ðr; XÞ ¼ yT ðXÞ h ðx; yÞvz T
ð8Þ
where f ðuÞ; u ¼ x; y; z is the vector of one dimensional orthogonal polynomials; g ðx; yÞ is a vector of the trial functions of iso-parameteric finite elements; hðx; yÞ is a vector of functions which equals to a constant in some region and equals to zero in
4
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227 þ
/ ¼ Vq Cðj j Þ qg ¼
X
Rs;g0 g /g0 þ
g 0 –g
BANK 2
þ
BANK 2
BANK 2
BANK 2
BANK 1
BANK 2
BANK 2
BANK 1
BANK 1
BANK 2
1 X v v Rf ;g0 /g0 k g g0
ð11Þ
where j and j stand for the vectors of the expansion coefficients of outgoing and incoming currents along the nodal surfaces. B, R, V and C are the nodal response matrices defined by geometry and material. Typical power iteration method is then applied to solve these three equations. The detailed derivation of this method can be found in our previous research paper (Zhang et al., 2017). HVNM solved the 3D transport equation directly without any coupling procedure. Thus, the main error source of HVNM is the spatial and angular discretization which has been analyzed in detail in reference (Zhang et al., 2017). The calculation parameters used in this paper is based on previous analysis.
BANK 2
BANK 1
ð10Þ
BANK 2
BANK 2
BANK 2
2.2. Formulations of 2D/1D-MOC BANK 2
The 2D/1D Method of Characteristic starts from the first-order neutron transport equation (the index of energy group is omitted):
X wðr; XÞ þ Rt ðr Þwðr; XÞ ¼ Q ðrÞ
ð12Þ
For a specific angle Xm , the equation can be written as: Fig. 5. Layout of NuScale problem.
nm other regions on the axial surface; yðXÞ is the vector of spherical harmonics for angular expansion. Then the expansion form of scalar flux and group source can be derived based on the expansion equations of even- and odd-parity flux. Submitting all the expansions into the nodal functional and employing the variational process finally results in the following set of three equations: þ
j ¼ Bq þ R j
ð9Þ
@wðr; Xm Þ @wðr; Xm Þ @wðr; Xm Þ þ gm þ lm þ Rt ðr Þwðr; Xm Þ @x @y @z
¼ Q ðr Þ
ð13Þ
where nm , gm and lm are the projection of Xm in the x, y and z direction respectively; Q is the isotropic source term including fission source and scattering source. Integrating the equation in z direction from zk-1/2 to zk+1/2 (the lower and upper surface of an axial mesh k) results in:
nm
@wm;k ðx;yÞ @x
þ gm
@wm;k ðx;yÞ @y
h i þ Dlzm wm;kþ1=2 ðx; yÞ wm;k1=2 ðx; yÞ k
þRt;k ðx; yÞwm;k ðx; yÞ ¼ Q k ðx; yÞ
ð14Þ
where
Table 3 Eigenvalue results of NuScale problem.
Dzk ¼ zk
Methods
HVNM
2D/1D-MOC
MC
keff Error /pcm
1.01337 90
1.01320 73
1.01247 ± 0.00003 /
þ 1=2
wm;k ðx; yÞ ¼
zk1=2
1 Dz k
Z
zkþ1=2
ð15Þ wðr; Xm Þdz
zk1=2
Fig. 6. Pin power error distribution of HVNM of NuScale problem /%.
ð16Þ
5
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Fig. 7. Pin power error distribution of 2D/1D-MOC of NuScale problem /%.
wm;k1=2 ðx; yÞ ¼ w x; y; zk1=2 ; Xm Q k ðx; yÞ ¼
1 Dz k
Z
zkþ1=2
ð17Þ
Q ðrÞdz
ð18Þ
zk1=2
TLAxial m;k ðx; yÞ ¼
lm h
Dzk
wm;kþ1=2 ðx; yÞ wm;k1=2 ðx; yÞ
i
ð19Þ
we can obtain:
@wm;k ðx; yÞ @wm;k ðx; yÞ þ gm þ Rt;k ðx; yÞwm;k ðx; yÞ nm @x @y
By defining the leakage term as:
¼ Q k ðx; yÞ TLAxial m;k ðx; yÞ
ð20Þ
0.31 -0.32 0.10 -0.12
0.11 -0.44
-0.29 -0.31
0.12 0.05
0.30 0.09
0.00
0.10 -0.50
0.34 0.13
0.09 -0.19
0.17 -0.35
0.26 -0.05
-0.06 -0.63
0.26 0.16
0.21 0.07
0.00
0.27 -0.36
0.27 -0.18
0.00
0.34 -0.02
0.18 -0.38
-0.06 -0.20
0.36 0.04
0.07 -0.23
0.36 0.17
0.20 -0.25
0.00 -0.49
-0.09 -0.31
-0.21 -0.51
-0.16 -0.36
0.13 -0.62
0.07 -0.07
0.07 -0.34
0.16 -0.23
-0.07 -0.63
-0.06 -0.43
-0.11 -0.73
0.24 -0.26
0.00
0.05 -0.32
0.10 -0.45
0.00
0.08 -0.56
0.12 -0.46
C B A
HVNM 2D/1D-MOC
0.00 Fig. 8. The locations of A, B and C assmblies.
Fig. 9. Pin power error distribution in assembly A /%.
6
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
This equation can be solved by 2D MOC by converting to this form:
dwðs; Xm Þ þ Rt ðsÞwðs; Xm Þ ¼ Q 0 ðsÞ ds
wm;y1=2 ðzÞ ¼
ð21Þ
wm;k ðzÞ ¼
where
Q 0 ðsÞ ¼ Q ðsÞ TLAxial m;k ðsÞ
ð22Þ
Q k ðzÞ ¼
Similarly, integrating Eq. (13) in x-y direction in the area xk1=2 to xkþ1=2 and yk1=2 to ykþ1=2 (the four surfaces of a radial mesh k) results in:
i g h i nm h wm;xþ1=2 ðzÞ wm;x1=2 ðzÞ þ m wm;yþ1=2 ðzÞ wm;y1=2 ðzÞ Dx Dy dwm;k ðzÞ þ Rt;k ðzÞwm;k ðzÞ þ lm dz ¼ Q k ðzÞ Z
ykþ1=2
w xk1=2 ; y; z; Xm dy
1 DxDy
TLRadial m;k ðzÞ ¼
Z
w x; yk1=2 ; z; Xm dx
xkþ1=2
ð25Þ
xk1=2
Z
Z
ykþ1=2
yk1=2
ykþ1=2 yk1=2
xkþ1=2
wðr; Xm Þdxdy
ð26Þ
xk1=2
Z
xkþ1=2
Q ðrÞdxdy
ð27Þ
xk1=2
i nm h wm;xþ1=2 ðzÞ wm;x1=2 ðzÞ Dx i g h þ m wm;yþ1=2 ðzÞ wm;y1=2 ðzÞ Dy
ð28Þ
the equation can be written as:
ð23Þ
lm 1 Dy
1 D xD y
Z
If we define the leakage term as:
where
wm;x1=2 ðzÞ ¼
1 Dx
dwm;k ðzÞ þ Rt;k ðzÞwm;k ðzÞ ¼ Q k ðzÞ TLRadial m;k ðzÞ dz
This is a one-dimensional transport equation which can be easily solved by any diffusion or transport methods. In our code, we choose Sn method as the 1D solver.
ð24Þ
yk1=2
0.80 0.71
0.45 0.20
0.64 0.19
0.07 -0.31
0.27 0.31
0.15 -0.04
0.00
0.37 -0.07
0.67 0.15
0.18 -0.19
0.12 0.01
0.51 0.12
0.58 0.37
0.35 0.28
0.18 -0.16
0.00
0.62 0.41
0.70 0.37
0.00
0.17 0.10
0.15 -0.20
0.14 0.23
0.23 -0.21
0.20 -0.28
0.44 0.15
0.43 -0.25
0.02 -0.47
0.30 0.17
0.24 0.26
0.18 -0.06
0.11 -0.02
0.01 0.07
0.36 0.07
0.47 0.18
0.24 0.18
0.26 -0.16
0.00
0.07 -0.60
0.37 0.17
0.00
0.23 -0.15
0.32 -0.23
0.00
-0.07 -0.42
0.28 -0.24
0.23 0.05
0.63 0.09
0.06 -0.32
-0.21 -0.23
0.16 -0.44
0.07 -0.28
0.14 0.11
0.42 -0.32
-0.21 -0.59
-0.02 0.08
0.09 0.15
0.01 -0.32
0.22 -0.54
0.00 0.08
0.02 -0.19
0.09 -0.28
0.16 0.03
0.21 0.29
0.37 -0.08
-0.05 -0.55
0.02 0.14
0.00
0.19 -0.06
-0.02 -0.44
0.00
0.31 -0.21
0.37 0.15
0.00
0.37 -0.30
0.48 0.28
0.00
0.23 -0.45
0.09 -0.23
0.01 -0.41
0.17 -0.31
-0.14 -0.15
-0.02 -0.17
0.37 -0.20
0.33 0.00
0.34 0.11
0.65 0.17
0.17 -0.01
0.30 0.11
0.21 -0.02
0.15 -0.37
0.24 -0.09
0.00
-0.02 -0.67
0.27 -0.36
0.12 -0.09
0.26 0.02
0.38 -0.29
0.47 0.04
0.18 -0.25
0.47 0.60
0.30 -0.34
0.00
0.40 -0.15
0.26 0.10
0.05 -0.28
-0.08 -0.36
0.15 -0.55
0.22 -0.04
0.00
0.23 -0.10
0.19 -0.33
0.00
0.27 -0.37
0.56 0.19
0.00
0.08 -0.47
0.23 -0.03
0.10 -0.34
0.20 0.04
0.29 0.21
0.10 -0.08
0.25 -0.06
0.35 -0.28
0.13 -0.02
0.25 -0.14
-0.04 -0.21
-0.06 -0.10
-0.01 -0.60
-0.02 -0.30
0.20 0.21
0.21 -0.01
0.09 0.06
0.23 -0.33
0.01 -0.32
0.35 0.10
0.41 0.19
0.28 0.10
0.19 0.01
0.13 -0.46
0.02 -0.19
0.07 -0.34
-0.06 -0.57
0.01 -0.50
0.35 0.15
0.17 -0.07
0.19 0.24
0.18 0.03
0.43 -0.05
0.44 0.25
0.20 -0.19
0.41 0.06
0.14 -0.49
HVNM 2D/1D-MOC
0.08 -0.36
ð29Þ
Fig. 10. Pin power error distribution in assembly B /%.
7
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
2D/1D-MOC respectively. Three pin-resolved problem are calculated, including the KASIT problem, the NuScale problem and the Beavrs problem. Both accuracy and efficiency are evaluated. The reference is generated by Monte Carlo (MC) calculation. All the calculations are run on the same machine and same convergence criteria are applied for fair comparisons.
We can find that the 2D and 1D equations are coupled by the leakage term. Therefore, an iteration between the 2D and 1D calculation is carried out to solve the three-dimensional equation, and we call this method as 2D/1D MOC. The derivation of MOC and the treatment of leakage terms can be detailed in Ref. (Chen et al., 2018). From above derivations, we can find that the 2D/1D MOC equations is obtained by integrating the three-dimensional neutron transport equation along the axial and radial direction, coupling by the leakage term. The derivation process does not introduce any approximation. Thus, if the integral region is small enough, the strict transport equation can be solved and very accurate results can be obtained. However, in practical calculations, considering the computational efficiency and memory requirements, approximations are usually applied which introduce some errors. In this paper, the source of error of the 2D/1D MOC is mainly from two aspects: first, the homogenization procedure in radial plane before doing the axial 1D calculation; second, the flat approximation of the leakage term on the surfaces of the coupling units (usually the pin-cell size in radial and several centimeters in axial).
3.1. Kaist problem This problem is based on the KAIST (Cho) benchmark problem. The quarter core configuration is shown in Fig. 1. The detailed description can be found in our previous research paper (Wang et al., 2019). In PANX, we treat each pin cell as one node (1.26 cm 1.26 cm) and divide the whole problem into 20 layer in axial evenly, with the height of 10 cm for each layer. Thus, the total number of nodes is 85 85 20. Within the pin-cell, the radius of the fuel rod (or guide tube) is 0.475 cm and 80 isoparametric finite elements are applied to explicitly describe the geometry of fuel rod and surroundered moderator. On the axial surface, those finite elements are divided into 8 zones for currents which are coupled between axial layers: 4 for the fuel and 4 for the moderator, as shown in Fig. 2. The expansion order for orthgonal ploynomials on the radial surfaces is 4th. In addition, the angular expansion order is P23
3. Results To compare the performance of HVNM and 2D/1D-MOC, two codes, PANX and NECP-X, were developed based on HVNM and
0.84 2.97
0.94 3.16
0.68 3.03
-0.06 1.56
-0.11 1.83
0.21 2.32
0.00
0.72 2.04
0.62 2.46
0.13 1.91
0.06 0.99
0.52 1.42
0.30 1.58
0.95 2.81
0.21 2.04
0.00
0.38 1.08
0.22 1.08
0.00
0.89 2.45
0.69 2.38
-0.01 0.60
0.23 0.68
0.08 0.95
0.04 0.97
1.25 2.25
0.95 2.56
0.62 2.49
0.09 0.41
-0.08 0.42
-0.14 0.10
-0.04 0.75
0.16 1.13
0.45 1.40
0.90 2.58
0.95 2.58
0.00
0.49 0.66
-0.29 -0.07
0.00
0.59 0.87
0.19 0.87
0.00
0.39 1.65
0.35 1.93
-0.09 0.15
-0.25 -0.19
-0.12 0.20
-0.01 0.29
0.56 0.85
0.17 0.86
0.22 1.13
0.17 1.06
-0.04 1.43
0.11 1.75
0.19 0.36
0.14 0.35
0.17 0.15
0.00 0.23
-0.08 0.22
0.31 0.44
0.12 0.64
0.57 1.46
0.04 0.81
0.31 1.74
0.42 1.97
0.00
0.42 0.28
0.42 0.28
0.00
0.44 0.51
-0.03 0.04
0.00
-0.28 0.03
0.36 0.77
0.00
0.19 1.25
0.48 1.91
0.22 0.33
0.52 0.36
0.28 0.37
0.29 0.40
0.40 0.39
0.19 0.37
0.01 0.25
-0.02 0.16
-0.34 0.18
0.10 0.44
-0.04 0.72
0.17 1.45
0.20 1.52
0.00
0.42 0.21
0.25 -0.01
0.18 0.33
0.28 0.48
0.45 0.37
0.19 0.42
0.32 0.59
0.16 0.15
-0.08 0.15
0.00
0.55 1.10
0.35 1.32
0.22 1.39
-0.16 -0.12
-0.09 -0.36
-0.01 -0.20
0.00
0.22 -0.01
0.43 0.16
0.00
0.40 0.25
0.27 0.25
0.00
0.31 0.45
0.13 0.39
-0.05 0.78
-0.16 0.88
0.18 1.22
0.45 0.37
0.14 0.05
0.44 0.37
0.01 0.03
0.14 -0.06
0.28 0.26
0.27 0.34
0.29 0.24
-0.05 0.04
0.16 0.34
0.12 0.16
-0.16 0.17
-0.06 0.35
-0.18 0.55
-0.13 0.76
0.02 0.79
0.07 -0.13
0.11 0.10
0.20 0.10
0.18 0.13
0.05 -0.02
0.32 0.36
0.23 0.32
0.30 0.39
0.02 0.09
0.13 0.26
0.31 0.58
0.11 0.41
0.14 0.56
-0.19 0.36
-0.08 0.61
0.01 0.78
HVNM 2D/1D-MOC
0.23 0.01
Fig. 11. Pin power error distribution in assembly C /%.
8
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
(Zhang et al., 2017) and P3 on the radial and axial surfaces, respectively. For NECP-X calculation, the problem is also divided into 20 layers in axial. The distance between the characteristic lines is 0.03 cm and the numbers of discreted polar and prone angle in each octant are 8 and 3 respectively. For Monte Carlo calculation, the partical number is 8 million in each cycle and 600 active cycles are carried out with first 300 cycles abandoned. Table 1 indicates both 2D/1D-MOC and HVNM obtain very accurate keff, the errors are around 50 pcm. Figs. 3. and 4 have shown the error of pin power distribution. We can find from Figs. 3 and 4 that: 1. the max. nagative error is 3.0% for HVNM and 6.0% for 2D/1D-MOC, which locates at the burnable rod in assembly UOX-2(BA16); 2. the max. positive error is 1.8% for HVNM and 6.0% for 2D/1D-MOC which locates at the interface of active core and reflector; 3. both HVNM and 2D/1D-MOC underestimate the pin power of MOX assemblies, while the assembly power error is about 1.0% and 2.0% for HVNM and 2D/1D-MOC respectively; 4. for UOX-2(CR) assemblies with control rods inserted, the max. error of is 2.0% and 1.0% for HVNM and 2D/1D-MOC at the center of the assemblies. Thus, we can conclude that the large errors of HVNM and 2D/1D-MOC usually occurs at strong absorbable rods, burnable rods and the outer surface of active core with sharp flux gradient; the error is less than 1.0% for other positions. To show the comparison more clearly, we summarize the results in Table 2. The error of the pin cell with max. power is 0.5% and the max. pin power error is 3.34% for HVNM, while the corresponding errors for 2D/1D-MOX is 1.1% and 6.44%. In addition, the AVG, RMS and MRE of HVNM is 0.50%, 0.01% and 0.59% respectively; the corresponding errors of 2D/1D-MOC is 1.18%, 0.02% and 1.17%. Therefore, we can find that HVNM achieved more accurate pin power distribution than 2D/1D-MOC. Additionally, the CPU time of HVNM is also 20% less than 2D/1DMOC.
1.6% UO2
2.4% UO2
3.1% UO2
Fig. 12. Layout of Beavrs problem.
Fuel pin cell Guide tube pin cell
3.2. NuScale core problem
Fission chamber pin cell
This problem is based on NuScale (Haugh and Mohamed, 2012) core design as shown in Fig. 5. The detailed description of the geometry and materila can be found in Ref. (Wang et al., 2019). The whole core is modeled in PANX. Each pin cell is treated as one node and the core is divided into 12 axial segaments. Thus, the total number of nodes is 239292. 96 isoparameteric finite elements are utilized to model the pin-resolved geometry. All other spatial and angular discretion is same as the KAIST problem. For NECP-X calculation, it employed same axial mesh size as PANX, the characteristic line distance and angular discretion is same as that in the KAIST problem. For MC calculation, same number of particals and active cycles is employed as before. Table 3 shows the results of eigenvalue. HVNM and 2D/1D-MOC acheves similar accuracy with the errors less than 100 pcm. The error of pin power distribution is shown in Figs. 6 and 7. We can find that: 1. most of the errors of HVNM are within 1.0% to 1.0%, a small number of errors close to the reflector are larger than 1.0% but still less than 1.5%; however, in 2D/1D-MOC results, a Table 4 Results of pin power and CPU time of NuScale problem. Methods Pin power error/%
CPU hour
Max. power Max. error AVG RMS MRE 120
Reflector
HVNM
2D/1D-MOC
MC
0.22 1.25 0.30 0.01 0.21 150
0.36 3.76 0.48 0.01 0.35 635.8
±0.08 ±0.17 ±0.07 ±0.01 ±0.06
Fig. 13. Configuration of fuel assembly of Beavrs problem.
larger number of errors adjacent to the reflector exceed 1.0% with 3.5% as the max. value. We compared the explicit pin power of three assemblies A, B, and C whose locations are marked in Fig. 8. The comparison is shown in Figs. 9–11. The upper and lower values in the squares are the normalized pin power obtained by HVNM and 2D/1DMOC, respectively. We can find the errors of assemblies A and B are all less than 1.0%. In assembly C which is adjacent to the reflector, there is only one pin cell has the error larger than 1.0% for HVNM. However, 41 pin power errors exceed 1.0% for 2D/1D-MOC calculation with 3.16% as the max. error. Obviously, HVNM obatin better pin power results than 2D/1D-MOC.
Table 5 Eigenvalue results of Beavrs problem. Methods
HVNM
2D/1D-MOC
MC
keff Error /pcm
1.26770 45
1.26766 41
1.26725 ± 0.00002 /
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
9
Fig. 14. Pin power distribution of Beavrs problem.
Table 4 summarizes the pin power comparison. The max. error, error of max. pin power, AVG, RMS and MRE of HVNM are 0.22%, 1.25%, 0.22%, 0.30% and 0.21%, and they are all less than those of 2D/1D-MOC. In terms of CPU time, the HVNM is 20% more efficient than 2D/1D-MOC. 3.3. Beavrs problem This problem is based on BEAVRS (Horelik et al., 2013) benchmark problem. It has 193 UO2 fuel assemblies which three enrichment: 1.6%, 2.4% and 3.1%. The active core is surrounded by the
reflector with the thick of one assembly. The configuration of the core is shown in Fig. 12. The layout of the assembly is typical 17 17 design as Fig. 13. The size of each pin cell is 1.26 cm and the radius of fuel rod (or guide tube) is 0.54 cm. The height of the core is 380 cm with 40 cm reflector on the top and at the bottom. The 7-group cross sections of the materials are listed in the appendix. Upper half core is calculated by PANX, NECP-X and MC code and the number of axial meshes is 12 for PANX and NECP-X. The total number of nodes for PANX is 891276. All other spatial and angular discretions are same as those of NuScale problem. As for MC
Fig. 15. Pin power error distribution of HVNM of Beavrs problem /%.
10
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Fig. 16. Pin power error distribution of 2D/1D-MOC of Beavrs problem /%.
calculation, the number of partical is 12,000,000 and 700 active cycles are carried out with first 300 cycles abandoned. Table 5 shows HVNM and 2D/1D-MOC obain accurate eigenvalue with the errors less than 50 pcm. The reference of pin power distribution calculated by MC code is shown in Fig. 14. It can be told that pin power in the inner core in relatively flat, while sharp power gradient occurs at the interface of the assemblies with high enrichment (3.1%) and low enrichment. The power peak is also locates at assemblies with high enrichment at the outer core. The error of pin power distribution of HVNM is shown in Fig. 15, compared with MC calculation. It can be found that the majority of errors are within 1.0%. The pin cells on the interface of high enrichment and low enrichment assemblies has relatively larger errors. The max. error occurs at the boundary of the active core which can be 2.0% ~ 2.5%. Overall, HVNM obtain accurate pin power, the average error is 0.45%. Fig. 16 shows the pin power error distribution of 2D/1D-MOC. Same as HVNM, most errors are less than 1.0%. However, there are some errors exceed 1.0% in the assemblies at the innermost core. Similarly, large negative errors occurs within the high enrichment assemblies. Table 6 lists the comparison of power distribution. It indicates that the error of the max. pin power calculated by 2D/1D-MOC is much larger than HVNM. They are 1.65% and 0.63% respectively. The max. error is similar for HVNM and 2D/1D-MOC, in the range of 2.0%–2.5%, both of them locate at the outermost
Table 6 Results of pin power and CPU time of Beavrs problem. Methods Pin power error /%
CPU hour
Max. power Max. error AVG RMS MRE 180
layer of the active core. The AVG and MRE of HVNM is lightly smaller than those of 2D/1D-MOC, which means HVNM acheves better accuracy. In terms of CPU time, the HVNM cost only half of the time of 2D/1D-MOC, but they are far less than MC. 4. Conclusions This paper compared the performance of HVNM and 2D/1DMOC in pin-resolved Neutron Transport Calculations of PWR. HVNM starts from even-parity form of neutron transport equation and employs iso-parameteric finite elements and spherical harmonics for spatial and angular discretization respectively. However, in 2D/1D-MOC, spatial and angular variables are discretized by characteristic lines in radial plane and 1D Sn calculation is employed in axial direction. PANX and NECP-X are developed based on the two methods. In numerical results section, the KASIT problem, the NuScale problem and the Beavrs problem are calculated for comparison. To summarize the numerical results, both HVNM and 2D/1DMOC can obtain accurate eigenvalue for pin-resolved problem with the error less than 100 pcm. However, HVNM obtains more accurate pin powers than 2D/1D-MOC. For the small reactor core KAIST problem, the error of the pin with maximum power, the maximum error and average error of 2D/1D-MOC are doubled compared with HVNM. Additionally, HVNM achieves relative higher computational efficiency than 2D/1D-MOC, especially for large reactor core like BEAVRS problem, the CPU time of HVNM is less than half of the time of 2D/1D-MOC. Acknowledgement
HVNM
2D/1D-MOC
MC
0.63 2.4 0.45 0.01 0.55 400
1.65 2.1 0.63 0.01 0.67 1250
±0.10 ±0.22 ±0.12 ±0.03 ±0.10
This work is financially supported by the National Natural Science Foundation of China (11775169 and 11735011). Appendix Tables A1–A6
11
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227 Table A1 Cross section of 1.6%UO2 in Beavrs problem (cm1). Group
1
2
3
4
5
6
7
Rgt Rgf
1.68490E1 1.39392E2
3.68396E1 8.50310E4
5.59531E1 9.78541E3
4.17033E1 1.74050E2
5.02485E1 7.66481E2
5.76455E1 1.28384E1
7.68239E1 2.58837E1
3.34539E2
2.04075E3
2.34850E2
4.17721E2
1.83956E1
3.08122E1
6.21208E1
5.92520E1 1.00870E1
4.07140E1 5.44760E2
3.31930E4 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
3.60500E1
2.36670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.86770E1
2.38670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.68140E1
1.56420E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.92050E3
3.54520E1
3.29690E2
7.63850E4
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.97670E2
3.20770E1
3.57300E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.44270E3
7.54980E2
3.33210E1
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
Table A2 Cross section of 2.4% UO2 in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt Rgf
1.68630E1 1.40792E2
3.68404E1 8.58856E4
5.59629E1 9.88375E3
4.17207E1 1.75800E2
5.03255E1 7.74185E2
5.77745E1 1.29674E1
7.70841E1 2.61439E1
3.37901E2
2.06126E3
2.37210E2
4.21919E2
1.85804E1
3.11218E1
6.27452E1
5.92520E1 1.00870E1
4.07140E1 5.44760E2
3.31930E4 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
3.60500E1
2.36670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.86770E1
2.38670E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.68140E1
1.56420E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.92050E3
3.54520E1
3.29690E2
7.63850E4
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.97670E2
3.20770E1
3.57300E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.44270E3
7.54980E2
3.33210E1
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
Table A3 Cross section of 3.1% UO2 in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt Rgf
1.68581E1 1.42282E2
3.67807E1 1.30784E3
5.68419E1 1.56254E2
4.32233E1 2.85611E2
5.58721E1 1.24018E1
6.73176E1 2.11689E1
9.59349E1 4.24566E1
3.41476E2
3.13881E3
3.75008E2
6.85466E2
2.97644E1
5.08052E1
1.01896E+0
5.92520E1 1.00780E1
4.07140E1 5.43980E2
3.31930E4 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
3.59480E1
2.33750E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.87360E1
2.30110E3
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
3.69670E1
1.53950E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.50720E3
3.58350E1
3.01000E2
6.92940E4
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.05850E2
3.22200E1
3.51460E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.47810E3
7.66610E2
3.33390E1
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
Table A4 Cross section of fission chamber in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt Rgf
1.26032E1 4.79002E9
2.93160E1 5.82564E9
2.84250E1 4.63719E7
2.81020E1 5.24406E6
3.34460E1 1.45390E7
5.65640E1 7.14972E7
1.17214E+0 2.08041E6
1.32340E8
1.43450E8
1.12860E6
1.27630E5
3.53850E7
1.74010E6
5.06330E6
5.87910E1 6.61659E2
4.11760E1 5.90700E2
3.39060E4 2.83340E4
1.17610E7 1.46220E6
0.00000E+0 2.06420E8
0.00000E+0 0.00000E+0
0.00000E+0 0.00000E+0
0.00000E+0
2.40377E1
5.24350E2
2.49900E4
1.92390E5
2.98750E6
4.21400E7
0.00000E+0
0.00000E+0
1.83425E1
9.22880E2
6.93650E3
1.07900E3
2.05430E4
0.00000E+0
0.00000E+0
0.00000E+0
7.90769E2
1.69990E1
2.58600E2
4.92560E3
0.00000E+0
0.00000E+0
0.00000E+0
3.73400E5
9.97570E2
2.06790E1
2.44780E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
9.17420E4
3.16774E1
2.38760E1
0.00000E+0
0.00000E+0
0.00000E + 0
0.00000E+0
0.00000E+0
4.97930E2
1.09910E+0
v Rgf vg
R1g s R2g s R3g s R4g s R5g s R6g s R7g s
12
Y. Wang et al. / Annals of Nuclear Energy 138 (2020) 107227
Table A5 Cross section of guide tube in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt R1g s R2g s R3g s R4g s R5g s R6g s R7g s
1.26032E1 6.61659E2
2.93160E1 5.90700E2
2.84240E1 2.83340E4
2.80960E1 1.46220E6
3.34440E1 2.06420E8
5.65640E1 0.00000E+0
1.17215E+0 0.00000E+0
0.00000E+0
2.40377E1
5.24350E2
2.49900E4
1.92390E5
2.98750E6
4.21400E7
0.00000E+0
0.00000E+0
1.83297E1
9.23970E2
6.94460E3
1.08030E3
2.05670E4
0.00000E+0
0.00000E+0
0.00000E+0
7.88511E2
1.70140E1
2.58810E2
4.92970E3
0.00000E+0
0.00000E+0
0.00000E+0
3.73330E5
9.97372E2
2.06790E1
2.44780E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
9.17260E4
3.16765E1
2.38770E1
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
4.97920E2
1.09912E+0
Table A6 Cross section of moderator and reflector in Beavrs problem (cm1). 能群g
1
2
3
4
5
6
7
Rgt R1g s R2g s R3g s R4g s R5g s R6g s R7g s
1.59206E1 4.44777E2
4.12970E1 1.13400E1
5.90310E1 7.23470E4
5.84350E1 3.74990E6
7.18000E1 5.31840E8
1.25445E+0 0.00000E+0
2.65038E+0 0.00000E+0
0.00000E+0
2.82334E1
1.29940E1
6.23400E4
4.80020E5
7.44860E6
1.04550E6
0.00000E+0
0.00000E+0
3.45256E1
2.24570E1
1.69990E2
2.64430E3
5.03440E4
0.00000E+0
0.00000E+0
0.00000E+0
9.10284E2
4.15510E1
6.37320E2
1.21390E2
0.00000E+0
0.00000E+0
0.00000E+0
7.14370E5
1.39138E1
5.11820E1
6.12290E2
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
2.21570E3
6.99913E1
5.37320E1
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
0.00000E+0
1.32440E1
2.48070E+0
References Cho NZ. Benchmark problems in reactor and particle transport physics[R/OL]. [2000-06-23]. http://nurapt.kaist.ac.kr/benchmark. Boyd, W., Shaner, S., Li, L., et al., 2014. The OpenMOC method of characteristics neutral partical transport code. Ann. Nucl. Energy 68, 43–52. Chen, J., Liu, Z.Y., Chen, Z., et al., 2018. A new high-fidelity neutronics code NECP-X. Ann. Nucl. Energy 116, 417–428. Chen, Q., Wu, H., Cao, L., 2008. Auto MOC—A 2D neutron transport code for arbitrary geometry based on the method of characteristics and customization of AutoCAD. Nucl. Eng. Design 238 (10), 2828–2833. Doriath, J.Y., Ruggieri, J.M., Buzzi, G., et al., 1994. Reactor analysis using the variational nodal method implemented in the ERANOS system. Physor, Knoxville, TN USA. Apirl 11-15 III:464-471. Fanning, T.H., Palmiotti, G., 1997. Variational nodal transport methods with heterogeneous nodes. Nucl. Sci. Eng. 127, 154–168. Haugh B, Mohamed A. 2012. Comparison of a NuScale SMR Conceptual Core Design Using CASMO5/SIMULATE5 and MCNP5[C]. Physor, Knoxville, Tennessee, USA, April 15–20. Horelik N, Herman B, Forget B, et al. 2013. Benchmark for evaluation and validation of reactor simulations (BEAVRS)[C]. International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Sun Valley, Idaho. Joo, H.G., Cho, Y., Kim, K.S., et al., 2004. Methods and performance of a threedimensional whole-core transport code DeCART. American Nuclear Society, Physor, USA. Kelley, B.W., Larsen, E.W., 2013. 2D/1D approximations to the 3D neutron transport equation. I: Theory May 5-9. In: Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering. American Nuclear Society, USA, pp. 1346–1357. Kochunas, B., Collins, B., Jabaay, D., et al., 2013. Overview of development and design of MPACT: Michigan parallel characteristics transport code. Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering. American Nuclear Society, USA, pp. 42–53. Li, Y., Wang, Y., Liang, B., et al., 2015. Partitioned-Matrix acceleration to the FissionSource iteration of the Variational Nodal Method. Progr. Nucl. Energy 85, 640– 647.
Liu Z. 2013. Research on Three-dimensional modular method of characteristics of neutron transport equation. Xi’an Jiaotong University, China. Marin-Lafleche A, Smith MA, and Lee CH. 2013. PROTEUS-MOC: A 3D Deterministic Solver Incorporating 2D Method of Characteristics. In: Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Sun Valley, ID, May 5–9. Palmiotti, G., Lewis, E.E., Carrico, C.B., 1995. VARIANT: VARIational Anisotropic Nodal Transport for multidimensional Cartesian and hexagonal geometry calculation. Argonne National Laboratory, Argonne, IL USA. ANL-95/40. Ryu, M., Jung, Y.S., Cho, H.H., et al., 2015. Solution of the BEAVRS benchmark using the nTRACER direct whole core calculation code. J. Nucl. Sci. Technol. 52, 961–969. Smith, K.S., 1986. Assembly homogenization techniques for light water reactor analysis. Progr. Nucl. Energy 17 (3), 303–335. Smith, M.A., Tsoulfanidis, N., Lewis, E.E., et al., 2003. A finite subelement generalization of the variational nodal method. Nucl. Sci. Eng. 144, 36–46. Talamo, A., 2013. Numerical solution of the time dependent neutron transport equation by the method of the characteristics. J. Comput. Phys. 240, 248–267. Tang, C.T., Zhang, S.H., 2009. Development and verification of an MOC code employing assembly modular ray tracing and efficient acceleration techniques. Ann. Nucl. Energy 36 (8), 1013–1020. Toppel, B.J., 1983. The Fuel Cycle Analysis Capability REBUS-3. Argonne National Laboratory, Argonne, IL USA. ANL-83-2. Wang YQ, Rabiti C, Palmiotti G. Krylov 2011. Solvers preconditioned with the loworder red-black algorithm for the PN hybrid FEM for the INSTANT code. In: International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Rio de Janeiro, RJ, Brazil,May 8–12. Wang, Y., Zhang, T., Lewis, E.E., et al., 2019. Three-dimensional variational nodal method parallelization for pin resolved neutron transport calculations. Progr. Nucl. Energy 117. Yuk, S., Cho, N.Z., 2015. Whole-core transport solutions with 2-D/1-D fusion kernel via p-CMFD acceleration and p-CMFD embedding of nonoverlapping local/global iterations. Nucl. Sci. Eng. 181 (1), 1–6. Zhang, T., Wang, Y., Lewis, E.E., et al., 2017. A three-dimensional variational nodal method for pin resolved neutron transport analysis of pressurized water reactors. Nucl. Sci. Eng. 188, 160. Zhang, T., Lewis, E.E., Smith, M.A., et al., 2017. A variational nodal approach to 2D/1D pin-resolved neutron transport for pressurized water reactors. Nucl. Sci. Eng. 186.