- Email: [email protected]
Local Iteration Method with SP3
Annals of Nuclear Energy 118 (2018) 49–60
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A Next Generation Method for Light Water Reactor core analysis by using Global/Local Iteration Method with SP3 Tzung-Yi Lin ⇑, Yen-Wan Hsueh Liu Institute of Nuclear Engineering and Science, National Tsing Hua University, 101, Section 2, Kuang-Fu Road, Hsinchu 30013, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 2 October 2017 Received in revised form 21 March 2018 Accepted 21 March 2018
Keywords: Next Generation Method (NGM) Simplified P3 method (SP3) Global/Local Iteration Method (GLIM) Even Parity Discontinuity Factor (EPDF)
a b s t r a c t The Global/Local Iteration Method (GLIM) capability was implemented into the in-house developed nodal code, NuCoT, which was based on the Hybrid Nodal Green’s Function Method (HNGFM) with the diffusion/SP3 approximations. A procedure was proposed for obtaining the Even Parity Discontinuity Factor (EPDF) in SP3 method. The accuracy of NuCoT with GLIM was verified by using two benchmark problems, KAIST and C5G7, with reference solutions calculated by NEWT and MCNP respectively. The results of NuCoT show that the accuracy of the GLIM is comparable to the whole core pin-wise calculation. The SP3 calculation with EPDF is better than the diffusion result. For the KAIST problem using GLIM, the Root-Mean-Square (RMS) of pin power error is 1%, the maximum pin power error is 6%, and the keff error is <83 pcm for All-Rod-In case and <135 pcm for All-Rod-Out case. For the C5G7 problem, the RMS of pin power is <0.8% and the maximum error is <4%. The keff error is <80 pcm in 2D problem and <200 pcm in 3D rodded problem. Ó 2018 Elsevier Ltd. All rights reserved.
1. Introduction The two-step method is efficient and widely used for the current Light Water Reactor (LWR) core analysis. Step 1 is single fuel assembly calculation, in which transport method and reflective boundary condition is used. The assembly-homogenized cross sections and assembly discontinuity factors are obtained in this step based on the Generalized Equivalence Theory (Smith, 1986) (GET). Step 2 is the whole core coarse mesh calculation. Normally two group nodal diffusion method is used. However, the reflective boundary used in the first step causes the homogenized parameters used in the second step incorrect for fuel assemblies near the core boundary, and for cores with highly heterogeneous fuel loading. The environment dependence of the assembly cross sections has been added into the homogenization parameters by several different approaches (Smith, 1994; Yamamoto and Tatsumi, 2006; Dall’Osso et al., 2010) with minor change of the two-step method. However, these approaches were challenged when the fuel and core design become more complicated such as the optimized fuel loading pattern, the mixed oxide fuel (MOX) and the diamond type water rod. The axially heterogeneous fuel and
⇑ Corresponding author. E-mail address: [email protected] (T.-Y. Lin). https://doi.org/10.1016/j.anucene.2018.03.030 0306-4549/Ó 2018 Elsevier Ltd. All rights reserved.
partially control rod insertion also affect the accuracy of the local peaking factor (Zhang et al., 2013) calculated by these approaches. There are several approaches proposed for the Next Generation Method of neutronic analysis of LWR to obtain more accuracy results. One is to perform the whole core 3D calculation by 2D(transport)/1D(diffusion) coupling method (Joo et al., 2004; Jung et al., 2013), which is accurate, but not practical considering current computing capability. One of more practical methods is the whole core multi-group 3D pin by pin calculation (pin-wise calculation) (Tatsumi et al., 2003), in which the homogenized cross sections and discontinuity factors of the pin cell rather than the assembly are generated in the first step. In the pin-wise calculation, the pin power reconstruction is not necessary. In order to obtain more accurate result, the Simplified P3 (SP3) (Gelbard, 1968) method is proposed for the pin-wise calculation to replace the diffusion method. However, without discontinuity factor, the solution of SP3 calculation is no better than that of the diffusion calculation. The discontinuity factor for the SP3 calculation has been investigated in many studies (Kozlowski et al., 2011; Yu et al., 2014; Yu and Chao, 2015; Yamamoto et al., 2016). The whole core pin-wise calculation using SP3 has been implemented into NuCoT and showed good accuracy (Lin and Liu, 2017) in the benchmark problems. Another practical method is the Global/Local Iteration Method (Kim, 1993) (GLIM). It is similar to the Iterative TransportDiffusion Methodology (Ivanov, 2007; Roberts, 2010; Colameco,
50
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
2012) (ITDM) or a hybrid method (Kozlowski et al., 2011; Lee and Downar, 2004; Lee et al., 2004; Bahadir et al., 2005; Knight et al., 2012; Yasseri and Rahnema, 2014; Kooreman and Rahnema, 2015). In GLIM, the assembly-homogenized cross sections of the assembly (or lattice) are generated by the transport calculation using the actual assembly surrounding environment. However, since most transport codes can only accept typical boundary conditions such as reflective, white or vacuum boundary, it is not easy to modify the existing transport code for used in the GLIM if one is not familiar with the code. One easier way is to apply the multigroup pin-wise calculation for the assembly calculation. This study presents the GLIM and its application to two benchmark problems, KAIST and C5G7. It includes the derivation of the Hybrid Nodal Green’s Function Method (Ougouag, 1981; Chang, 2011) (HNGFM) with the diffusion and SP3 method for the global and local calculation, the discontinuity factor calculation, and the interface between the global and local calculation.
/gx ðxÞ ¼
The GLIM consists of two parts of calculations using different nodal size. The global calculation is a whole core calculation using the assembly-wise size nodal method, which is usually a fewgroup nodal diffusion method. The local calculation is an assembly calculation using the pin-cell size nodal method, which in this study can be a fine-group diffusion or SP3 nodal method. The local calculation is applied to every assembly nodes and the reflector nodes. The iteration process is shown in Fig. 1. The interface between the global and local calculation called the submesh interface is described in Section 2.4. 2.1. Global calculation using HNGFM with diffusion approximation The steady-state three-dimensional diffusion equation within one node on the Cartesian system for energy group g is typically written in the following form.
X
2 X /gxl Pl ð^xÞ
where Pl ð^xÞ is the l-th order Legendre polynomial and ^x ¼ x=ax . And the transverse leakage profiles were assumed to be flat. The weighted residual procedure was used to solve /gxl using Eq. (2). Then the unsolved term of Eq. (2) at the surface, hQGgx ðax Þi, can be obtained.
hQGgx ðax Þi ¼
2 X
hQ gxl Pl Ggx ðax Þi
ð4Þ
l¼0
where
hPl Ggx ðax Þi
Q gxl
1 2ax
Z
ax
ax
dx Ggx ðax jxÞP l ð^xÞ;
8 S Lgx ; > l¼0 > < gx G G X X ¼ vg f s > > : k 0 mRg0 /g0 xl þ 0 Rgg0 /g0 xl ; l ¼ 1; 2 g –g
g
The continuity of the surface current and discontinuity of the surface flux at the interface of two adjacent nodes provide the rela by using tionship between the current J and nodal average flux / Eq. (2) and 1D nodal balance equation. This relationship can be substituted into the 3D nodal balance equation of node N to obtain Eq. (5).
N BNg / g
Xh
i ðN;m1Þ ðN;m1Þ ðN;mþ1Þ ¼ RN /g / Bgu þ BðN;mþ1Þ gu g g
ð5Þ
u¼x;y;z
where
m ¼ node index related to the u direction; ( BNg
¼R
r;N g
Xh
1þ
ANgu aNu
þ
ANgu aNu
i
) ;
u¼x;y;z
ð1Þ
where
Rrg ¼ Rtg Rsgg ;
ðN;m1Þ Bgu ¼
m1 2am1 Rr;m1 xm1 au j am1 u gu u g ; 2aNu DxNgu ðaNu Þ Xh
RNg ¼ SNg
i C Ngu ðaNu Þ þ C Ngu ðaNu Þ ;
u¼x;y;z
and
Sg ðrÞ ¼
ð3Þ
l¼0
and
2. Global/Local Iteration Method
@2 Dg 2 /g ðrÞ þ Rrg /g ðrÞ ¼ Sg ðrÞ for ðrÞ @u u¼x;y;z 2 ½ax ; ax ; ay ; ay ; ½az ; az
To evaluate hQGgx ðxÞi at the surface, the quadratic polynomial was used to approximate the spatial flux distribution.
SN ¼ vg g k
G G X vg X mRgf 0 /g0 ðrÞ þ Rsgg0 /g0 ðrÞ:
k
g0
g 0 –g
After deriving one-dimensional equation from Eq. (1) and applying the Green’s function with Neumann (or second-type) boundary condition, the 1-D flux in x-direction within a node can be written as:
/gx ðxÞ ¼ Ggx ðxj ax ÞJ gx ðax Þ Ggx ðxjax ÞJ gx ðax Þ þ 2ax hQGgx ðxÞi
ð2Þ
Ggx ðxjx0 Þ = Green’s function at x from x0 with Neumann boundary condition, Q gx ðxÞ ¼ Sgx ðxÞ Lgx ðxÞ; P v P Sgx ðxÞ ¼ kg Gg0 mRgf 0 /g0 x ðxÞ þ Gg0 –g Rsgg 0 /g0 x ðxÞ; and Lgx ðxÞ is the transverse leakage profile.
G X
g0
g 0 –g
mRgf 0 / Ng0 þ
N0 ; Rsgg0 / g
xNgu aNu j aNu ANgu aNu ¼ ; DxNgu ðaNu Þ C Ngu ðaNu Þ ¼
where d J gx ¼ Dg dx /gx ; Ra hQGgx ðxÞi 2a1x ax x dx0 Ggx ðxjx0 ÞQ gx ðx0 Þ;
G X
N ðaN Þ DhQGiNgu ðaNu Þ DQ gu u ; N N 2au Dxgu ðaNu Þ
xNgu ðAjBÞ ¼ f Ngu ðAÞGNgu ðAjBÞ; N
f gu ðAÞ ¼ Discontinuity factor at the location A in u direction; h i DxNgu ðaNu Þ ¼ xNgu aNu j aNu xNgu aNu j aNu h m1 m1 i xm1 au j am1 au j am1 xm1 ; gu u gu u
51
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Fig. 1. The flow of the Global/Local Iteration Method.
N N N ðaN Þ ¼ 2aN Q N N DQ gu u u gu xgu au j au m1 xm1 am1 j am1 ; 2am1 Q u
gu
gu
u
2.2. Local calculation using HNGFM with SP3 approximation
u
The typical SP3 equation of group g of a node on the Cartesian system is given below:
N ¼ SN LN ; Q gu g gu
LN ¼ gu
X
@2 F ng ðrÞ þ Rrng F ng ðrÞ ¼ T ng ðrÞ; 2 @u u¼x;y;z for r 2 ½ax ; ax ; ½ay ; ay ; ½az ; az
i X 1 hN J gh ðah Þ J Ngh ðah Þ ; 2ah h–u
Dng
ð6Þ
where
D1g ¼
and N
n ¼ 1; 2
m1
hDQGgu iN ðaNu Þ ¼ 2aNu f gu ðaNu ÞhQGNgu ðaNu Þi 2am1 f gu ðam1 Þ u u m1 Þi: hQGm1 gu ðau
The available boundary conditions are albedo, zero flux and known incoming partial current. The surface current and surface flux can be related by the albedo. After using Eq. (2) and 1D nodal balance equation, the surface current on the boundary can be related to its nodal average flux to obtain a modified Eq. (5) for the node on the boundary. If the incoming partial current is known, it can be treated as an external source with albedo set to zero.
1 ; 3Rtrg
D2g ¼
1 ; 7Rtg
Rr1g ¼ Rtg Rsgg ;
2 T 1g ðrÞ ¼ Sg ðrÞ þ Rr1g F 2g ðrÞ; 3 Sg ðrÞ ¼
2 2 T 2g ðrÞ ¼ Sg ðrÞ þ Rr1g F 1g ðrÞ; 3 3
G G X vg X mRgf 0 /g0 0 ðrÞ þ Rsgg0 /g0 0 ðrÞ
k
g0
4 9
Rr2g ¼ Rtg Rsgg ;
g 0 –g
F 1g ðrÞ ¼ /g0 ðrÞ þ 2/g2 ðrÞ; F 2g ðrÞ ¼ 3/g2 ðrÞ; and /g0 ðrÞ and /g2 ðrÞ are the flux moments.
52
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Similar to the diffusion approximation, the one-dimensional pseudo fluxes in x-direction within a node can be written as
C Nngu ðaNu Þ ¼
N ðaN Þ DhQGiNngu ðaNu Þ DQ ngu u : 2aNu DxNngu ðaNu Þ
F ngx ðxÞ ¼ Gngx ðxj ax ÞJ ngx ðax Þ Gngx ðxjax ÞJ ngx ðax Þ þ 2ax hQGngx ðxÞi n ¼ 1;2 ð7Þ
where
xNngu ðAjBÞ ¼ f Nngu ðAÞGNngu ðAjBÞ; N
f ngu ðAÞ ¼ Discontinuity factor at the location A in u direction;
d J ngx ¼ Dng F ngx ; dx
m1 DxNngu ðaNu Þ ¼ ½xNngu ðaNu j aNu Þ xNngu ðaNu j aNu Þ ½xm1 j ngu ðau
Gngx ðxjx0 Þ ¼ Green’s function at x0 with Neumann boundary condition;
hQGngx ðax Þi
1 2ax
Z
ax
ax
m1 Þ xm1 j am1 Þ; am1 u ngu ðau u N N N m1 m1 m1 N ðaN Þ ¼ 2aN Q N Q ngu xngu ðam1 DQ j ngu u u ngu xngu ðau j au Þ 2au u
dx0 Gngx ðax jx0 ÞQ ngx ðx0 Þ
Þ; am1 u
Q ngx ðxÞ ¼ T ngx ðxÞ Lngx ðxÞ; N ¼ T N LN ; Q ngu ng ngu
2 T 1g ðxÞ ¼ Sgx ðxÞ þ Rr1g F 2gx ðxÞ; 3
LN ¼ ngu
2 2 T 2g ðxÞ ¼ Sgx ðxÞ þ Rr1g F 1gx ðxÞ; 3 3 Sgx ðxÞ ¼
vg
G X
G X
k
g0
g 0 –g
mRgf 0 /g0 0x ðxÞ þ
and N
Xh
i ðN;m1Þ þ BðN;mþ1Þ F ðN;mþ1Þ ¼ RN ; BðN;m1Þ F ngu ngu ng ng ng
m1 Þi: hQGm1 ngu ðau
n ¼ 1; 2
u¼x;y;z
ð8Þ where
m ¼ node index related to the u direction; ( BNng ¼ Rr;N 1þ ng
Xh
) i ; ANngu ðaNu Þ þ ANngu aNu
u¼x;y;z
BðN;m1Þ ¼ ngu
m1 2am1 Rr;m1 xm1 j am1 Þ u ngu ðau u ng ; N N N 2au nxngu ðau Þ
Xh
RNng ¼ T Nng
i C Nngu ðaNu Þ þ C Nngu ðaNu Þ ;
u¼x;y;z
2 N ; T N1g ¼ SNg þ Rr1g F 2g 3 SN ¼ vg g k
2 2 N ; T N2g ¼ SNg þ Rr1g F 1g 3 3
G X
G X
g0
g 0 –g
mRgf 0 / Ng0 0 þ
ANngu ðaNu Þ ¼
x
N0 ; Rsgg0 / g0
N N N ngu ðau j au Þ ; N N D ngu ðau Þ
x
m1
DhQGiNngu ðaNu Þ ¼ 2aNu f ngu ðaNu ÞhQGNngu ðaNu Þi 2am1 f ngu ðam1 Þ u u
Rsgg0 /g0 0x ðxÞ;
and Lngx ðxÞ is the transverse leakage profile which was assumed to be flat. The unsolved term of Eq. (7) at the surface, hQGngx ðax Þi, can be obtained by the same procedure as in the diffusion method. Also, the relationship between the surface current J n and nodal average n was used to obtain Eq. (8) by utilizing the continupseudo flux F ity of the surface current ðJ 1 and J 2 Þ and discontinuity of the surface pseudo flux ðF 1 and F 2 Þ, Eq. (7), at the interface of two adjacent nodes, following the same procedure for obtaining Eq. (5) in the diffusion method.
N BNng F ng
i X 1 h J Nngh ðah Þ J Nngh ðah Þ ; 2a h h–u
The boundary conditions available in SP3 calculation are albedo, zero pseudo flux and known incoming partial current. The partial current moments (J 1 and J 2 ), based on the Marshak coefficient (Hébert, 2009) are:
(
1 J 1gu ðuÞ ¼ 14 F 1gu ðuÞ 16 F 2gu ðuÞ 12 J 1gu ðuÞ 1 7 F 1gu ðuÞ þ 48 F 2gu ðuÞ 12 J 2gu ðuÞ J 2gu ðuÞ ¼ 16
ð9Þ
Therefore, for the albedo boundary condition, the surface currents (J 1 and J 2 ) and the surface pseudo fluxes (F 1 and F 2 ) can be related by the albedo. Then a modified Eq. (8) for the boundary node can be obtained by using the similar procedure as in the diffusion method. 2.3. Even Parity Discontinuity Factor (EPDF) In SP3 calculation, the discontinuity factors of two pseudo fluxes cannot be easily generated because the 2nd and 3rd angular flux moments were not available from the widely used transport codes. EPDF was developed for the transport code by Yamamoto et al. (2011). It was unified in the diffusion/SP3 calculation by Yu and Chao (2015). The idea is to make the homogeneous even parity surface flux Uhom (not the scalar flux /hom ) between two adjacent 0 nodes continuous after applying EPDF. The Uhom on the left- and right-hand side of the interface are continuous by introducing the factor d. hom F Right Uhom Right ¼ F Left ULeft
ð10Þ
where
F Right ¼ 1 d;
F Left ¼ 1 þ d;
and
d¼
hom Uhom Right ULeft hom Uhom Right þ ULeft
:
The even parity surface flux can be rewritten as the average value of the partial incoming and outgoing currents.
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
U¼
J in þ J out 2
In diffusion approximation, the even parity flux was:
1 4
U ffi /0
ð11aÞ
In SP3 approximation:
1 4
U ffi /0 þ
5 1 1 / ¼ F1 F2 16 2 4 16
or Left
¼
U
U
hom Right or Left
2.4. Submesh interface program
ð13Þ
where
U ¼
hom Uhom Right þ ULeft
2
:
After applying EPDF, the homogeneous even parity surface fluxes are made continuous. In the single fuel assembly transport calculation using reflective boundary condition implies that the neighboring assemblies are the same. Hence, EPDF on the assembly boundary is one. However, the neighboring fuel assemblies in the core may not be the same, in which case EPDF will not be equal to one. Yu and Chao (2015) used the same idea of the conventional discontinuity factor to set U of Eq. (13) to the reference solution Uhet , named Renormalized EPDF (REPDF). The original EPDF definition was named Prime EPDF (PEPDF). The difficulty was to obtain the reference 2nd flux moment (F 2 or /2 ) and current (J 2 or /3 ) for the homogeneous even parity surface flux of the SP3 approximation of HNGFM. In this study, the procedure for obtaining the EPDF for SP3 calculation could have followed Yu et. al.’s idea to iterate EPDF utilizing the available information in the reference solutions, such as surface current J 1 , nodal 0 and keff. However, an alternative method using less average flux / iteration parameters is adopted. The two equations of nodal balance equations from Eq. (6) can be reduced to Eq. (14) by substituting the reference solution of surface current J 1 and nodal average 0 . flux /
X
1 2g ¼ 2 Hg J 2gu ðau Þ J 2gu ðau Þ þ Rt2g F 2a 3 u u¼x;y;z
ð14Þ
where
5 9
Rt2g ¼ Rtg ; Hg ¼
Eq. (14) also appeared in Kozlowski et al. (2011) (Eqs. (28)– (30)), where it was used in a different way but with approximations to calculate the discontinuity factors. It is easy to obtain the EPDF for the diffusion approximation. When the surface net current, nodal average flux and keff are from the reference solution, the only one unknown of Eq. (2) for calculating the homogeneous surface flux is QGgx ðax Þ which needs an iteration to obtain the moment terms of Eq. (4).
ð12bÞ
In this study, Eq. (10) was rewritten to the form similar to the conventional discontinuity factor:
EPDF Right
53
X
1 J 1gu ðau Þ J 1gu ðau Þ 2a u u¼x;y;z
Therefore, a similar form of Eq. (8) can be easily obtained and only 2 matrix is necessarily to be solved. F The iterative procedure for obtaining EPDF is described below: Assuming EPDF is one. 2 matrix from the similar form of Eq. (8) and (1) Solve the F obtain J 2 . (2) Obtain two pseudo surface fluxes by solving Eq. (7), which needs another iteration for QGngx ðau Þ. (3) Update EPDF by using Eq. (13). (4) Repeat step (1)–(3) until converged.
The submesh interface program connects the global and local calculation. The data requirements of the global calculation are the assembly-homogenized cross sections and the assembly discontinuity factors. In the traditional two-step procedure of the core analysis, the transport calculation generates these parameters by using reflective boundary condition. In the GLIM of this study, these parameters were generated by the submesh interface program from the results of the local calculation using boundary conditions based on the global calculation. This local calculation result can be regarded as the heterogeneous solution for the assemblyhomogenized cross sections and assembly discontinuity factors. The assembly discontinuity factors were calculated by the method described at the end of Section 2.3 for the global calculation using the diffusion method. Based on the theory of EPDF, the homogenous even parity surface fluxes on the interface between the left- and right-side nodes (Eq. (10)) were needed at the same time. However, the homogeneous even parity surface flux of the neighboring node at the assembly boundary was unknown because the local calculation was a single assembly calculation. Therefore, REPDF was applied for the global calculation, not PEPDF. The homogenized cross sections and the discontinuity factors of the pin cell used for the local calculation were input data. The pin cell homogenized cross sections were obtained from the single assembly transport calculation with the reflective boundary condition. The pin cell discontinuity factors for the diffusion or SP3 method were calculated in advance by the method described in Section 2.3. The keff for the local calculation was given by the global calculation, instead of being calculated by using source iteration. In other words, the local calculation was a fixed-k calculation. Since the GLIM was an iteration process between the global and local calculation, the boundary condition for the local calculation should be updated after the global calculation using new cross sections and discontinuity factors. The realistic boundary condition (Ivanov, 2007; Roberts, 2010; Colameco, 2012) for the local calculation was the incoming neutron current distribution which was just the outgoing neutron current distribution of the adjacent assembly from the local calculation. In order to apply new global calculation result (step N) for the local calculation, the outgoing current distribution of the adjacent assembly from the previous local calculation (step N-1) was normalized to one. Then the newly averaged outgoing current of the adjacent assembly of the global calculation (step N) multiplied by this normalized outgoing current distribution became the new incoming current distribution of the calculated assembly. In the local calculation using SP3 method, since the second kind of averaged outgoing current did not exist in the global calculation in which diffusion method is used, the second kind of outgoing current distribution was normalized to the averaged outgoing current (first kind) by assuming the ratios of these two kinds of averaged outgoing current were the same. One of the advantages in the GLIM is the size of pin cell of the two neighboring assemblies can be different. By adding additional continuity conditions on the surfaces of the two neighboring assemblies with different pin cell size, the nodal method can be used without further modification. The submesh interface program
54
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
is used to generate the mapped incoming current from the neighboring assembly. In Tada et. al.’s study (Tada et al., 2009), the incoming current can be approximated by a flat, linear or quadratic function in 2D problem. For 3D problems, the incoming current is on a 2D plane, not a 1D line. The linear or quadratic function becomes complicated. The simplest way is to use the areaweighted calculation in this study, i.e. the flat approximation. The process of GLIM can be terminated by examining convergence of parameters such as the surface current or averaged nodal flux. In this study, the submesh interface program used the convergence of the relative nodal power fraction and keff of the global calculation as the stopping criteria. 3. Application of the GLIM The GLIM was implemented into our in-house code, NuCoT. Two well-known benchmark problems, KAIST (Cho, 2000) and C5G7 (Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation, 2003, 2005), were selected to verify its accuracy. The source iteration used in NuCoT was accelerated by the Two-Mode Asymptotic Source Extrapolation (Ougouag, 1981). Root-Mean-Square (RMS) error was used to compare the pin power result of NuCoT with the reference solution.
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P ðP Pref Þ2 =N P RMSð%Þ ¼ 100 Pref =N
ð15Þ
where N = total number of the fuel pin, P = pin power, and Pref=pin power of reference solution. Pin-wise results are included for comparison. In case ‘‘PEPDF”, PEPDF was applied to all surfaces of the pin cells of the local calculation. In case ‘‘REPDF”, REPDF was applied to the surfaces of the pin cells on the assembly/reflector boundary and PEPDF was applied to the rest of the surfaces in the local calculation. The accuracy of the diffusion method and SP3 method in the local calculation was compared. 3.1. KAIST benchmark problem KAIST benchmark problem is a small 2D quarter PWR core consists of three types of UO2 and two types of zoned MOX 17 17 fuels. There are two core conditions, all rods out (ARO) and all rods in (ARI), as shown in Fig. 2. The seven-group macro cross section set including anisotropic scattering cross section was provided by the KAIST 2A benchmark problem (Cho, 2000). The KAIST benchmark problem consists of six types of fuel assemblies, 2.0% UO2 fuel (UOX1), 3.3% UO2 fuel with or without control rod insertion (UOX2 or UOX2cr), UOX2 fuel containing 16 gadolinium (Gd) rods (UOX2ba16), zoned MOX fuel (MOX1) and a MOX1 fuel with 8 Gd rods (MOX1ba8). The reference solution for obtaining EPDF and homogenized pin-cell cross sections of the single fuel assembly were performed by using NEWT of SCALE 6.1 (Scale: A Comprehensive Modeling and Simulation Suite for Nuclear Safety Analysis and Design, 2011). The angular quadrature was S16, and the convergence criteria of keff, flux and source were all 105. The accuracy of NuCoT results using diffusion and SP3 approximation these six fuel types of assembly are shown in Table 1. The convergence criteria for the keff, source and flux were all 105. As expected, without applying discontinuity factors, the error in kinf was large. After applying
discontinuity factors, the NuCoT results were very close to the reference results. The reference solution of the whole core calculation was also obtained by using NEWT. The angular quadrature and the convergence criteria used were the same as in the single assembly calculation. In the NuCoT whole core pin-wise calculation, two pin-cell discontinuity factor types, PEPDF and REPDF were applied. The convergence criteria were all 105. The PEPDF and REPDF of pin cells in the fuel assembly were obtained by the single assembly calculation. The reflector including the baffle was divided into 17 17 pin cells as the fuel assembly. Since the pin cells of the reflector were homogeneous, the pin-cell discontinuity factors were set to one, except the one at the reflector boundary adjacent to the heterogeneous fuel assembly. The accuracy of the pin-wise and the GLIM capabilities of NuCoT are shown in Table 2. In the ARO case, the results are insensitive to the choice of discontinuity factors. In the ARI case, using PEPDF gives better results than using REPDF. The SP3 results are better than the diffusion results. In SP3 calculation using PEPDF, the maximum pin power error is <6% and the RMS of pin power error is 1%. The location of the maximum pin power error is near the reflector and its power is much smaller than one. The keff error is <70 pcm for ARI case and <130 pcm for ARO case. PEPDF was selected for the following GLIM calculation. In the NuCoT GLIM calculation, the convergence criteria for the keff, source and flux of the global and local calculations were all 105. The convergence criteria for the keff and relative power fraction in the submesh interface were 104 and 5x103, respectively. The energy groups used in the local calculation was seven. In the global calculation, both two and seven energy groups are used (named GLIM-2G or GLIM-7G). As shown in Table 2, the accuracy of GLIM calculation is comparable to the pin-wise calculation. The RMS of pin power and error of keff are slightly larger than the results of direct pin-by-pin calculations. The GLIM-7G result was not better than GLIM-2G result. It may be caused by the inconsistency of the average surface outgoing partial current (or net current) between the local calculation and the global calculation since it is not one of the convergent criteria. The computing time of GLIM was 2–3 times of the pin-wise calculation as shown in Table 3. Since the local calculations of each assembly-size node are independent, parallelism can be applied to the GLIM to reduce computing time.
3.2. Results of C5G7 benchmark problem C5G7 benchmark problem is a mini quarter PWR core containing UO2 and zoned MOX 17 17 fuel assemblies, surrounded by the reflector. This benchmark problem includes one 2D case and three 3D cases, Unrodded, RoddedA and RoddedB, as shown in Fig. 3. The fuel patterns of 3D cases are the same as the 2D case. The reference solution for EPDF iteration and homogenized pincell cross sections of the fuel assembly were obtained by the single assembly calculation using NEWT using S16 angular quadrature set. The convergence criteria for keff, flux and source were 105. The seven group macro cross section set, including the transportcorrected total and scattering cross sections, were provided by the benchmark problem (Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation, 2003, 2005). Table 4. shows the results of NuCoT pin-wise single assembly calculation compared with the reference solution for the fuel types, UO2, MOX and the rodded fuels (UO2cr and MOXcr), using both the diffusion and SP3 approximation. NuCoT results obtained by using discontinuity factors agree well with the reference results.
55
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Fig. 2. KAIST benchmark problem (2A).
Table 1 Results of KAIST single assembly benchmark problem by NuCoT pin-wise calculation with and without using discontinuity factors in diffusion and SP3 calculation. Fuel type
a b c d e
Diffusion
SP3
No DFa
EPDF
No DF
EPDF
MOX1 (1.176597)b
69.2c 0.67d (1.17)e
0.7 0.00 (0.01)
44.5 0.36 (0.83)
0.5 0.03 (0.07)
MOX1ba8 (1.136757)
333 0.78 (4.07)
0.3 0.00 (0.01)
240 0.48 (2.88)
0.5 0.04 (0.10)
UOX1 (1.100852)
72.9 0.30 (0.72)
4.0 0.00 (0.01)
34.4 0.25 (0.59)
3.7 0.01 (0.02)
UOX2 (1.252353)
84.4 0.36 (0.86)
4.3 0.00 (0.00)
40.9 0.30 (0.68)
4.6 0.01 (0.03)
UOX2ba16 (1.043023)
2459 1.47 (8.24)
4.2 0.00 (0.00)
1426 1.23 (4.57)
7.7 0.09 (0.22)
UOX2cr (0.818892)
4927 2.37 (4.66)
1.7 0.00 (0.00)
2596 1.90 (3.52)
5.9 0.14 (0.42)
Without using discontinuity factor. Reference kinf. kinf error (pcm). RMS error (%) of pin power. Maximum pin power error (%).
3.2.1. 2D problem Table 5 shows the accuracy of NuCoT results of C5G7 2D problem using the pin-wise and GLIM calculations. The reflector of the 2D problem was divided into 17 17 cells as the fuel assembly. Since the reflector is homogeneous, the pin-cell discontinuity factor of the reflector was applied by the same approach as for the KAIST benchmark problem. The reference solution of C5G7 2D
problem was obtained from MCNP (MNCP, 2008) calculation using the seven-group macro cross section set provided by the benchmark problem (Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation, 2003, 2005). The statistic error of keff was 7 pcm and of pin power was <0.5%. The convergence criteria of NuCoT for keff, flux and power were all 105.
56
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Table 2 Results of 2D whole core KAIST benchmark problem using the pin-wise calculation and GLIM capability of NuCoT. Case
Pin-wise
ARO Reference keff = 1.136577
Diffusion SP3
GLIMc
Diffusion SP3
a b c d e
ARI Reference keff = 0.983737
keff Error (pcm)
Pin Power Error(%) RMS/Max
keff Error (pcm)
Pin Power Error(%) RMS/Max
REPDFa PEPDFb REPDF PEPDF
190.8 196.0 142.7 129.8
1.15/5.46 0.97/5.07 1.05/4.75 0.78/5.46
314.1 229.2 140.3 67.8
2.49/5.28 1.73/4.34 1.47/5.03 1.10/5.56
GLIM-7Gd GLIM-2Ge GLIM-7G GLIM-2G
197.0 194.7 134.3 129.9
0.98/5.01 0.97/4.49 0.92/6.02 0.88/5.55
247.0 242.7 76.7 82.6
1.82/4.92 1.79/4.52 1.23/6.02 1.21/5.54
Renormalized EPDF. Prime EPDF. Using prime EPDF. 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
Table 3 CPU time comparison of 2D whole core calculations of KAIST benchmark problem between using the pin-wise calculation and GLIM capability of NuCoT. Case
a b c d
ARO
ARI
CPUa time (sec)
CPU time ratio
CPU time (sec)
CPU time ratio
Diffusion
Pin-wiseb GLIM-7Gc GLIM-2Gd
96.1 172.0 287.7
1 1.8 3.0
73.3 150.9 474.4
1 2.1 6.5
SP3
Pin-wise GLIM-7G GLIM-2G
264.4 893.3 885.0
1 3.4 3.3
273.3 677.0 568.4
1 2.5 2.1
Intel Core i7 930 2.8 GHz. 7 energy groups for the pin-wise calculation with PEPDF (diffusion or SP3). 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
In the NuCoT pin-wise calculation, the keff errors of the diffusion and SP3 method with either PEPDF or REPDF discontinuity factor types are all <100 pcm. The pin power result of the SP3 method is more accurate than the diffusion method. Using PEPDF in SP3 calculation gives slightly more accurate results than using REPDF, in which the maximum error and RMS of pin power are 2.42% and 0.38%, respectively. The maximum pin power error is located in the MOX fuel or outer UO2 fuel which both are near reflector. PEPDF was therefore selected for the GLIM calculation. Using SP3 method in the local calculation gives more accurate pin power than the diffusion method. Using 7-group or 2-group for the global calculation does not affect the results of SP3 calculation. The results of GLIM-7G and GLIM-2G are similar. The maximum error and RMS error of pin power are 3.33% and 0.46% respectively, slightly larger than the pin-wise results. The CPU time of C5G7 2D problem using NuCoT is shown in Table 6. SP3 calculation took 3–7 times than the diffusion calculation. The GLIM required 30–80% more computing time than the pin-wise calculation using SP3. However, parallelism can easily be applied to the GLIM to make it more efficient in the future development. 3.2.2. 3D problem In C5G7 3D problem, for simplicity, the top reflector was assumed to be pure reflector without control rods or holes. The discontinuity factor of the top reflector were set to one. There are three different control rod insertion cases in 3D problems. In ‘‘Unrodded” case, all rods are out. In ‘‘RoddedA” case, one control rod is inserted into the inner UO2 fuel. In RoddedB case, three con-
trol rods are inserted into the inner UO2 fuel and MOX fuels. PEPDF was used for all 3D cases. The axial nodal size was 21.42 cm. The reference solutions of C5G7 3D problems were the MCNP results using seven group macro cross section set. The statistic error of keff was 3 pcm and of nodal pin power was <1%. The convergence criteria of NuCoT calculation for 3D problems were the same as 2D problem, except the relative power fraction convergence criterion of GLIM was set to 103 due to the convergence difficulty. Because of the rodded conditions, the axial nodal size was shorten from 21.42 cm to 7.14 cm to unify the axial levels of assemblies for the whole core pin-wise calculation in all 3D cases. For GLIM, the nodal height remained to be 21.42 cm for the global calculation. The nodal height used in the local calculation of GLIM in NuCoT is calculated automatically in order to obtain the axial flux distribution within the node considering the homogenized cross sections and axial flux discontinuity factors used in the global calculation. NuCoT sets one third of the nodal height of the global calculation as the maximum nodal height in the local calculation to accommodate the axial heterogeneous lattice distribution such as the partial control rod insertion. In this study, the nodal height in the local calculation is 7.14 cm for C5G7 problems the same as for the pin-wise calculation. Table 7 shows the results of C5G7 3D problem using NuCoT. The comparison of the pin power includes whole pin (integrated pin) and nodal pin (axial length is 21.42 cm). The pin power result of SP3 calculation is more accurate than the diffusion calculation. The accuracy of the GLIM results and pin-wise results are comparable. Although the keff error of GLIM with local calculation using the diffusion method is smaller than using SP3, it is believed to
57
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Fig. 3. Modified C5G7 benchmark problem.
Table 4 Results of C5G7 single assembly benchmark problem by NuCoT pin-wise calculation with and without using discontinuity factors in diffusion and SP3 calculation. Fuel type
Diffusion No DFa
a b c d e
c
SP3 EPDF
No DF
EPDF
UO2 (1.333463)b
80.3 0.61d (1.47)e
1.0 0.00 (0.01)
28.2 0.54 (1.3)
0.6 0.02 (0.05)
UO2cr (0.998971)
2780 1.20 (3.04)
0.5 0.00 (0.01)
1365 1.14 (2.64)
0.2 0.06 (0.15)
MOX (1.186470)
98.8 1.03 (1.86)
0.6 0.00 (0.00)
78.6 0.70 (1.33)
0.1 0.03 (0.06)
MOXcr (0.990750)
876 0.57 (1.49)
0.4 0.00 (0.01)
426 0.40 (1.02)
0.3 0.02 (0.06)
Without using discontinuity factor. Reference kinf. kinf error (pcm). RMS error (%) of pin power. Maximum pin power error (%).
58
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Table 5 Results of 2D whole core C5G7 benchmark problem using the pin-wise calculation and GLIM capability of NuCoT. Reference keff = 1.18645 (1r = 7 105) Pin-wise
Diffusion SP3
GLIMc
Diffusion SP3
a b c d e
keff Error (pcm)
Pin Power Error(%) RMS/Max MOX
Inner UO2
Outer UO2
All
REPDFa PEPDFb REPDF PEPDF
45.9 75.3 51.1 54.8
0.88/4.62 1.28/6.35 0.68/3.01 0.49/2.41
0.36/0.67 0.77/1.70 0.54/1.02 0.19/0.91
1.10/4.01 1.33/4.90 0.78/3.76 0.72/2.42
0.66/4.62 1.08/6.35 0.67/3.76 0.38/2.42
GLIM-7Gd GLIM-2Ge GLIM-7G GLIM-2G
53.9 59.9 79.8 75.1
1.59/6.76 1.46/5.01 0.59/2.81 0.53/2.18
1.60/2.20 1.48/2.18 0.25/1.82 0.27/1.35
2.25/6.10 1.90/4.52 0.82/3.33 0.83/3.44
1.84/6.76 1.69/5.01 0.46/3.33 0.45/3.44
Renormalized EPDF. Prime EPDF. Using prime EPDF. 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
Table 6 CPU time comparison of 2D whole core calculations of C5G7 benchmark problem between using the pin-wise calculation and GLIM capability of NuCoT. Case
CPUa time (sec)
CPU time ratio
Diffusion
Pin-wiseb GLIM-7Gc GLIM-2Gd
8.5 28.3 18.8
1 3.3 2.2
SP3
Pin-wise GLIM-7G GLIM-2G
63.8 85.3 115.0
1 1.3 1.8
a
Intel Core i7 930 2.8 GHz. 7 energy groups for the pin-wise calculation with PEPDF (diffusion or SP3). c 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). d 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion). b
be due to the cancellation of errors since the pin power error is large. More computational time is needed for GLIM calculation as shown in Table 8.
4. Conclusions The Global/Local Iteration Method incorporated in NuCoT using SP3 method together with PEPDF for the local calculation provides quite satisfactory results for the LWR benchmark problems, KAIST and C5G7, more accurate pin power than using the diffusion method. The results of C5G7 and KAIST benchmark problems show that both GLIM and pin-wise calculations in NuCoT give comparable accuracy. For C5G7 problem, compared with MCNP result, the keff error is <200 pcm, the RMS of pin power is < 0.8% and the maximum pin power error is <4% using GLIM. In the pin-wise calculation, the keff error is <160 pcm, the RMS and maximum error of pin power are <0.9% and 2.5% respectively. For KAIST problem, compared with NEWT results, the keff error (<135 pcm), the RMS of pin power (1%) and maximum pin power error (6%) of GLIM calculation are close to the pin-wise calculation results. Currently, in GLIM method, the incoming partial current is used as the boundary condition of the local calculation. It is observed that GLIM-7G result was not better than GLIM-2G result. It may
Table 7 Results of 3D whole core C5G7 benchmark problem using the pin-wise calculation and GLIM capability of NuCoT. Casea
Unrodded 1.14474(1r = 3 105)
keff Error (pcm)
Diffusion
SP3
RoddedA 1.12871(1r = 3 105)
Diffusion
SP3
RoddedB 1.07806(1r = 3 105)
Diffusion
SP3
a b c d e f
Power Error(%) RMS/Max Nodal Pine
Whole Pinf
Pin-wiseb GLIM-7Gc GLIM-2Gd Pin-wise GLIM-7G GLIM-2G
172.3 58.7 57.6 63.1 79.7 85.0
1.17/6.86 1.75/7.24 1.68/5.31 0.51/2.47 0.56/3.65 0.53/3.88
1.43/6.39 1.62/5.97 1.56/4.18 0.59/2.24 0.35/2.51 0.32/2.64
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
125.9 2.5 8.3 116.0 136.6 141.4
1.08/6.32 1.68/6.59 1.57/4.81 0.54/2.39 0.53/3.69 0.51/3.79
0.98/6.17 1.46/5.74 1.37/4.03 0.41/2.28 0.37/2.44 0.35/2.48
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
-129.4 49.6 45.8 156.3 194.8 191.9
1.12/6.43 1.69/5.87 1.45/3.91 0.86/2.22 0.79/3.93 0.76/3.75
1.02/5.94 1.45/5.62 1.24/3.70 0.46/2.10 0.42/2.62 0.40/2.46
Using prime EPDF. 7 energy groups for the pin-wise calculation (diffusion or SP3). 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion). The pin in the node. Nodal pin is axially integrated.
59
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60 Table 8 CPU time comparison of 3D whole core calculations of C5G7 benchmark problem between using the pin-wise calculation and GLIM capability of NuCoT. CPUa time(sec)
CPU time ratio
Pin-wise GLIM-7Gc GLIM-2Gd Pin-wise GLIM-7G GLIM-2G
163.7 551.8 650.6 785.8 1609.2 2135.8
1 3.4 4.0 1 2.0 2.7
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
138.7 565.9 659.6 795.6 1624.2 2191.1
1 4.1 4.8 1 2.0 2.8
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
142.3 669.6 616.3 813.8 1756.2 2239.8
1 4.7 4.3 1 2.2 2.8
Case Unrodded
Diffusion
SP3
RoddedA
Diffusion
SP3
RoddedB
Diffusion
SP3
a b c d
b
Intel Core i7 930 2.8 GHz. 7 energy groups for the pin-wise calculation with PEPDF (diffusion or SP3). 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
be due to the inconsistency of the average surface outgoing partial current (or net current) between the local calculation and the global calculation. Further study is needed to resolve this issue. Although the GLIM calculation takes more computing time than the pin-wise calculation, it is more useful for the current complicated reactor core, such as core loading of different fuel assembly type, partial control rod insertion, and heterogeneity of fuel assembly due to existence of spacers. Parallelism can easily be implemented to the GLIM in the future development to make it more efficient. The SP3 method can be accelerated by the coarse-mesh finite difference method (Yamamoto et al., 2016) to shorten the computing time.
Acknowledgments The authors would like to thank Yung-An Chao for helpful discussions on the even parity flux discontinuity factor. This work was supported by Atomic Energy Council of Taiwan (ROC) of Republic of China.
References Bahadir, T., Lindahl, S., Palmtag, S.P., 2005, ‘‘SIMULATE-4 Multigroup Nodal Code with Microscopic Depletion Model,” M&C 2005, Palais des Papes, Avignon, France, September 12-15, 2005, on CD-ROM, American Nuclear Society, LaGrange Park, IL. ‘‘Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation: A 2-D/3-D MOX Fuel Assembly Benchmark,” NEA/NSC/DOC (2003)16. ‘‘Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation: MOX Fuel Assembly 3-D Extension Case,” NEA/NSC/DOC (2005)16. Chang, C.W., 2011. Exploratory Investigations of Time and Space Scales of Heat Deposition and Conduction for Modeling Doppler Effect in High Temperature Reactors Master thesis. National Tsing-Hua University, Taiwan. Cho, N.Z., 2000, ‘‘Benchmark Problems in Reactor and Particle Transport Physics,” http://nurapt.kaist.ac.kr/benchmark/. Colameco, D.V., 2012. Next Generation Iterative Transport-Diffusion Methodology (ITDM) For LWR Core Analysis Ph.D thesis. The Pennsylvania State University, Pennsylvania. Dall’Osso, A., Tomatis, D., Du Yun, 2010, Improving Cross Sections via Spectral Rehomogenization, PHYSOR 2010, Pittsburgh, Pennsylvania, USA, May 9-14, 2010, on CD-ROM, American Nuclear Society, LaGrange Park, IL.
Gelbard, E.M., 1968. Spherical Harmonics Methods: PL and Double-PL Approximations. In: Greenspan, H., Kelber, C.N., Okrent, D. (Eds.), Computing Methods in Reactor Physics. Gordon and Breach Science Publishers, New York. Hébert Alain, Applied Reactor Physics, Presses inter Polytechnique, pp.101, 2009. Ivanov, B.D., 2007. Methodology for Embedded Transport Core Calculation Ph.D thesis. The Pennsylvania State University, Pennsylvania. Joo, H.G., Cho, J.Y., Kim, K.S., Lee, C.C., Zee, S.Q., 2004, ‘‘Methods and Performance of a Three-Dimensional Whole-Core Transport Code DeCART,” PHYSOR 2004, Chicago, Illinois, April 25-29, 2004, on CD-ROM, American Nuclear Society, Lagrange Park, IL. Jung, Y.S., Shim, C.B., Lim, C.H., Joo, H.G., 2013. Practical numerical reactor employing direct whole core neutron transport and subchannel thermal/ hydraulic solvers. Ann. Nucl. Energy 62, 357–374. Kim, H.R., 1993. Global/Local Iterative Homogenization Methods for Neutron Diffusion Nodal Theory Ph.D thesis. Korea Advanced Institute of Science and Technology, Korea. Knight, M., Bryce, P., Hall, S., 2012, ‘‘WIMS/PANTHER Analysis of UO2/MOX Cores Using Embedded Supercells,” PHYSOR 2012, Knoxville, Tennessee, USA, April 15-20, 2012, on CD-ROM, American Nuclear Society, LaGrange Park, IL. Kooreman, G., Rahnema, F., 2015. Enhanced hybrid diffusion-transport spatial homogenization method. Nucl. Technol. 192, 264–277. Kozlowski, T., Xu, Y., Downar, T.J., Lee, D., 2011. Cell Homogenization Method for Pin-by-Pin Neutron Transport Calculations. Nucl. Sci. Eng. 169, 1–18. Lee, C.H., Downar, T.J., 2004. A hybrid nodal diffusion/SP3 method using one-node coarse-mesh finite difference formulation. Nucl. Sci. Eng. 146, 176–187. Lee, D., Downar, T.J., Kim, Y., 2004. A nodal and finite difference hybrid method for Pin-by-Pin heterogeneous three-dimensional light water reactor diffusion calculations. Nucl. Sci. Eng. 146, 319–339. Lin Tzung-Yi, Liu Yen-Wan Hsueh, 2017, ‘‘Hybrid Nodal Green’s Function Method with SP3 for Pin-by-pin Calculation,” M&C 2017 - International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering, Jeju, Korea, April 16-20, 2017. ‘‘MCNP — A General Monte Carlo N-particle transport code, Version 5,” LA-UR-031987, 2008. Ougouag, A.M., 1981. Coarse-Mesh Nodal Method for Multigroup Multidimensional Neutron Diffusion Computations Master thesis. University of Illinois, Urbana. Roberts, D.R., 2010. Development of Iterative Transport–diffusion Methodology for LWR Analysis Master thesis. The Pennsylvania State University, Pennsylvania. ‘‘Scale: A Comprehensive Modeling and Simulation Suite for Nuclear Safety Analysis and Design,” ORNL/TM-2005/39, Version 6.1, 2011. Smith, K.S., 1986. Assembly homogenization techniques for light water reactor analysis. Prog. Nucl. Energy 17 (3), 303–335. Smith, K.S., 1994. Practical and efficient iterative method for LWR fuel assembly homogenization. Trans. Am. Nucl. Soc. 71, 238–241. Tada, K., Yamamoto, A., Yamane, Y., 2009. Treatment of staggered mesh for BWR pin by pin core analysis. J. Nucl. Sci. Technol. 46 (2), 163–174. Tatsumi, M., Yamamoto, A., Nagano, H., Sengoku, K., 2003. Advanced PWR core calculation based on multi-group nodal-transport method in three-dimensional pin-by-pin geometry. J. Nucl. Sci. Technol. 40 (6), 376–387. Yamamoto, A., Endo, T., Chao, Y.A., 2011, A derivation of discontinuity factor for angular flux in integro-differential transport equation, Transactions of the American Nuclear Society, vol. 104, Hollywood, Florida, June 26–30, 2011.
60
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Yamamoto, A., Tatsumi, M., 2006. Improvement of spatial discretization error on the semi-analytic nodal method using the scattered source subtraction method. J. Nucl. Sci. Technol. 43 (12), 1481–1489. Yamamoto, A., Sakamoto, T., Endo, T., 2016. Discontinuity factors for Simplified P3 theory. Nucl. Sci. Eng. 183, 39–51. Yamamoto, A., Sakamoto, T., Endo, T., 2016. A CMFD acceleration method for SP3 advanced nodal method. Nucl. Sci. Eng. 184, 168–173. Yasseri, S., Rahnema, F., 2014. Consistent spatial homogenization in transport theory. Nucl. Sci. Eng. 176, 292–311.
Yu, L., Chao, Y.A., 2015. A Unified Generic Theory on Discontinuity Factors in Diffusion, SP3 and Transport Calculations. Ann. Nucl. Energy 75, 239–248. Yu, L., Lu, D., Chao, Y.A., 2014. The Calculation Method for SP3 Discontinuity Factor and Its Application. Ann. Nucl. Energy 69, 14–24. Zhang, B., Mayhue, L., Huria, H., Ivanov, B., 2013. Development of a threedimensional pseudo pin-by-pin calculation methodology in ANC. Nucl. Technol. 183, 527–534.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A Next Generation Method for Light Water Reactor core analysis by using Global/Local Iteration Method with SP3 Tzung-Yi Lin ⇑, Yen-Wan Hsueh Liu Institute of Nuclear Engineering and Science, National Tsing Hua University, 101, Section 2, Kuang-Fu Road, Hsinchu 30013, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 2 October 2017 Received in revised form 21 March 2018 Accepted 21 March 2018
Keywords: Next Generation Method (NGM) Simplified P3 method (SP3) Global/Local Iteration Method (GLIM) Even Parity Discontinuity Factor (EPDF)
a b s t r a c t The Global/Local Iteration Method (GLIM) capability was implemented into the in-house developed nodal code, NuCoT, which was based on the Hybrid Nodal Green’s Function Method (HNGFM) with the diffusion/SP3 approximations. A procedure was proposed for obtaining the Even Parity Discontinuity Factor (EPDF) in SP3 method. The accuracy of NuCoT with GLIM was verified by using two benchmark problems, KAIST and C5G7, with reference solutions calculated by NEWT and MCNP respectively. The results of NuCoT show that the accuracy of the GLIM is comparable to the whole core pin-wise calculation. The SP3 calculation with EPDF is better than the diffusion result. For the KAIST problem using GLIM, the Root-Mean-Square (RMS) of pin power error is 1%, the maximum pin power error is 6%, and the keff error is <83 pcm for All-Rod-In case and <135 pcm for All-Rod-Out case. For the C5G7 problem, the RMS of pin power is <0.8% and the maximum error is <4%. The keff error is <80 pcm in 2D problem and <200 pcm in 3D rodded problem. Ó 2018 Elsevier Ltd. All rights reserved.
1. Introduction The two-step method is efficient and widely used for the current Light Water Reactor (LWR) core analysis. Step 1 is single fuel assembly calculation, in which transport method and reflective boundary condition is used. The assembly-homogenized cross sections and assembly discontinuity factors are obtained in this step based on the Generalized Equivalence Theory (Smith, 1986) (GET). Step 2 is the whole core coarse mesh calculation. Normally two group nodal diffusion method is used. However, the reflective boundary used in the first step causes the homogenized parameters used in the second step incorrect for fuel assemblies near the core boundary, and for cores with highly heterogeneous fuel loading. The environment dependence of the assembly cross sections has been added into the homogenization parameters by several different approaches (Smith, 1994; Yamamoto and Tatsumi, 2006; Dall’Osso et al., 2010) with minor change of the two-step method. However, these approaches were challenged when the fuel and core design become more complicated such as the optimized fuel loading pattern, the mixed oxide fuel (MOX) and the diamond type water rod. The axially heterogeneous fuel and
⇑ Corresponding author. E-mail address: [email protected] (T.-Y. Lin). https://doi.org/10.1016/j.anucene.2018.03.030 0306-4549/Ó 2018 Elsevier Ltd. All rights reserved.
partially control rod insertion also affect the accuracy of the local peaking factor (Zhang et al., 2013) calculated by these approaches. There are several approaches proposed for the Next Generation Method of neutronic analysis of LWR to obtain more accuracy results. One is to perform the whole core 3D calculation by 2D(transport)/1D(diffusion) coupling method (Joo et al., 2004; Jung et al., 2013), which is accurate, but not practical considering current computing capability. One of more practical methods is the whole core multi-group 3D pin by pin calculation (pin-wise calculation) (Tatsumi et al., 2003), in which the homogenized cross sections and discontinuity factors of the pin cell rather than the assembly are generated in the first step. In the pin-wise calculation, the pin power reconstruction is not necessary. In order to obtain more accurate result, the Simplified P3 (SP3) (Gelbard, 1968) method is proposed for the pin-wise calculation to replace the diffusion method. However, without discontinuity factor, the solution of SP3 calculation is no better than that of the diffusion calculation. The discontinuity factor for the SP3 calculation has been investigated in many studies (Kozlowski et al., 2011; Yu et al., 2014; Yu and Chao, 2015; Yamamoto et al., 2016). The whole core pin-wise calculation using SP3 has been implemented into NuCoT and showed good accuracy (Lin and Liu, 2017) in the benchmark problems. Another practical method is the Global/Local Iteration Method (Kim, 1993) (GLIM). It is similar to the Iterative TransportDiffusion Methodology (Ivanov, 2007; Roberts, 2010; Colameco,
50
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
2012) (ITDM) or a hybrid method (Kozlowski et al., 2011; Lee and Downar, 2004; Lee et al., 2004; Bahadir et al., 2005; Knight et al., 2012; Yasseri and Rahnema, 2014; Kooreman and Rahnema, 2015). In GLIM, the assembly-homogenized cross sections of the assembly (or lattice) are generated by the transport calculation using the actual assembly surrounding environment. However, since most transport codes can only accept typical boundary conditions such as reflective, white or vacuum boundary, it is not easy to modify the existing transport code for used in the GLIM if one is not familiar with the code. One easier way is to apply the multigroup pin-wise calculation for the assembly calculation. This study presents the GLIM and its application to two benchmark problems, KAIST and C5G7. It includes the derivation of the Hybrid Nodal Green’s Function Method (Ougouag, 1981; Chang, 2011) (HNGFM) with the diffusion and SP3 method for the global and local calculation, the discontinuity factor calculation, and the interface between the global and local calculation.
/gx ðxÞ ¼
The GLIM consists of two parts of calculations using different nodal size. The global calculation is a whole core calculation using the assembly-wise size nodal method, which is usually a fewgroup nodal diffusion method. The local calculation is an assembly calculation using the pin-cell size nodal method, which in this study can be a fine-group diffusion or SP3 nodal method. The local calculation is applied to every assembly nodes and the reflector nodes. The iteration process is shown in Fig. 1. The interface between the global and local calculation called the submesh interface is described in Section 2.4. 2.1. Global calculation using HNGFM with diffusion approximation The steady-state three-dimensional diffusion equation within one node on the Cartesian system for energy group g is typically written in the following form.
X
2 X /gxl Pl ð^xÞ
where Pl ð^xÞ is the l-th order Legendre polynomial and ^x ¼ x=ax . And the transverse leakage profiles were assumed to be flat. The weighted residual procedure was used to solve /gxl using Eq. (2). Then the unsolved term of Eq. (2) at the surface, hQGgx ðax Þi, can be obtained.
hQGgx ðax Þi ¼
2 X
hQ gxl Pl Ggx ðax Þi
ð4Þ
l¼0
where
hPl Ggx ðax Þi
Q gxl
1 2ax
Z
ax
ax
dx Ggx ðax jxÞP l ð^xÞ;
8 S Lgx ; > l¼0 > < gx G G X X ¼ vg f s > > : k 0 mRg0 /g0 xl þ 0 Rgg0 /g0 xl ; l ¼ 1; 2 g –g
g
The continuity of the surface current and discontinuity of the surface flux at the interface of two adjacent nodes provide the rela by using tionship between the current J and nodal average flux / Eq. (2) and 1D nodal balance equation. This relationship can be substituted into the 3D nodal balance equation of node N to obtain Eq. (5).
N BNg / g
Xh
i ðN;m1Þ ðN;m1Þ ðN;mþ1Þ ¼ RN /g / Bgu þ BðN;mþ1Þ gu g g
ð5Þ
u¼x;y;z
where
m ¼ node index related to the u direction; ( BNg
¼R
r;N g
Xh
1þ
ANgu aNu
þ
ANgu aNu
i
) ;
u¼x;y;z
ð1Þ
where
Rrg ¼ Rtg Rsgg ;
ðN;m1Þ Bgu ¼
m1 2am1 Rr;m1 xm1 au j am1 u gu u g ; 2aNu DxNgu ðaNu Þ Xh
RNg ¼ SNg
i C Ngu ðaNu Þ þ C Ngu ðaNu Þ ;
u¼x;y;z
and
Sg ðrÞ ¼
ð3Þ
l¼0
and
2. Global/Local Iteration Method
@2 Dg 2 /g ðrÞ þ Rrg /g ðrÞ ¼ Sg ðrÞ for ðrÞ @u u¼x;y;z 2 ½ax ; ax ; ay ; ay ; ½az ; az
To evaluate hQGgx ðxÞi at the surface, the quadratic polynomial was used to approximate the spatial flux distribution.
SN ¼ vg g k
G G X vg X mRgf 0 /g0 ðrÞ þ Rsgg0 /g0 ðrÞ:
k
g0
g 0 –g
After deriving one-dimensional equation from Eq. (1) and applying the Green’s function with Neumann (or second-type) boundary condition, the 1-D flux in x-direction within a node can be written as:
/gx ðxÞ ¼ Ggx ðxj ax ÞJ gx ðax Þ Ggx ðxjax ÞJ gx ðax Þ þ 2ax hQGgx ðxÞi
ð2Þ
Ggx ðxjx0 Þ = Green’s function at x from x0 with Neumann boundary condition, Q gx ðxÞ ¼ Sgx ðxÞ Lgx ðxÞ; P v P Sgx ðxÞ ¼ kg Gg0 mRgf 0 /g0 x ðxÞ þ Gg0 –g Rsgg 0 /g0 x ðxÞ; and Lgx ðxÞ is the transverse leakage profile.
G X
g0
g 0 –g
mRgf 0 / Ng0 þ
N0 ; Rsgg0 / g
xNgu aNu j aNu ANgu aNu ¼ ; DxNgu ðaNu Þ C Ngu ðaNu Þ ¼
where d J gx ¼ Dg dx /gx ; Ra hQGgx ðxÞi 2a1x ax x dx0 Ggx ðxjx0 ÞQ gx ðx0 Þ;
G X
N ðaN Þ DhQGiNgu ðaNu Þ DQ gu u ; N N 2au Dxgu ðaNu Þ
xNgu ðAjBÞ ¼ f Ngu ðAÞGNgu ðAjBÞ; N
f gu ðAÞ ¼ Discontinuity factor at the location A in u direction; h i DxNgu ðaNu Þ ¼ xNgu aNu j aNu xNgu aNu j aNu h m1 m1 i xm1 au j am1 au j am1 xm1 ; gu u gu u
51
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Fig. 1. The flow of the Global/Local Iteration Method.
N N N ðaN Þ ¼ 2aN Q N N DQ gu u u gu xgu au j au m1 xm1 am1 j am1 ; 2am1 Q u
gu
gu
u
2.2. Local calculation using HNGFM with SP3 approximation
u
The typical SP3 equation of group g of a node on the Cartesian system is given below:
N ¼ SN LN ; Q gu g gu
LN ¼ gu
X
@2 F ng ðrÞ þ Rrng F ng ðrÞ ¼ T ng ðrÞ; 2 @u u¼x;y;z for r 2 ½ax ; ax ; ½ay ; ay ; ½az ; az
i X 1 hN J gh ðah Þ J Ngh ðah Þ ; 2ah h–u
Dng
ð6Þ
where
D1g ¼
and N
n ¼ 1; 2
m1
hDQGgu iN ðaNu Þ ¼ 2aNu f gu ðaNu ÞhQGNgu ðaNu Þi 2am1 f gu ðam1 Þ u u m1 Þi: hQGm1 gu ðau
The available boundary conditions are albedo, zero flux and known incoming partial current. The surface current and surface flux can be related by the albedo. After using Eq. (2) and 1D nodal balance equation, the surface current on the boundary can be related to its nodal average flux to obtain a modified Eq. (5) for the node on the boundary. If the incoming partial current is known, it can be treated as an external source with albedo set to zero.
1 ; 3Rtrg
D2g ¼
1 ; 7Rtg
Rr1g ¼ Rtg Rsgg ;
2 T 1g ðrÞ ¼ Sg ðrÞ þ Rr1g F 2g ðrÞ; 3 Sg ðrÞ ¼
2 2 T 2g ðrÞ ¼ Sg ðrÞ þ Rr1g F 1g ðrÞ; 3 3
G G X vg X mRgf 0 /g0 0 ðrÞ þ Rsgg0 /g0 0 ðrÞ
k
g0
4 9
Rr2g ¼ Rtg Rsgg ;
g 0 –g
F 1g ðrÞ ¼ /g0 ðrÞ þ 2/g2 ðrÞ; F 2g ðrÞ ¼ 3/g2 ðrÞ; and /g0 ðrÞ and /g2 ðrÞ are the flux moments.
52
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Similar to the diffusion approximation, the one-dimensional pseudo fluxes in x-direction within a node can be written as
C Nngu ðaNu Þ ¼
N ðaN Þ DhQGiNngu ðaNu Þ DQ ngu u : 2aNu DxNngu ðaNu Þ
F ngx ðxÞ ¼ Gngx ðxj ax ÞJ ngx ðax Þ Gngx ðxjax ÞJ ngx ðax Þ þ 2ax hQGngx ðxÞi n ¼ 1;2 ð7Þ
where
xNngu ðAjBÞ ¼ f Nngu ðAÞGNngu ðAjBÞ; N
f ngu ðAÞ ¼ Discontinuity factor at the location A in u direction;
d J ngx ¼ Dng F ngx ; dx
m1 DxNngu ðaNu Þ ¼ ½xNngu ðaNu j aNu Þ xNngu ðaNu j aNu Þ ½xm1 j ngu ðau
Gngx ðxjx0 Þ ¼ Green’s function at x0 with Neumann boundary condition;
hQGngx ðax Þi
1 2ax
Z
ax
ax
m1 Þ xm1 j am1 Þ; am1 u ngu ðau u N N N m1 m1 m1 N ðaN Þ ¼ 2aN Q N Q ngu xngu ðam1 DQ j ngu u u ngu xngu ðau j au Þ 2au u
dx0 Gngx ðax jx0 ÞQ ngx ðx0 Þ
Þ; am1 u
Q ngx ðxÞ ¼ T ngx ðxÞ Lngx ðxÞ; N ¼ T N LN ; Q ngu ng ngu
2 T 1g ðxÞ ¼ Sgx ðxÞ þ Rr1g F 2gx ðxÞ; 3
LN ¼ ngu
2 2 T 2g ðxÞ ¼ Sgx ðxÞ þ Rr1g F 1gx ðxÞ; 3 3 Sgx ðxÞ ¼
vg
G X
G X
k
g0
g 0 –g
mRgf 0 /g0 0x ðxÞ þ
and N
Xh
i ðN;m1Þ þ BðN;mþ1Þ F ðN;mþ1Þ ¼ RN ; BðN;m1Þ F ngu ngu ng ng ng
m1 Þi: hQGm1 ngu ðau
n ¼ 1; 2
u¼x;y;z
ð8Þ where
m ¼ node index related to the u direction; ( BNng ¼ Rr;N 1þ ng
Xh
) i ; ANngu ðaNu Þ þ ANngu aNu
u¼x;y;z
BðN;m1Þ ¼ ngu
m1 2am1 Rr;m1 xm1 j am1 Þ u ngu ðau u ng ; N N N 2au nxngu ðau Þ
Xh
RNng ¼ T Nng
i C Nngu ðaNu Þ þ C Nngu ðaNu Þ ;
u¼x;y;z
2 N ; T N1g ¼ SNg þ Rr1g F 2g 3 SN ¼ vg g k
2 2 N ; T N2g ¼ SNg þ Rr1g F 1g 3 3
G X
G X
g0
g 0 –g
mRgf 0 / Ng0 0 þ
ANngu ðaNu Þ ¼
x
N0 ; Rsgg0 / g0
N N N ngu ðau j au Þ ; N N D ngu ðau Þ
x
m1
DhQGiNngu ðaNu Þ ¼ 2aNu f ngu ðaNu ÞhQGNngu ðaNu Þi 2am1 f ngu ðam1 Þ u u
Rsgg0 /g0 0x ðxÞ;
and Lngx ðxÞ is the transverse leakage profile which was assumed to be flat. The unsolved term of Eq. (7) at the surface, hQGngx ðax Þi, can be obtained by the same procedure as in the diffusion method. Also, the relationship between the surface current J n and nodal average n was used to obtain Eq. (8) by utilizing the continupseudo flux F ity of the surface current ðJ 1 and J 2 Þ and discontinuity of the surface pseudo flux ðF 1 and F 2 Þ, Eq. (7), at the interface of two adjacent nodes, following the same procedure for obtaining Eq. (5) in the diffusion method.
N BNng F ng
i X 1 h J Nngh ðah Þ J Nngh ðah Þ ; 2a h h–u
The boundary conditions available in SP3 calculation are albedo, zero pseudo flux and known incoming partial current. The partial current moments (J 1 and J 2 ), based on the Marshak coefficient (Hébert, 2009) are:
(
1 J 1gu ðuÞ ¼ 14 F 1gu ðuÞ 16 F 2gu ðuÞ 12 J 1gu ðuÞ 1 7 F 1gu ðuÞ þ 48 F 2gu ðuÞ 12 J 2gu ðuÞ J 2gu ðuÞ ¼ 16
ð9Þ
Therefore, for the albedo boundary condition, the surface currents (J 1 and J 2 ) and the surface pseudo fluxes (F 1 and F 2 ) can be related by the albedo. Then a modified Eq. (8) for the boundary node can be obtained by using the similar procedure as in the diffusion method. 2.3. Even Parity Discontinuity Factor (EPDF) In SP3 calculation, the discontinuity factors of two pseudo fluxes cannot be easily generated because the 2nd and 3rd angular flux moments were not available from the widely used transport codes. EPDF was developed for the transport code by Yamamoto et al. (2011). It was unified in the diffusion/SP3 calculation by Yu and Chao (2015). The idea is to make the homogeneous even parity surface flux Uhom (not the scalar flux /hom ) between two adjacent 0 nodes continuous after applying EPDF. The Uhom on the left- and right-hand side of the interface are continuous by introducing the factor d. hom F Right Uhom Right ¼ F Left ULeft
ð10Þ
where
F Right ¼ 1 d;
F Left ¼ 1 þ d;
and
d¼
hom Uhom Right ULeft hom Uhom Right þ ULeft
:
The even parity surface flux can be rewritten as the average value of the partial incoming and outgoing currents.
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
U¼
J in þ J out 2
In diffusion approximation, the even parity flux was:
1 4
U ffi /0
ð11aÞ
In SP3 approximation:
1 4
U ffi /0 þ
5 1 1 / ¼ F1 F2 16 2 4 16
or Left
¼
U
U
hom Right or Left
2.4. Submesh interface program
ð13Þ
where
U ¼
hom Uhom Right þ ULeft
2
:
After applying EPDF, the homogeneous even parity surface fluxes are made continuous. In the single fuel assembly transport calculation using reflective boundary condition implies that the neighboring assemblies are the same. Hence, EPDF on the assembly boundary is one. However, the neighboring fuel assemblies in the core may not be the same, in which case EPDF will not be equal to one. Yu and Chao (2015) used the same idea of the conventional discontinuity factor to set U of Eq. (13) to the reference solution Uhet , named Renormalized EPDF (REPDF). The original EPDF definition was named Prime EPDF (PEPDF). The difficulty was to obtain the reference 2nd flux moment (F 2 or /2 ) and current (J 2 or /3 ) for the homogeneous even parity surface flux of the SP3 approximation of HNGFM. In this study, the procedure for obtaining the EPDF for SP3 calculation could have followed Yu et. al.’s idea to iterate EPDF utilizing the available information in the reference solutions, such as surface current J 1 , nodal 0 and keff. However, an alternative method using less average flux / iteration parameters is adopted. The two equations of nodal balance equations from Eq. (6) can be reduced to Eq. (14) by substituting the reference solution of surface current J 1 and nodal average 0 . flux /
X
1 2g ¼ 2 Hg J 2gu ðau Þ J 2gu ðau Þ þ Rt2g F 2a 3 u u¼x;y;z
ð14Þ
where
5 9
Rt2g ¼ Rtg ; Hg ¼
Eq. (14) also appeared in Kozlowski et al. (2011) (Eqs. (28)– (30)), where it was used in a different way but with approximations to calculate the discontinuity factors. It is easy to obtain the EPDF for the diffusion approximation. When the surface net current, nodal average flux and keff are from the reference solution, the only one unknown of Eq. (2) for calculating the homogeneous surface flux is QGgx ðax Þ which needs an iteration to obtain the moment terms of Eq. (4).
ð12bÞ
In this study, Eq. (10) was rewritten to the form similar to the conventional discontinuity factor:
EPDF Right
53
X
1 J 1gu ðau Þ J 1gu ðau Þ 2a u u¼x;y;z
Therefore, a similar form of Eq. (8) can be easily obtained and only 2 matrix is necessarily to be solved. F The iterative procedure for obtaining EPDF is described below: Assuming EPDF is one. 2 matrix from the similar form of Eq. (8) and (1) Solve the F obtain J 2 . (2) Obtain two pseudo surface fluxes by solving Eq. (7), which needs another iteration for QGngx ðau Þ. (3) Update EPDF by using Eq. (13). (4) Repeat step (1)–(3) until converged.
The submesh interface program connects the global and local calculation. The data requirements of the global calculation are the assembly-homogenized cross sections and the assembly discontinuity factors. In the traditional two-step procedure of the core analysis, the transport calculation generates these parameters by using reflective boundary condition. In the GLIM of this study, these parameters were generated by the submesh interface program from the results of the local calculation using boundary conditions based on the global calculation. This local calculation result can be regarded as the heterogeneous solution for the assemblyhomogenized cross sections and assembly discontinuity factors. The assembly discontinuity factors were calculated by the method described at the end of Section 2.3 for the global calculation using the diffusion method. Based on the theory of EPDF, the homogenous even parity surface fluxes on the interface between the left- and right-side nodes (Eq. (10)) were needed at the same time. However, the homogeneous even parity surface flux of the neighboring node at the assembly boundary was unknown because the local calculation was a single assembly calculation. Therefore, REPDF was applied for the global calculation, not PEPDF. The homogenized cross sections and the discontinuity factors of the pin cell used for the local calculation were input data. The pin cell homogenized cross sections were obtained from the single assembly transport calculation with the reflective boundary condition. The pin cell discontinuity factors for the diffusion or SP3 method were calculated in advance by the method described in Section 2.3. The keff for the local calculation was given by the global calculation, instead of being calculated by using source iteration. In other words, the local calculation was a fixed-k calculation. Since the GLIM was an iteration process between the global and local calculation, the boundary condition for the local calculation should be updated after the global calculation using new cross sections and discontinuity factors. The realistic boundary condition (Ivanov, 2007; Roberts, 2010; Colameco, 2012) for the local calculation was the incoming neutron current distribution which was just the outgoing neutron current distribution of the adjacent assembly from the local calculation. In order to apply new global calculation result (step N) for the local calculation, the outgoing current distribution of the adjacent assembly from the previous local calculation (step N-1) was normalized to one. Then the newly averaged outgoing current of the adjacent assembly of the global calculation (step N) multiplied by this normalized outgoing current distribution became the new incoming current distribution of the calculated assembly. In the local calculation using SP3 method, since the second kind of averaged outgoing current did not exist in the global calculation in which diffusion method is used, the second kind of outgoing current distribution was normalized to the averaged outgoing current (first kind) by assuming the ratios of these two kinds of averaged outgoing current were the same. One of the advantages in the GLIM is the size of pin cell of the two neighboring assemblies can be different. By adding additional continuity conditions on the surfaces of the two neighboring assemblies with different pin cell size, the nodal method can be used without further modification. The submesh interface program
54
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
is used to generate the mapped incoming current from the neighboring assembly. In Tada et. al.’s study (Tada et al., 2009), the incoming current can be approximated by a flat, linear or quadratic function in 2D problem. For 3D problems, the incoming current is on a 2D plane, not a 1D line. The linear or quadratic function becomes complicated. The simplest way is to use the areaweighted calculation in this study, i.e. the flat approximation. The process of GLIM can be terminated by examining convergence of parameters such as the surface current or averaged nodal flux. In this study, the submesh interface program used the convergence of the relative nodal power fraction and keff of the global calculation as the stopping criteria. 3. Application of the GLIM The GLIM was implemented into our in-house code, NuCoT. Two well-known benchmark problems, KAIST (Cho, 2000) and C5G7 (Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation, 2003, 2005), were selected to verify its accuracy. The source iteration used in NuCoT was accelerated by the Two-Mode Asymptotic Source Extrapolation (Ougouag, 1981). Root-Mean-Square (RMS) error was used to compare the pin power result of NuCoT with the reference solution.
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P ðP Pref Þ2 =N P RMSð%Þ ¼ 100 Pref =N
ð15Þ
where N = total number of the fuel pin, P = pin power, and Pref=pin power of reference solution. Pin-wise results are included for comparison. In case ‘‘PEPDF”, PEPDF was applied to all surfaces of the pin cells of the local calculation. In case ‘‘REPDF”, REPDF was applied to the surfaces of the pin cells on the assembly/reflector boundary and PEPDF was applied to the rest of the surfaces in the local calculation. The accuracy of the diffusion method and SP3 method in the local calculation was compared. 3.1. KAIST benchmark problem KAIST benchmark problem is a small 2D quarter PWR core consists of three types of UO2 and two types of zoned MOX 17 17 fuels. There are two core conditions, all rods out (ARO) and all rods in (ARI), as shown in Fig. 2. The seven-group macro cross section set including anisotropic scattering cross section was provided by the KAIST 2A benchmark problem (Cho, 2000). The KAIST benchmark problem consists of six types of fuel assemblies, 2.0% UO2 fuel (UOX1), 3.3% UO2 fuel with or without control rod insertion (UOX2 or UOX2cr), UOX2 fuel containing 16 gadolinium (Gd) rods (UOX2ba16), zoned MOX fuel (MOX1) and a MOX1 fuel with 8 Gd rods (MOX1ba8). The reference solution for obtaining EPDF and homogenized pin-cell cross sections of the single fuel assembly were performed by using NEWT of SCALE 6.1 (Scale: A Comprehensive Modeling and Simulation Suite for Nuclear Safety Analysis and Design, 2011). The angular quadrature was S16, and the convergence criteria of keff, flux and source were all 105. The accuracy of NuCoT results using diffusion and SP3 approximation these six fuel types of assembly are shown in Table 1. The convergence criteria for the keff, source and flux were all 105. As expected, without applying discontinuity factors, the error in kinf was large. After applying
discontinuity factors, the NuCoT results were very close to the reference results. The reference solution of the whole core calculation was also obtained by using NEWT. The angular quadrature and the convergence criteria used were the same as in the single assembly calculation. In the NuCoT whole core pin-wise calculation, two pin-cell discontinuity factor types, PEPDF and REPDF were applied. The convergence criteria were all 105. The PEPDF and REPDF of pin cells in the fuel assembly were obtained by the single assembly calculation. The reflector including the baffle was divided into 17 17 pin cells as the fuel assembly. Since the pin cells of the reflector were homogeneous, the pin-cell discontinuity factors were set to one, except the one at the reflector boundary adjacent to the heterogeneous fuel assembly. The accuracy of the pin-wise and the GLIM capabilities of NuCoT are shown in Table 2. In the ARO case, the results are insensitive to the choice of discontinuity factors. In the ARI case, using PEPDF gives better results than using REPDF. The SP3 results are better than the diffusion results. In SP3 calculation using PEPDF, the maximum pin power error is <6% and the RMS of pin power error is 1%. The location of the maximum pin power error is near the reflector and its power is much smaller than one. The keff error is <70 pcm for ARI case and <130 pcm for ARO case. PEPDF was selected for the following GLIM calculation. In the NuCoT GLIM calculation, the convergence criteria for the keff, source and flux of the global and local calculations were all 105. The convergence criteria for the keff and relative power fraction in the submesh interface were 104 and 5x103, respectively. The energy groups used in the local calculation was seven. In the global calculation, both two and seven energy groups are used (named GLIM-2G or GLIM-7G). As shown in Table 2, the accuracy of GLIM calculation is comparable to the pin-wise calculation. The RMS of pin power and error of keff are slightly larger than the results of direct pin-by-pin calculations. The GLIM-7G result was not better than GLIM-2G result. It may be caused by the inconsistency of the average surface outgoing partial current (or net current) between the local calculation and the global calculation since it is not one of the convergent criteria. The computing time of GLIM was 2–3 times of the pin-wise calculation as shown in Table 3. Since the local calculations of each assembly-size node are independent, parallelism can be applied to the GLIM to reduce computing time.
3.2. Results of C5G7 benchmark problem C5G7 benchmark problem is a mini quarter PWR core containing UO2 and zoned MOX 17 17 fuel assemblies, surrounded by the reflector. This benchmark problem includes one 2D case and three 3D cases, Unrodded, RoddedA and RoddedB, as shown in Fig. 3. The fuel patterns of 3D cases are the same as the 2D case. The reference solution for EPDF iteration and homogenized pincell cross sections of the fuel assembly were obtained by the single assembly calculation using NEWT using S16 angular quadrature set. The convergence criteria for keff, flux and source were 105. The seven group macro cross section set, including the transportcorrected total and scattering cross sections, were provided by the benchmark problem (Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation, 2003, 2005). Table 4. shows the results of NuCoT pin-wise single assembly calculation compared with the reference solution for the fuel types, UO2, MOX and the rodded fuels (UO2cr and MOXcr), using both the diffusion and SP3 approximation. NuCoT results obtained by using discontinuity factors agree well with the reference results.
55
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Fig. 2. KAIST benchmark problem (2A).
Table 1 Results of KAIST single assembly benchmark problem by NuCoT pin-wise calculation with and without using discontinuity factors in diffusion and SP3 calculation. Fuel type
a b c d e
Diffusion
SP3
No DFa
EPDF
No DF
EPDF
MOX1 (1.176597)b
69.2c 0.67d (1.17)e
0.7 0.00 (0.01)
44.5 0.36 (0.83)
0.5 0.03 (0.07)
MOX1ba8 (1.136757)
333 0.78 (4.07)
0.3 0.00 (0.01)
240 0.48 (2.88)
0.5 0.04 (0.10)
UOX1 (1.100852)
72.9 0.30 (0.72)
4.0 0.00 (0.01)
34.4 0.25 (0.59)
3.7 0.01 (0.02)
UOX2 (1.252353)
84.4 0.36 (0.86)
4.3 0.00 (0.00)
40.9 0.30 (0.68)
4.6 0.01 (0.03)
UOX2ba16 (1.043023)
2459 1.47 (8.24)
4.2 0.00 (0.00)
1426 1.23 (4.57)
7.7 0.09 (0.22)
UOX2cr (0.818892)
4927 2.37 (4.66)
1.7 0.00 (0.00)
2596 1.90 (3.52)
5.9 0.14 (0.42)
Without using discontinuity factor. Reference kinf. kinf error (pcm). RMS error (%) of pin power. Maximum pin power error (%).
3.2.1. 2D problem Table 5 shows the accuracy of NuCoT results of C5G7 2D problem using the pin-wise and GLIM calculations. The reflector of the 2D problem was divided into 17 17 cells as the fuel assembly. Since the reflector is homogeneous, the pin-cell discontinuity factor of the reflector was applied by the same approach as for the KAIST benchmark problem. The reference solution of C5G7 2D
problem was obtained from MCNP (MNCP, 2008) calculation using the seven-group macro cross section set provided by the benchmark problem (Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation, 2003, 2005). The statistic error of keff was 7 pcm and of pin power was <0.5%. The convergence criteria of NuCoT for keff, flux and power were all 105.
56
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Table 2 Results of 2D whole core KAIST benchmark problem using the pin-wise calculation and GLIM capability of NuCoT. Case
Pin-wise
ARO Reference keff = 1.136577
Diffusion SP3
GLIMc
Diffusion SP3
a b c d e
ARI Reference keff = 0.983737
keff Error (pcm)
Pin Power Error(%) RMS/Max
keff Error (pcm)
Pin Power Error(%) RMS/Max
REPDFa PEPDFb REPDF PEPDF
190.8 196.0 142.7 129.8
1.15/5.46 0.97/5.07 1.05/4.75 0.78/5.46
314.1 229.2 140.3 67.8
2.49/5.28 1.73/4.34 1.47/5.03 1.10/5.56
GLIM-7Gd GLIM-2Ge GLIM-7G GLIM-2G
197.0 194.7 134.3 129.9
0.98/5.01 0.97/4.49 0.92/6.02 0.88/5.55
247.0 242.7 76.7 82.6
1.82/4.92 1.79/4.52 1.23/6.02 1.21/5.54
Renormalized EPDF. Prime EPDF. Using prime EPDF. 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
Table 3 CPU time comparison of 2D whole core calculations of KAIST benchmark problem between using the pin-wise calculation and GLIM capability of NuCoT. Case
a b c d
ARO
ARI
CPUa time (sec)
CPU time ratio
CPU time (sec)
CPU time ratio
Diffusion
Pin-wiseb GLIM-7Gc GLIM-2Gd
96.1 172.0 287.7
1 1.8 3.0
73.3 150.9 474.4
1 2.1 6.5
SP3
Pin-wise GLIM-7G GLIM-2G
264.4 893.3 885.0
1 3.4 3.3
273.3 677.0 568.4
1 2.5 2.1
Intel Core i7 930 2.8 GHz. 7 energy groups for the pin-wise calculation with PEPDF (diffusion or SP3). 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
In the NuCoT pin-wise calculation, the keff errors of the diffusion and SP3 method with either PEPDF or REPDF discontinuity factor types are all <100 pcm. The pin power result of the SP3 method is more accurate than the diffusion method. Using PEPDF in SP3 calculation gives slightly more accurate results than using REPDF, in which the maximum error and RMS of pin power are 2.42% and 0.38%, respectively. The maximum pin power error is located in the MOX fuel or outer UO2 fuel which both are near reflector. PEPDF was therefore selected for the GLIM calculation. Using SP3 method in the local calculation gives more accurate pin power than the diffusion method. Using 7-group or 2-group for the global calculation does not affect the results of SP3 calculation. The results of GLIM-7G and GLIM-2G are similar. The maximum error and RMS error of pin power are 3.33% and 0.46% respectively, slightly larger than the pin-wise results. The CPU time of C5G7 2D problem using NuCoT is shown in Table 6. SP3 calculation took 3–7 times than the diffusion calculation. The GLIM required 30–80% more computing time than the pin-wise calculation using SP3. However, parallelism can easily be applied to the GLIM to make it more efficient in the future development. 3.2.2. 3D problem In C5G7 3D problem, for simplicity, the top reflector was assumed to be pure reflector without control rods or holes. The discontinuity factor of the top reflector were set to one. There are three different control rod insertion cases in 3D problems. In ‘‘Unrodded” case, all rods are out. In ‘‘RoddedA” case, one control rod is inserted into the inner UO2 fuel. In RoddedB case, three con-
trol rods are inserted into the inner UO2 fuel and MOX fuels. PEPDF was used for all 3D cases. The axial nodal size was 21.42 cm. The reference solutions of C5G7 3D problems were the MCNP results using seven group macro cross section set. The statistic error of keff was 3 pcm and of nodal pin power was <1%. The convergence criteria of NuCoT calculation for 3D problems were the same as 2D problem, except the relative power fraction convergence criterion of GLIM was set to 103 due to the convergence difficulty. Because of the rodded conditions, the axial nodal size was shorten from 21.42 cm to 7.14 cm to unify the axial levels of assemblies for the whole core pin-wise calculation in all 3D cases. For GLIM, the nodal height remained to be 21.42 cm for the global calculation. The nodal height used in the local calculation of GLIM in NuCoT is calculated automatically in order to obtain the axial flux distribution within the node considering the homogenized cross sections and axial flux discontinuity factors used in the global calculation. NuCoT sets one third of the nodal height of the global calculation as the maximum nodal height in the local calculation to accommodate the axial heterogeneous lattice distribution such as the partial control rod insertion. In this study, the nodal height in the local calculation is 7.14 cm for C5G7 problems the same as for the pin-wise calculation. Table 7 shows the results of C5G7 3D problem using NuCoT. The comparison of the pin power includes whole pin (integrated pin) and nodal pin (axial length is 21.42 cm). The pin power result of SP3 calculation is more accurate than the diffusion calculation. The accuracy of the GLIM results and pin-wise results are comparable. Although the keff error of GLIM with local calculation using the diffusion method is smaller than using SP3, it is believed to
57
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Fig. 3. Modified C5G7 benchmark problem.
Table 4 Results of C5G7 single assembly benchmark problem by NuCoT pin-wise calculation with and without using discontinuity factors in diffusion and SP3 calculation. Fuel type
Diffusion No DFa
a b c d e
c
SP3 EPDF
No DF
EPDF
UO2 (1.333463)b
80.3 0.61d (1.47)e
1.0 0.00 (0.01)
28.2 0.54 (1.3)
0.6 0.02 (0.05)
UO2cr (0.998971)
2780 1.20 (3.04)
0.5 0.00 (0.01)
1365 1.14 (2.64)
0.2 0.06 (0.15)
MOX (1.186470)
98.8 1.03 (1.86)
0.6 0.00 (0.00)
78.6 0.70 (1.33)
0.1 0.03 (0.06)
MOXcr (0.990750)
876 0.57 (1.49)
0.4 0.00 (0.01)
426 0.40 (1.02)
0.3 0.02 (0.06)
Without using discontinuity factor. Reference kinf. kinf error (pcm). RMS error (%) of pin power. Maximum pin power error (%).
58
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Table 5 Results of 2D whole core C5G7 benchmark problem using the pin-wise calculation and GLIM capability of NuCoT. Reference keff = 1.18645 (1r = 7 105) Pin-wise
Diffusion SP3
GLIMc
Diffusion SP3
a b c d e
keff Error (pcm)
Pin Power Error(%) RMS/Max MOX
Inner UO2
Outer UO2
All
REPDFa PEPDFb REPDF PEPDF
45.9 75.3 51.1 54.8
0.88/4.62 1.28/6.35 0.68/3.01 0.49/2.41
0.36/0.67 0.77/1.70 0.54/1.02 0.19/0.91
1.10/4.01 1.33/4.90 0.78/3.76 0.72/2.42
0.66/4.62 1.08/6.35 0.67/3.76 0.38/2.42
GLIM-7Gd GLIM-2Ge GLIM-7G GLIM-2G
53.9 59.9 79.8 75.1
1.59/6.76 1.46/5.01 0.59/2.81 0.53/2.18
1.60/2.20 1.48/2.18 0.25/1.82 0.27/1.35
2.25/6.10 1.90/4.52 0.82/3.33 0.83/3.44
1.84/6.76 1.69/5.01 0.46/3.33 0.45/3.44
Renormalized EPDF. Prime EPDF. Using prime EPDF. 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
Table 6 CPU time comparison of 2D whole core calculations of C5G7 benchmark problem between using the pin-wise calculation and GLIM capability of NuCoT. Case
CPUa time (sec)
CPU time ratio
Diffusion
Pin-wiseb GLIM-7Gc GLIM-2Gd
8.5 28.3 18.8
1 3.3 2.2
SP3
Pin-wise GLIM-7G GLIM-2G
63.8 85.3 115.0
1 1.3 1.8
a
Intel Core i7 930 2.8 GHz. 7 energy groups for the pin-wise calculation with PEPDF (diffusion or SP3). c 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). d 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion). b
be due to the cancellation of errors since the pin power error is large. More computational time is needed for GLIM calculation as shown in Table 8.
4. Conclusions The Global/Local Iteration Method incorporated in NuCoT using SP3 method together with PEPDF for the local calculation provides quite satisfactory results for the LWR benchmark problems, KAIST and C5G7, more accurate pin power than using the diffusion method. The results of C5G7 and KAIST benchmark problems show that both GLIM and pin-wise calculations in NuCoT give comparable accuracy. For C5G7 problem, compared with MCNP result, the keff error is <200 pcm, the RMS of pin power is < 0.8% and the maximum pin power error is <4% using GLIM. In the pin-wise calculation, the keff error is <160 pcm, the RMS and maximum error of pin power are <0.9% and 2.5% respectively. For KAIST problem, compared with NEWT results, the keff error (<135 pcm), the RMS of pin power (1%) and maximum pin power error (6%) of GLIM calculation are close to the pin-wise calculation results. Currently, in GLIM method, the incoming partial current is used as the boundary condition of the local calculation. It is observed that GLIM-7G result was not better than GLIM-2G result. It may
Table 7 Results of 3D whole core C5G7 benchmark problem using the pin-wise calculation and GLIM capability of NuCoT. Casea
Unrodded 1.14474(1r = 3 105)
keff Error (pcm)
Diffusion
SP3
RoddedA 1.12871(1r = 3 105)
Diffusion
SP3
RoddedB 1.07806(1r = 3 105)
Diffusion
SP3
a b c d e f
Power Error(%) RMS/Max Nodal Pine
Whole Pinf
Pin-wiseb GLIM-7Gc GLIM-2Gd Pin-wise GLIM-7G GLIM-2G
172.3 58.7 57.6 63.1 79.7 85.0
1.17/6.86 1.75/7.24 1.68/5.31 0.51/2.47 0.56/3.65 0.53/3.88
1.43/6.39 1.62/5.97 1.56/4.18 0.59/2.24 0.35/2.51 0.32/2.64
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
125.9 2.5 8.3 116.0 136.6 141.4
1.08/6.32 1.68/6.59 1.57/4.81 0.54/2.39 0.53/3.69 0.51/3.79
0.98/6.17 1.46/5.74 1.37/4.03 0.41/2.28 0.37/2.44 0.35/2.48
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
-129.4 49.6 45.8 156.3 194.8 191.9
1.12/6.43 1.69/5.87 1.45/3.91 0.86/2.22 0.79/3.93 0.76/3.75
1.02/5.94 1.45/5.62 1.24/3.70 0.46/2.10 0.42/2.62 0.40/2.46
Using prime EPDF. 7 energy groups for the pin-wise calculation (diffusion or SP3). 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion). The pin in the node. Nodal pin is axially integrated.
59
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60 Table 8 CPU time comparison of 3D whole core calculations of C5G7 benchmark problem between using the pin-wise calculation and GLIM capability of NuCoT. CPUa time(sec)
CPU time ratio
Pin-wise GLIM-7Gc GLIM-2Gd Pin-wise GLIM-7G GLIM-2G
163.7 551.8 650.6 785.8 1609.2 2135.8
1 3.4 4.0 1 2.0 2.7
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
138.7 565.9 659.6 795.6 1624.2 2191.1
1 4.1 4.8 1 2.0 2.8
Pin-wise GLIM-7G GLIM-2G Pin-wise GLIM-7G GLIM-2G
142.3 669.6 616.3 813.8 1756.2 2239.8
1 4.7 4.3 1 2.2 2.8
Case Unrodded
Diffusion
SP3
RoddedA
Diffusion
SP3
RoddedB
Diffusion
SP3
a b c d
b
Intel Core i7 930 2.8 GHz. 7 energy groups for the pin-wise calculation with PEPDF (diffusion or SP3). 7 energy groups for the local calculation (diffusion or SP3) and the global calculation (diffusion). 7 energy groups for the local calculation (diffusion or SP3) and 2 energy groups for the global calculation (diffusion).
be due to the inconsistency of the average surface outgoing partial current (or net current) between the local calculation and the global calculation. Further study is needed to resolve this issue. Although the GLIM calculation takes more computing time than the pin-wise calculation, it is more useful for the current complicated reactor core, such as core loading of different fuel assembly type, partial control rod insertion, and heterogeneity of fuel assembly due to existence of spacers. Parallelism can easily be implemented to the GLIM in the future development to make it more efficient. The SP3 method can be accelerated by the coarse-mesh finite difference method (Yamamoto et al., 2016) to shorten the computing time.
Acknowledgments The authors would like to thank Yung-An Chao for helpful discussions on the even parity flux discontinuity factor. This work was supported by Atomic Energy Council of Taiwan (ROC) of Republic of China.
References Bahadir, T., Lindahl, S., Palmtag, S.P., 2005, ‘‘SIMULATE-4 Multigroup Nodal Code with Microscopic Depletion Model,” M&C 2005, Palais des Papes, Avignon, France, September 12-15, 2005, on CD-ROM, American Nuclear Society, LaGrange Park, IL. ‘‘Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation: A 2-D/3-D MOX Fuel Assembly Benchmark,” NEA/NSC/DOC (2003)16. ‘‘Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation: MOX Fuel Assembly 3-D Extension Case,” NEA/NSC/DOC (2005)16. Chang, C.W., 2011. Exploratory Investigations of Time and Space Scales of Heat Deposition and Conduction for Modeling Doppler Effect in High Temperature Reactors Master thesis. National Tsing-Hua University, Taiwan. Cho, N.Z., 2000, ‘‘Benchmark Problems in Reactor and Particle Transport Physics,” http://nurapt.kaist.ac.kr/benchmark/. Colameco, D.V., 2012. Next Generation Iterative Transport-Diffusion Methodology (ITDM) For LWR Core Analysis Ph.D thesis. The Pennsylvania State University, Pennsylvania. Dall’Osso, A., Tomatis, D., Du Yun, 2010, Improving Cross Sections via Spectral Rehomogenization, PHYSOR 2010, Pittsburgh, Pennsylvania, USA, May 9-14, 2010, on CD-ROM, American Nuclear Society, LaGrange Park, IL.
Gelbard, E.M., 1968. Spherical Harmonics Methods: PL and Double-PL Approximations. In: Greenspan, H., Kelber, C.N., Okrent, D. (Eds.), Computing Methods in Reactor Physics. Gordon and Breach Science Publishers, New York. Hébert Alain, Applied Reactor Physics, Presses inter Polytechnique, pp.101, 2009. Ivanov, B.D., 2007. Methodology for Embedded Transport Core Calculation Ph.D thesis. The Pennsylvania State University, Pennsylvania. Joo, H.G., Cho, J.Y., Kim, K.S., Lee, C.C., Zee, S.Q., 2004, ‘‘Methods and Performance of a Three-Dimensional Whole-Core Transport Code DeCART,” PHYSOR 2004, Chicago, Illinois, April 25-29, 2004, on CD-ROM, American Nuclear Society, Lagrange Park, IL. Jung, Y.S., Shim, C.B., Lim, C.H., Joo, H.G., 2013. Practical numerical reactor employing direct whole core neutron transport and subchannel thermal/ hydraulic solvers. Ann. Nucl. Energy 62, 357–374. Kim, H.R., 1993. Global/Local Iterative Homogenization Methods for Neutron Diffusion Nodal Theory Ph.D thesis. Korea Advanced Institute of Science and Technology, Korea. Knight, M., Bryce, P., Hall, S., 2012, ‘‘WIMS/PANTHER Analysis of UO2/MOX Cores Using Embedded Supercells,” PHYSOR 2012, Knoxville, Tennessee, USA, April 15-20, 2012, on CD-ROM, American Nuclear Society, LaGrange Park, IL. Kooreman, G., Rahnema, F., 2015. Enhanced hybrid diffusion-transport spatial homogenization method. Nucl. Technol. 192, 264–277. Kozlowski, T., Xu, Y., Downar, T.J., Lee, D., 2011. Cell Homogenization Method for Pin-by-Pin Neutron Transport Calculations. Nucl. Sci. Eng. 169, 1–18. Lee, C.H., Downar, T.J., 2004. A hybrid nodal diffusion/SP3 method using one-node coarse-mesh finite difference formulation. Nucl. Sci. Eng. 146, 176–187. Lee, D., Downar, T.J., Kim, Y., 2004. A nodal and finite difference hybrid method for Pin-by-Pin heterogeneous three-dimensional light water reactor diffusion calculations. Nucl. Sci. Eng. 146, 319–339. Lin Tzung-Yi, Liu Yen-Wan Hsueh, 2017, ‘‘Hybrid Nodal Green’s Function Method with SP3 for Pin-by-pin Calculation,” M&C 2017 - International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering, Jeju, Korea, April 16-20, 2017. ‘‘MCNP — A General Monte Carlo N-particle transport code, Version 5,” LA-UR-031987, 2008. Ougouag, A.M., 1981. Coarse-Mesh Nodal Method for Multigroup Multidimensional Neutron Diffusion Computations Master thesis. University of Illinois, Urbana. Roberts, D.R., 2010. Development of Iterative Transport–diffusion Methodology for LWR Analysis Master thesis. The Pennsylvania State University, Pennsylvania. ‘‘Scale: A Comprehensive Modeling and Simulation Suite for Nuclear Safety Analysis and Design,” ORNL/TM-2005/39, Version 6.1, 2011. Smith, K.S., 1986. Assembly homogenization techniques for light water reactor analysis. Prog. Nucl. Energy 17 (3), 303–335. Smith, K.S., 1994. Practical and efficient iterative method for LWR fuel assembly homogenization. Trans. Am. Nucl. Soc. 71, 238–241. Tada, K., Yamamoto, A., Yamane, Y., 2009. Treatment of staggered mesh for BWR pin by pin core analysis. J. Nucl. Sci. Technol. 46 (2), 163–174. Tatsumi, M., Yamamoto, A., Nagano, H., Sengoku, K., 2003. Advanced PWR core calculation based on multi-group nodal-transport method in three-dimensional pin-by-pin geometry. J. Nucl. Sci. Technol. 40 (6), 376–387. Yamamoto, A., Endo, T., Chao, Y.A., 2011, A derivation of discontinuity factor for angular flux in integro-differential transport equation, Transactions of the American Nuclear Society, vol. 104, Hollywood, Florida, June 26–30, 2011.
60
T.-Y. Lin, Y.-W.H Liu / Annals of Nuclear Energy 118 (2018) 49–60
Yamamoto, A., Tatsumi, M., 2006. Improvement of spatial discretization error on the semi-analytic nodal method using the scattered source subtraction method. J. Nucl. Sci. Technol. 43 (12), 1481–1489. Yamamoto, A., Sakamoto, T., Endo, T., 2016. Discontinuity factors for Simplified P3 theory. Nucl. Sci. Eng. 183, 39–51. Yamamoto, A., Sakamoto, T., Endo, T., 2016. A CMFD acceleration method for SP3 advanced nodal method. Nucl. Sci. Eng. 184, 168–173. Yasseri, S., Rahnema, F., 2014. Consistent spatial homogenization in transport theory. Nucl. Sci. Eng. 176, 292–311.
Yu, L., Chao, Y.A., 2015. A Unified Generic Theory on Discontinuity Factors in Diffusion, SP3 and Transport Calculations. Ann. Nucl. Energy 75, 239–248. Yu, L., Lu, D., Chao, Y.A., 2014. The Calculation Method for SP3 Discontinuity Factor and Its Application. Ann. Nucl. Energy 69, 14–24. Zhang, B., Mayhue, L., Huria, H., Ivanov, B., 2013. Development of a threedimensional pseudo pin-by-pin calculation methodology in ANC. Nucl. Technol. 183, 527–534.