Reactivity Effect Evaluation of Fast Reactor Based on Angular-Dependent Few-Group Cross Sections Generation

: Among all the possible occurring reactivity effects of a fast reactor, the situations whereby the control rod was inserted, or the coolant was voided could lead to strong anisotropy of neutron ﬂux distribution, therefore the angular dependence on neutron ﬂux should be considered during the few-group cross-sections generation. Therefore, the purpose of this paper is to compare the inﬂuence whether the angular dependence on neutron ﬂux is considered in the calculation of few-group cross sections for the reactivity effect calculation. In the study, the 1-D SN ﬁnite difference neutron transport equation solver was implemented in the TULIP of SARAX code system so that the high-order neutron ﬂux could be obtained. Meanwhile, the improved Tone’s method was also applied. The numerical results were obtained based on three experimental FR cores, the JOYO MK-I core, ZPPR-9 core, and ZPPR-10B core. Both control rod worth and sodium void reactivity were calculated and compared with the measurement data. By summarizing and comparing the results of 46 cases, signiﬁcant differences were found between different consideration of the neutronic analysis. The consideration of angular dependence on neutron ﬂux distribution in the few-group cross-sections generation was beneﬁcial to the neutronic design analysis of FR, especially for the reactivity effect calculation.


Introduction
Fast reactors (FR), including the sodium-cooled fast reactor, the gas-cooled fast reactor, and the lead-cooled fast reactor, are the main candidates in developing the next generation of advanced nuclear reactor. The research and development of FR are based on the neutronic numerical analysis. Since the two-step method is still the mainstream for the design work, noticeable progress of code development has been made around the world. In the two-step method, the few-group cross sections of each material were obtained through assembly-wise transport calculation. Then the power distribution of 3-D core was calculated by solving the 3-D neutron transport/diffusion equation with obtained few-group cross sections.
Considering that the neutron behavior of FR is quite different from the light water reactor, many codes for FR calculation specifically have been developed. Owing to the long mean free path of fast neutron, the equivalent 1-D assembly model is applied to take the heterogeneity effect into account during the calculation of few-group cross sections. As a result, the energy group could be divided into thousands of groups so that the resonance could be well treated. In order to balance the efficiency, the collision probability method (CPM) is widely used, such as ECCO [1], MC 2 -3 [2], and SLAROM-UF [3].
In the two-step calculation of FR, the neutron flux is used as weight function in order to get the few-group cross sections. According to the homogenization theory, the anisotropy of neutron flux could be well treated by using the angular neutron flux weighted total cross sections and scattering matrices [4,5]. However, in the derivation of CPM, the isotropic assumption of the scattering and fission source term was applied so that only the neutron Energies 2021, 14, 4042 2 of 12 scalar flux could be obtained. When facing problems with obvious anisotropic neutron flux distribution, such as the coolant void case or control rod insertion case, such an assumption of CPM could introduce significant differences. In order to solve this problem, the neutron transport equation with anisotropic scattering should be solved during the few-group cross sections calculation. The utilization of another neutron transport method could be an effective way to solve this problem, such as the SN method or the MOC method.
The SARAX code system [6], which was developed by Xi'an Jiaotong University, considers the requirements in modeling different kinds of reactors designed in China, including the critical and subcritical systems with sodium or lead as the coolant. The SARAX code consists of four parts: The cross-sections generation module TULIP [7], the reactor steady-state analysis module LAVENDER [8], the reactor transient analysis module DAISY [9], and the neutron spectrum modification and shielding design module HYDRA [10]. Many verification and validation works have been conducted during the code development stage, such as the OECD/NEA sodium-cooled fast reactors benchmark, the JOYO MK-I zero power experiments, and several assembly layouts of zero power physics reactor (ZPPR) series experiments [11][12][13].
In the original TULIP code, the CPM was also used to obtain the neutron flux. Therefore, the motivation of current study was to improve the accuracy of reactivity effect calculation, especially for the sodium void reactivity and control rod worth. In order to be compatible with current SARAX code system, the 1-D SN finite difference solver of Hydra code was implemented in the TULIP code. At the same time, the background cross sections of resonance treatment were calculated by using 1-D SN solver based on the improved Tone's method [14], which is similar to APOLLO3 [15] or EXUS-F [16] code.
The remainder of this paper is organized as follows. Section 2 briefly describes the new development about latest TULIP code. The numerical results based on JOYO MK-I [17], ZPPR-9 [18], and ZPPR-10B [19] reactors are summarized in Section 3. The conclusion of current study is drawn in Section 4.

