- Email: [email protected]
IFDF acceleration method with adaptive diffusion coefficients for SN nodal calculation in SARAX code system
Annals of Nuclear Energy 136 (2020) 107056
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
IFDF acceleration method with adaptive diffusion coefficients for SN nodal calculation in SARAX code system Zhitao Xu, Youqi Zheng, Yongping Wang ⇑, Hongchun Wu School of Nuclear Science and Technology, Xi’an Jiaotong University, 28 West Xianning Road, Xi’an, Shaanxi 710049, China
a r t i c l e
i n f o
Article history: Received 10 May 2019 Received in revised form 9 September 2019 Accepted 11 September 2019 Available online 19 September 2019 Keywords: IFDF CMFD SN nodal method Fast reactor
a b s t r a c t SN nodal method has been successfully applied to fast reactor core neutronics analysis. The computational efficiency can be enhanced efficiently by the popular coarse-mesh finite difference (CMFD) method. However, the CMFD method suffers from the problem of instability, especially for optically thick problems. In this paper, another efficient nonlinear acceleration method, named interface flux discontinuity factor (IFDF) acceleration, was studied and improved for the SN nodal calculation in hexagonal prism and cuboid meshes. With adaptive diffusion coefficients, the IFDF is guaranteed to be positive and within a reasonable region, which makes the acceleration stable. It was found that the IFDF and CMFD methods obtained similar speedups for the tested fast reactor problems. However, the stability of the IFDF method was much better than the CMFD method even with artificial diffusion improvement. Additionally, the rebalance accelerated IFDF method was applied to the external-source multiplication problem besides the common eigenvalue problem. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Fast reactor core calculation is critical for fast reactor nuclear design analysis. Due to the complex energy spectrum (Yang, 2012) and the strong neutronics anisotropy in fast reactor core, multigroup (for example, 33 energy groups (Ruggieri et al., 2006) transport calculation is essential. However, the computational cost of multigroup transport increases by about three orders of magnitude compared with the traditional few-group diffusion calculation. Thus, it is necessary to study efficient numerical methods and acceleration techniques to obtain higher efficiency. The existing researches show that the discrete ordinates (SN) nodal method well balances the calculation accuracy, stability, memory and efficiency (Xu et al., 2018). Meanwhile, in fast reactors, the mean free path of neutrons is comparable to the size of one assembly and homogeneous cross sections are usually used in the core calculation, so it is not necessary to subdivide the mesh within the hexagonal or rectangular assembly. What’s more, the ray-effect of the SN method are negligible for the reactor core calculation because the fission and scattering sources are distributed all over the core. Therefore, it is of great significance and application value to study efficient SN nodal method in hexagonal prism and cuboid meshes.
⇑ Corresponding author. E-mail address: [email protected] (Y. Wang). https://doi.org/10.1016/j.anucene.2019.107056 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
Nonlinear acceleration methods which employ the strategy of nonlinear iteration between accurate high-order problem and efficient low-order problem are one kind of efficient acceleration method and they usually have better performance than the Krylov methods (Willert et al., 2014). The CMFD acceleration method (Full-core, 2002) is one of the most widely applied nonlinear methods employed in the method of characteristics (MOC) or SN calculation. It can provide a speedup from 10 to 50 for most cases. However, for some cases, especially for optically thick problems, CMFD could cause iterative divergence due to the inconsistency between the high-order problem and the low-order problem. Many innovative schemes for improving the stability and acceleration performance of the CMFD acceleration have been proposed (Cho, 2017; Jarrett et al., 2016; Zhu et al., 2016; Yuk and Cho, 2017; Xiao et al., 2018; Yamamoto et al., 2018; Masiello, 2018; Masiello et al., 2009; Xu et al., 2012) in recent years. These methods focus on improving the consistency of high-order problem and low-order problem by modifying coupling formulation or coupling meshes. The CMFD acceleration with stability techniques including the artificial diffusion and multiple transport sweeps mentioned in reference (Jarrett et al., 2016) has been implemented in the SN nodal methods (Xu et al., 2018; Lu and Wu, 2007; Xu et al., 2019) for hexagonal-z or triangular-z geometries in SARAX (Zheng et al., 2018; Zheng et al., 2018) code system. However, for a large MOX-3600 fast reactor with large assembly size of 21 cm, the CMFD acceleration with the artificial diffusion and multiple
2
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
transport sweeps techniques for the SN nodal calculation with hexagonal prism meshes still failed. Therefore, another nonlinear acceleration method (Aragones and Ahnert, 1986) employing the interface flux discontinuity factor (IFDF) was studied in this paper. For convenience, it is called IFDF acceleration. IFDF acceleration has been studied for short characteristics transport calculation (Badruzzaman, 1985) but the benefits was not clear, because the original reference (Aragones and Ahnert, 1986) did not clearly illustrated how to control the transport correction factor. In this paper, the IFDF acceleration is modified by giving a detailed formulation of adaptive diffusion coefficients to adjust the transport correction factor for the stability of acceleration and it is implemented in the SN nodal calculation with hexagonal prism and cuboid meshes. The paper is organized as follows. In Section 2, the introduction of the SN nodal formulations, the IFDF formulations, the adaptive diffusion coefficients, and the iteration schemes for both the eigenvalue and the external-source multiplication problems are presented. In Section 3, numerical results for a MOX-3600 fast reactor, a subcritical small fast reactor with external-source, and a ZPPR-10B fast reactor are given and discussed. Finally, Section 4 gives the conclusions of this paper.
X
lmd @
wm ðx; y; zÞ þ Rt wm ðx; y; zÞ ¼ Q m ðx; y; zÞ
d m m m y ðxÞwm x ðxÞ þ px ys ðxÞwx ðxÞ ¼ qx ðxÞ dx s
ð2Þ
where
pffiffiffi ys ðxÞ ¼ ð1 jxjÞ= 3 1 ¼ 2ys ðxÞ
wm x ð xÞ
pm x ¼
Zze
ð3Þ
yZ s ðx Þ
wm ðx; y; zÞdydz
ð5Þ
lmx
hr
m m ys ðxÞ Q m x ðxÞ Lax ðxÞ Lr x ðxÞ
qm x ðxÞ ¼
lmx
Qm x ðxÞ ¼
1 2ys ðxÞ
Lam x ð xÞ ¼
wm z ðxÞ
ð4Þ
zb ys ðxÞ
Rt hr
2.1. Introduction of the SN nodal formulations To conduct the multigroup neutronics transport calculation, the spatial domain is discretized into hexagonal prism or cuboid nodes (meshes). With the SN discretization for angle, nodal response relations which couple the nodal outgoing and incoming surface fluxes, along with the nodal interior fluxes and source, can be obtained for every node. The nodal response relations are employed for transport sweeps to get the final solution. The SN nodal response relations can be written as the SN nodal formulations which are introduced as follows. First let’s consider the hexagonal prism meshes. For a hexagonal node with the flat distance hd (unit: cm) along d (d ¼ v ; x; u; z) channel, the discrete neutron transport equation within one node can be written in a standard hexagonal prism illustrated in Fig. 1 as
ð1Þ
where lm d is direction cosine relative to d coordinate axis direction (channel), hd is the actual side-length (unit: cm) of the node along d channel, hv ¼ hx ¼ hy ¼ hr is satisfied and hr is the radial flat distance of the node, wm is the neutron flux (unit: cm2s1) for discrete angle m, Rt is the total cross section (unit: cm1) of this node, Q m is the neutron source (unit: cm3s1) including scattering source, fission source and external source. With the geometric symmetry of the node and the selected SN angular quadrature, we only need to consider the case when lmd > 0. The transverse integral technique is employed to decouple the 3-D equation, take x direction as an example,
2. Theory The introduction of the SN nodal formulations with hexagonal prism and cuboid meshes is given in subsection 2.1, which is the foundation of this work. In subsection 2.2., the IFDF formulations are derived in detail. The equations of the transport correction factor will be given, and the iteration schemes for both the eigenvalue and the external-source multiplication problems will be presented. Then, the formulations of adaptive diffusion coefficients are derived in subsection 2.3.
hd @d
d¼v ;x;u;z
lmz hz
Zze
yZ s ðx Þ
Q m ðx; y; zÞdydz
ð7Þ
zb ys ðxÞ
m wm zþ ðxÞ wz ðxÞ
1 ¼ 2ys ðxÞ
ð6Þ
ð8Þ
yZ s ðx Þ
wm ðx; y; 1=2Þdy
ð9Þ
ys ðxÞ
sgnðxÞ m m lu wusgnðxÞ x ðxÞ þ lmv wmv sgnðxÞ x ðxÞ Lrm x ðxÞ ¼ pffiffiffi 3hr
ð10Þ
Zze wm u x ðxÞ
¼
w½x; ys ðxÞ; zdz
ð11Þ
wm ½x; ys ðxÞ; zdz
ð12Þ
zb
Zze wm v x ðxÞ ¼ zb
pffiffiffi
l ¼ m u
lmx þ 3lmy 2
ð13Þ
pffiffiffi
lv ¼ m
Fig. 1. Illustration of a standard hexagonal prism.
lmx 3lmy 2
ð14Þ
For other directions, similar transverse integral process can be conducted. We expand the 1-D nodal interior flux and source with the orthogonal polynomials. Take x direction as an example,
3
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056 3 X
wm x ð xÞ ¼
wm xi hi ðxÞ
ð15Þ
Qm xi hi ðxÞ
ð16Þ
i¼0
Qm x ð xÞ ¼
3 X i¼0
where hi ðxÞ is a selected orthogonal polynomials for radial channels. And we expand the 1-D nodal surface fluxes as: m m wm z ðxÞ ¼ wz þ wz x1 h1 ðxÞ
m where em z and F z0 are also nodal integral coefficients which can be calculated before the transport sweeps. These additional highorder nodal surface flux formulations are employed to update the m m transverse terms Lam d1 , Lzd1 and Lz1 . Then let’s consider the cuboid meshes. For a cuboid node with the flat distance hs (cm) along d (d ¼ x; y; z) direction, the discrete neutron transport equation of one node (shown as Fig. 2) can be written as
ð17Þ
X
d ¼ v; u
ð18Þ
d¼x;y;z
ðzÞ ¼ wm þ wm f ðzÞ; d ¼ v ; x; u wm d z d d z1 1
ð19Þ
After similar transverse integral process and applying polynomial expansion as the z direction of the hexagonal prism case, the nodal formulations for cuboid meshes can be obtained. Take the z direction as an example,
ð xÞ wm d
¼
wm ; d
where f 1 ðzÞ is a selected orthogonal polynomials for z direction. Through some complex mathematical derivations, we can obtain the formulations coupling the outgoing nodal surface flux m , with the incoming nodal surface wm and the nodal average flux w dþ
m m flux wm ¼ Qm d0 ): d and the nodal interior source Q dj (Q
m m m Am þ wþ ¼ A w þ d
wm zþ
m
ð20Þ
m Am þ wm 1 Bm þ Dm ¼w z z z z
ð21Þ
lmd
¼ w m
@ m w ðx; y; zÞ þ Rt wm ðx; y; zÞ ¼ Q m ðx; y; zÞ @d
m þ P wm Bm Dm lm =hd Q d d d d d¼x;y;z
Rt þ
P d¼x;y;z
m Am d ld =hd
m m m m wm þ Dm zþ ¼ w Az þ wz 1 Bz z
where
wm þ ¼ wm ¼
wm vþ
wm xþ
wm uþ
m T w
wm v
wm x
wm u
wm z
T
3 3 P m m m Qm v j C v j Lav 1 C v 1 7 6 j¼1 7 6 7 6 3 7 6 P m m m m 7 6 Q C La C x1 x1 xj xj 7 6 j¼1 7 6 m 7 d ¼6 7 6 3 P m m m 7 6 m Q C La C 7 6 u1 u1 uj uj 7 6 j¼1 7 6 ! 7 6 3 5 4 m lmz P m m m m Q hz Q zj C zj Lz1 C z1
ð22Þ ð23Þ
2
ð24Þ
Dm z ¼
3
X m m Qm zj Lzj W zij
m m Qm zj Lzj C zj
ð25Þ
j¼1 m m m m and Am þ , A , Ad , Bd and C d are nodal integral coefficients which can be calculated before the transport sweeps. Then, we can obtain the formulations coupling the nodal interior flux wm di with the m incoming nodal surface flux wm d and nodal interior source Q dj ,
ð32Þ
ð33Þ
j¼1
Different from the hexagonal prism nodal formulations, nonlinear terms (Xu et al., 2019) are not appeared in these formulations. Therefore these formulations obey the standard nodal-equivalent finite difference (NEFD) (Azmy, 1997) scheme. However, for the additional high-order nodal surface flux formulations, this paper extends the simplified method to cuboid mesh to update the high-order transverse leakage term Lm d1 ,
ð34Þ
3
X m m m m m m Sm xj Laxj W xij Lnx0 Wni0 Lpx0 Wpi0
With the nodal response relations presented in Section 2.1, the transport solutions can be obtained by iterations with assumed initial solution. Since the critical characteristic is important for reactor physics analysis, the eigenvalue (k-eigenvalue) problem is a common topic. While for the analysis of reactor neutron dynamics, subcritical reactor neutronics and uncertainty quantification, the external-source multiplication problem is encountered. To unify the two kinds of problem, the steady-state neutron transport equation can be formulated as follows,
j¼1
ð26Þ m m m m wm zi ¼ wz0 U zi þ wz V zi þ
ð31Þ
2.2. IFDF formulations
3 X
m m m m wm xi ¼ wx0 U xi þ wx V xi þ
m m m m wm zi ¼ wz0 U zi þ wz V zi þ
m m m m wm x1zþ wx1z ez þ Q x1 F z0
j¼1
ð30Þ
3
X m m Qm zj Lzj W zij
ð27Þ
j¼1 m m m m where U m di , V di , W di , Wndi , Wpdi are also nodal integral coefficients which can be calculated before the transport sweeps. Additional high-order nodal surface flux formulations are obtained by a simplified method (Xu et al., 2019):
m m m m wm x1zþ wx1z ez þ Q x1 F z0
ð28Þ
m m wm xþ z1 Q z1 W z11
ð29Þ
Fig. 2. Illustration of a standard cuboid.
4
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
1 Lw ¼ SU þ FU þ Q ex k
ð35Þ
where L is an operator including neutron leakage and total reaction rate; w is the angular flux vector; S is the scattering operator; U is the scalar flux vector. k equals to the effective multiplication factor keff for the eigenvalue problem and it equals to constant 1.0 for the external-source multiplication problem. F is the fission operator and Q ex is the external-source vector which only exists for the external-source multiplication problem. For external-source multiplication problem, to evaluate the multiplication performance of the system, a source multiplication factor ks is defined as follows,
ks ¼
hFUi hFUi hQ ex i
ð36Þ
Lw
ðn1Þ
¼ SU
þ
1 k
ðn1Þ
ðn1Þ
FU
þ Q ex
Uðn1=2Þ ¼ Twðn1=2Þ D ðn1=2Þ
k
ðn1Þ
¼k
ð37Þ ð38Þ
E
FUðn1=2Þ D E FUðn1Þ
ð39Þ
where T is the operator to transfer the angular flux vector to the scalar flux vector, n is the iteration index, ðn 1=2Þ means the value has just been updated by the high-order (SN nodal) operators from the ðn 1Þ-th value, and Eq. (39) is only valid for the eigenvalue problem. However, the source iteration converges very slowly, especially for problems with large dominance ratio. For large fast reactor, the dominance ratio is close to 1.0 and the number of meshes is also very large. Therefore, an efficient acceleration method is important. Here we consider the IFDF acceleration method. The IFDF acceleration employs the nonlinear iteration between the highorder (SN nodal) problem and the low-order (diffusion finitedifference) problem. Since hexagonal meshes cannot be subdivided or combined, we use the same meshes for the low-order problem as the high-order meshes (nodes). In addition, the high-order and the low-order problem utilize the same number of energy groups to avoid additional inconsistency caused by cross sections homogenization. The low-order problem to be solved can be written as below,
ML U00 ¼
1 FL U00 þ Q L;ex k
ð40Þ
where ML is an operator including diffusion, absorption and scattering; U00 is the nodal average scalar flux vector of all groups; F is the fission operator for U00 ; Q L;ex is the external-source vector for low-order problem. Eq. (40) is much less expensive to be solved compared with Eq. (35), but the accuracy is lower than that of Eq. (35). A proper nonlinear iteration between the high-order problem and the low-order problem can combine their advantages. The nonlinear iteration scheme for the IFDF acceleration is similar with that of the CMFD acceleration, while the distinct difference is the formulations of the low-order problem, which are presented as follows. For hexagonal prism meshes, the diffusion finite difference formulation for a single group can be written as below,
X J K d þ1=2 J K d 1=2 AK d þ Rr;K /K V K ¼ Q K V K
d¼v ;x;u;z
X J K d þ1=2 J K d 1=2 AK d þ Rr;K /K V K ¼ Q K V K
ð42Þ
d¼x;y;z
where h i is the operator for integral over the full phase space. The final solutions can be obtained by the source iteration (SI). In this work, a outer fission-source iteration and an inner scattering-source iteration for the nodal average flux of every group are employed. Let’s focus on the outer iteration, which can be expressed as below, ðn1=2Þ
where K is the nodal index; JK d 1=2 is the neutron current (unit: cm2s1) in the normal direction d along mesh surface K d 1=2; K d means the value is related with the direction d of mesh K; K d 1=2 means the value is related with the mesh surface; AK d is the mesh surface area (unit: cm2) along direction d of mesh K; Rr;K is the removal cross section (unit: cm1) of mesh K; /K is the scalar flux (unit: cm2s1) of mesh K; V K is the volume (unit: cm3) of mesh K; Q K is the neutron source (unit: cm3s1) of mesh K, which includes scattering source from other groups, fission source and external source. For cuboid meshes, the diffusion finite difference formulation is similar,
ð41Þ
The Fick’s law gives the relation between neutron current and neutron flux. For an inner mesh K, the Fick’s law can be written as below,
J K d þ1=2 ¼
DK d /K d þ1=2 /K d
ð43Þ
hK d =2
where hK d is the mesh size along direction d of mesh K. However, the Fick’s law gives large error for transport problem. A discontinuK þ1=2
ity factor f K dd
J K d þ1=2 ¼
is employed:
K þ1=2 DK d f K dd /K d þ1=2 /K d hK d =2
ð44Þ
For mesh K þ 1,
J K d þ1=2 ¼
K þ1=2 DK d þ1 /K d þ1 f K dd þ1 /K d þ1=2 hK d þ1 =2
ð45Þ
Combine Eqs. (44) and (45), we obtain
J K d þ1=2 ¼
DK d þ1=2 rK d þ1=2 /K d þ1 /K d hK d þ1=2
ð46Þ
where
DK d þ1=2 ¼
2DK d DK d þ1 hK d þ1=2 DK d r K d þ1=2 hK d þ1 þ DK d þ1 hK d hK
/K d J K d þ1=2 2DKd
K þ1=2
rK d þ1=2 ¼
hK d þ1=2 ¼
f K dd
K þ1=2
f K dd þ1
ð47Þ
¼
d
hK
þ1
/K d þ1 þ J K d þ1=2 2DKd þ1
ð48Þ
d
hK d þ hK d þ1 2
ð49Þ
The ratio rK d þ1=2 of two discontinuity factors are called as IFDF by reference (Aragones and Ahnert, 1986), and it is also called as transport correction factor in this paper. In the nonlinear acceleration process, the transport correction factor is updated with Eq. (48) by substituting the flux and current with the corresponding values from high-order calculation, hK
rK d þ1=2 ¼
SN d /SN K d J K d þ1=2 2DK
d
hK
þ1
SN d /SN K d þ1 þ J K d þ1=2 2DK þ1
ð50Þ
d
where the superscript SN means the value is obtained from the high-order SN nodal calculation. For the case that mesh K exists and mesh K d þ 1 does not exists, the modified Fick’s law is written as below,
5
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
J þK d þ1=2
/K þ1=2 DK d ¼ d 4 2 /K þ1=2 DK d ¼ d þ 4 2
K þ1=2 f K dd /K d þ1=2
¼
K 1
ð51Þ
/K d
Jþ K d þ1=2
and
J K d þ1=2
is expressed as below,
bK d þ1=2 J þK d þ1=2
ð53Þ ð54Þ
where bK d þ1=2 is the boundary albedo of mesh surface K d þ 1=2. Combine Eqs. (51)–(54), and substitute the flux and current with the corresponding values from high-order calculation, we obtain
DK d þ1=2 /K d hK d þ1=2
ð55Þ
where
ð56Þ
hK d aK d þ1=2 þ 2DK d r K dd
K þ1=2
hK
¼
K þ1=2 f K dd
¼
SN d /SN K d J K d þ1=2 2DK
ð57Þ
d
JSN K þ1=2 d
aK d þ1=2
1 bK d þ1=2
¼ 2 1 þ bK d þ1=2
aK d þ1=2
ð58Þ
hK d þ1=2 ¼ hK d
ð59Þ
For the case that mesh K exists and mesh K d 1 does not exists, similar formulations are obtained:
J K d 1=2 ¼
DK d 1=2 rK d 1=2 /K d hK d 1=2
ð60Þ
where
dK ¼ Rr;K þ
DK d 1=2
K þ1
¼ cK dd
K 1
¼ cK dd
dK dd
K 1=2
¼
dK dd
K 1=2
aK d 1=2
1
¼
þ
JSN K 1=2 hK d
2DK
ð62Þ d
d
1 bK d 1=2
¼ 2 1 þ bK d 1=2
ð63Þ
ð64Þ
Combing Eq. (41), (46), (55), and (60), for hexagonal prism meshes we obtain
X
d¼v ;x;u;z
where
dK ¼ Rr;K þ K þ1
dK dd
K þ1
/K d þ1 dK dd
X
d¼v ;x;u;z
K 1
/K d 1 dK dd
¼ QK
X K þ1=2 K 1=2 cK dd DK d þ1=2 þ cK dd DK d 1=2 r K d 1=2
d¼x;u;v ;z
K þ1=2
¼ cK dd
DK d þ1=2 r K d þ1=2
ð65Þ
J SN K d þ1=2 ¼
ð67Þ
ð71Þ ð72Þ
DK d 1=2
ð73Þ
1 hK d hK d 1=2
ð74Þ
X
es;K d þ1=2
s¼x;y;z
X m
lms xm wmKd þ1=2
ð75Þ
where es;K d þ1=2 is the projection value of the normal vector of mesh surface K d þ 1=2 at s direction, and xm is the quadrature weight of the discrete angle m. With the solution of Eqs. (65) or (70) equation for every group, the solution of Eq. (40) can be obtained by fission-source iteration. The high-order values obtained by Eqs. (38) and (39) are employed as the initial values of the low-order fission-source iteration, which is expressed as below,
Uð00n1;lÞ ¼
1 k
ðn1;l1Þ
¼k
ðn1;l1Þ
ðn1;l1Þ
D D
FL U00
ðn1;lÞ
FL U00
ð76Þ
E
ðn1;l1Þ
FL U00
þ Q L;ex
E
ð77Þ
where the low-order iteration index l is employed, and Eq. (77) is only valid for the eigenvalue problem. For the external-source multiplication problem, the fission-source iteration converges slowly, the rebalance iteration (Xu et al., 2017) is employed to further accelerate the low-order fission-source iteration. The rebalance iteration employs a global rebalance factor f by assuming that the spatial shape of the flux does not change quickly during iterations. The original iteration (without eigenvalue) is expressed as: ðn1Þ
ML
Uð00n1;l1=2Þ ¼ FL Uð00n1;l1Þ þ Q L;ex
ð78Þ
The rebalance method introduces a factor to enforce the updated flux to satisfy the neutron balance equation:
f ð66Þ
ð70Þ
The linear system constituted by equation or equation can be solved by the red–black algorithm (Zhou et al., 2018). After the high-order fission-source iteration as expressed by Eqs. (37)–(39), the updated high-order flux and current are used to calculate IFDF with Eqs. (50), (57), and (62). The current is calculated with the angular flux. Take the current at mesh surface K d þ 1=2 as an example,
ðn1;lÞ
aK d 1=2
/SN Kd
¼ QK
DK d þ1=2 rK d þ1=2
ð61Þ
hK d 1=2 ¼ hK d
/K dK
¼
d
K 1=2
K 1
/K d 1 dK dd
X K þ1=2 K 1=2 cK dd DK d þ1=2 þ cK dd DK d 1=2 r K d 1=2
J SN K 1=2
f K dd
X
d¼x;y;z
K þ1=2
k r K dd
ð69Þ
d¼x;y;z
ML
2DK d aK d 1=2 hK d 1=2 ¼ hK d aK d 1=2 r K d 1=2 þ 2DK d
K þ1
/K d þ1 dK dd
d¼x;y;z
ðn1Þ
where
2=ð3hr hr Þ; d ¼ v ; x; u 1= hK d hK d 1=2 ; d ¼ z
X
/K dK
cK dd
2DK d aK d þ1=2 hK d þ1=2
¼
ð68Þ
DK d 1=2
For cuboid meshes we obtain
K 1=2
DK d þ1=2 ¼
K þ1=2 r K dd
K 1=2
¼ cK dd
(
ð52Þ
J K d þ1=2 ¼ J þK d þ1=2 J K d þ1=2
J K d þ1=2 ¼
dK dd
K 1=2 cK dd
hK d =2
The relation between
J K d þ1=2
/K d þ1=2 /K d hK d =2
J K d þ1=2
K þ1=2
f K dd
ðn1;lÞ
D
ðn1Þ
ML
E
D
E
Uð00n1;l1=2Þ ¼ f ðn1;lÞ FL Uð00n1;l1=2Þ þ Q L;ex
ð79Þ
Thus, combine Eqs. (78) and (79) to obtain the expression of the factor:
f
ðn1;lÞ
¼
Q L;ex
E
D E D ðn1;l1=2Þ ðn1;l1Þ Q L;ex FL U00 þ FL U00
ð80Þ
6
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
Then the flux is modified by multiplying the factor: ðn1;lÞ 00
U
¼f
ðn1;lÞ
ðn1;l1=2Þ 00
U
ð81Þ
Through above process, the fission-source iteration can be well accelerated. After the low-order problem is converged, we update the highorder values as follows, ðnÞ
ðn1=2Þ
/i;K ¼ /i;K ðn Þ
k
/DIF K ðn1=2Þ /0;K
;i P 0
ð82Þ
DIF
¼k
ð83Þ
where the superscript DIF means the value is from the converged low-order eigenvalue problem, i is the spatial moment index of the nodal polynomial expansion, and Eq. (83) is only valid for the eigenvalue problem. The overall nonlinear iteration scheme of the IFDF iteration is illustrated in the flowchart of Fig. 3. The loworder problem is calculated for certain high-order fission-source iterations to avoid unnecessary low-order calculation efforts and improve acceleration stability. The diffusion coefficients are updated to ensure the IFDF is positive but not too large, which will be explained in subection 2.3. 2.3. Adaptive diffusion coefficients
J SN K d þ1=2
> 0, the fol-
lowing condition must be satisfied,
hK d
ð84Þ
2/SN Kd
If J SN K d þ1=2 < 0, the following condition must be satisfied,
DK d þ1 > J SN K d þ1=2
hK d þ1
ð85Þ
2/SN K d þ1
and the following condition must also be satisfied for adjacent mesh K þ 1,
DK d > J SN K d 1=2
hK d
ð86Þ
2/SN Kd
For the case that mesh K exists and mesh K d þ 1 does not exist, IFDF is defined by Eq. (57). For the reflective boundary condition, DK d þ1=2 equals to zero and IFDF makes no sense. Thus, we only need to consider the vacuum boundary condition, which means the condition that J SN K d þ1=2 > 0 is always satisfied. The same condition of Eq. (84) must be satisfied. For the case that mesh K exists and mesh K d 1 does not exist, IFDF is defined by Eq. (57). For the same reason as above, the condition that J SN K d 1=2 < 0 is always satisfied. The same condition as Eq. (86) must be satisfied. Combine Eqs. (84) and (86), the condition is written as below,
DK d >
As IFDF is defined as the ratio of the discontinuity factors along one mesh surface, it should always be positive to be physically significant. In addition, it should stay within a reasonable region (positive but not too large). This subsection will give the condition to bound IFDF or the transport correction factor. First let’s find the condition that makes IFDF to be positive. For an inner mesh K, IFDF is defined by Eq. (50). If
DK d > J SN K d þ1=2
1 2/SN Kd
SN max hK d J SN K d 1=2 ; hK d J K d þ1=2
ð87Þ
Eq. (87) must be satisfied for every channels to ensure that all the IFDFs around a mesh are positive, so the condition is written as below,
DK > D0K
ð88Þ
for hexagonal meshes,
D0K ¼
1 2/SN K
SN SN SN max hK v J SN K v 1=2 ; hK v J K v þ1=2 ; hK x J K x 1=2 ; hK x J K x þ1=2 ;
SN SN SN hK u J SN K u 1=2 ; hK u J K u þ1=2 ; hK z J K z 1=2 ; hK z J K z þ1=2
ð89Þ
for cuboid meshes,
D0K ¼
1 2/SN K
SN SN SN max hK x J SN K x 1=2 ; hK x J K x þ1=2 ; hK y J K y 1=2 ; hK y J K y þ1=2 ;
SN hK z J SN K z 1=2 ; hK z J K z þ1=2
ð90Þ
To ensure IFDF is positive but not too large, this paper gives the following formulation for adaptive determination of the diffusion coefficients of every mesh before update IFDF,
1 DK ¼ max ; cD0K 3Rt;K
ð91Þ
where Rt;K is the total cross section of mesh K, c is a damping parameter obtained by experience. The damping parameter is used to keep r in a reasonable region (as is shown by Eq. (94)) and to adjust the ‘step size’ of the iteration towards convergence. It can be regarded as an adjuster between the speedup and stability of this acceleration method. Here the diffusion coefficient is selected for the adaptive value because it can help IFDF improve the interface current consistency between the high-order problem and the loworder problem without affecting the reaction rates. To understand the effect of the damping parameter c, let’s consider a simple model that all the mesh surface currents equal to jJ j and jJ j is large enough that cD0K > 3R1t;K and all the mesh sizes Fig. 3. Nonlinear iteration scheme of the IFDF acceleration.
equal to h, the following equations can be obtained,
7
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
DK ¼ cD0K ¼ c
hjJ j
ð92Þ
2/SN K
DKþ1 ¼ cD0K ¼ c
hjJ j
ð93Þ
2/SN Kþ1
Substitute Eqs. (92) and (93) into (50), we obtain
r K d þ1=2 ¼
/SN K /SN Kþ1
1þ
2 c1
ð94Þ
From Eq. (94) it’s seen that if c near 1, rK d þ1=2 will be very large and even exceed the computer floating number limit; if c is very SN large, rK d þ1=2 will be bounded by /SN K =/Kþ1 . For this work we found that if c is very large, the acceleration speedup is low, so a proper value c ¼ 5 is employed after many numerical tests. From the numerical tests below it will be seen that c ¼ 5 is conservative to guarantee iteration stability. Smaller c is recommended if one expects higher speedup and larger c is recommended if one expects higher iteration stability.
3. Results The IFDF acceleration method described above has been implemented in the SARAX-DNTH (Discrete Nodal Transport Code for Hexagonal-Z Meshes) code and SARAX-DNTC (Discrete Nodal Transport Code for Cuboid Meshes) code in NECP-SARAX (SARAX) code system. In this section, a MOX-3600 fast reactor with large assemblies, a small subcritical fast reactor and a ZPPR-10B large fast reactor with square assemblies were calculated to evaluate the modified IFDF acceleration method.
Fig. 4. Radial core layout of 3600 MWth oxide core (Takeda and Ikeda, 1991).
For all the calculations, S4 Legendre-Chebyshev angular quadrature was used; the nodal interior variables were expanded with 2nd-order polynomials; the nodal surface variables were expanded with 1st-order polynomials; the fission-source iteration criterion of the high-order problem was 1 105 ; the scattering-source iteration criterion of the high-order problem was 5 106 ; an optimized iteration limit of 16 (Bernnat et al., 2015) was employed if only source iteration was carried out (no acceleration); an iteration limit of 5 was employed if nonlinear acceleration was applied; the fissionsource iteration criterion of the low-order problem was 2:5 106 ; one low-order problem was solved before four high-order fissionsource iterations were performed. All the following calculations were performed on a 3.4 GHz Intel Xeon Gold 6134 CPU core. 3.1. The MOX-3600 fast reactor problem This is a benchmark problem of fast reactor (Takeda and Ikeda, 1991) with large assembly size. The radial flat distance of the hexagonal assembly is 21.2205 cm and the height of the reactor core is 301.6974 cm. This core is consisted of 225 inner fuel assemblies, 228 outer fuel assemblies, 331 radial reflector assemblies and 33 control assemblies. The radial layout of the core is illustrated in Fig. 4. 33-group homogenized cross sections for every assembly were obtained by the TULIP code in SARAX. With these cross sections, we have conducted the whole-core adjoint transport calculation with SARAX-DNTH (DNTH). The core was divided by 21 axial meshes and 817 radial hexagonal meshes. The overview calculation results are presented in Table 1, in which SI means that only source iteration was employed (no acceleration); CMFD/adCMFD means that the traditional CMFD or artificial diffusion (g ¼ 1=4) CMFD (Jarrett et al., 2016) acceleration was employed; IFDF (c = 0) means the IFDF acceleration without adaptive diffusion coefficients was employed; IFDF (c = 1.001), IFDF (c = 1.1), IFDF (c = 5) and IFDF (c = 1000) mean that the improved IFDF acceleration with different c in Eq. (91) was employed. In addition, one transport sweep means conducting flux calculation for all angles and meshes for one group. As shown in Table 1, CMFD/ adCMFD, IFDF (c = 0) and IFDF (c = 1.001) could not converge for this problem, while IFDF (c = 1.1), IFDF (c = 5) and IFDF (c = 1000) gave effective results. IFDF (c = 0) failed due to the adaptive diffusion coefficients were not employed, while IFDF (c = 1.0001) failed due to IFDF was too large. The calculated keff of IFDF (c = 1.1) was 1.00671, which was the same as the 1.00671 of SI, while the calculated keff of IFDF (c = 5) and IFDF (c = 1000) was slightly different from that of SI, which may be caused by the possible inconsistency between the high-order problem and low-order problem. The total CPU time of IFDF (c = 1.1) was significantly reduced from 5364.7 s to 345.1 s, which was mainly due to the obvious reduction of the number of transport sweep. The low-order calculation time of IFDF (c = 1.1) was 41.3 s, which was the overhead of the IFDF acceleration. With the increase of c, the speedup was decreased but the divergence risk of large IFDF was also decreased. It was remarkable that the low-order calculation time of IFDF (c = 1000) was smaller than that of IFDF (c = 5). Fig. 5 presents the convergence process
Table 1 Results for the MOX-3600 problem. Item
SI
CMFD /adCMFD
IFDF (c = 0)
IFDF (c = 1.001)
IFDF (c = 1.1)
IFDF (c = 5)
IFDF (c = 1000)
keff Number of fission-source iteration Number of transport sweep Total CPU time/s Speedup Low-order calculation time/s
1.00671 392 48,408 5364.7 1.0
Not converged
Not converged
Not converged
1.00671 22 2732 345.1 15.5 41.3
1.00670 26 3197 402.0 13.3 42.7
1.00673 84 9954 1147.0 4.7 22.3
8
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
of SI and IFDF (c = 5). It clearly shows that IFDF (c = 5) converged much faster than SI. Table 2 presents the region average adjoint flux results, which shows IFDF acceleration caused small derivation to the distributed solutions. 3.2. The small subcritical fast reactor problem
Fig. 5. Iteration errors of different methods.
Table 2 Region average adjoint flux results for the MOX-3600 problem. Material region
Inner fuel 1 Inner fuel 2 Inner fuel 3 Inner fuel 4 Inner fuel 5 Outer fuel 1 Outer fuel 2 Outer fuel 3 Outer fuel 4 Outer fuel 5 Lower gas Lower reflector Lower gas Upper reflector Radial reflector Primary control Secondary control Duct
Normalized region average adjoint flux SI
IFDF(c = 5)
Deviation/%
3.911E08 5.614E08 6.102E08 5.290E08 3.370E08 2.566E08 3.723E08 4.079E08 3.584E08 2.331E08 2.148E09 1.111E08 1.526E08 2.985E09 1.379E09 3.157E09 2.645E09 2.344E08
3.915E08 5.619E08 6.106E08 5.291E08 3.368E08 2.556E08 3.709E08 4.063E08 3.569E08 2.320E08 2.163E09 1.114E08 1.515E08 2.997E09 1.393E09 3.162E09 2.650E09 2.345E08
0.1 0.1 0.1 0.0 0.0 0.4 0.4 0.4 0.4 0.5 0.7 0.3 0.7 0.4 1.0 0.1 0.2 0.0
This is a small subcritical fast reactor problem modified from the KNK-II (Lesage, 2001) problem. The radial flat distance of the hexagonal assembly is 12.9904 cm and the reactor core has a height of 190 cm and is consisted of 169 hexagonal assemblies. Fig. 6 indicates the radial layout of the reactor core. The core layout and 4-group cross sections are totally the same as the control rod half-inserted case described in reference (Lesage, 2001). However, in this case, an external source is added to the central assembly in the height range of 90–100 cm. The external-source density from the first group to the last group are 1.0, 0.1, 0.001, 0.0001 (unit: cm3s1) respectively. This problem has been calculated with DNTH by different iteration schemes, including SI, IFDF, IFDF-RE (rebalance accelerated IFDF), CMFD-RE (rebalance accelerated CMFD), and adCMFD-RE (rebalance accelerated adCMFD, g ¼ 1=4). The parameter c in Eq. (91) was set as 5. 33 axial meshes were used in the calculation. We have normalized the distribution values by requiring the total neutron production rate to be 1 s1 (Lesage, 2001). The overview results are presented in Table 3. Small derivations of the ks were found, which were mainly caused by different convergence level of the inner iteration. When IFDF was employed, Table 3 Results for the small subcritical fast reactor problem. Item
SI
IFDF
IFDF-RE
CMFDRE
adCMFDRE
ks Number of fissionsource iteration Number of transport sweep Total CPU time/s Speedup Low-order calculation time/s
0.98729 575
0.98740 37
0.98738 37
0.98738 37
0.98738 48
29,850
621
618
585
844
705.9 1.0
23.9 29.5 5.6
20.0 35.3 1.5
19.1 37.0 1.5
29.0 24.3 0.8
Fig. 6. Radial core layout of the small subcritical fast reactor (Lesage, 2001).
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056 Table 4 Region average flux results for the small subcritical fast reactor problem. Group
SI
IFDF-RE
Derivation/%
Axial Blanket 1 2 3 4
2.221E05 2.778E05 1.581E05 1.140E05
2.221E05 2.778E05 1.582E05 1.141E05
0.00 0.01 0.03 0.12
Test Zone 1 2 3 4
1.558E04 1.163E04 2.677E05 2.405E06
1.557E04 1.163E04 2.676E05 2.406E06
0.06 0.04 0.02 0.02
Driver without moderator 1 1.042E04 2 7.915E05 3 2.530E05 4 5.220E06
1.042E04 7.914E05 2.530E05 5.221E06
0.01 0.01 0.01 0.02
Control rod 1 2 3 4
4.712E05 3.785E05 8.045E06 6.201E07
4.709E05 3.783E05 8.042E06 6.201E07
0.06 0.06 0.04 0.01
Control rod follower 1 1.160E04 2 9.800E05 3 3.597E05 4 1.049E05
1.160E04 9.801E05 3.598E05 1.049E05
0.01 0.01 0.02 0.03
9
the total calculation time was reduced from 705.9 s to 23.9 s. When the low-order problem was further accelerated by the rebalance acceleration method, the low-order calculation time of IFDF was further reduced from 5.6 s to 1.5 s of IFDF-RE which achieved a speedup of 35 compared with SI. The traditional CMFD acceleration with rebalance accelerated low-order calculation could also converge and obtained similar speedup with IFDF-RE, while adCMFD-RE gave significantly lower speedup than IFDF-RE. Table 4 compares the region average flux derivations between SI and IFDFRE. The derivations were insignificant, which further proved the effectiveness of the IFDF acceleration method for the externalsource multiplication problem. 3.3. The ZPPR-10B fast reactor problem This is a large fast reactor problem of the ZPPR-10, phase B (Lesage, 2001) (ZPPR-10B) fast critical experiments of Argonne National Laboratory. The radial flat distance of the square assembly is 5.5245 cm. The reactor core has a height of 209.56 cm. Fig. 7 gives the radial layout of the reactor core. 33-group homogenized cross sections for every assembly were obtained by the TULIP code in SARAX. With these cross sections, we have conducted the whole-core transport calculation with NECP-DNTC (DNTC). The core was divided by 24 axial meshes and 59 61 radial square meshes. We have normalized the distribution values by requiring the total power to be 1 MW.
Fig. 7. Radial core layout of the ZPPR-10B fast reactor (Lesage, 2001).
10
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
Table 5 Results for the ZPPR-10B problem.
Acknowledgements
Item
SI
IFDF
keff Number of fission-source iteration Number of transport sweep Total CPU time/s Speedup Low-order calculation time/s
0.99619 152 24,100 9377.0 1.0
0.99620 13 1993 898.3 10.4 89.4
This work was partially supported by the National Natural Science Foundation of China (11775170) and partially supported by the National Key R&D Program of China (2017YFE0302200). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.107056.
Table 6 Normalized neutron production rate results for the ZPPR-10B problem. Material region
Inner fuel 1 Inner fuel 2 Inner fuel 3 Outer fuel Axial blanket Radial blanket
Normalized region average adjoint flux DNTH-SI
DNTH-IFDF
Deviation/%
4.202E+10 3.911E+08 2.256E+08 4.223E+10 1.148E+09 1.518E+09
4.203E+10 3.915E+08 2.259E+08 4.222E+10 1.148E+09 1.517E+09
0.0 0.1 0.1 0.0 0.0 0.1
The overview results are presented in Table 5. For this problem, only SI and IFDF were considered. The parameter c in Eq. (91) was set as 5. SI gave a keff of 0.99619, while IFDF gave a close keff of 0.99620. The total calculation time of SI was 9377.0 s, while that of IFDF was only 898.3 s, which meant that a speedup of 10.4 was achieved when the low-order calculation overhead was 89.4 s. Table 6 compares the normalized neutron production rates of SI and IFDF. It is seen that the deviations were within 0.1%, which further proved the effectiveness of the IFDF acceleration method in this paper for cuboid meshes. 4. Conclusion Firstly, the SN nodal formulations in hexagonal prism and cuboid meshes have been introduced. Then, the improved interface flux discontinuity factor (IFDF) acceleration method was proposed. The detailed IFDF formulations were derived and applied to the acceleration of the SN nodal calculation for the eigenvalue problem and external-source multiplication problem. The IFDF acceleration method was improved by proposing a formulation to obtain the adaptive diffusion coefficients to ensure IFDF to be positive and keep it in a reasonable region, which could effectively improve the convergence stability. Finally, the theory has been implemented in the NECP-DNTH and NECP-DNTC code and three fast reactor problems were tested. For the large assembly MOX-3600 problem, the IFDF acceleration method obtained a speedup of 13, while the traditional CMFD or the artificial diffusion (g ¼ 1=4) CMFD (adCMFD) could not converge, which proved that the IFDF acceleration method was efficient and stable. For the small subcritical fast reactor problem, the rebalance accelerated IFDF acceleration achieved a speedup of 35, which was similar with that of the rebalance accelerated traditional CMFD and significantly higher than that of rebalance accelerated adCMFD. For the large ZPPR-10B problem with square assemblies, IFDF obtained a speedup of 10, which proved the effectiveness of the IFDF acceleration method for cuboid meshes. In summary, the proposed adaptive diffusion coefficients formulation with a proper damping factor (c ¼ 5) was valid for all the three tested cases and effectively improved the convergence stability for the MOX-3600 problem with large assembly size.
References Yang, W.S., 2012. Fast reactor physics and computational methods. Nucl. Eng. Technol. 44 (44), 177–198. Ruggieri JM, Tommasi J, Lebrat JF, et al. ERANOS 2.1: International code system for GEN IV fast reactor analysis[C]. Proceedings of ICAPP ’06, Reno, NV USA, June 48, 2006. Xu ZT, Wu HC, Zheng YQ. Development of an optimized transport solver in SARAX for fast reactor analysis[C]. In Proceedings of the 26th International Conference on Nuclear Engineering. London, UK, 2018. Willert, J., Park, H., Knoll, D.A., 2014. A comparison of acceleration methods for solving the neutron transport k-eigenvalue problem. J. Comput. Phys. 274, 681– 694. Smith KS. Full-core, 2-D, LWR core calculations with CASMO-4E[C]. PHYSOR 2002, Seoul, Korea, 2002. Cho, N.Z., 2017. An overview of p-CMFD acceleration and its applications to reactor physics transport calculation. Trans. Am. Nucl. Soc. 117, 1247–1250. Jarrett, M., Kochunas, B., Zhu, A., et al., 2016. Analysis of stabilization techniques for CMFD acceleration of neutron transport problems. Nucl. Sci. Eng. 184 (2), 208– 227. Zhu, A., Jarrett, M., Xu, Y., et al., 2016. An optimally diffusive coarse mesh finite difference method to accelerate neutron transport calculations. Ann. Nucl. Energy 95, 116–124. Yuk, S., Cho, N.Z., 2017. Two-level convergence speedup schemes for p-CMFD acceleration in neutron transport calculation. Nucl. Sci. Eng. 188 (1), 1–14. Xiao, S., Ren, K., Wang, D., 2018. A local adaptive coarse-mesh nonlinear diffusion acceleration scheme for neutron transport calculations. Nucl. Sci. Eng. 189 (3), 272–281. Yamamoto, A., Endo, T., Giho, A., 2018. Transport Consistent Diffusion Coefficient for CMFD Acceleration. Trans. Am. Nucl. Soc. 119, 1179–1181. Masiello, E., 2018. ‘‘On the fly” stabilization of the coarse-mesh finite difference acceleration for multidimensional discrete-ordinates transport calculations. J. Comput. Phys. 373, 1–27. Masiello, E., Sanchez, R., Zmijarevic, I., 2009. New numerical solution with the method of short characteristics for 2-d heterogeneous cartesian cells in the apollo2 code: Numerical analysis and tests. Nucl. Sci. Eng. 161 (3), 257–278. Xu YL, Downar T. Convergence analysis of a CMFD method based on generalized equivalence theory[C]. PHYSOR 2012, Knoxville, Tennessee, USA, 2012. Lu, H.L., Wu, H.C., 2007. A nodal SN transport method for three-dimensional triangular-z geometry. Nucl. Eng. Des. 237 (8), 830–839. Xu, Z.T., Wu, H.C., Zheng, Y.Q., 2019. An improved hexagonal-z nodal SN method in SARAX code system. Prog. Nucl. Energy. submitted for publication. Zheng, Y.Q., Du, X.N., Xu, Z.T., et al., 2018. SARAX: a new code for fast reactor analysis part I: methods. Nucl. Eng. Des., 421–430 Zheng, Y.Q., Qiao, L., Zhai, Z.A., et al., 2018. SARAX: A new code for fast reactor analysis part II: verification, validation and uncertainty quantification. Nucl. Eng. Des. 331, 41–53. Aragones, J.M., Ahnert, C., 1986. A linear discontinuous finite difference formulation for synthetic coarse-mesh few-group diffusion calculations. Nucl. Sci. Eng. 94 (4), 309–322. Badruzzaman, A., 1985. An efficient algorithm for nodal-transport solutions in multidimensional geometry. Nucl. Sci. Eng. 89 (3), 281–290. Azmy, Y.Y., 1997. Multiprocessing for neutron diffusion and deterministic transport methods. Prog. Nucl. Energy 31 (3), 317–368. Zhou SC, Wu HC, Cao LZ, et al. Comparison of three acceleration techniques for the fixed source neutron transport calculation[C]. In Proceedings of the PHYSOR 2018, Cancun, Mexico, 2018. Xu, Z.T., Wu, H.C., Zheng, Y.Q., 2017. Hexagonal CMFD accelerated discrete nodal transport method in triangular-z geometry. Trans. Am. Nucl. Soc. 117, 718–720. Bernnat W, Blanchet D, Brun E, et al. Benchmark for neutronic analysis of sodiumcooled fast reactor cores with various fuel types and core sizes[R]. Nuclear Energy Agency Technical Report, NEA/NSC, 2015. Takeda T, Ikeda H. 3-D neutron transport benchmarks[R]. NEACRP-L-330, OECD/ NEA Committee on Reactor Physics, 1991. Lesage LG. An overview of the Argonne National Laboratory fast critical experiments [R]. ANL-NT-175, Argonne National Laboratory, 2001.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
IFDF acceleration method with adaptive diffusion coefficients for SN nodal calculation in SARAX code system Zhitao Xu, Youqi Zheng, Yongping Wang ⇑, Hongchun Wu School of Nuclear Science and Technology, Xi’an Jiaotong University, 28 West Xianning Road, Xi’an, Shaanxi 710049, China
a r t i c l e
i n f o
Article history: Received 10 May 2019 Received in revised form 9 September 2019 Accepted 11 September 2019 Available online 19 September 2019 Keywords: IFDF CMFD SN nodal method Fast reactor
a b s t r a c t SN nodal method has been successfully applied to fast reactor core neutronics analysis. The computational efficiency can be enhanced efficiently by the popular coarse-mesh finite difference (CMFD) method. However, the CMFD method suffers from the problem of instability, especially for optically thick problems. In this paper, another efficient nonlinear acceleration method, named interface flux discontinuity factor (IFDF) acceleration, was studied and improved for the SN nodal calculation in hexagonal prism and cuboid meshes. With adaptive diffusion coefficients, the IFDF is guaranteed to be positive and within a reasonable region, which makes the acceleration stable. It was found that the IFDF and CMFD methods obtained similar speedups for the tested fast reactor problems. However, the stability of the IFDF method was much better than the CMFD method even with artificial diffusion improvement. Additionally, the rebalance accelerated IFDF method was applied to the external-source multiplication problem besides the common eigenvalue problem. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Fast reactor core calculation is critical for fast reactor nuclear design analysis. Due to the complex energy spectrum (Yang, 2012) and the strong neutronics anisotropy in fast reactor core, multigroup (for example, 33 energy groups (Ruggieri et al., 2006) transport calculation is essential. However, the computational cost of multigroup transport increases by about three orders of magnitude compared with the traditional few-group diffusion calculation. Thus, it is necessary to study efficient numerical methods and acceleration techniques to obtain higher efficiency. The existing researches show that the discrete ordinates (SN) nodal method well balances the calculation accuracy, stability, memory and efficiency (Xu et al., 2018). Meanwhile, in fast reactors, the mean free path of neutrons is comparable to the size of one assembly and homogeneous cross sections are usually used in the core calculation, so it is not necessary to subdivide the mesh within the hexagonal or rectangular assembly. What’s more, the ray-effect of the SN method are negligible for the reactor core calculation because the fission and scattering sources are distributed all over the core. Therefore, it is of great significance and application value to study efficient SN nodal method in hexagonal prism and cuboid meshes.
⇑ Corresponding author. E-mail address: [email protected] (Y. Wang). https://doi.org/10.1016/j.anucene.2019.107056 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
Nonlinear acceleration methods which employ the strategy of nonlinear iteration between accurate high-order problem and efficient low-order problem are one kind of efficient acceleration method and they usually have better performance than the Krylov methods (Willert et al., 2014). The CMFD acceleration method (Full-core, 2002) is one of the most widely applied nonlinear methods employed in the method of characteristics (MOC) or SN calculation. It can provide a speedup from 10 to 50 for most cases. However, for some cases, especially for optically thick problems, CMFD could cause iterative divergence due to the inconsistency between the high-order problem and the low-order problem. Many innovative schemes for improving the stability and acceleration performance of the CMFD acceleration have been proposed (Cho, 2017; Jarrett et al., 2016; Zhu et al., 2016; Yuk and Cho, 2017; Xiao et al., 2018; Yamamoto et al., 2018; Masiello, 2018; Masiello et al., 2009; Xu et al., 2012) in recent years. These methods focus on improving the consistency of high-order problem and low-order problem by modifying coupling formulation or coupling meshes. The CMFD acceleration with stability techniques including the artificial diffusion and multiple transport sweeps mentioned in reference (Jarrett et al., 2016) has been implemented in the SN nodal methods (Xu et al., 2018; Lu and Wu, 2007; Xu et al., 2019) for hexagonal-z or triangular-z geometries in SARAX (Zheng et al., 2018; Zheng et al., 2018) code system. However, for a large MOX-3600 fast reactor with large assembly size of 21 cm, the CMFD acceleration with the artificial diffusion and multiple
2
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
transport sweeps techniques for the SN nodal calculation with hexagonal prism meshes still failed. Therefore, another nonlinear acceleration method (Aragones and Ahnert, 1986) employing the interface flux discontinuity factor (IFDF) was studied in this paper. For convenience, it is called IFDF acceleration. IFDF acceleration has been studied for short characteristics transport calculation (Badruzzaman, 1985) but the benefits was not clear, because the original reference (Aragones and Ahnert, 1986) did not clearly illustrated how to control the transport correction factor. In this paper, the IFDF acceleration is modified by giving a detailed formulation of adaptive diffusion coefficients to adjust the transport correction factor for the stability of acceleration and it is implemented in the SN nodal calculation with hexagonal prism and cuboid meshes. The paper is organized as follows. In Section 2, the introduction of the SN nodal formulations, the IFDF formulations, the adaptive diffusion coefficients, and the iteration schemes for both the eigenvalue and the external-source multiplication problems are presented. In Section 3, numerical results for a MOX-3600 fast reactor, a subcritical small fast reactor with external-source, and a ZPPR-10B fast reactor are given and discussed. Finally, Section 4 gives the conclusions of this paper.
X
lmd @
wm ðx; y; zÞ þ Rt wm ðx; y; zÞ ¼ Q m ðx; y; zÞ
d m m m y ðxÞwm x ðxÞ þ px ys ðxÞwx ðxÞ ¼ qx ðxÞ dx s
ð2Þ
where
pffiffiffi ys ðxÞ ¼ ð1 jxjÞ= 3 1 ¼ 2ys ðxÞ
wm x ð xÞ
pm x ¼
Zze
ð3Þ
yZ s ðx Þ
wm ðx; y; zÞdydz
ð5Þ
lmx
hr
m m ys ðxÞ Q m x ðxÞ Lax ðxÞ Lr x ðxÞ
qm x ðxÞ ¼
lmx
Qm x ðxÞ ¼
1 2ys ðxÞ
Lam x ð xÞ ¼
wm z ðxÞ
ð4Þ
zb ys ðxÞ
Rt hr
2.1. Introduction of the SN nodal formulations To conduct the multigroup neutronics transport calculation, the spatial domain is discretized into hexagonal prism or cuboid nodes (meshes). With the SN discretization for angle, nodal response relations which couple the nodal outgoing and incoming surface fluxes, along with the nodal interior fluxes and source, can be obtained for every node. The nodal response relations are employed for transport sweeps to get the final solution. The SN nodal response relations can be written as the SN nodal formulations which are introduced as follows. First let’s consider the hexagonal prism meshes. For a hexagonal node with the flat distance hd (unit: cm) along d (d ¼ v ; x; u; z) channel, the discrete neutron transport equation within one node can be written in a standard hexagonal prism illustrated in Fig. 1 as
ð1Þ
where lm d is direction cosine relative to d coordinate axis direction (channel), hd is the actual side-length (unit: cm) of the node along d channel, hv ¼ hx ¼ hy ¼ hr is satisfied and hr is the radial flat distance of the node, wm is the neutron flux (unit: cm2s1) for discrete angle m, Rt is the total cross section (unit: cm1) of this node, Q m is the neutron source (unit: cm3s1) including scattering source, fission source and external source. With the geometric symmetry of the node and the selected SN angular quadrature, we only need to consider the case when lmd > 0. The transverse integral technique is employed to decouple the 3-D equation, take x direction as an example,
2. Theory The introduction of the SN nodal formulations with hexagonal prism and cuboid meshes is given in subsection 2.1, which is the foundation of this work. In subsection 2.2., the IFDF formulations are derived in detail. The equations of the transport correction factor will be given, and the iteration schemes for both the eigenvalue and the external-source multiplication problems will be presented. Then, the formulations of adaptive diffusion coefficients are derived in subsection 2.3.
hd @d
d¼v ;x;u;z
lmz hz
Zze
yZ s ðx Þ
Q m ðx; y; zÞdydz
ð7Þ
zb ys ðxÞ
m wm zþ ðxÞ wz ðxÞ
1 ¼ 2ys ðxÞ
ð6Þ
ð8Þ
yZ s ðx Þ
wm ðx; y; 1=2Þdy
ð9Þ
ys ðxÞ
sgnðxÞ m m lu wusgnðxÞ x ðxÞ þ lmv wmv sgnðxÞ x ðxÞ Lrm x ðxÞ ¼ pffiffiffi 3hr
ð10Þ
Zze wm u x ðxÞ
¼
w½x; ys ðxÞ; zdz
ð11Þ
wm ½x; ys ðxÞ; zdz
ð12Þ
zb
Zze wm v x ðxÞ ¼ zb
pffiffiffi
l ¼ m u
lmx þ 3lmy 2
ð13Þ
pffiffiffi
lv ¼ m
Fig. 1. Illustration of a standard hexagonal prism.
lmx 3lmy 2
ð14Þ
For other directions, similar transverse integral process can be conducted. We expand the 1-D nodal interior flux and source with the orthogonal polynomials. Take x direction as an example,
3
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056 3 X
wm x ð xÞ ¼
wm xi hi ðxÞ
ð15Þ
Qm xi hi ðxÞ
ð16Þ
i¼0
Qm x ð xÞ ¼
3 X i¼0
where hi ðxÞ is a selected orthogonal polynomials for radial channels. And we expand the 1-D nodal surface fluxes as: m m wm z ðxÞ ¼ wz þ wz x1 h1 ðxÞ
m where em z and F z0 are also nodal integral coefficients which can be calculated before the transport sweeps. These additional highorder nodal surface flux formulations are employed to update the m m transverse terms Lam d1 , Lzd1 and Lz1 . Then let’s consider the cuboid meshes. For a cuboid node with the flat distance hs (cm) along d (d ¼ x; y; z) direction, the discrete neutron transport equation of one node (shown as Fig. 2) can be written as
ð17Þ
X
d ¼ v; u
ð18Þ
d¼x;y;z
ðzÞ ¼ wm þ wm f ðzÞ; d ¼ v ; x; u wm d z d d z1 1
ð19Þ
After similar transverse integral process and applying polynomial expansion as the z direction of the hexagonal prism case, the nodal formulations for cuboid meshes can be obtained. Take the z direction as an example,
ð xÞ wm d
¼
wm ; d
where f 1 ðzÞ is a selected orthogonal polynomials for z direction. Through some complex mathematical derivations, we can obtain the formulations coupling the outgoing nodal surface flux m , with the incoming nodal surface wm and the nodal average flux w dþ
m m flux wm ¼ Qm d0 ): d and the nodal interior source Q dj (Q
m m m Am þ wþ ¼ A w þ d
wm zþ
m
ð20Þ
m Am þ wm 1 Bm þ Dm ¼w z z z z
ð21Þ
lmd
¼ w m
@ m w ðx; y; zÞ þ Rt wm ðx; y; zÞ ¼ Q m ðx; y; zÞ @d
m þ P wm Bm Dm lm =hd Q d d d d d¼x;y;z
Rt þ
P d¼x;y;z
m Am d ld =hd
m m m m wm þ Dm zþ ¼ w Az þ wz 1 Bz z
where
wm þ ¼ wm ¼
wm vþ
wm xþ
wm uþ
m T w
wm v
wm x
wm u
wm z
T
3 3 P m m m Qm v j C v j Lav 1 C v 1 7 6 j¼1 7 6 7 6 3 7 6 P m m m m 7 6 Q C La C x1 x1 xj xj 7 6 j¼1 7 6 m 7 d ¼6 7 6 3 P m m m 7 6 m Q C La C 7 6 u1 u1 uj uj 7 6 j¼1 7 6 ! 7 6 3 5 4 m lmz P m m m m Q hz Q zj C zj Lz1 C z1
ð22Þ ð23Þ
2
ð24Þ
Dm z ¼
3
X m m Qm zj Lzj W zij
m m Qm zj Lzj C zj
ð25Þ
j¼1 m m m m and Am þ , A , Ad , Bd and C d are nodal integral coefficients which can be calculated before the transport sweeps. Then, we can obtain the formulations coupling the nodal interior flux wm di with the m incoming nodal surface flux wm d and nodal interior source Q dj ,
ð32Þ
ð33Þ
j¼1
Different from the hexagonal prism nodal formulations, nonlinear terms (Xu et al., 2019) are not appeared in these formulations. Therefore these formulations obey the standard nodal-equivalent finite difference (NEFD) (Azmy, 1997) scheme. However, for the additional high-order nodal surface flux formulations, this paper extends the simplified method to cuboid mesh to update the high-order transverse leakage term Lm d1 ,
ð34Þ
3
X m m m m m m Sm xj Laxj W xij Lnx0 Wni0 Lpx0 Wpi0
With the nodal response relations presented in Section 2.1, the transport solutions can be obtained by iterations with assumed initial solution. Since the critical characteristic is important for reactor physics analysis, the eigenvalue (k-eigenvalue) problem is a common topic. While for the analysis of reactor neutron dynamics, subcritical reactor neutronics and uncertainty quantification, the external-source multiplication problem is encountered. To unify the two kinds of problem, the steady-state neutron transport equation can be formulated as follows,
j¼1
ð26Þ m m m m wm zi ¼ wz0 U zi þ wz V zi þ
ð31Þ
2.2. IFDF formulations
3 X
m m m m wm xi ¼ wx0 U xi þ wx V xi þ
m m m m wm zi ¼ wz0 U zi þ wz V zi þ
m m m m wm x1zþ wx1z ez þ Q x1 F z0
j¼1
ð30Þ
3
X m m Qm zj Lzj W zij
ð27Þ
j¼1 m m m m where U m di , V di , W di , Wndi , Wpdi are also nodal integral coefficients which can be calculated before the transport sweeps. Additional high-order nodal surface flux formulations are obtained by a simplified method (Xu et al., 2019):
m m m m wm x1zþ wx1z ez þ Q x1 F z0
ð28Þ
m m wm xþ z1 Q z1 W z11
ð29Þ
Fig. 2. Illustration of a standard cuboid.
4
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
1 Lw ¼ SU þ FU þ Q ex k
ð35Þ
where L is an operator including neutron leakage and total reaction rate; w is the angular flux vector; S is the scattering operator; U is the scalar flux vector. k equals to the effective multiplication factor keff for the eigenvalue problem and it equals to constant 1.0 for the external-source multiplication problem. F is the fission operator and Q ex is the external-source vector which only exists for the external-source multiplication problem. For external-source multiplication problem, to evaluate the multiplication performance of the system, a source multiplication factor ks is defined as follows,
ks ¼
hFUi hFUi hQ ex i
ð36Þ
Lw
ðn1Þ
¼ SU
þ
1 k
ðn1Þ
ðn1Þ
FU
þ Q ex
Uðn1=2Þ ¼ Twðn1=2Þ D ðn1=2Þ
k
ðn1Þ
¼k
ð37Þ ð38Þ
E
FUðn1=2Þ D E FUðn1Þ
ð39Þ
where T is the operator to transfer the angular flux vector to the scalar flux vector, n is the iteration index, ðn 1=2Þ means the value has just been updated by the high-order (SN nodal) operators from the ðn 1Þ-th value, and Eq. (39) is only valid for the eigenvalue problem. However, the source iteration converges very slowly, especially for problems with large dominance ratio. For large fast reactor, the dominance ratio is close to 1.0 and the number of meshes is also very large. Therefore, an efficient acceleration method is important. Here we consider the IFDF acceleration method. The IFDF acceleration employs the nonlinear iteration between the highorder (SN nodal) problem and the low-order (diffusion finitedifference) problem. Since hexagonal meshes cannot be subdivided or combined, we use the same meshes for the low-order problem as the high-order meshes (nodes). In addition, the high-order and the low-order problem utilize the same number of energy groups to avoid additional inconsistency caused by cross sections homogenization. The low-order problem to be solved can be written as below,
ML U00 ¼
1 FL U00 þ Q L;ex k
ð40Þ
where ML is an operator including diffusion, absorption and scattering; U00 is the nodal average scalar flux vector of all groups; F is the fission operator for U00 ; Q L;ex is the external-source vector for low-order problem. Eq. (40) is much less expensive to be solved compared with Eq. (35), but the accuracy is lower than that of Eq. (35). A proper nonlinear iteration between the high-order problem and the low-order problem can combine their advantages. The nonlinear iteration scheme for the IFDF acceleration is similar with that of the CMFD acceleration, while the distinct difference is the formulations of the low-order problem, which are presented as follows. For hexagonal prism meshes, the diffusion finite difference formulation for a single group can be written as below,
X J K d þ1=2 J K d 1=2 AK d þ Rr;K /K V K ¼ Q K V K
d¼v ;x;u;z
X J K d þ1=2 J K d 1=2 AK d þ Rr;K /K V K ¼ Q K V K
ð42Þ
d¼x;y;z
where h i is the operator for integral over the full phase space. The final solutions can be obtained by the source iteration (SI). In this work, a outer fission-source iteration and an inner scattering-source iteration for the nodal average flux of every group are employed. Let’s focus on the outer iteration, which can be expressed as below, ðn1=2Þ
where K is the nodal index; JK d 1=2 is the neutron current (unit: cm2s1) in the normal direction d along mesh surface K d 1=2; K d means the value is related with the direction d of mesh K; K d 1=2 means the value is related with the mesh surface; AK d is the mesh surface area (unit: cm2) along direction d of mesh K; Rr;K is the removal cross section (unit: cm1) of mesh K; /K is the scalar flux (unit: cm2s1) of mesh K; V K is the volume (unit: cm3) of mesh K; Q K is the neutron source (unit: cm3s1) of mesh K, which includes scattering source from other groups, fission source and external source. For cuboid meshes, the diffusion finite difference formulation is similar,
ð41Þ
The Fick’s law gives the relation between neutron current and neutron flux. For an inner mesh K, the Fick’s law can be written as below,
J K d þ1=2 ¼
DK d /K d þ1=2 /K d
ð43Þ
hK d =2
where hK d is the mesh size along direction d of mesh K. However, the Fick’s law gives large error for transport problem. A discontinuK þ1=2
ity factor f K dd
J K d þ1=2 ¼
is employed:
K þ1=2 DK d f K dd /K d þ1=2 /K d hK d =2
ð44Þ
For mesh K þ 1,
J K d þ1=2 ¼
K þ1=2 DK d þ1 /K d þ1 f K dd þ1 /K d þ1=2 hK d þ1 =2
ð45Þ
Combine Eqs. (44) and (45), we obtain
J K d þ1=2 ¼
DK d þ1=2 rK d þ1=2 /K d þ1 /K d hK d þ1=2
ð46Þ
where
DK d þ1=2 ¼
2DK d DK d þ1 hK d þ1=2 DK d r K d þ1=2 hK d þ1 þ DK d þ1 hK d hK
/K d J K d þ1=2 2DKd
K þ1=2
rK d þ1=2 ¼
hK d þ1=2 ¼
f K dd
K þ1=2
f K dd þ1
ð47Þ
¼
d
hK
þ1
/K d þ1 þ J K d þ1=2 2DKd þ1
ð48Þ
d
hK d þ hK d þ1 2
ð49Þ
The ratio rK d þ1=2 of two discontinuity factors are called as IFDF by reference (Aragones and Ahnert, 1986), and it is also called as transport correction factor in this paper. In the nonlinear acceleration process, the transport correction factor is updated with Eq. (48) by substituting the flux and current with the corresponding values from high-order calculation, hK
rK d þ1=2 ¼
SN d /SN K d J K d þ1=2 2DK
d
hK
þ1
SN d /SN K d þ1 þ J K d þ1=2 2DK þ1
ð50Þ
d
where the superscript SN means the value is obtained from the high-order SN nodal calculation. For the case that mesh K exists and mesh K d þ 1 does not exists, the modified Fick’s law is written as below,
5
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
J þK d þ1=2
/K þ1=2 DK d ¼ d 4 2 /K þ1=2 DK d ¼ d þ 4 2
K þ1=2 f K dd /K d þ1=2
¼
K 1
ð51Þ
/K d
Jþ K d þ1=2
and
J K d þ1=2
is expressed as below,
bK d þ1=2 J þK d þ1=2
ð53Þ ð54Þ
where bK d þ1=2 is the boundary albedo of mesh surface K d þ 1=2. Combine Eqs. (51)–(54), and substitute the flux and current with the corresponding values from high-order calculation, we obtain
DK d þ1=2 /K d hK d þ1=2
ð55Þ
where
ð56Þ
hK d aK d þ1=2 þ 2DK d r K dd
K þ1=2
hK
¼
K þ1=2 f K dd
¼
SN d /SN K d J K d þ1=2 2DK
ð57Þ
d
JSN K þ1=2 d
aK d þ1=2
1 bK d þ1=2
¼ 2 1 þ bK d þ1=2
aK d þ1=2
ð58Þ
hK d þ1=2 ¼ hK d
ð59Þ
For the case that mesh K exists and mesh K d 1 does not exists, similar formulations are obtained:
J K d 1=2 ¼
DK d 1=2 rK d 1=2 /K d hK d 1=2
ð60Þ
where
dK ¼ Rr;K þ
DK d 1=2
K þ1
¼ cK dd
K 1
¼ cK dd
dK dd
K 1=2
¼
dK dd
K 1=2
aK d 1=2
1
¼
þ
JSN K 1=2 hK d
2DK
ð62Þ d
d
1 bK d 1=2
¼ 2 1 þ bK d 1=2
ð63Þ
ð64Þ
Combing Eq. (41), (46), (55), and (60), for hexagonal prism meshes we obtain
X
d¼v ;x;u;z
where
dK ¼ Rr;K þ K þ1
dK dd
K þ1
/K d þ1 dK dd
X
d¼v ;x;u;z
K 1
/K d 1 dK dd
¼ QK
X K þ1=2 K 1=2 cK dd DK d þ1=2 þ cK dd DK d 1=2 r K d 1=2
d¼x;u;v ;z
K þ1=2
¼ cK dd
DK d þ1=2 r K d þ1=2
ð65Þ
J SN K d þ1=2 ¼
ð67Þ
ð71Þ ð72Þ
DK d 1=2
ð73Þ
1 hK d hK d 1=2
ð74Þ
X
es;K d þ1=2
s¼x;y;z
X m
lms xm wmKd þ1=2
ð75Þ
where es;K d þ1=2 is the projection value of the normal vector of mesh surface K d þ 1=2 at s direction, and xm is the quadrature weight of the discrete angle m. With the solution of Eqs. (65) or (70) equation for every group, the solution of Eq. (40) can be obtained by fission-source iteration. The high-order values obtained by Eqs. (38) and (39) are employed as the initial values of the low-order fission-source iteration, which is expressed as below,
Uð00n1;lÞ ¼
1 k
ðn1;l1Þ
¼k
ðn1;l1Þ
ðn1;l1Þ
D D
FL U00
ðn1;lÞ
FL U00
ð76Þ
E
ðn1;l1Þ
FL U00
þ Q L;ex
E
ð77Þ
where the low-order iteration index l is employed, and Eq. (77) is only valid for the eigenvalue problem. For the external-source multiplication problem, the fission-source iteration converges slowly, the rebalance iteration (Xu et al., 2017) is employed to further accelerate the low-order fission-source iteration. The rebalance iteration employs a global rebalance factor f by assuming that the spatial shape of the flux does not change quickly during iterations. The original iteration (without eigenvalue) is expressed as: ðn1Þ
ML
Uð00n1;l1=2Þ ¼ FL Uð00n1;l1Þ þ Q L;ex
ð78Þ
The rebalance method introduces a factor to enforce the updated flux to satisfy the neutron balance equation:
f ð66Þ
ð70Þ
The linear system constituted by equation or equation can be solved by the red–black algorithm (Zhou et al., 2018). After the high-order fission-source iteration as expressed by Eqs. (37)–(39), the updated high-order flux and current are used to calculate IFDF with Eqs. (50), (57), and (62). The current is calculated with the angular flux. Take the current at mesh surface K d þ 1=2 as an example,
ðn1;lÞ
aK d 1=2
/SN Kd
¼ QK
DK d þ1=2 rK d þ1=2
ð61Þ
hK d 1=2 ¼ hK d
/K dK
¼
d
K 1=2
K 1
/K d 1 dK dd
X K þ1=2 K 1=2 cK dd DK d þ1=2 þ cK dd DK d 1=2 r K d 1=2
J SN K 1=2
f K dd
X
d¼x;y;z
K þ1=2
k r K dd
ð69Þ
d¼x;y;z
ML
2DK d aK d 1=2 hK d 1=2 ¼ hK d aK d 1=2 r K d 1=2 þ 2DK d
K þ1
/K d þ1 dK dd
d¼x;y;z
ðn1Þ
where
2=ð3hr hr Þ; d ¼ v ; x; u 1= hK d hK d 1=2 ; d ¼ z
X
/K dK
cK dd
2DK d aK d þ1=2 hK d þ1=2
¼
ð68Þ
DK d 1=2
For cuboid meshes we obtain
K 1=2
DK d þ1=2 ¼
K þ1=2 r K dd
K 1=2
¼ cK dd
(
ð52Þ
J K d þ1=2 ¼ J þK d þ1=2 J K d þ1=2
J K d þ1=2 ¼
dK dd
K 1=2 cK dd
hK d =2
The relation between
J K d þ1=2
/K d þ1=2 /K d hK d =2
J K d þ1=2
K þ1=2
f K dd
ðn1;lÞ
D
ðn1Þ
ML
E
D
E
Uð00n1;l1=2Þ ¼ f ðn1;lÞ FL Uð00n1;l1=2Þ þ Q L;ex
ð79Þ
Thus, combine Eqs. (78) and (79) to obtain the expression of the factor:
f
ðn1;lÞ
¼
Q L;ex
E
D E D ðn1;l1=2Þ ðn1;l1Þ Q L;ex FL U00 þ FL U00
ð80Þ
6
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
Then the flux is modified by multiplying the factor: ðn1;lÞ 00
U
¼f
ðn1;lÞ
ðn1;l1=2Þ 00
U
ð81Þ
Through above process, the fission-source iteration can be well accelerated. After the low-order problem is converged, we update the highorder values as follows, ðnÞ
ðn1=2Þ
/i;K ¼ /i;K ðn Þ
k
/DIF K ðn1=2Þ /0;K
;i P 0
ð82Þ
DIF
¼k
ð83Þ
where the superscript DIF means the value is from the converged low-order eigenvalue problem, i is the spatial moment index of the nodal polynomial expansion, and Eq. (83) is only valid for the eigenvalue problem. The overall nonlinear iteration scheme of the IFDF iteration is illustrated in the flowchart of Fig. 3. The loworder problem is calculated for certain high-order fission-source iterations to avoid unnecessary low-order calculation efforts and improve acceleration stability. The diffusion coefficients are updated to ensure the IFDF is positive but not too large, which will be explained in subection 2.3. 2.3. Adaptive diffusion coefficients
J SN K d þ1=2
> 0, the fol-
lowing condition must be satisfied,
hK d
ð84Þ
2/SN Kd
If J SN K d þ1=2 < 0, the following condition must be satisfied,
DK d þ1 > J SN K d þ1=2
hK d þ1
ð85Þ
2/SN K d þ1
and the following condition must also be satisfied for adjacent mesh K þ 1,
DK d > J SN K d 1=2
hK d
ð86Þ
2/SN Kd
For the case that mesh K exists and mesh K d þ 1 does not exist, IFDF is defined by Eq. (57). For the reflective boundary condition, DK d þ1=2 equals to zero and IFDF makes no sense. Thus, we only need to consider the vacuum boundary condition, which means the condition that J SN K d þ1=2 > 0 is always satisfied. The same condition of Eq. (84) must be satisfied. For the case that mesh K exists and mesh K d 1 does not exist, IFDF is defined by Eq. (57). For the same reason as above, the condition that J SN K d 1=2 < 0 is always satisfied. The same condition as Eq. (86) must be satisfied. Combine Eqs. (84) and (86), the condition is written as below,
DK d >
As IFDF is defined as the ratio of the discontinuity factors along one mesh surface, it should always be positive to be physically significant. In addition, it should stay within a reasonable region (positive but not too large). This subsection will give the condition to bound IFDF or the transport correction factor. First let’s find the condition that makes IFDF to be positive. For an inner mesh K, IFDF is defined by Eq. (50). If
DK d > J SN K d þ1=2
1 2/SN Kd
SN max hK d J SN K d 1=2 ; hK d J K d þ1=2
ð87Þ
Eq. (87) must be satisfied for every channels to ensure that all the IFDFs around a mesh are positive, so the condition is written as below,
DK > D0K
ð88Þ
for hexagonal meshes,
D0K ¼
1 2/SN K
SN SN SN max hK v J SN K v 1=2 ; hK v J K v þ1=2 ; hK x J K x 1=2 ; hK x J K x þ1=2 ;
SN SN SN hK u J SN K u 1=2 ; hK u J K u þ1=2 ; hK z J K z 1=2 ; hK z J K z þ1=2
ð89Þ
for cuboid meshes,
D0K ¼
1 2/SN K
SN SN SN max hK x J SN K x 1=2 ; hK x J K x þ1=2 ; hK y J K y 1=2 ; hK y J K y þ1=2 ;
SN hK z J SN K z 1=2 ; hK z J K z þ1=2
ð90Þ
To ensure IFDF is positive but not too large, this paper gives the following formulation for adaptive determination of the diffusion coefficients of every mesh before update IFDF,
1 DK ¼ max ; cD0K 3Rt;K
ð91Þ
where Rt;K is the total cross section of mesh K, c is a damping parameter obtained by experience. The damping parameter is used to keep r in a reasonable region (as is shown by Eq. (94)) and to adjust the ‘step size’ of the iteration towards convergence. It can be regarded as an adjuster between the speedup and stability of this acceleration method. Here the diffusion coefficient is selected for the adaptive value because it can help IFDF improve the interface current consistency between the high-order problem and the loworder problem without affecting the reaction rates. To understand the effect of the damping parameter c, let’s consider a simple model that all the mesh surface currents equal to jJ j and jJ j is large enough that cD0K > 3R1t;K and all the mesh sizes Fig. 3. Nonlinear iteration scheme of the IFDF acceleration.
equal to h, the following equations can be obtained,
7
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
DK ¼ cD0K ¼ c
hjJ j
ð92Þ
2/SN K
DKþ1 ¼ cD0K ¼ c
hjJ j
ð93Þ
2/SN Kþ1
Substitute Eqs. (92) and (93) into (50), we obtain
r K d þ1=2 ¼
/SN K /SN Kþ1
1þ
2 c1
ð94Þ
From Eq. (94) it’s seen that if c near 1, rK d þ1=2 will be very large and even exceed the computer floating number limit; if c is very SN large, rK d þ1=2 will be bounded by /SN K =/Kþ1 . For this work we found that if c is very large, the acceleration speedup is low, so a proper value c ¼ 5 is employed after many numerical tests. From the numerical tests below it will be seen that c ¼ 5 is conservative to guarantee iteration stability. Smaller c is recommended if one expects higher speedup and larger c is recommended if one expects higher iteration stability.
3. Results The IFDF acceleration method described above has been implemented in the SARAX-DNTH (Discrete Nodal Transport Code for Hexagonal-Z Meshes) code and SARAX-DNTC (Discrete Nodal Transport Code for Cuboid Meshes) code in NECP-SARAX (SARAX) code system. In this section, a MOX-3600 fast reactor with large assemblies, a small subcritical fast reactor and a ZPPR-10B large fast reactor with square assemblies were calculated to evaluate the modified IFDF acceleration method.
Fig. 4. Radial core layout of 3600 MWth oxide core (Takeda and Ikeda, 1991).
For all the calculations, S4 Legendre-Chebyshev angular quadrature was used; the nodal interior variables were expanded with 2nd-order polynomials; the nodal surface variables were expanded with 1st-order polynomials; the fission-source iteration criterion of the high-order problem was 1 105 ; the scattering-source iteration criterion of the high-order problem was 5 106 ; an optimized iteration limit of 16 (Bernnat et al., 2015) was employed if only source iteration was carried out (no acceleration); an iteration limit of 5 was employed if nonlinear acceleration was applied; the fissionsource iteration criterion of the low-order problem was 2:5 106 ; one low-order problem was solved before four high-order fissionsource iterations were performed. All the following calculations were performed on a 3.4 GHz Intel Xeon Gold 6134 CPU core. 3.1. The MOX-3600 fast reactor problem This is a benchmark problem of fast reactor (Takeda and Ikeda, 1991) with large assembly size. The radial flat distance of the hexagonal assembly is 21.2205 cm and the height of the reactor core is 301.6974 cm. This core is consisted of 225 inner fuel assemblies, 228 outer fuel assemblies, 331 radial reflector assemblies and 33 control assemblies. The radial layout of the core is illustrated in Fig. 4. 33-group homogenized cross sections for every assembly were obtained by the TULIP code in SARAX. With these cross sections, we have conducted the whole-core adjoint transport calculation with SARAX-DNTH (DNTH). The core was divided by 21 axial meshes and 817 radial hexagonal meshes. The overview calculation results are presented in Table 1, in which SI means that only source iteration was employed (no acceleration); CMFD/adCMFD means that the traditional CMFD or artificial diffusion (g ¼ 1=4) CMFD (Jarrett et al., 2016) acceleration was employed; IFDF (c = 0) means the IFDF acceleration without adaptive diffusion coefficients was employed; IFDF (c = 1.001), IFDF (c = 1.1), IFDF (c = 5) and IFDF (c = 1000) mean that the improved IFDF acceleration with different c in Eq. (91) was employed. In addition, one transport sweep means conducting flux calculation for all angles and meshes for one group. As shown in Table 1, CMFD/ adCMFD, IFDF (c = 0) and IFDF (c = 1.001) could not converge for this problem, while IFDF (c = 1.1), IFDF (c = 5) and IFDF (c = 1000) gave effective results. IFDF (c = 0) failed due to the adaptive diffusion coefficients were not employed, while IFDF (c = 1.0001) failed due to IFDF was too large. The calculated keff of IFDF (c = 1.1) was 1.00671, which was the same as the 1.00671 of SI, while the calculated keff of IFDF (c = 5) and IFDF (c = 1000) was slightly different from that of SI, which may be caused by the possible inconsistency between the high-order problem and low-order problem. The total CPU time of IFDF (c = 1.1) was significantly reduced from 5364.7 s to 345.1 s, which was mainly due to the obvious reduction of the number of transport sweep. The low-order calculation time of IFDF (c = 1.1) was 41.3 s, which was the overhead of the IFDF acceleration. With the increase of c, the speedup was decreased but the divergence risk of large IFDF was also decreased. It was remarkable that the low-order calculation time of IFDF (c = 1000) was smaller than that of IFDF (c = 5). Fig. 5 presents the convergence process
Table 1 Results for the MOX-3600 problem. Item
SI
CMFD /adCMFD
IFDF (c = 0)
IFDF (c = 1.001)
IFDF (c = 1.1)
IFDF (c = 5)
IFDF (c = 1000)
keff Number of fission-source iteration Number of transport sweep Total CPU time/s Speedup Low-order calculation time/s
1.00671 392 48,408 5364.7 1.0
Not converged
Not converged
Not converged
1.00671 22 2732 345.1 15.5 41.3
1.00670 26 3197 402.0 13.3 42.7
1.00673 84 9954 1147.0 4.7 22.3
8
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
of SI and IFDF (c = 5). It clearly shows that IFDF (c = 5) converged much faster than SI. Table 2 presents the region average adjoint flux results, which shows IFDF acceleration caused small derivation to the distributed solutions. 3.2. The small subcritical fast reactor problem
Fig. 5. Iteration errors of different methods.
Table 2 Region average adjoint flux results for the MOX-3600 problem. Material region
Inner fuel 1 Inner fuel 2 Inner fuel 3 Inner fuel 4 Inner fuel 5 Outer fuel 1 Outer fuel 2 Outer fuel 3 Outer fuel 4 Outer fuel 5 Lower gas Lower reflector Lower gas Upper reflector Radial reflector Primary control Secondary control Duct
Normalized region average adjoint flux SI
IFDF(c = 5)
Deviation/%
3.911E08 5.614E08 6.102E08 5.290E08 3.370E08 2.566E08 3.723E08 4.079E08 3.584E08 2.331E08 2.148E09 1.111E08 1.526E08 2.985E09 1.379E09 3.157E09 2.645E09 2.344E08
3.915E08 5.619E08 6.106E08 5.291E08 3.368E08 2.556E08 3.709E08 4.063E08 3.569E08 2.320E08 2.163E09 1.114E08 1.515E08 2.997E09 1.393E09 3.162E09 2.650E09 2.345E08
0.1 0.1 0.1 0.0 0.0 0.4 0.4 0.4 0.4 0.5 0.7 0.3 0.7 0.4 1.0 0.1 0.2 0.0
This is a small subcritical fast reactor problem modified from the KNK-II (Lesage, 2001) problem. The radial flat distance of the hexagonal assembly is 12.9904 cm and the reactor core has a height of 190 cm and is consisted of 169 hexagonal assemblies. Fig. 6 indicates the radial layout of the reactor core. The core layout and 4-group cross sections are totally the same as the control rod half-inserted case described in reference (Lesage, 2001). However, in this case, an external source is added to the central assembly in the height range of 90–100 cm. The external-source density from the first group to the last group are 1.0, 0.1, 0.001, 0.0001 (unit: cm3s1) respectively. This problem has been calculated with DNTH by different iteration schemes, including SI, IFDF, IFDF-RE (rebalance accelerated IFDF), CMFD-RE (rebalance accelerated CMFD), and adCMFD-RE (rebalance accelerated adCMFD, g ¼ 1=4). The parameter c in Eq. (91) was set as 5. 33 axial meshes were used in the calculation. We have normalized the distribution values by requiring the total neutron production rate to be 1 s1 (Lesage, 2001). The overview results are presented in Table 3. Small derivations of the ks were found, which were mainly caused by different convergence level of the inner iteration. When IFDF was employed, Table 3 Results for the small subcritical fast reactor problem. Item
SI
IFDF
IFDF-RE
CMFDRE
adCMFDRE
ks Number of fissionsource iteration Number of transport sweep Total CPU time/s Speedup Low-order calculation time/s
0.98729 575
0.98740 37
0.98738 37
0.98738 37
0.98738 48
29,850
621
618
585
844
705.9 1.0
23.9 29.5 5.6
20.0 35.3 1.5
19.1 37.0 1.5
29.0 24.3 0.8
Fig. 6. Radial core layout of the small subcritical fast reactor (Lesage, 2001).
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056 Table 4 Region average flux results for the small subcritical fast reactor problem. Group
SI
IFDF-RE
Derivation/%
Axial Blanket 1 2 3 4
2.221E05 2.778E05 1.581E05 1.140E05
2.221E05 2.778E05 1.582E05 1.141E05
0.00 0.01 0.03 0.12
Test Zone 1 2 3 4
1.558E04 1.163E04 2.677E05 2.405E06
1.557E04 1.163E04 2.676E05 2.406E06
0.06 0.04 0.02 0.02
Driver without moderator 1 1.042E04 2 7.915E05 3 2.530E05 4 5.220E06
1.042E04 7.914E05 2.530E05 5.221E06
0.01 0.01 0.01 0.02
Control rod 1 2 3 4
4.712E05 3.785E05 8.045E06 6.201E07
4.709E05 3.783E05 8.042E06 6.201E07
0.06 0.06 0.04 0.01
Control rod follower 1 1.160E04 2 9.800E05 3 3.597E05 4 1.049E05
1.160E04 9.801E05 3.598E05 1.049E05
0.01 0.01 0.02 0.03
9
the total calculation time was reduced from 705.9 s to 23.9 s. When the low-order problem was further accelerated by the rebalance acceleration method, the low-order calculation time of IFDF was further reduced from 5.6 s to 1.5 s of IFDF-RE which achieved a speedup of 35 compared with SI. The traditional CMFD acceleration with rebalance accelerated low-order calculation could also converge and obtained similar speedup with IFDF-RE, while adCMFD-RE gave significantly lower speedup than IFDF-RE. Table 4 compares the region average flux derivations between SI and IFDFRE. The derivations were insignificant, which further proved the effectiveness of the IFDF acceleration method for the externalsource multiplication problem. 3.3. The ZPPR-10B fast reactor problem This is a large fast reactor problem of the ZPPR-10, phase B (Lesage, 2001) (ZPPR-10B) fast critical experiments of Argonne National Laboratory. The radial flat distance of the square assembly is 5.5245 cm. The reactor core has a height of 209.56 cm. Fig. 7 gives the radial layout of the reactor core. 33-group homogenized cross sections for every assembly were obtained by the TULIP code in SARAX. With these cross sections, we have conducted the whole-core transport calculation with NECP-DNTC (DNTC). The core was divided by 24 axial meshes and 59 61 radial square meshes. We have normalized the distribution values by requiring the total power to be 1 MW.
Fig. 7. Radial core layout of the ZPPR-10B fast reactor (Lesage, 2001).
10
Z. Xu et al. / Annals of Nuclear Energy 136 (2020) 107056
Table 5 Results for the ZPPR-10B problem.
Acknowledgements
Item
SI
IFDF
keff Number of fission-source iteration Number of transport sweep Total CPU time/s Speedup Low-order calculation time/s
0.99619 152 24,100 9377.0 1.0
0.99620 13 1993 898.3 10.4 89.4
This work was partially supported by the National Natural Science Foundation of China (11775170) and partially supported by the National Key R&D Program of China (2017YFE0302200). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.107056.
Table 6 Normalized neutron production rate results for the ZPPR-10B problem. Material region
Inner fuel 1 Inner fuel 2 Inner fuel 3 Outer fuel Axial blanket Radial blanket
Normalized region average adjoint flux DNTH-SI
DNTH-IFDF
Deviation/%
4.202E+10 3.911E+08 2.256E+08 4.223E+10 1.148E+09 1.518E+09
4.203E+10 3.915E+08 2.259E+08 4.222E+10 1.148E+09 1.517E+09
0.0 0.1 0.1 0.0 0.0 0.1
The overview results are presented in Table 5. For this problem, only SI and IFDF were considered. The parameter c in Eq. (91) was set as 5. SI gave a keff of 0.99619, while IFDF gave a close keff of 0.99620. The total calculation time of SI was 9377.0 s, while that of IFDF was only 898.3 s, which meant that a speedup of 10.4 was achieved when the low-order calculation overhead was 89.4 s. Table 6 compares the normalized neutron production rates of SI and IFDF. It is seen that the deviations were within 0.1%, which further proved the effectiveness of the IFDF acceleration method in this paper for cuboid meshes. 4. Conclusion Firstly, the SN nodal formulations in hexagonal prism and cuboid meshes have been introduced. Then, the improved interface flux discontinuity factor (IFDF) acceleration method was proposed. The detailed IFDF formulations were derived and applied to the acceleration of the SN nodal calculation for the eigenvalue problem and external-source multiplication problem. The IFDF acceleration method was improved by proposing a formulation to obtain the adaptive diffusion coefficients to ensure IFDF to be positive and keep it in a reasonable region, which could effectively improve the convergence stability. Finally, the theory has been implemented in the NECP-DNTH and NECP-DNTC code and three fast reactor problems were tested. For the large assembly MOX-3600 problem, the IFDF acceleration method obtained a speedup of 13, while the traditional CMFD or the artificial diffusion (g ¼ 1=4) CMFD (adCMFD) could not converge, which proved that the IFDF acceleration method was efficient and stable. For the small subcritical fast reactor problem, the rebalance accelerated IFDF acceleration achieved a speedup of 35, which was similar with that of the rebalance accelerated traditional CMFD and significantly higher than that of rebalance accelerated adCMFD. For the large ZPPR-10B problem with square assemblies, IFDF obtained a speedup of 10, which proved the effectiveness of the IFDF acceleration method for cuboid meshes. In summary, the proposed adaptive diffusion coefficients formulation with a proper damping factor (c ¼ 5) was valid for all the three tested cases and effectively improved the convergence stability for the MOX-3600 problem with large assembly size.
References Yang, W.S., 2012. Fast reactor physics and computational methods. Nucl. Eng. Technol. 44 (44), 177–198. Ruggieri JM, Tommasi J, Lebrat JF, et al. ERANOS 2.1: International code system for GEN IV fast reactor analysis[C]. Proceedings of ICAPP ’06, Reno, NV USA, June 48, 2006. Xu ZT, Wu HC, Zheng YQ. Development of an optimized transport solver in SARAX for fast reactor analysis[C]. In Proceedings of the 26th International Conference on Nuclear Engineering. London, UK, 2018. Willert, J., Park, H., Knoll, D.A., 2014. A comparison of acceleration methods for solving the neutron transport k-eigenvalue problem. J. Comput. Phys. 274, 681– 694. Smith KS. Full-core, 2-D, LWR core calculations with CASMO-4E[C]. PHYSOR 2002, Seoul, Korea, 2002. Cho, N.Z., 2017. An overview of p-CMFD acceleration and its applications to reactor physics transport calculation. Trans. Am. Nucl. Soc. 117, 1247–1250. Jarrett, M., Kochunas, B., Zhu, A., et al., 2016. Analysis of stabilization techniques for CMFD acceleration of neutron transport problems. Nucl. Sci. Eng. 184 (2), 208– 227. Zhu, A., Jarrett, M., Xu, Y., et al., 2016. An optimally diffusive coarse mesh finite difference method to accelerate neutron transport calculations. Ann. Nucl. Energy 95, 116–124. Yuk, S., Cho, N.Z., 2017. Two-level convergence speedup schemes for p-CMFD acceleration in neutron transport calculation. Nucl. Sci. Eng. 188 (1), 1–14. Xiao, S., Ren, K., Wang, D., 2018. A local adaptive coarse-mesh nonlinear diffusion acceleration scheme for neutron transport calculations. Nucl. Sci. Eng. 189 (3), 272–281. Yamamoto, A., Endo, T., Giho, A., 2018. Transport Consistent Diffusion Coefficient for CMFD Acceleration. Trans. Am. Nucl. Soc. 119, 1179–1181. Masiello, E., 2018. ‘‘On the fly” stabilization of the coarse-mesh finite difference acceleration for multidimensional discrete-ordinates transport calculations. J. Comput. Phys. 373, 1–27. Masiello, E., Sanchez, R., Zmijarevic, I., 2009. New numerical solution with the method of short characteristics for 2-d heterogeneous cartesian cells in the apollo2 code: Numerical analysis and tests. Nucl. Sci. Eng. 161 (3), 257–278. Xu YL, Downar T. Convergence analysis of a CMFD method based on generalized equivalence theory[C]. PHYSOR 2012, Knoxville, Tennessee, USA, 2012. Lu, H.L., Wu, H.C., 2007. A nodal SN transport method for three-dimensional triangular-z geometry. Nucl. Eng. Des. 237 (8), 830–839. Xu, Z.T., Wu, H.C., Zheng, Y.Q., 2019. An improved hexagonal-z nodal SN method in SARAX code system. Prog. Nucl. Energy. submitted for publication. Zheng, Y.Q., Du, X.N., Xu, Z.T., et al., 2018. SARAX: a new code for fast reactor analysis part I: methods. Nucl. Eng. Des., 421–430 Zheng, Y.Q., Qiao, L., Zhai, Z.A., et al., 2018. SARAX: A new code for fast reactor analysis part II: verification, validation and uncertainty quantification. Nucl. Eng. Des. 331, 41–53. Aragones, J.M., Ahnert, C., 1986. A linear discontinuous finite difference formulation for synthetic coarse-mesh few-group diffusion calculations. Nucl. Sci. Eng. 94 (4), 309–322. Badruzzaman, A., 1985. An efficient algorithm for nodal-transport solutions in multidimensional geometry. Nucl. Sci. Eng. 89 (3), 281–290. Azmy, Y.Y., 1997. Multiprocessing for neutron diffusion and deterministic transport methods. Prog. Nucl. Energy 31 (3), 317–368. Zhou SC, Wu HC, Cao LZ, et al. Comparison of three acceleration techniques for the fixed source neutron transport calculation[C]. In Proceedings of the PHYSOR 2018, Cancun, Mexico, 2018. Xu, Z.T., Wu, H.C., Zheng, Y.Q., 2017. Hexagonal CMFD accelerated discrete nodal transport method in triangular-z geometry. Trans. Am. Nucl. Soc. 117, 718–720. Bernnat W, Blanchet D, Brun E, et al. Benchmark for neutronic analysis of sodiumcooled fast reactor cores with various fuel types and core sizes[R]. Nuclear Energy Agency Technical Report, NEA/NSC, 2015. Takeda T, Ikeda H. 3-D neutron transport benchmarks[R]. NEACRP-L-330, OECD/ NEA Committee on Reactor Physics, 1991. Lesage LG. An overview of the Argonne National Laboratory fast critical experiments [R]. ANL-NT-175, Argonne National Laboratory, 2001.