- Email: [email protected]
Whole-core forward-adjoint neutron transport solutions with coupled 2-D MOC and 1-D SN and kinetics parameter calculation
Progress in Nuclear Energy 108 (2018) 310–318
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
Whole-core forward-adjoint neutron transport solutions with coupled 2-D MOC and 1-D SN and kinetics parameter calculation
T
Qu Wua,b,∗, Xingjie Pengb, Xiao Tanga,b, Yingrui Yub, Qing Lib, Kan Wanga a b
Department of Engineering Physics, Tsinghua University, Beijing, 100084, China Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, 610041, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Forward-adjoint neutron transport 2-D/1D coupling method Kinetics parameter C5G7
The 2-D/1D whole-core transport method has been widely studied for 3-D fine flux or power distributions and implemented in many deterministic codes. KYCORE, a 2-D/1-D transport code, developed a new iteration strategy to calculate forward flux. In order to perform kinetics parameter computation and sensitivity analysis, an adjoint flux solution is required in modern reactor physics analysis. In this study, a new adjoint neutron transport solver, named KYADJ, was developed by utilizing the 2-D MOC and 1-D SN coupling method. A 3 × 3 lattice test problem and C5G7 OECD/NEA 3-D benchmarks were used to verify the forward-adjoint neutron transport calculation of KYADJ. Forward-adjoint multiplication factors, forward flux, adjoint flux and kinetics parameters were compared with reference results generated by the Reactor Monte Carlo code RMC. Results showed that KYADJ agreed well with RMC and KYADJ had the ability to provide fine pin-by-pin forward and adjoint flux distributions and accurately compute kinetics parameters.
1. Introduction An adjoint flux distribution physically represents neutron importance that measures the contribution of neutrons to the steady power in a critical reactor core (Cacuci, 2010). It has two aspects of applications: the adjoint sensitivity analysis based on the perturbation theory Cacuci, 2003) and the quasi-static method in point kinetics equations (Stacey, 2007). The point kinetics equation solutions require the accurate calculation of point reactor kinetics parameters which use the adjoint flux as their weight function. Therefore, the adjoint neutron transport equation has been extensively studied in many advanced codes (Han et al., 2015; Peng et al., 2017; Pusa, 2012; Shokueifar et al., 2016; Zhu et al., 2015a,b). The “two-step” approach dominates the conventional industrial application and reactor analysis due to its mature technology and convenience. Firstly, a lattice code provides homogenized parameters such as few-group cross-sections and discontinuous factors by conducting the transport calculation of an isolated assembly. Secondly, nodal codes apply these parameters to the whole-core diffusion calculation. The accuracy of this approach remains a serious problem as a result of approximations in models and numerical methods. Compared with the conventional “two-step” approach, the whole-core transport calculation can provide more precise physical quantities at the expense of time and memory (Yuk and Cho, 2015). At the present stage, three-
∗
dimensional (3-D) method of characteristics (MOC) encounters great obstacles due to the huge computational overhead. Recent studies have focused on the two-dimensional (2-D)/one-dimensional (1-D) coupled transport calculation (Joo et al., 2004; Tang et al., 2017; Yuk and Cho, 2015; Zhu et al., 2015a,b). Generally, 2-D MOC is used in the radial direction where most inhomogeneity occurs, and different varieties of 1-D difference methods such as diffusion and SN are adopted in the axial direction (Zhu et al., 2015a,b). The 2-D adjoint transport solutions have been investigated by many researchers in lattice codes (Han et al., 2015; Pusa, 2012; Shokueifar et al., 2016). The MOC and 3-D coarse mesh finite difference (CMFD) adjoint methods have been achieved in MPACT (Zhu et al., 2015a,b). However, few researchers attempt to solve the 2-D/1-D coupled adjoint transport equations. It would be of interest to figure out the fine distribution of 3-D pin-by-pin adjoint flux. We aimed, therefore, to apply the radial MOC and axial SN to solve the forward-adjoint Boltzmann equation, and then verify the validity via kinetics parameter solutions. KYCORE (Tang et al., 2017), developed by Nuclear Power Institute of China, is a modern 2-D/1-D whole-core transport code designed for the direct transport calculation without any assembly homogenization steps. KYCORE uses the 2-D MOC and 1-D SN coupled by the radial and axial leakage and accelerates the iteration convergence with 3-D CMFD. An optimized procedure has been explored in KYCORE to decrease the computational time and stabilize the convergence of the iteration (Tang
Corresponding author. Department of Engineering Physics, Tsinghua University, Beijing, 100084, China. E-mail address: [email protected] (Q. Wu).
https://doi.org/10.1016/j.pnucene.2018.06.006 Received 23 September 2017; Received in revised form 21 May 2018; Accepted 10 June 2018 0149-1970/ © 2018 Elsevier Ltd. All rights reserved.
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
et al., 2017). On the basis of KYCORE, we have developed a new code, called KYADJ, to perform the forward-adjoint neutron transport calculation. The paper mainly clarifies the theory and method of the forward-adjoint neutron transport in detail, and then shows numerical results for a 3 × 3 lattice test problem and the C5G7 OECD/NEA 3-D benchmarks. Eventually, the paper compares the point reactor kinetics parameters with those calculated from the Reactor Monte Carlo (RMC) code (Wang et al., 2015) for verification. 2. The forward-adjoint neutron transport solutions In KYADJ, the steady 3-D neutron forward transport equations are solved by the angular direction and energy discretization assuming isotropic scattering. In the cylindrical coordinate system, the equations can be written as
ξm
∂ψg (r, Ωm) ∂r
Σt , g (r) ψg (r, Ωm) =
1 4π
+ ηm
(
∂ψg (r, Ωm) r ∂θ
−
G ∑g′= 1 Σs, g′→ g ϕ g′ (r)
∂ψg (r, Ωm)
+
)+μ
m
r ∂ω
1 χ 4πk eff g
∂ψg (r, Ωm) ∂z
+
G ∑g′= 1 υΣf , g′ ϕ g′ (r)
(1) where Ωm refers to the discrete angular direction, ω denotes the angle between the projection of Ωm on the plane (er , eθ ) and er , ξm = Ωm⋅er , ηm = Ωm⋅eθ , μm = Ωm⋅ez , r is the spatial position, g is the energy group index, ψg is the angular flux, keff is the effective multiplication factor, Σt , g Σs , χg , υ and Σf , g′ are the total cross section, the scattering cross section, the fission neutron energy spectrum, the average number of fission neutrons and the fission cross section, respectively. ϕ g′ represents the scalar flux approximated by quadrature formula:
ϕg (r) =
∑ ωm ψg (r, Ωm),
Fig. 1. The code modules and efficient iteration strategy for forward-adjoint calculation.
(2)
m
LgAxial , m, l (r ) ≡
where ωm is the quadrature weight. The corresponding steady 3-D neutron adjoint transport equations can be shown as
− ξm
∂ψg∗ (r, Ωm)
Σt , g (r) ψg∗ (r, Ωm) =
∂r 1 4π
∂ψg∗ (r, Ωm)
− ηm ⎛ r ∂θ − ⎝ G ∑g′= 1 Σs, g → g′ ϕ g∗′ (r) +
∂ψg∗ (r, Ωm)
Once Qg, m, l and are provided, Eq. (4) can be solved by the 2-D MOC which is discussed in detail in Chai's paper (Chai et al., 2016). The corresponding 2-D MOC adjoint equations can be derived:
∂ψ∗ (r, Ω )
m ⎞ − μm g + ∂z ⎠ G υΣf , g ∑g′= 1 χ g′ ϕ g∗′ (r)
− ξm (3)
Qg∗, m, l (r ) ≡
+ Σt , g ψg, m, l (r ) = Qg, m, l (r ) − LgAxial , m, l (r ),
1 Qg, m, l (r ) ≡ 4π
1 ∑ Σs,g′→g ϕ g′,l (r) + 4πk χg eff g ′= 1
∂ψ (r, Ωm)
∑ υΣf ,g′ ϕ g′,l (r),
g ′= 1
1 υΣf , g ∗ 4πkeff
G
∑ χ g′ ϕ g∗′ (r), g ′= 1
μm ∗ ⎡ψ (r ) − ψ∗ (r ) ⎤. g , m, l − 1 2 Δz l ⎣ g, m, l + 1 2 ⎦
(8)
(9)
In a similar way, integrate both sides of Eq. (1) over the radial mesh p and divide the area Δs . A serial of radial 1-D SN forward equations can be obtained (Tang et al., 2017):
μm
G g ′= 1
G
∑ Σs,g →g′ ϕ g∗′ (r) +
2.2. 1-D SN adjoint equations
(4)
direction is considered. The differential terms g and g are r ∂θ r ∂ω ignored in the 2-D MOC. The total source term and the axial leakage term are expressed as G
(7)
In comparison with the 2-D forward equations, the adjoint equations exchange υΣf , g with χg and transpose the P0 scattering matrix. The first term in Eq. (7) and the leakage term have the inverse direction with those in Eq. (4). In other words, if the total source term Qg, m, l is replaced merely by the adjoint source term Qg∗, m, l in Eq. (4), the 2-D MOC solver can offer the adjoint angular flux ψg∗, m′, l , where the subscript m′ signifies the inverse direction with m .
where Qg, m, l and LgAxial , m, l are defined as the total source term and the axial leakage term, respectively. Note that only the change of flux along the r ∂ψ (r, Ωm)
+ Σt , g ψg∗, m, l (r ) = Qg∗, m, l (r ) − Lg∗,Axial m, l (r ), Lg∗,Axial m, l
1 4π
Lg∗,Axial m, l (r ) ≡ −
Suppose all cross sections remain unchanged within the axial region [z l − 1/2, z l + 1/2]. Integrating both sides of Eq. (1) over the axial coordinate z and dividing the axial mesh span Δz l , a serial of radial 2-D forward equations can be obtained:
∂r
∂r
where and are defined as the total adjoint source term and the axial adjoint leakage term is expressed as
2.1. 2-D MOC adjoint equations
∂ψg, m, l (r )
∂ψg∗, m, l (r ) Qg∗, m, l
where ψg∗ and ϕg∗ are the adjoint angular flux and adjoint scalar flux, and ∗ is the effective multiplication factor. keff The procedure of the forward neutron transport equation solution has been clarified exhaustively in the previous paper (Tang et al., 2017). Hence we focus on the solving of the adjoint neutron transport equation, especially its differences from the forward neutron transport equation in this paper.
ξm
(6)
LgAxial , m, l
r ∂ω 1 ∗ 4πkeff
μm ⎡ψ 1 (r ) − ψg , m, l − 1 (r ) ⎤ . 2 ⎦ Δz l ⎣ g, m, l + 2
∂ψg, m, p (z ) ∂z
+ Σt , g ψg, m, p (z ) = Qg, m, p (z ) − LgRadial , m, p (z ),
(10)
where Qg, m, p is the total source term which has the same format as Eq. (5). LgRadial , m, p is the radial leakage provided by the 2-D MOC sweep that
(5) 311
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 2. 3 × 3 lattice benchmark model.
Table 1 Multiplication factors for the test problem.
Σx =
Method
k eff
Std/Diff (pcm)
RMC Forward Adjoint
1.09742 1.09804 1.09747
5 62 5
1 − μm2 Δdk
∑
Δs
k
⎡ψg, m, p + 1 (z ) − ψg, m, p − 1 (z ) ⎤, 2 2 ⎣ ⎦
− μm
∂z
+
Σt , g ψg∗, m, p (z )
=
Qg∗, m, p (z )
−
Radial Lg∗, m , p (z ),
,
ϕi =
∑fi ∈ i ϕfi Vfi Vfi
∼i, j i, j Jgi, j = −D (ϕg, j − ϕg, i ) + Dˆ (ϕg, j + ϕg, i ),
(11)
2Di Dj
∼i, j D =
where Δdk represents the ray segment width in the 2-D MOC. In KYADJ, Eq. (10) is solved by 1-D SN sweep. The 1-D SN adjoint equations have a similar express as
∂ψg∗, m, p (z )
∑fi ∈ i ϕfi Vfi
(15) ∼i, j i, j ˆ The net current is estimated by Eq. (16), where D and D are the ∼i, j coupling coefficients between the coarse mesh i and j . D is derived by Fick's law (Stacey, 2007) and expressed as Eq. (17), where D is the diffusion coefficient (Zhu et al., 2016) and h is the width of the coarse mesh. The aim is to conform the net current in the low order diffusion calculation to the net current in the high order transport calculation. i, j Dˆ is computed by Eq. (18), where Jghigh, i, j means the net current provided by the 2-D MOC and 1-D SN sweep.
satisfies
LgRadial , m, p (z ) ≡
∑fi ∈ i Σx , fi ϕfi Vfi
Di hj + Dj hi
,
∼i, j Jghigh, i, j + D (ϕg, j − ϕg, i )
i, j Dˆ =
ϕg, i + ϕg, j
(12)
(16)
(17)
. (18)
For simplicity, Eq. (14) is rewritten into the matrix format: Radial Lg∗, m , p (z )
≡ −∑
1 − μm2 Δdk
k
Δs
∗ ⎡ψ∗ ⎤ 1 (z ) − ψg , m, p − 1 (z ) . 2 ⎣ g , m, p + 2 ⎦
Aϕ =
(13)
As mentioned before, the adjoint angular flux can be obtained if the P0 scattering matrix is transposed and υΣf , g is exchanged with χg in the identical 1-D SN sweep, where the subscript m′ signifies the inverse direction with m .
ψg∗, m′, p
1 Bϕ. keff
(19)
The adjoint format of Eq. (14) is written as
A∗ϕ∗ =
1 ∗ ∗ Bϕ. ∗ keff
(20)
A∗
Research (bell) shows that is the transposed matrix of A in the diffusion theory and B is converted into B∗ when υΣf , g is exchanged with χg . Eqs. (19) and (20) can both be solved by the generalized minimal residual method (GMRES) in the Intel MKL database as long as the corresponding matrices are transmitted to the GMRES functions input parameter.
2.3. CMFD adjoint equations In KYADJ, CMFD is used to accelerate the source iteration. The CMFD equation is considered as the neutron balance equation in the space coarse mesh i and is shown as Eq. (14)
∑ Jgi,j S i,j + Vi Σt ϕg,i − j
Vi 4π
G
χg
∑ Σs,g′→g ϕ g′,i = Vi 4πk
g ′= 1
G
3. The kinetics parameter calculations
∑ υΣf ,g′ ϕ g′,i .
eff g ′= 1
The spatially dependent point kinetics equations (Stacey, 2007) are solved by the quasi-static method used in the reactor transient calculation. The main idea is that the flux is written as a product of an amplitude function which changes acutely with time and a shape function which changes slowly with time. The point kinetics equations can be derived by weighting the time-dependence neutron equation with the adjoint flux and integral over the whole space. It is significant to figure out the kinetics parameters before solving these equations.
(14) where Jgi, j is the net current in the surface between the coarse mesh i and j , S i, j is the area of the surface and Vi is the volume of the space coarse mesh i . The cross sections Σx where the subscript x means different types of cross section and the scalar flux ϕg, i are homogenized in the coarse mesh i . As is shown in Eq. (15), fi is the fine mesh index contained in the coarse mesh i . 312
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 3. Geometries of C5G7 OECD/NEA benchmarks (a) Benchmark fuel pin compositions, (b) core configuration and fuel pin layout, (c) Geometry for the unrodded and the rodded A configurations. 1
Table 2 Multiplication factors for 3-D C5G7 benchmarks.
Λ=
Method
k eff (UR)
Std/Diff (pcm)
Time (hour)
k eff (RA)
Std/Diff (pcm)
Time (hour)
RMC Forward Adjoint
1.14322 1.14326 1.14301
3 4 21
126.8 6.3 9.5
1.12830 1.12822 1.12794
4 8 36
72.3 5.5 6.5
∫V ∫E ∫Ω ψ∗ (r , Ω, E ) v ψ (r , Ω, E ) dΩdEdV ∫V ∫E ∫Ω ψ∗ (r , Ω, E )
χ (r , E ) ∫E′ ∫Ω′ νΣf 4π
,
(r , E′) ψ (r , Ω′, E′) dΩ′dE′dΩdEdV
(22)
D
where βi is the i′th delayed neutron fraction, χi is the i′th delayed neutron group spectra and v is the speed of neutron. In KYADJ, the angular flux and the adjoint angular flux are provided in the form of multi-group and angular discretization. Eqs. (21) and (22), hence, can be approximated by
Therefore, this study focuses on the result of the effective delayed neutron fraction βeff and the neutron generation time Λ , which are denoted as
βeff , i =
∑fi ∑g ∑m wm ψfi∗, g, m
βi χiD 4π χg
∑g′ νΣf ϕfi, g′
∑fi ∑g ∑m wm ψfi∗, g, m 4π ∑g′ νΣf , g′ ϕfi, g′
, (23)
1
βeff , i =
βi χiD (r , E ) ∫E′ ∫Ω′ νΣf (r , E′) ψ (r , Ω′, E′) dΩ′dE′dΩdEdV 4π χ (r , E ) E) νΣf (r , E′) ψ (r , Ω′, E′) dΩ′dE′dΩdEdV ∫ E ′ ∫Ω′ 4π
∫V ∫E ∫Ω ψ∗ (r , Ω, E ) ∫V ∫E ∫Ω ψ∗ (r , Ω,
Λ=
,
(21)
∑fi ∑g ∑m wm ψfi∗, g, m v ψfi, g, m χg
∑fi ∑g ∑m wm ψfi∗, g, m 4π ∑g′ νΣf , g′ ϕfi, g′
, (24)
where the subscripts and physical quantities have the same meanings as those mentioned in Section 2. 313
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 4. The convergence curve of the RA case using the forward-adjoint calculation. (a) The convergence curve of the multiplication factor. (b) The convergence curve of the flux error.
314
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 5. The seven group forward flux distributions of 3-D C5G7 benchmarks.
4. Implementations
4.2. Iteration strategy
4.1. Boundary condition and incident flux
A new iteration strategy has been developed for the forward flux calculation in KYCORE. The iteration strategy is directly applied to the adjoint flux calculation in KYADJ. Fig. 1 shows the code modules and efficient iteration strategy for the forward-adjoint calculation. Both the forward and adjoint calculations begin with 3-D CMFD to offer the initial effective multiply factor and the scalar flux. Next, 1-D SN sweep provides axial leakage term for 2-D MOC source term. It is worth noting that 1-D SN sweep is executed again between 2-D MOC and 3-D CMFD calculation, which accelerates the convergence of the iteration. In the former code version, the iteration between 1-D SN and 2-D MOC spends considerable time to converge. Such an improvement greatly accelerates the convergence of the leakage source. Although unit inner iteration time increases, the total time of the convergence decreases because outer iteration speeds up (Tang et al., 2017). The convergence acceleration has been testified by the Fourier analysis method.
In KYADJ, the vacuum boundary condition or the reflective boundary condition can be used. The vacuum boundary conditions of Eqs. (1) and (3) can be expressed as
ψg (Γ , Ωm) = 0, Ωm⋅n < 0,
(25)
ψg∗ (Γ , Ωm) = 0, Ωm⋅n > 0,
(26)
where n represents the outward normal direction of the boundary face Γ . The reflective boundary condition of Eqs. (1) and (3) can be expressed as
ψg (Γ , Ωm) = ψg (Γ , Ω′m), Ωm⋅n < 0,
(27)
ψg∗ (Γ , Ωm) = ψg∗ (Γ , Ω′m), Ωm⋅n > 0,
(28)
5. Results and discussion 5.1. Test problem: 3 × 3 lattice benchmark
where Ω′m is the reflective direction of Ωm . According to Eq. (25)–(28), the forward solver can provide the adjoint flux with the opposite direction. The incident flux of every ray is required in advance in the 2-D MOC sweep. With regard to the vacuum boundary condition, the cyclic trajectory layout is adopted and the incident flux can be solved without iteration. The ray prolongation method is applied to the complex geometry and solves the incident flux combined with the reflective boundary condition. The ray layout methods and the solving of incident flux are exhaustively discussed in Chai's paper (Chai et al., 2016).
A 3 × 3 lattice problem for test and verification is described by Fig. 2. The radius of the cylinder is 0.54 cm which is filled with UO2 and surrounded by water moderator. The side length of the lattice is 1.26 cm, and the axial height of the model is 20 cm. The top boundary is a vacuum boundary condition and the rest is a reflection boundary condition. The right of Fig. 2 shows the division of the radial meshes. In axial direction, one layer is used for the 2-D MOC sweep as a result of axial homogeneity. Ten layers are divided for CMFD calculation in each MOC layer and ten layers are divided for SN sweep in each CMFD layer. 315
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 6. The seven group adjoint flux distributions of 3-D C5G7 benchmarks.
reference value comes from the RMC calculation using the same geometry and multi-group data file (OECD, 2005). The multiplication factors for 3-D C5G7 benchmarks (UR and RA cases) calculated by KYADJ are shown in Table 2. To some extent, the results confirm that the multiplication factors for the adjoint transport equation and forward transport equation are the same. KYADJ has high accuracy to compute the integral parameter keff , which agrees with RMC results. Besides, the calculation time is compared in Table 2. Note that RMC uses the parallel computing and KYADJ uses the serial computing. The calculation time is represented by the actual computing time multiplied by threads for comparison. The calculation time of KYADJ is far less than that of RMC. The KYADJ calculation can converge fast due to the new iteration strategy mentioned in Section 4.2, particularly in the multiplication factor calculation. Fig. 4a shows that the multiplication factors fully converge in the approximately tenth outer iteration using both forward and adjoint methods. In order to ensure the convergence of the forward and adjoint calculations, the convergence curve of the flux error is presented in Fig. 4b. The ordinate uses the logarithmic coordinates and the convergence criterion is set as 0.0001. The maximum error of flux denotes the maximum value of the relative errors between the coarse fluxes in current iteration and those in last iteration. The forward calculation converges after the 34th iteration and the adjoint calculation converges after the 40th iteration. Further on, the fine distribution of the local quantity can insure an accurate core code calculation. KYADJ can be unambiguously evaluated by fine forward and adjoint flux distributions. In this paper, the fuel assembly regions in the third MOC sweep lay from the bottom of C5G7 benchmarks are selected to describe the forward and adjoint flux distributions. The seven group (Group 1 represents the highest energy group) forward flux distributions are presented in Fig. 5. The figure
The model is computed by KYADJ and RMC for comparison. Note that all RMC calculations are based on the iterated fission probability (IFP) method to obtain the adjoint flux (Qiu et al., 2017). Table 1 shows the results of the multiplication factor (keff ) calculated by the forward-adjoint transport and Monte Carlo method. The statistical standard deviations (Std) for RMC results and the differences (Diff) of the forward-adjoint results compared with RMC results are given in form of pcm, i.e., 10−5. It can be seen that the multiplication factor of the forward-adjoint calculation are both quite precise. It should be pointed out that the multiplication factor by the adjoint calculation identifies with that by RMC calculation. 5.2. 3-D C5G7 benchmarks The 3-D C5G7 problems (OECD, 2005) are generated by OECD/NEA to identify the capacity of deterministic codes without homogenization. The unrodded (UR) and the rodded A (RA) cases are selected to test the ability of KYADJ in the forward-adjoint calculation. As Fig. 3 shows, the benchmarks consist of four assemblies which are made up of 17 × 17 UO2 and MOX pin cells and are surrounded by reflectors. Axially, a total height of the problem is 64.26 cm and the height for the reflector at the top of the fuel region is 21.42 cm. Control rods are inserted into the top reflector region for the unrodded case and the 1/3 height of one UO2 assembly for the rodded A case. The calculation details are as follows. The axial region is divided into four layers for the MOC sweep, which are subdivided into seven, seven, seven and ten layers for CMFD calculation respectively. Each CMFD layer is divided into ten SN layers. The radial region uses 82156 rays for the MOC. Six polar angles are used for the axial direction and twenty azimuthal angles for the radial direction. The 3-D C5G7 problems use seven group macroscopic cross sections for calculation. The 316
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 7. The relative error distributions of the seven group forward flux between KYADJ and RMC. Table 3 Kinetics parameters for and the test problem and 3-D C5G7 benchmarks. Cases
Lattice benchmark C5G7-UR C5G7-RA
Λ (us)
βeff
RMC
KYADJ
Re (%)
RMC
KYADJ
Re (%)
1.6151E+01 1.5107E+01 1.4906E+01
1.5945E+01 1.5022E+01 1.4893E+01
−1.3 −0.6 −0.1
7.0342E-03 5.8044E-03 5.6954E-03
7.0912E-03 5.8680E-03 5.7937E-03
0.8 1.1 1.7
benchmark and C5G7 benchmarks. The computational parameters, such as χiD , βi and v in Eqs. (23) and (24), are from another report of C5G7 benchmarks (Boyarinov et al., 2016). The reference value comes from the RMC calculation (Qiu et al., 2017) using the same delayed neutron data. Kinetics parameters are presented in Table 3. The relative errors in the table show that the results computed by KYADJ agree with RMC. Therefore, the consistency of the kinetics parameters verifies the reliability of the adjoint transport calculation.
presents that the distribution becomes smoother with the energy group decreasing. The maximum value of flux appears in the center of the reactor cote for fast groups, but not for thermal groups. The seven group adjoint flux distributions are presented in Fig. 6. The adjoint flux shows an entirely different distribution from the forward flux. For both the fast groups and the thermal groups, the adjoint flux in the central regions is larger than that in the peripheral regions. Considering the physical meaning of the adjoint flux, neutrons in the central regions have greater contribution to fission power, i.e. neutron value, compared with neutrons in the peripheral regions. RMC can provide the pin-by-pin flux distribution. Fig. 7 shows the relative error distributions with RMC results as reference. The relative errors are below 2% in majority regions. In minority peripheral regions, the relative errors do not exceed 5%. It is worth noting that the statistical errors of local quantities are large in peripheral regions for the Monte Carlo method. Therefore, the forward flux distributions provided by KYADJ have high accuracy and fidelity.
6. Conclusions In this work, the 2-D MOC and 1-D SN coupling forward-adjoint transport calculation were implemented in the deterministic code KYADJ on the basis of KYCORE. The new iteration strategy was applied to the adjoint transport solution. Our results show that iteration can converge fast and steadily and the code can provide the pin-by-pin forward and adjoint flux distributions. The multiplication factors and the forward flux distributions agree with RMC results. We fail to compare the adjoint flux distributions with RMC because RMC lacks the adjoint flux tally at present. However, both KYADJ and RMC can calculate kinetics parameters using the adjoint flux as the weight function.
5.3. Verification of kinetics parameters The kinetics parameters βeff and Λ are calculated for the 3 × 3 lattice 317
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
The adjoint flux distributions can be verified by comparing the kinetics parameters. Our study guarantees the precision of the adjoint flux. Comparing with previous work, we calculated the 3-D fine adjoint flux rather than 2-D adjoint flux using advanced 2-D MOC and 1-D SN method. Moreover, KYADJ outperforms the Monte Carlo method in terms of memory and time consumption. At present, aforementioned calculations are in serial mode. In future work, parallelization will be realized in KYADJ. The forward and the adjoint flux will be applied to compute sensitivity coefficients and uncertainties.
Han, T.Y., Lee, H.C., Noh, J.M., 2015. Development of a sensitivity and uncertainty analysis code for high temperature gas-cooled reactor physics based on the generalized perturbation theory. Ann. Nucl. Energy 85, 501–511. Joo, H.G., Cho, J.Y., Kim, K.S., 2004. Methods and Performance of a Three-dimensional Whole-core Transport Code DeCART. PHYSOR 2004, Chicago, Illinois, April 25–29. OECD, 2005. Benchmark on Deterministic Transport Calculations without Spatial Homogenisation, OECD Report ISBN 92-64-01069-6. MOX Fuel Assembly 3DExtension Case, NEA No. 5420. Peng, X., Liang, J., Forget, B., 2017. A New Generalized Sensitivity Calculation Approach Based on Continuous-energy Monte Carlo Method. RPHA17, Chengdu, China, August 24–25. Pusa, M., 2012. Incorporating sensitivity and uncertainty analysis to a lattice physics code with application to CASMO-4. Ann. Nucl. Energy 40, 153–162. Qiu, Y., Wang, Z., Li, K., 2017. Calculation of adjoint-weighted kinetic parameters with the Reactor Monte Carlo code RMC. Prog. Nucl. Energy xxx, 1–11. Shokueifar, V., Zangian, M., Minuchehr, A., et al., 2016. Synchronized forward-adjoint neutron transport using method of characteristics. Prog. Nucl. Energy 92, 211–219. Stacey, W.M., 2007. Nuclear Reactor Physics. WILEY- VCH Verlag GmbH & Co. KGaA. Tang, X., Li, Q., Chai, X., 2017. Efficient Procedure for Radial MOC and Axial SN Coupled 3D Neutron Transport Calculation. M&C2017, Jeju, Korea, April 16–20. Wang, K., Li, Z., She, D., et al., 2015. RMC - a Monte Carlo code for reactor core analysis. Ann. Nucl. Energy 82, 121–129. Yuk, S., Cho, N.Z., 2015. Whole-Core transport solutions with 2-D/1-D fusion kernel via pCMFD acceleration and p-CMFD embedding of nonoverlapping local/global iteration. Nucl. Sci. Eng. 181, 1–16. Zhu, A., Jarrett, M., Xu, Y., 2016. An optimally diffusive coarse mesh finite difference method to accelerate neutron transport calculations. Ann. Nucl. Energy 95, 116–124. Zhu, A., Xu, Y., Graham, A., et al., 2015a. Transient Methods for Pin-resolved Whole Core Transport Using the 2D-1D Methodology in MPACT. ANS MC2015, Nashville, TN, April 19–23. Zhu, A., Xu, Y., Saller, T., et al., 2015b. The Implementation and Analysis of the MOC and CMFD Adjoint Capabilities in the 2D-1D Code MPACT. ANS MC2015, Nashville, TN, April 19–23.
Acknowledgements The work is supported by the Innovation Center of Nuclear Power Technology for National Defense Industry NO. HDLCXZX-2017-HD018. The authors are grateful to Mr. Xiaotong Shang and Mr. Guanlin Shi of Tsinghua University for his advice on RMC and providing the data for comparison in this study. References Boyarinov, V.F., Fomichenko, P.A., Hou, J., 2016. Deterministic Time-dependent Neutron Transport Benchmark without Spatial Homogenization (C5G7-td). NEA/NSC/ DOC(2016). Version 1.7. Cacuci, D.G., 2003. Sensitivity and Uncertainty Analysis, vol. 1 Chapman & Hall/CRC, Boca Raton, Fla, USA. Cacuci, D.G., 2010. Handbook of Nuclear Engineering. Springer Science Publications. Chai, X., Tu, X., Lu, W., et al., 2016. Development and V&V of MOC module in advanced neutronics lattice code KYLIN-2. Nucl. Power Eng. 37 (4), 154–159.
318
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
Whole-core forward-adjoint neutron transport solutions with coupled 2-D MOC and 1-D SN and kinetics parameter calculation
T
Qu Wua,b,∗, Xingjie Pengb, Xiao Tanga,b, Yingrui Yub, Qing Lib, Kan Wanga a b
Department of Engineering Physics, Tsinghua University, Beijing, 100084, China Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, 610041, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Forward-adjoint neutron transport 2-D/1D coupling method Kinetics parameter C5G7
The 2-D/1D whole-core transport method has been widely studied for 3-D fine flux or power distributions and implemented in many deterministic codes. KYCORE, a 2-D/1-D transport code, developed a new iteration strategy to calculate forward flux. In order to perform kinetics parameter computation and sensitivity analysis, an adjoint flux solution is required in modern reactor physics analysis. In this study, a new adjoint neutron transport solver, named KYADJ, was developed by utilizing the 2-D MOC and 1-D SN coupling method. A 3 × 3 lattice test problem and C5G7 OECD/NEA 3-D benchmarks were used to verify the forward-adjoint neutron transport calculation of KYADJ. Forward-adjoint multiplication factors, forward flux, adjoint flux and kinetics parameters were compared with reference results generated by the Reactor Monte Carlo code RMC. Results showed that KYADJ agreed well with RMC and KYADJ had the ability to provide fine pin-by-pin forward and adjoint flux distributions and accurately compute kinetics parameters.
1. Introduction An adjoint flux distribution physically represents neutron importance that measures the contribution of neutrons to the steady power in a critical reactor core (Cacuci, 2010). It has two aspects of applications: the adjoint sensitivity analysis based on the perturbation theory Cacuci, 2003) and the quasi-static method in point kinetics equations (Stacey, 2007). The point kinetics equation solutions require the accurate calculation of point reactor kinetics parameters which use the adjoint flux as their weight function. Therefore, the adjoint neutron transport equation has been extensively studied in many advanced codes (Han et al., 2015; Peng et al., 2017; Pusa, 2012; Shokueifar et al., 2016; Zhu et al., 2015a,b). The “two-step” approach dominates the conventional industrial application and reactor analysis due to its mature technology and convenience. Firstly, a lattice code provides homogenized parameters such as few-group cross-sections and discontinuous factors by conducting the transport calculation of an isolated assembly. Secondly, nodal codes apply these parameters to the whole-core diffusion calculation. The accuracy of this approach remains a serious problem as a result of approximations in models and numerical methods. Compared with the conventional “two-step” approach, the whole-core transport calculation can provide more precise physical quantities at the expense of time and memory (Yuk and Cho, 2015). At the present stage, three-
∗
dimensional (3-D) method of characteristics (MOC) encounters great obstacles due to the huge computational overhead. Recent studies have focused on the two-dimensional (2-D)/one-dimensional (1-D) coupled transport calculation (Joo et al., 2004; Tang et al., 2017; Yuk and Cho, 2015; Zhu et al., 2015a,b). Generally, 2-D MOC is used in the radial direction where most inhomogeneity occurs, and different varieties of 1-D difference methods such as diffusion and SN are adopted in the axial direction (Zhu et al., 2015a,b). The 2-D adjoint transport solutions have been investigated by many researchers in lattice codes (Han et al., 2015; Pusa, 2012; Shokueifar et al., 2016). The MOC and 3-D coarse mesh finite difference (CMFD) adjoint methods have been achieved in MPACT (Zhu et al., 2015a,b). However, few researchers attempt to solve the 2-D/1-D coupled adjoint transport equations. It would be of interest to figure out the fine distribution of 3-D pin-by-pin adjoint flux. We aimed, therefore, to apply the radial MOC and axial SN to solve the forward-adjoint Boltzmann equation, and then verify the validity via kinetics parameter solutions. KYCORE (Tang et al., 2017), developed by Nuclear Power Institute of China, is a modern 2-D/1-D whole-core transport code designed for the direct transport calculation without any assembly homogenization steps. KYCORE uses the 2-D MOC and 1-D SN coupled by the radial and axial leakage and accelerates the iteration convergence with 3-D CMFD. An optimized procedure has been explored in KYCORE to decrease the computational time and stabilize the convergence of the iteration (Tang
Corresponding author. Department of Engineering Physics, Tsinghua University, Beijing, 100084, China. E-mail address: [email protected] (Q. Wu).
https://doi.org/10.1016/j.pnucene.2018.06.006 Received 23 September 2017; Received in revised form 21 May 2018; Accepted 10 June 2018 0149-1970/ © 2018 Elsevier Ltd. All rights reserved.
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
et al., 2017). On the basis of KYCORE, we have developed a new code, called KYADJ, to perform the forward-adjoint neutron transport calculation. The paper mainly clarifies the theory and method of the forward-adjoint neutron transport in detail, and then shows numerical results for a 3 × 3 lattice test problem and the C5G7 OECD/NEA 3-D benchmarks. Eventually, the paper compares the point reactor kinetics parameters with those calculated from the Reactor Monte Carlo (RMC) code (Wang et al., 2015) for verification. 2. The forward-adjoint neutron transport solutions In KYADJ, the steady 3-D neutron forward transport equations are solved by the angular direction and energy discretization assuming isotropic scattering. In the cylindrical coordinate system, the equations can be written as
ξm
∂ψg (r, Ωm) ∂r
Σt , g (r) ψg (r, Ωm) =
1 4π
+ ηm
(
∂ψg (r, Ωm) r ∂θ
−
G ∑g′= 1 Σs, g′→ g ϕ g′ (r)
∂ψg (r, Ωm)
+
)+μ
m
r ∂ω
1 χ 4πk eff g
∂ψg (r, Ωm) ∂z
+
G ∑g′= 1 υΣf , g′ ϕ g′ (r)
(1) where Ωm refers to the discrete angular direction, ω denotes the angle between the projection of Ωm on the plane (er , eθ ) and er , ξm = Ωm⋅er , ηm = Ωm⋅eθ , μm = Ωm⋅ez , r is the spatial position, g is the energy group index, ψg is the angular flux, keff is the effective multiplication factor, Σt , g Σs , χg , υ and Σf , g′ are the total cross section, the scattering cross section, the fission neutron energy spectrum, the average number of fission neutrons and the fission cross section, respectively. ϕ g′ represents the scalar flux approximated by quadrature formula:
ϕg (r) =
∑ ωm ψg (r, Ωm),
Fig. 1. The code modules and efficient iteration strategy for forward-adjoint calculation.
(2)
m
LgAxial , m, l (r ) ≡
where ωm is the quadrature weight. The corresponding steady 3-D neutron adjoint transport equations can be shown as
− ξm
∂ψg∗ (r, Ωm)
Σt , g (r) ψg∗ (r, Ωm) =
∂r 1 4π
∂ψg∗ (r, Ωm)
− ηm ⎛ r ∂θ − ⎝ G ∑g′= 1 Σs, g → g′ ϕ g∗′ (r) +
∂ψg∗ (r, Ωm)
Once Qg, m, l and are provided, Eq. (4) can be solved by the 2-D MOC which is discussed in detail in Chai's paper (Chai et al., 2016). The corresponding 2-D MOC adjoint equations can be derived:
∂ψ∗ (r, Ω )
m ⎞ − μm g + ∂z ⎠ G υΣf , g ∑g′= 1 χ g′ ϕ g∗′ (r)
− ξm (3)
Qg∗, m, l (r ) ≡
+ Σt , g ψg, m, l (r ) = Qg, m, l (r ) − LgAxial , m, l (r ),
1 Qg, m, l (r ) ≡ 4π
1 ∑ Σs,g′→g ϕ g′,l (r) + 4πk χg eff g ′= 1
∂ψ (r, Ωm)
∑ υΣf ,g′ ϕ g′,l (r),
g ′= 1
1 υΣf , g ∗ 4πkeff
G
∑ χ g′ ϕ g∗′ (r), g ′= 1
μm ∗ ⎡ψ (r ) − ψ∗ (r ) ⎤. g , m, l − 1 2 Δz l ⎣ g, m, l + 1 2 ⎦
(8)
(9)
In a similar way, integrate both sides of Eq. (1) over the radial mesh p and divide the area Δs . A serial of radial 1-D SN forward equations can be obtained (Tang et al., 2017):
μm
G g ′= 1
G
∑ Σs,g →g′ ϕ g∗′ (r) +
2.2. 1-D SN adjoint equations
(4)
direction is considered. The differential terms g and g are r ∂θ r ∂ω ignored in the 2-D MOC. The total source term and the axial leakage term are expressed as G
(7)
In comparison with the 2-D forward equations, the adjoint equations exchange υΣf , g with χg and transpose the P0 scattering matrix. The first term in Eq. (7) and the leakage term have the inverse direction with those in Eq. (4). In other words, if the total source term Qg, m, l is replaced merely by the adjoint source term Qg∗, m, l in Eq. (4), the 2-D MOC solver can offer the adjoint angular flux ψg∗, m′, l , where the subscript m′ signifies the inverse direction with m .
where Qg, m, l and LgAxial , m, l are defined as the total source term and the axial leakage term, respectively. Note that only the change of flux along the r ∂ψ (r, Ωm)
+ Σt , g ψg∗, m, l (r ) = Qg∗, m, l (r ) − Lg∗,Axial m, l (r ), Lg∗,Axial m, l
1 4π
Lg∗,Axial m, l (r ) ≡ −
Suppose all cross sections remain unchanged within the axial region [z l − 1/2, z l + 1/2]. Integrating both sides of Eq. (1) over the axial coordinate z and dividing the axial mesh span Δz l , a serial of radial 2-D forward equations can be obtained:
∂r
∂r
where and are defined as the total adjoint source term and the axial adjoint leakage term is expressed as
2.1. 2-D MOC adjoint equations
∂ψg, m, l (r )
∂ψg∗, m, l (r ) Qg∗, m, l
where ψg∗ and ϕg∗ are the adjoint angular flux and adjoint scalar flux, and ∗ is the effective multiplication factor. keff The procedure of the forward neutron transport equation solution has been clarified exhaustively in the previous paper (Tang et al., 2017). Hence we focus on the solving of the adjoint neutron transport equation, especially its differences from the forward neutron transport equation in this paper.
ξm
(6)
LgAxial , m, l
r ∂ω 1 ∗ 4πkeff
μm ⎡ψ 1 (r ) − ψg , m, l − 1 (r ) ⎤ . 2 ⎦ Δz l ⎣ g, m, l + 2
∂ψg, m, p (z ) ∂z
+ Σt , g ψg, m, p (z ) = Qg, m, p (z ) − LgRadial , m, p (z ),
(10)
where Qg, m, p is the total source term which has the same format as Eq. (5). LgRadial , m, p is the radial leakage provided by the 2-D MOC sweep that
(5) 311
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 2. 3 × 3 lattice benchmark model.
Table 1 Multiplication factors for the test problem.
Σx =
Method
k eff
Std/Diff (pcm)
RMC Forward Adjoint
1.09742 1.09804 1.09747
5 62 5
1 − μm2 Δdk
∑
Δs
k
⎡ψg, m, p + 1 (z ) − ψg, m, p − 1 (z ) ⎤, 2 2 ⎣ ⎦
− μm
∂z
+
Σt , g ψg∗, m, p (z )
=
Qg∗, m, p (z )
−
Radial Lg∗, m , p (z ),
,
ϕi =
∑fi ∈ i ϕfi Vfi Vfi
∼i, j i, j Jgi, j = −D (ϕg, j − ϕg, i ) + Dˆ (ϕg, j + ϕg, i ),
(11)
2Di Dj
∼i, j D =
where Δdk represents the ray segment width in the 2-D MOC. In KYADJ, Eq. (10) is solved by 1-D SN sweep. The 1-D SN adjoint equations have a similar express as
∂ψg∗, m, p (z )
∑fi ∈ i ϕfi Vfi
(15) ∼i, j i, j ˆ The net current is estimated by Eq. (16), where D and D are the ∼i, j coupling coefficients between the coarse mesh i and j . D is derived by Fick's law (Stacey, 2007) and expressed as Eq. (17), where D is the diffusion coefficient (Zhu et al., 2016) and h is the width of the coarse mesh. The aim is to conform the net current in the low order diffusion calculation to the net current in the high order transport calculation. i, j Dˆ is computed by Eq. (18), where Jghigh, i, j means the net current provided by the 2-D MOC and 1-D SN sweep.
satisfies
LgRadial , m, p (z ) ≡
∑fi ∈ i Σx , fi ϕfi Vfi
Di hj + Dj hi
,
∼i, j Jghigh, i, j + D (ϕg, j − ϕg, i )
i, j Dˆ =
ϕg, i + ϕg, j
(12)
(16)
(17)
. (18)
For simplicity, Eq. (14) is rewritten into the matrix format: Radial Lg∗, m , p (z )
≡ −∑
1 − μm2 Δdk
k
Δs
∗ ⎡ψ∗ ⎤ 1 (z ) − ψg , m, p − 1 (z ) . 2 ⎣ g , m, p + 2 ⎦
Aϕ =
(13)
As mentioned before, the adjoint angular flux can be obtained if the P0 scattering matrix is transposed and υΣf , g is exchanged with χg in the identical 1-D SN sweep, where the subscript m′ signifies the inverse direction with m .
ψg∗, m′, p
1 Bϕ. keff
(19)
The adjoint format of Eq. (14) is written as
A∗ϕ∗ =
1 ∗ ∗ Bϕ. ∗ keff
(20)
A∗
Research (bell) shows that is the transposed matrix of A in the diffusion theory and B is converted into B∗ when υΣf , g is exchanged with χg . Eqs. (19) and (20) can both be solved by the generalized minimal residual method (GMRES) in the Intel MKL database as long as the corresponding matrices are transmitted to the GMRES functions input parameter.
2.3. CMFD adjoint equations In KYADJ, CMFD is used to accelerate the source iteration. The CMFD equation is considered as the neutron balance equation in the space coarse mesh i and is shown as Eq. (14)
∑ Jgi,j S i,j + Vi Σt ϕg,i − j
Vi 4π
G
χg
∑ Σs,g′→g ϕ g′,i = Vi 4πk
g ′= 1
G
3. The kinetics parameter calculations
∑ υΣf ,g′ ϕ g′,i .
eff g ′= 1
The spatially dependent point kinetics equations (Stacey, 2007) are solved by the quasi-static method used in the reactor transient calculation. The main idea is that the flux is written as a product of an amplitude function which changes acutely with time and a shape function which changes slowly with time. The point kinetics equations can be derived by weighting the time-dependence neutron equation with the adjoint flux and integral over the whole space. It is significant to figure out the kinetics parameters before solving these equations.
(14) where Jgi, j is the net current in the surface between the coarse mesh i and j , S i, j is the area of the surface and Vi is the volume of the space coarse mesh i . The cross sections Σx where the subscript x means different types of cross section and the scalar flux ϕg, i are homogenized in the coarse mesh i . As is shown in Eq. (15), fi is the fine mesh index contained in the coarse mesh i . 312
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 3. Geometries of C5G7 OECD/NEA benchmarks (a) Benchmark fuel pin compositions, (b) core configuration and fuel pin layout, (c) Geometry for the unrodded and the rodded A configurations. 1
Table 2 Multiplication factors for 3-D C5G7 benchmarks.
Λ=
Method
k eff (UR)
Std/Diff (pcm)
Time (hour)
k eff (RA)
Std/Diff (pcm)
Time (hour)
RMC Forward Adjoint
1.14322 1.14326 1.14301
3 4 21
126.8 6.3 9.5
1.12830 1.12822 1.12794
4 8 36
72.3 5.5 6.5
∫V ∫E ∫Ω ψ∗ (r , Ω, E ) v ψ (r , Ω, E ) dΩdEdV ∫V ∫E ∫Ω ψ∗ (r , Ω, E )
χ (r , E ) ∫E′ ∫Ω′ νΣf 4π
,
(r , E′) ψ (r , Ω′, E′) dΩ′dE′dΩdEdV
(22)
D
where βi is the i′th delayed neutron fraction, χi is the i′th delayed neutron group spectra and v is the speed of neutron. In KYADJ, the angular flux and the adjoint angular flux are provided in the form of multi-group and angular discretization. Eqs. (21) and (22), hence, can be approximated by
Therefore, this study focuses on the result of the effective delayed neutron fraction βeff and the neutron generation time Λ , which are denoted as
βeff , i =
∑fi ∑g ∑m wm ψfi∗, g, m
βi χiD 4π χg
∑g′ νΣf ϕfi, g′
∑fi ∑g ∑m wm ψfi∗, g, m 4π ∑g′ νΣf , g′ ϕfi, g′
, (23)
1
βeff , i =
βi χiD (r , E ) ∫E′ ∫Ω′ νΣf (r , E′) ψ (r , Ω′, E′) dΩ′dE′dΩdEdV 4π χ (r , E ) E) νΣf (r , E′) ψ (r , Ω′, E′) dΩ′dE′dΩdEdV ∫ E ′ ∫Ω′ 4π
∫V ∫E ∫Ω ψ∗ (r , Ω, E ) ∫V ∫E ∫Ω ψ∗ (r , Ω,
Λ=
,
(21)
∑fi ∑g ∑m wm ψfi∗, g, m v ψfi, g, m χg
∑fi ∑g ∑m wm ψfi∗, g, m 4π ∑g′ νΣf , g′ ϕfi, g′
, (24)
where the subscripts and physical quantities have the same meanings as those mentioned in Section 2. 313
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 4. The convergence curve of the RA case using the forward-adjoint calculation. (a) The convergence curve of the multiplication factor. (b) The convergence curve of the flux error.
314
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 5. The seven group forward flux distributions of 3-D C5G7 benchmarks.
4. Implementations
4.2. Iteration strategy
4.1. Boundary condition and incident flux
A new iteration strategy has been developed for the forward flux calculation in KYCORE. The iteration strategy is directly applied to the adjoint flux calculation in KYADJ. Fig. 1 shows the code modules and efficient iteration strategy for the forward-adjoint calculation. Both the forward and adjoint calculations begin with 3-D CMFD to offer the initial effective multiply factor and the scalar flux. Next, 1-D SN sweep provides axial leakage term for 2-D MOC source term. It is worth noting that 1-D SN sweep is executed again between 2-D MOC and 3-D CMFD calculation, which accelerates the convergence of the iteration. In the former code version, the iteration between 1-D SN and 2-D MOC spends considerable time to converge. Such an improvement greatly accelerates the convergence of the leakage source. Although unit inner iteration time increases, the total time of the convergence decreases because outer iteration speeds up (Tang et al., 2017). The convergence acceleration has been testified by the Fourier analysis method.
In KYADJ, the vacuum boundary condition or the reflective boundary condition can be used. The vacuum boundary conditions of Eqs. (1) and (3) can be expressed as
ψg (Γ , Ωm) = 0, Ωm⋅n < 0,
(25)
ψg∗ (Γ , Ωm) = 0, Ωm⋅n > 0,
(26)
where n represents the outward normal direction of the boundary face Γ . The reflective boundary condition of Eqs. (1) and (3) can be expressed as
ψg (Γ , Ωm) = ψg (Γ , Ω′m), Ωm⋅n < 0,
(27)
ψg∗ (Γ , Ωm) = ψg∗ (Γ , Ω′m), Ωm⋅n > 0,
(28)
5. Results and discussion 5.1. Test problem: 3 × 3 lattice benchmark
where Ω′m is the reflective direction of Ωm . According to Eq. (25)–(28), the forward solver can provide the adjoint flux with the opposite direction. The incident flux of every ray is required in advance in the 2-D MOC sweep. With regard to the vacuum boundary condition, the cyclic trajectory layout is adopted and the incident flux can be solved without iteration. The ray prolongation method is applied to the complex geometry and solves the incident flux combined with the reflective boundary condition. The ray layout methods and the solving of incident flux are exhaustively discussed in Chai's paper (Chai et al., 2016).
A 3 × 3 lattice problem for test and verification is described by Fig. 2. The radius of the cylinder is 0.54 cm which is filled with UO2 and surrounded by water moderator. The side length of the lattice is 1.26 cm, and the axial height of the model is 20 cm. The top boundary is a vacuum boundary condition and the rest is a reflection boundary condition. The right of Fig. 2 shows the division of the radial meshes. In axial direction, one layer is used for the 2-D MOC sweep as a result of axial homogeneity. Ten layers are divided for CMFD calculation in each MOC layer and ten layers are divided for SN sweep in each CMFD layer. 315
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 6. The seven group adjoint flux distributions of 3-D C5G7 benchmarks.
reference value comes from the RMC calculation using the same geometry and multi-group data file (OECD, 2005). The multiplication factors for 3-D C5G7 benchmarks (UR and RA cases) calculated by KYADJ are shown in Table 2. To some extent, the results confirm that the multiplication factors for the adjoint transport equation and forward transport equation are the same. KYADJ has high accuracy to compute the integral parameter keff , which agrees with RMC results. Besides, the calculation time is compared in Table 2. Note that RMC uses the parallel computing and KYADJ uses the serial computing. The calculation time is represented by the actual computing time multiplied by threads for comparison. The calculation time of KYADJ is far less than that of RMC. The KYADJ calculation can converge fast due to the new iteration strategy mentioned in Section 4.2, particularly in the multiplication factor calculation. Fig. 4a shows that the multiplication factors fully converge in the approximately tenth outer iteration using both forward and adjoint methods. In order to ensure the convergence of the forward and adjoint calculations, the convergence curve of the flux error is presented in Fig. 4b. The ordinate uses the logarithmic coordinates and the convergence criterion is set as 0.0001. The maximum error of flux denotes the maximum value of the relative errors between the coarse fluxes in current iteration and those in last iteration. The forward calculation converges after the 34th iteration and the adjoint calculation converges after the 40th iteration. Further on, the fine distribution of the local quantity can insure an accurate core code calculation. KYADJ can be unambiguously evaluated by fine forward and adjoint flux distributions. In this paper, the fuel assembly regions in the third MOC sweep lay from the bottom of C5G7 benchmarks are selected to describe the forward and adjoint flux distributions. The seven group (Group 1 represents the highest energy group) forward flux distributions are presented in Fig. 5. The figure
The model is computed by KYADJ and RMC for comparison. Note that all RMC calculations are based on the iterated fission probability (IFP) method to obtain the adjoint flux (Qiu et al., 2017). Table 1 shows the results of the multiplication factor (keff ) calculated by the forward-adjoint transport and Monte Carlo method. The statistical standard deviations (Std) for RMC results and the differences (Diff) of the forward-adjoint results compared with RMC results are given in form of pcm, i.e., 10−5. It can be seen that the multiplication factor of the forward-adjoint calculation are both quite precise. It should be pointed out that the multiplication factor by the adjoint calculation identifies with that by RMC calculation. 5.2. 3-D C5G7 benchmarks The 3-D C5G7 problems (OECD, 2005) are generated by OECD/NEA to identify the capacity of deterministic codes without homogenization. The unrodded (UR) and the rodded A (RA) cases are selected to test the ability of KYADJ in the forward-adjoint calculation. As Fig. 3 shows, the benchmarks consist of four assemblies which are made up of 17 × 17 UO2 and MOX pin cells and are surrounded by reflectors. Axially, a total height of the problem is 64.26 cm and the height for the reflector at the top of the fuel region is 21.42 cm. Control rods are inserted into the top reflector region for the unrodded case and the 1/3 height of one UO2 assembly for the rodded A case. The calculation details are as follows. The axial region is divided into four layers for the MOC sweep, which are subdivided into seven, seven, seven and ten layers for CMFD calculation respectively. Each CMFD layer is divided into ten SN layers. The radial region uses 82156 rays for the MOC. Six polar angles are used for the axial direction and twenty azimuthal angles for the radial direction. The 3-D C5G7 problems use seven group macroscopic cross sections for calculation. The 316
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
Fig. 7. The relative error distributions of the seven group forward flux between KYADJ and RMC. Table 3 Kinetics parameters for and the test problem and 3-D C5G7 benchmarks. Cases
Lattice benchmark C5G7-UR C5G7-RA
Λ (us)
βeff
RMC
KYADJ
Re (%)
RMC
KYADJ
Re (%)
1.6151E+01 1.5107E+01 1.4906E+01
1.5945E+01 1.5022E+01 1.4893E+01
−1.3 −0.6 −0.1
7.0342E-03 5.8044E-03 5.6954E-03
7.0912E-03 5.8680E-03 5.7937E-03
0.8 1.1 1.7
benchmark and C5G7 benchmarks. The computational parameters, such as χiD , βi and v in Eqs. (23) and (24), are from another report of C5G7 benchmarks (Boyarinov et al., 2016). The reference value comes from the RMC calculation (Qiu et al., 2017) using the same delayed neutron data. Kinetics parameters are presented in Table 3. The relative errors in the table show that the results computed by KYADJ agree with RMC. Therefore, the consistency of the kinetics parameters verifies the reliability of the adjoint transport calculation.
presents that the distribution becomes smoother with the energy group decreasing. The maximum value of flux appears in the center of the reactor cote for fast groups, but not for thermal groups. The seven group adjoint flux distributions are presented in Fig. 6. The adjoint flux shows an entirely different distribution from the forward flux. For both the fast groups and the thermal groups, the adjoint flux in the central regions is larger than that in the peripheral regions. Considering the physical meaning of the adjoint flux, neutrons in the central regions have greater contribution to fission power, i.e. neutron value, compared with neutrons in the peripheral regions. RMC can provide the pin-by-pin flux distribution. Fig. 7 shows the relative error distributions with RMC results as reference. The relative errors are below 2% in majority regions. In minority peripheral regions, the relative errors do not exceed 5%. It is worth noting that the statistical errors of local quantities are large in peripheral regions for the Monte Carlo method. Therefore, the forward flux distributions provided by KYADJ have high accuracy and fidelity.
6. Conclusions In this work, the 2-D MOC and 1-D SN coupling forward-adjoint transport calculation were implemented in the deterministic code KYADJ on the basis of KYCORE. The new iteration strategy was applied to the adjoint transport solution. Our results show that iteration can converge fast and steadily and the code can provide the pin-by-pin forward and adjoint flux distributions. The multiplication factors and the forward flux distributions agree with RMC results. We fail to compare the adjoint flux distributions with RMC because RMC lacks the adjoint flux tally at present. However, both KYADJ and RMC can calculate kinetics parameters using the adjoint flux as the weight function.
5.3. Verification of kinetics parameters The kinetics parameters βeff and Λ are calculated for the 3 × 3 lattice 317
Progress in Nuclear Energy 108 (2018) 310–318
Q. Wu et al.
The adjoint flux distributions can be verified by comparing the kinetics parameters. Our study guarantees the precision of the adjoint flux. Comparing with previous work, we calculated the 3-D fine adjoint flux rather than 2-D adjoint flux using advanced 2-D MOC and 1-D SN method. Moreover, KYADJ outperforms the Monte Carlo method in terms of memory and time consumption. At present, aforementioned calculations are in serial mode. In future work, parallelization will be realized in KYADJ. The forward and the adjoint flux will be applied to compute sensitivity coefficients and uncertainties.
Han, T.Y., Lee, H.C., Noh, J.M., 2015. Development of a sensitivity and uncertainty analysis code for high temperature gas-cooled reactor physics based on the generalized perturbation theory. Ann. Nucl. Energy 85, 501–511. Joo, H.G., Cho, J.Y., Kim, K.S., 2004. Methods and Performance of a Three-dimensional Whole-core Transport Code DeCART. PHYSOR 2004, Chicago, Illinois, April 25–29. OECD, 2005. Benchmark on Deterministic Transport Calculations without Spatial Homogenisation, OECD Report ISBN 92-64-01069-6. MOX Fuel Assembly 3DExtension Case, NEA No. 5420. Peng, X., Liang, J., Forget, B., 2017. A New Generalized Sensitivity Calculation Approach Based on Continuous-energy Monte Carlo Method. RPHA17, Chengdu, China, August 24–25. Pusa, M., 2012. Incorporating sensitivity and uncertainty analysis to a lattice physics code with application to CASMO-4. Ann. Nucl. Energy 40, 153–162. Qiu, Y., Wang, Z., Li, K., 2017. Calculation of adjoint-weighted kinetic parameters with the Reactor Monte Carlo code RMC. Prog. Nucl. Energy xxx, 1–11. Shokueifar, V., Zangian, M., Minuchehr, A., et al., 2016. Synchronized forward-adjoint neutron transport using method of characteristics. Prog. Nucl. Energy 92, 211–219. Stacey, W.M., 2007. Nuclear Reactor Physics. WILEY- VCH Verlag GmbH & Co. KGaA. Tang, X., Li, Q., Chai, X., 2017. Efficient Procedure for Radial MOC and Axial SN Coupled 3D Neutron Transport Calculation. M&C2017, Jeju, Korea, April 16–20. Wang, K., Li, Z., She, D., et al., 2015. RMC - a Monte Carlo code for reactor core analysis. Ann. Nucl. Energy 82, 121–129. Yuk, S., Cho, N.Z., 2015. Whole-Core transport solutions with 2-D/1-D fusion kernel via pCMFD acceleration and p-CMFD embedding of nonoverlapping local/global iteration. Nucl. Sci. Eng. 181, 1–16. Zhu, A., Jarrett, M., Xu, Y., 2016. An optimally diffusive coarse mesh finite difference method to accelerate neutron transport calculations. Ann. Nucl. Energy 95, 116–124. Zhu, A., Xu, Y., Graham, A., et al., 2015a. Transient Methods for Pin-resolved Whole Core Transport Using the 2D-1D Methodology in MPACT. ANS MC2015, Nashville, TN, April 19–23. Zhu, A., Xu, Y., Saller, T., et al., 2015b. The Implementation and Analysis of the MOC and CMFD Adjoint Capabilities in the 2D-1D Code MPACT. ANS MC2015, Nashville, TN, April 19–23.
Acknowledgements The work is supported by the Innovation Center of Nuclear Power Technology for National Defense Industry NO. HDLCXZX-2017-HD018. The authors are grateful to Mr. Xiaotong Shang and Mr. Guanlin Shi of Tsinghua University for his advice on RMC and providing the data for comparison in this study. References Boyarinov, V.F., Fomichenko, P.A., Hou, J., 2016. Deterministic Time-dependent Neutron Transport Benchmark without Spatial Homogenization (C5G7-td). NEA/NSC/ DOC(2016). Version 1.7. Cacuci, D.G., 2003. Sensitivity and Uncertainty Analysis, vol. 1 Chapman & Hall/CRC, Boca Raton, Fla, USA. Cacuci, D.G., 2010. Handbook of Nuclear Engineering. Springer Science Publications. Chai, X., Tu, X., Lu, W., et al., 2016. Development and V&V of MOC module in advanced neutronics lattice code KYLIN-2. Nucl. Power Eng. 37 (4), 154–159.
318