New Development in TULIP Code
Originally, the CPM was used in the TULIP code for the resonance and transport calculation. In the resonance calculation, Tone's method [20] was applied in order to obtain region, nuclide, and energy group-wise background cross sections, as shown in Equation (1) where σ 0,i,r,g is the background cross section of nuclide i, region r, and group g. P j→r,g is the collision probability from region j to region r of group g. N is the nuclide density and V is volume, respectively. The subscript t stands for the total cross sections. Once the background cross sections are obtained, the group-wise escape cross sections Σ e,g,i,r and self-shielded cross sections σ x,r,i,g are calculated by using following equations under narrow resonance approximation: After finishing the resonance calculation, the neutron flux is calculated by using 1-D CPM solver. As a consequence, the few-group cross sections are calculated with following equations: In the above equations, the subscript n indicates the order number of Legendre polynomials, and the summation of volume is ignored. In the right-hand-side of Equations (4) and (6), the nth order expansions of the total and scattering cross sections are calculated based on the following equation.
In order to take the anisotropic neutron flux distribution into account, the 1-D SN finite difference solver was implemented in the TULIP code. Therefore, according to the improved Tone's method, the background cross sections can be alternatively obtained with two fixed-source equations, which could be solved by using SN solver: Ω · ∇φ 1,i,g (r, Ω) + Σ t,g (r)φ 1,i,g (r, Ω) = ∑ k =i Σ t,g,k (r) Ω · ∇φ 2,i,g (r, Ω) + Σ t,g (r)φ 2,i,g (r, Ω) = N i (r) (8) Furthermore, the background cross sections are determined by using the solution of above equations: By solving neutron transport equation with 1-D SN solver, the angular neutron flux could be obtained, and the Equations (4) and (6) change into the following form: Based on this, the 1-D SN solver and improved Tone's method were implemented into original TULIP code, and the computational flow chart of updated TULIP code is shown in Figure 1. In this figure, the OP1 and OP2 are the CPM and SN methods discussed above. The biggest difference between OP1 and OP2 is the consideration of angular dependence on neutron flux. The SARAX code system has a three-step standard procedure for steadystate analysis, which is TULIP-HYDRA-LAVENDER. It should be noted that during the HYDRA calculation, the global-wise space-and angular-dependent neutron flux is obtained. Therefore, in the previous validation studies, the angular dependence is actually considered. In order to discuss the effect of angular dependence on neutron flux, the TULIP-LAVENDER sequence is applied to the following section's calculation and the few-group cross sections are obtained directly from the TULIP code. Considering that the CPM method and SN method has the different convergence properties of spatial discretization, in the following calculation, a different mesh was used in the TULIP code to make sure the neutron flux solved by OP1 or OP2 was converged.  (11) Based on this, the 1-D SN solver and improved Tone's method were implemented into original TULIP code, and the computational flow chart of updated TULIP code is shown in Figure 1. In this figure, the OP1 and OP2 are the CPM and SN methods discussed above. The biggest difference between OP1 and OP2 is the consideration of angular dependence on neutron flux. The SARAX code system has a three-step standard procedure for steady-state analysis, which is TULIP-HYDRA-LAVENDER. It should be noted that during the HYDRA calculation, the global-wise space-and angular-dependent neutron flux is obtained. Therefore, in the previous validation studies, the angular dependence is actually considered. In order to discuss the effect of angular dependence on neutron flux, the TULIP-LAVENDER sequence is applied to the following section's calculation and the few-group cross sections are obtained directly from the TULIP code. Considering that the CPM method and SN method has the different convergence properties of spatial discretization, in the following calculation, a different mesh was used in the TULIP code to make sure the neutron flux solved by OP1 or OP2 was converged.

Numerical Results
In this section, several experimental reactors with measurement data were calculated and the difference between different options was compared.

JOYO MK-I Experimental Tests
JOYO is Japan's first experimental fast reactor that was constructed at the O-arai Engineering Center of the Power Reactor and Nuclear Fuel Development Corporation to acquire various aspects of data relating to actual fast reactor performances. The core was loaded with MOX fuel and surrounded by radial and axial blanket of depleted UO 2 . Several experiments were performed based on JOYO MK-I core, including critical measurements and reactivity effects measurements. Since the control rod and sodium void introduce high angular dependence on neutron flux, in the following calculation, the control rod worth and the sodium void reactivity of JOYO MK-I reactor were calculated and analyzed.

Control Rod Worth Calculation of JOYO MK-I Core
In the JOYO MK-I core, the six control rod assemblies were divided into two groups, regulating control rod (RR) and safety control rod (SR). These two groups and the radial position of each control rod assembly are shown in Figure 2. The control rod worth was obtained from two SARAX criticality calculations. The first one was obtained by full withdrawal of all six control rods, and the second one was obtained by full insertion of a specific control rod. The Equation (12) gives the definition of final control rod worth (CRW) results. The effective delayed neutron fraction was 521 pcm from benchmark report.
quire various aspects of data relating to actual fast reactor performances. The core was loaded with MOX fuel and surrounded by radial and axial blanket of depleted UO2. Several experiments were performed based on JOYO MK-I core, including critical measurements and reactivity effects measurements. Since the control rod and sodium void introduce high angular dependence on neutron flux, in the following calculation, the control rod worth and the sodium void reactivity of JOYO MK-I reactor were calculated and analyzed.

Control Rod Worth Calculation of JOYO MK-I Core
In the JOYO MK-I core, the six control rod assemblies were divided into two groups, regulating control rod (RR) and safety control rod (SR). These two groups and the radial position of each control rod assembly are shown in Figure 2. The control rod worth was obtained from two SARAX criticality calculations. The first one was obtained by full withdrawal of all six control rods, and the second one was obtained by full insertion of a specific control rod. The Equation (12) gives the definition of final control rod worth (CRW) results. The effective delayed neutron fraction was 521 pcm from benchmark report.
In the equation, subscript 1 indicates the full withdrawal case and 2 indicates the insertion of one specific rod. The effective delayed neutron fraction was 521 pcm according to the benchmark report.
The calculated CRW values of each control rod assembly are summarized in Table 1. Compared with measurement, the SR1 has the biggest relative difference in −5.91% of OP1 and −3.00% of OP2, respectively. It can be seen from the results that the relative difference In the equation, subscript 1 indicates the full withdrawal case and 2 indicates the insertion of one specific rod. The effective delayed neutron fraction was 521 pcm according to the benchmark report.
The calculated CRW values of each control rod assembly are summarized in Table 1. Compared with measurement, the SR1 has the biggest relative difference in −5.91% of OP1 and −3.00% of OP2, respectively. It can be seen from the results that the relative difference calculated by OP2 is generally less than that calculated by OP1, and all the results calculated by OP2 are within the uncertainty of the experimental results of each control rod worth measurement.   Table 2.  Figure 3 shows the position of voided assemblies during the measurement. In total, six assemblies were voided in the experiment.  (13) was used to obtain the sodium void reactivity (SVR) value. Subscript 1 and 2 indicate the voided case and normal case. During the measurement, two RRs were inserted into the core and all the SRs were out of the core. The position of each control rod assembly is shown in  By using different options of few-group cross sections generation, the calculated SVR value and the comparison are listed in Table 3. For such a small reactor core, the neutron By using different options of few-group cross sections generation, the calculated SVR value and the comparison are listed in Table 3. For such a small reactor core, the neutron leakage effect is the main contributor to the SVR in fuel region so that the SVR values are negative. In the blanket region, the positive SVR value is caused by the harden of neutron spectrum. The calculated SVR value of both OP1 and OP2 are within the uncertainty of the experimental results.

ZPPR Series Experimental Tests
Experiments at the ZPPR aimed to provide an initial data base for assessment of physics issues in large liquid metal fast breeder reactor cores. These experiments belong to a joint program, called the Japanese-United States Program of Integral Tests and Experimental Research (JUPITER). The JUPITER program consisted of twenty-one experimental assemblies. In this section, one clean benchmark core (ZPPR-9) and one engineering Energies 2021, 14, 4042 7 of 12 mockup core (ZPPR-10B) based on JUPITER-I were selected. The specifications of ZPPR-9 and ZPPR-10B are illustrated in Figures 4 and 5. According to the configuration of ZPPR-9 and ZPPR-10B, more types of drawers were loaded in the ZPPR-10B reactor core. 9 and ZPPR-10B are illustrated in Figures 4 and 5. According to the configuration of ZPPR-9 and ZPPR-10B, more types of drawers were loaded in the ZPPR-10B reactor core.
Unlike the JOYO MK-I core, the ZPPR cores were assembled from several types of small plates into stainless steel drawers, such as depleted uranium, Pu-U-Mo fuel, sodium, and stainless steel. In the calculation, the as-built drawers were simplified into onedimensional plate-stretch models, as shown in Figure 6. In addition, the canned sodium plates in drawers in the affected zone were replaced by void cans during the sodium void reactivity experiment, and the affected zone contained at least 9 drawers for each case. Therefore, the special design and specific experiment of ZPPR made the neutron flux of core more anisotropic. Consequently, the control rod worth and sodium void reactivity experiment of two cores were calculated in the following sections.    Unlike the JOYO MK-I core, the ZPPR cores were assembled from several types of small plates into stainless steel drawers, such as depleted uranium, Pu-U-Mo fuel, sodium, and stainless steel. In the calculation, the as-built drawers were simplified into onedimensional plate-stretch models, as shown in Figure 6. In addition, the canned sodium plates in drawers in the affected zone were replaced by void cans during the sodium void reactivity experiment, and the affected zone contained at least 9 drawers for each case. Therefore, the special design and specific experiment of ZPPR made the neutron flux of core more anisotropic. Consequently, the control rod worth and sodium void reactivity experiment of two cores were calculated in the following sections.

Control Rod Worth Calculation of ZPPR Cores
There were 14 and 6 cases of control rod worth experiment for ZPPR-9 and ZPPR-10B, respectively. The worth of absorber and follower at different position were measured. Tables 4 and 5 show the calculated CRW values for two cores. Results of the CRW calculation were normalized to a unit dollar, of which the normalized factor is 343.6p cm in the benchmark report. The definition of each case could be found in the benchmark report [18,19].
By comparing with measurement of ZPPR-9 core, the results of OP1 have relative differences ranging from 3% to 11%, while the results of OP2 have relative differences

Control Rod Worth Calculation of ZPPR Cores
There were 14 and 6 cases of control rod worth experiment for ZPPR-9 and ZPPR-10B, respectively. The worth of absorber and follower at different position were measured. Tables 4 and 5 show the calculated CRW values for two cores. Results of the CRW calculation were normalized to a unit dollar, of which the normalized factor is 343.6p cm in the benchmark report. The definition of each case could be found in the benchmark report [18,19].
By comparing with measurement of ZPPR-9 core, the results of OP1 have relative differences ranging from 3% to 11%, while the results of OP2 have relative differences ranging from −1% to 7%. In the calculation of CRW, the relative difference of the method used in OP2 is generally smaller than that of OP1. In this case, only the follower part is inserted into the core, such as Cr-4, Cr-5, and Cr-13, and the differences between OP1 and OP2 are around 1%. On the contrary, the differences between two options become larger when absorber is inserted. Therefore, OP2 has a better performance on the treatment of the stronger absorber.
For the ZPPR-10B core, the results obtained by using OP2 have also improved around 7%. Compared with ZPPR-9 core, the fuel drawer distribution tends to be more nonuniform. Therefore, the environment effect could lead to the discrepancy, especially for control rod drawer in outside fuel region. In the current 1-D assembly calculation, it is hard to rigorously take the environment effect into account. However, it still could be concluded that the utilization of OP2 could obtain better results.

Sodium Void Reactivity Calculation of ZPPR Core
During the experiment, the voided region was expanded step-by-step, as illustrated in Figures 7 and 8. The results of each step are summarized in Tables 6 and 7.
For the case of small sodium void reactivity of ZPPR-9 core, i.e., Step-1, Step-2, and Step-3, the relative difference of SVR value of OP1 is 7%, and that of OP2 is 16%. It should be noted that the difference between OP1 and OP2 of Step-1 is only 0.3 cent, which is less than 1 pcm and is smaller than the k eff convergence criteria. Therefore, for the case of small sodium void reactivity, it could be concluded that OP1 and OP2 have comparable results. For the case of large value, i.e., Step-4, Step-5, and Step-6, the maximum relative difference of the method based on OP1 is up to 40%. At that time, the maximum relative difference of the method based on OP2 is only 8%. Since the last two steps of this experiment have a larger void region, the OP2 gives accurate results because of the treatment of anisotropic neutron flux distribution. In addition, for all the cases of this experiment, the absolute values of relative difference based on OP1 and OP2 are in the range of 1%~40% and 5%~16%, respectively. It can be seen from the results that the treatment of anisotropic neutron flux distribution in few-group cross-sections generation could receive stable results.
From the results of ZPPR-10B, the maximum relative difference of SVR of the method based on OP1 is 38% while the maximum relative difference of the method based on OP2 is 6%. The comparison of these two options shows that the use of OP2 greatly improves the calculation results of sodium void reactivity. From the results of ZPPR-10B, both ZPPR-10B and ZPPR-9 could reach the same conclusion, which is for the case of strong anisotropy, where the relative difference of the method based on angular-dependent neutron flux is smaller and the accuracy is higher than that based on scalar neutron flux.
From the results of ZPPR-10B, the maximum relative difference of SVR of the method based on OP1 is 38% while the maximum relative difference of the method based on OP2 is 6%. The comparison of these two options shows that the use of OP2 greatly improves the calculation results of sodium void reactivity. From the results of ZPPR-10B, both ZPPR-10B and ZPPR-9 could reach the same conclusion, which is for the case of strong anisotropy, where the relative difference of the method based on angular-dependent neutron flux is smaller and the accuracy is higher than that based on scalar neutron flux.    sults. From the results of ZPPR-10B, the maximum relative difference of SVR of the method based on OP1 is 38% while the maximum relative difference of the method based on OP2 is 6%. The comparison of these two options shows that the use of OP2 greatly improves the calculation results of sodium void reactivity. From the results of ZPPR-10B, both ZPPR-10B and ZPPR-9 could reach the same conclusion, which is for the case of strong anisotropy, where the relative difference of the method based on angular-dependent neutron flux is smaller and the accuracy is higher than that based on scalar neutron flux.

Summary and Conclusions
The evaluation of reactivity effect is important to FR neutronic design. Among all the possible occurring reactivity effects, the control rod insertion or coolant void could lead to strong anisotropy of neutron flux distribution. In order to simulate this anisotropy accurately, the angular-dependent neutron flux should be used in the process of the generation of few-group cross sections based on homogenization theory. Therefore, the purpose of this paper is to compare the influence regarding whether the angular dependence on neutron flux is considered in the calculation of few-group cross sections for the reactivity effect calculation. In the study, the 1-D SN finite difference neutron transport equation solver was implemented in the TULIP code so that the highorder neutron flux could be obtained during the assembly calculation. Meanwhile, the improved Tone's method was applied in terms of implementation of SN solver.
The numerical results were obtained based on three experimental FR cores, JOYO MK-I core, ZPPR-9 core, and ZPPR-10B core. Both control rod worth and sodium void reactivity were calculated and compared with the measurement. Figure 9 summarizes all the numerical results and illustrated as the statistical histogram in terms of relative difference. In total, 46 cases were calculated and compared in this study. When the fewgroup cross sections are generated based on OP1, only six cases have small deviation compared with measured data within the range from 0% to 3%. At the same time, the relative differences of 15 cases are larger than 7%. Considering the angular dependence by using the OP2, nearly half of these cases have small deviation. At the same time, the number of cases of largest difference (>15%) has decreased obviously. We should note that the OP1 and OP2 obtained comparable results in a few cases, which only happened when the control rod follower part was inserted into the active core. That is to say, OP1 works well only at the condition that neutron flux distribution is not that anisotropic.

Summary and Conclusions
The evaluation of reactivity effect is important to FR neutronic design. Among all the possible occurring reactivity effects, the control rod insertion or coolant void could lead to strong anisotropy of neutron flux distribution. In order to simulate this anisotropy accurately, the angular-dependent neutron flux should be used in the process of the generation of few-group cross sections based on homogenization theory. Therefore, the purpose of this paper is to compare the influence regarding whether the angular dependence on neutron flux is considered in the calculation of few-group cross sections for the reactivity effect calculation. In the study, the 1-D SN finite difference neutron transport equation solver was implemented in the TULIP code so that the highorder neutron flux could be obtained during the assembly calculation. Meanwhile, the improved Tone's method was applied in terms of implementation of SN solver.
The numerical results were obtained based on three experimental FR cores, JOYO MK-I core, ZPPR-9 core, and ZPPR-10B core. Both control rod worth and sodium void reactivity were calculated and compared with the measurement. Figure 9 summarizes all the numerical results and illustrated as the statistical histogram in terms of relative difference. In total, 46 cases were calculated and compared in this study. When the few-group cross sections are generated based on OP1, only six cases have small deviation compared with measured data within the range from 0% to 3%. At the same time, the relative differences of 15 cases are larger than 7%. Considering the angular dependence by using the OP2, nearly half of these cases have small deviation. At the same time, the number of cases of largest difference (>15%) has decreased obviously. We should note that the OP1 and OP2 obtained comparable results in a few cases, which only happened when the control rod follower part was inserted into the active core. That is to say, OP1 works well only at the condition that neutron flux distribution is not that anisotropic.  To sum up, the consideration of angular dependence on neutron flux distribution during the few-group cross sections generation is beneficial to the neutronic design analysis of FR, especially for the reactivity effect calculation.