Resonance self-shielding methodology of new neutron transport code STREAM

Abstract

This paper reports on the development and verification of three new resonance self-shielding methods. The verifications were performed using the new neutron transport code, STREAM. The new methodologies encompass the extension of energy range for resonance treatment, the development of optimum rational approximation, and the application of resonance treatment to isotopes in the cladding region. (1) The extended resonance energy range treatment has been developed to treat the resonances below 4 eV of three resonance isotopes and shows significant improvements in the accuracy of effective cross sections (XSs) in that energy range. (2) The optimum rational approximation can eliminate the geometric limitations of the conventional approach of equivalence theory and can also improve the accuracy of fuel escape probability. (3) The cladding resonance treatment method makes it possible to treat resonances in cladding material which have not been treated explicitly in the conventional methods. These three new methods have been implemented in the new lattice physics code STREAM and the improvement in the accuracy of effective XSs is demonstrated through detailed verification calculations.

Keywords:

• resonance
• resonance integral
• cross section
• PWR type reactor
• equivalence theory
• self-shielding
• resonance treatment
• fuel escape probability
• Dancoff factor

1. Introduction

One of the most important phenomena in nuclear reactor analysis is resonance self-shielding and its treatment is still not an easy problem even with modern lattice physics codes. Equivalence theory [1] accounts for the self-shielding effect in heterogeneous systems and is widely used in production calculations to compute effective multi-group (MG) XSs for commercial reactor analysis. It provides a theoretical basis for treating the following two different systems as an equivalent system: the first system is a homogeneous problem with pre-generated XS library corresponding to selected background XSs; and the second system is a heterogeneous problem given at hand. The background XS which represents the magnitude of the self-shielding effect, plays the role of a bridge connecting the two systems.

Modern lattice physics codes such as CASMO, NEWT, HELIOS, and DeCART implemented resonance treatments based on equivalence theory [2–5]. CASMO uses Dancoff factors and rational approximation to represent the fuel escape probability for an infinite array of fuel pins. The Dancoff factors are adjusted considering surrounding pincells. CASMO is known as a fast and accurate lattice code designed for commercial pressurized water reactors (PWRs) and boiling water reactors (BWRs) and it does not support a general geometry like HELIOS. NEWT uses a point-wise (PW) flux spectrum to compute MG XSs, which can reduce errors caused by resonance interferences. However, the Wigner-Seitz approximation [6], adopted for the calculation of PW flux spectrum, can cause significant errors for heterogeneous problems. HELIOS and DeCART use the subgroup method to calculate the fuel escape probability, which enables both lattice codes to handle general geometries. However, there are difficulties in generating accurate subgroup parameters and the subgroup method takes a relatively long time to generate effective XSs because of multiple fixed source transport calculations which are needed to compute background XSs for each subgroup level in each resonance group. It is noteworthy that all contemporary lattice physics codes are based on the equivalence theory for resonance treatment. They have the following drawbacks.

First, the resonance treatment energy range has been set from 4 eV to several 10s of keV which is optimized for ²³⁸U [6,7]. However, some resonance isotopes such as Pu, Am, and Er have resonances even below 4 eV. The lattice codes treat these low lying resonances by adopting dozens of fine energy groups, which is not sufficient for high accuracy [6,8]. In order to treat the resonances directly by fine energy groups, an ultra-fine-group structure of typically more than 10,000 groups is required [9]. In this work, a new approach has been adopted to accurately treat the low lying resonances below 4 eV, i.e., the resonance integral table has been extended down to 0.3 eV and the background XSs are evaluated for that energy range to compute more accurate effective XSs. The new approach shows significant improvements in accuracy.

Second, the conventional equivalence theory uses rational approximations such as Wigner's one-term or Carlvik's two-term rational approximations for the fuel escape probability or fuel-to-fuel collision probability [1]. Since the analytic expressions are based on black pin assumption, there are uncertainties in applying the rational form to real world problems which are not black pins. Also, the conventional usage of rational approximation does not support general geometries since they are derived for cylindrical and slab geometry. In this paper, a new optimum rational approximation has been developed where a numerical method is used to directly evaluate fuel escape probability. Since the new method does not depend on analytic expression, general geometries can be supported. There are several differences in computing coefficients of rational approximation between the new method in this paper and the gray resonance treatment method [10]. The optimum rational approximation computes the fuel escape probability directly by counting the actual number of neutrons.

Third, most contemporary lattice physics codes consider resonance isotopes only in the fuel region, control rod region, and burnable absorbers. However, zirconium isotopes in the cladding material have resonances too which can cause 200–300 pcm errors for PWRs if the resonances are ignored [6]. Traditionally, the resonances in cladding material have been ignored or the XSs are generated for a typical background XS (i.e., 300 barn) [1,6]. The coating materials in integral fuel burnable absorbers [11] also contain resonance isotopes and this can cause error if the resonances are not treated correctly. In this work, a new resonance treatment method for the cladding material and coated burnable absorber has been developed. In the new method, the annular ring is approximated to infinite slab geometry, and then Roman's two-term approximation [1] is applied to find the fuel/cladding escape probabilities. In a comparison of reaction rates and XSs to those of Monte Carlo calculations by Monte Carlo N-Particle (MCNP) [12], it was verified that the new method can reduce difference caused by ignoring the resonance in cladding materials.

In addition to the three new methods, processes to generate a resonance upscattering correction factor are presented. Recently, it has been reported that a neutron upscattering event has significant impact on reactivity and the fuel temperature coefficient (FTC) [13], and several methods have been developed by Monte Carlo methods [13,14]. It is not a simple task to model the resonance upscattering event by a deterministic method. Therefore, the resonance upscattering correction factor is generated by the Monte Carlo method and is applied to MG XS library in order to consider an equivalent effect of the resonance upscattering event.

The newly developed methods are described in this paper, as well as the conventional methods. The new methods are implemented into the lattice physics code STREAM which is being developed at UNIST [15,16]. The impacts of each method are evaluated and the overall accuracy of STREAM with the new resonance treatment methods is compared to that of Monte Carlo calculations. It is noted that the order of difference of STREAM code calculations is around ∼100 pcm for PWR analysis.

2. Nuclear data library generation

2.1. Multi-group cross section library

The XS processing program NJOY [17] is used to generate MG XSs as a function of a neutron energy group, background XS, and temperature. In order to generate XSs in the resonance energy range, NJOY solves the slowing down equation for the homogeneous medium with target isotopes whose XS will be computed for variations of background XS. The XSs of resonance isotopes are generated for 10 background XSs. The XSs are generated for several tens of or hundreds of energy groups. The energy groups are divided into three ranges: fast energy group; resonance energy group; and thermal energy group. The resonance energy range of many lattice codes is from 4 to 10 keV. The lower boundary 4 eV is selected from the upper boundary of thermal scattering treatment.

2.2. Extended resonance integral

In generating the MG library of STREAM, NJOY was used to generate self-shielded XSs from ENDF/B-VII.0. The XS library was generated for PWR and contains 391 isotopes. In the NJOY/GROUPR slowing down calculation, hydrogen was used as a moderator material.

The energy group structure for STREAM is similar to the 70 energy group structures of CASMO [2] except for some modifications of the fast and resonance energy groups. The energy group structure of STREAM is indicated in Table 1. The upper boundary of the resonance energy group is extended from about 10 to 25 keV so that it contains the whole resolved resonance energy range of ²³⁸U in ENDF/B-VII.0. The resonance treatment is performed in the resonance energy groups for most of the resonance isotopes. Some resonance isotopes (i.e., Pu, Am, etc.) have resonances below 4 eV. There are two drawbacks in treating resonance below 4 eV. First, many isotopes including Pu, Am have wide resonance in that energy range. Generally, energy group boundaries are adjusted in order to put important resonance close to the center of groups. However, it is difficult to decide energy boundaries because of wide resonances. Second, the resonance treatment does not consider up-scattering effect. For most of isotopes, 4 eV is assumed that cutoff of for upscattering. For these reasons, the resonances below 4 eV have been modeled with several 10s of fine energy groups. However, the number of energy groups is still not enough to treat those complicated resonances accurately. For those isotopes, resonance energy ranges have been extended and resonance integral is generated for those energy ranges. Figure 1 shows the capture XSs and the extended resonance integral energy range of U, Pu and Am. Although five resonance isotopes are indicated in the figure, STREAM uses the extended library for wider variety of isotopes (i.e., Er and Gd). The kinds of isotopes can be adjusted depending on the energy group. For example, a lattice code using hundreds of energy groups may not need the extended resonance integral because of its sufficiently fine energy group.

Table 1. Energy group structure.

Description	Group index	Energy range [eV]
Fast groups	1–12	24780–10^7
Resonance groups	13–27	4–24780
Extended resonance groups	28–51	0.3–4

Figure 1. Resonances below 4 eV and the extended resonance energy range for each isotope.

The extended resonance integral (or XS) has been generated by assuming downscattering only. In spite of the assumption, the extended resonance integral can reduce an error in the effective XS. This will be discussed in the verification section. In the extended resonance energy range, intermediate resonance parameters and the fuel escape probability are assumed to be the same as those of group 27 (4–9.877 eV) because the resonance is usually wide in low energy.

In addition to the extension of the energy range, the number of background XSs in the XS library is extended. Ten background XSs are not enough to compute the accurate effective XSs. The number of background XSs in the library is extended up to 19 for the dominant resonance isotopes (i.e., U and Pu), which can reduce interpolation errors associated with background XS intervals. Other resonance isotopes (i.e., Zr) are generated for 10 background XSs.

2.3. IR parameter generation by Monte Carlo code

The IR parameters for the IR approximation [1] adopted in STREAM are generated by comparing the slowing down capability of the isotope of interest with that of hydrogen. The solution of slowing down equations can be obtained with either ultra-fine group or continuous energy XSs [18]. For STREAM, a continuous energy Monte Carlo code MCS [19] is used in the following way:

1. The absorption XSs of ²³⁸U, σ^U238_{a, g}, are calculated for a range of variations in dilution, i.e., for background XSs σ₀ of 30 barns, 50 barns, 70 barns, etc. and the results are made into a σ^U238_{a, g} versus σ₀ table.

The slowing down spectrum calculations for this table generation are done for an infinite homogeneous mixture of ²³⁸U/¹H. The background XS σ₀ of ²³⁸U is defined as (1) $${\sigma }_0=\frac{1}{N^R}\sum _{i\ne R}{N}^i{\sigma }_p^i,$$(1) where N^R is the number density of ²³⁸U and σⁱ_p is the scattering potential XS of isotope i in the mixture.

The IR parameter of hydrogen is assumed to be 1 because the MG XS library is generated based on a hydrogen moderator. In the MCS calculations, all isotopes except ²³⁸U are set to have constant scattering XS only (no absorption).

1. After the σ^U238_{a, g} versus σ₀ table is generated, another neutron slowing down problem with a mixture of ²³⁸U/¹H/X is solved where X is the isotope whose λ^X_g is to be determined.

2. The IR parameter λ^X_g for the isotope X is calculated as follows: (2) $${\lambda }_g^X=\frac{N^R{\sigma }_{0,g}^{\mathrm{eq}}-N^{H1}{\sigma }_p^{H1}}{N^X{\sigma }_p^X},$$(2)

where σ^eq_{0, g} is an equivalent background XS of ²³⁸U determined by a table interpolation using the σ^U238_{a, g}calculated in step 2 as a parameter, N^Xσ^X_p is a macroscopic potential XS of isotope X, and N^H1σ^H1_p is the macroscopic potential XS of hydrogen.

Equation (2)(2) $${\lambda }_g^X=\frac{N^R{\sigma }_{0,g}^{\mathrm{eq}}-N^{H1}{\sigma }_p^{H1}}{N^X{\sigma }_p^X},$$(2) can be rearranged as Equation (3)(3) $${\sigma }_{0,g}^{\mathrm{eq}}=\frac{N^{H1}{\sigma }_p^{H1}+{\lambda }_g^X{N}^X{\sigma }_p^X}{N^R}.$$(3) . In Equation (3)(3) $${\sigma }_{0,g}^{\mathrm{eq}}=\frac{N^{H1}{\sigma }_p^{H1}+{\lambda }_g^X{N}^X{\sigma }_p^X}{N^R}.$$(3) , the IR parameter adjusts the potential XS of nuclide X to fit the equivalent background XS compared to the σ^U238_{a, g} versus σ₀ table which is generated with a hydrogen moderator. (3) $${\sigma }_{0,g}^{\mathrm{eq}}=\frac{N^{H1}{\sigma }_p^{H1}+{\lambda }_g^X{N}^X{\sigma }_p^X}{N^R}.$$(3)

The IR parameters are calculated and stored in advance for the 10 selected isotopes which covers the range of atomic masses and tabularized such that STREAM can interpolate λ^X_g as a function of the mass of given isotopes. The STREAM IR parameters are shown in Figure 2. The use of Monte Carlo code in the IR parameter calculation automatically takes into account the neutron upscattering in the epi-thermal energy range which is not a simple problem in conventional methods. A detailed description for this issue is presented in the next section.

Figure 2. IR parameter for 10 isotopes.

2.4. Resonance upscattering correction

Recently, it has been reported that a correct modeling of neutron upscattering events in the low lying resonance energy range can lead up to 10% difference in Doppler coefficients [13]. The discrepancy comes from an assumption in the asymptotic elastic scattering model of NJOY, which assumes that a target nuclide is at rest. The asymptotic elastic scattering model only considers the downscattering event depending on the mass of the target nuclide. Many of the deterministic codes using NJOY cannot consider the resonance upscattering event.

Monte Carlo codes [12,20] also have a similar problem. Monte Carlo codes consider the thermal motion of target nuclides if the energy of an incident neutron is comparable to the thermal motion of the target nuclide. The energy of a target nuclide can be sampled with Maxwellian distribution by assuming a constant XS of the target nuclide. This method is referred to as sampling the velocity of the target nucleus (SVT) [12]. This assumption is valid for the thermal energy region. However, it can make problems for the resonance energy region because XS of resonance dramatically changes. Therefore, additional rejection techniques are necessary to correct the collision probability depending on XS. This method is known as Doppler broadening rejection correction (DBRC) [21].

In order to consider the resonance upscattering event in STREAM, MG XS of ²³⁸U in the STREAM library is modified with results by the Monte Carlo method. The asymptotic scattering kernel and DBRC method are implemented in MCS. Two sets of resonance XS libraries of ²³⁸U are generated by using the two kinds of scattering kernel as a function of the background XSs, energy group, and temperatures. The upscattering correction factor is calculated as a ratio of the ²³⁸U XSs, and it is used to correct the XS from NJOY as follows: (4) $${\sigma }_{b,g,t}^{\text{'}}={\sigma }_{b,g,t}^{\mathrm{NJOY}}{f}_{b,g,t}^{\mathrm{RUC}}={\sigma }_{b,g,t}^{\mathrm{NJOY}}\frac{{\sigma }_{b,g,t}^{\mathrm{DBRC}}}{{\sigma }_{b,g,t}^{\mathrm{Asy}.}},$$(4) where, b, g and t are indexes of the background XS, energy group and temperature; σ^NJOY_{b, g, t} is XS of ²³⁸U generated by NJOY; f^RUC_{b, g, t} is the resonance upscattering correction factor; σ^DBRC_{b, g, t} and σ^Asy._{b, g, t} are XSs of ²³⁸U generated by using the DBRC method and the asymptotic scattering kernel of MCS, respectively.

When the upscattering correction factor is generated, MCS performs fixed source calculations in a homogeneous medium. In STREAM, the correction factor is saved as a library, and used when the upscattering correction is performed. A similar work was already accomplished in CASMO5 [13].

The exact scattering kernel and the DBRC method are widely known as accurate methods to treat the resonance scattering event. However, still many of the Monte Carlo codes do not have the features officially. Therefore, an additional correction factor is generated to consider the SVT method as follows: (5) $${\sigma }_{b,g,t}^{\text{'}\text{'}}={\sigma }_{b,g,t}^{\mathrm{NJOY}}{f}_{b,g,t}^{\mathrm{SVT}}={\sigma }_{b,g,t}^{\mathrm{NJOY}}\frac{{\sigma }_{b,g,t}^{\mathrm{SVT}}}{{\sigma }_{b,g,t}^{\mathrm{Asy}.}},$$(5) where f^SVT_{b, g, t} is the SVT correction factor; and σ^SVT_{b, g, t} is XSs of ²³⁸U generated by using the SVT method of MCS.

The SVT correction factor is used optionally when a result of STREAM is compared to that of the Monte Carlo code. The SVT correction factor is also generated for only ²³⁸U. In STREAM code, XS library from NJOY is used as the default option. The impacts of the resonance upscattering correction factor and the SVT correction factor are presented in Section 5.

The IR parameter is also generated by the three kinds of scattering kernels which are the asymptotic scattering kernel, SVT method and DBRC method. Each IR parameter is used with corresponding options.

3. Rational approximation for fuel escape probability

3.1. Conventional rational approximation with Dancoff factor

In the resonance treatment methods using equivalence theory, the prediction of fuel escape probability is key for high accuracy of effective XSs in heterogeneous problems. Traditionally, the fuel escape probability has been approximated by Carlvik's two-term rational approximation [1] with the Dancoff factor for cylindrical geometry. In general, the N-term rational approximation for fuel-to-fuel collision probability P_ff in a lattice configuration can be written as (6) $$P_{ff}=\sum _{n=1}^N\frac{{\beta }_{n}X}{X+{\alpha }_n},$$(6) where coefficients α_nand β_n are calculated with the Dancoff factor, and optical length X is defined as $X={\Sigma }_t/\phantom{{\Sigma }_t{\Sigma }_e}{\Sigma }_e$. The escape XS Σ_e varies depending on geometry.

The current practice of using Dancoff factors and rational approximations has sources of inaccuracy.

1. A collision probability method with the Wigner–Seitz approximation is used in computing Dancoff factors, which can cause errors due to the cylindrical geometry assumption rather than the real square lattice. This problem can be mitigated by adopting MOC with a square lattice when computing the Dancoff factors [2,22].

2. Even with the MOC Dancoff factors, it is still difficult to apply the Dancoff factor to general geometry since the analytical rational approximations are derived for cylindrical or slab fuel.

3. The Dancoff factor and rational approximations are based on the black pin assumption, which can cause errors for real fuel pins which are gray.

3.2. Methodology of optimum rational approximation

A new optimum rational approximation has been developed which uses a numerical approach to find a more accurate fuel escape probability. The optimum rational approximation is described as follows. The transport equation by the collision probability method for the flux in the fuel is expressed as (7) $${\Sigma }_t^f{\phi }_f{V}_f=q_f{V}_f{P}_{ff}+\sum _{K\ne f}{q}_K{V}_K{P}_{Kf},$$(7) where Σ^f_t is the macroscopic total XS; φ_f is the flux in the fuel; V_f and V_K are the volume of fuel and region K, respectively; q_f and q_K are neutron sources born in the fuel and region K, respectively; P_ff is the fuel-to-fuel collision probability; and P_Kf is the K-to-fuel collision probability.

From Equation (7)(7) $${\Sigma }_t^f{\phi }_f{V}_f=q_f{V}_f{P}_{ff}+\sum _{K\ne f}{q}_K{V}_K{P}_{Kf},$$(7) , the fuel-to-fuel collision probability can be expressed as (8) $$P_{ff}=\frac{{\Sigma }_t^f{\phi }_f}{q_f}-\frac{{\sum }_{K\ne f}{q}_K{V}_K{P}_{Kf}}{q_f{V}_f}.$$(8) When q_K is set to be zero, the fuel-to-fuel collision probability is computed as (9) $$P_{ff}=\frac{{\Sigma }_t^f{\phi }_f}{q_f}.$$(9)

The fuel-to-fuel collision probability in Equation (9)(9) $$P_{ff}=\frac{{\Sigma }_t^f{\phi }_f}{q_f}.$$(9) is computed with solution of MOC fixed source calculation. In the MOC calculation, the macroscopic total XSs of coolant and cladding regions are defined as (10) $${\Sigma }_t^K\approx \lambda {\Sigma }_p^K=\sum _{i\in K}{N}^i{\lambda }^i{\sigma }_p^i,$$(10) where Nⁱ is number density, λⁱ is the IR parameter, and σⁱ_p is the potential XS of isotope i.

For the total XS of the fuel, several values of total XS are chosen in order to include various optical lengths. For example, the optical lengths can be 10⁻², 10⁻¹ and 10¹. After the total XSs are assigned in all regions, fixed source calculations are performed with the MOC module. In the fixed source calculations, neutron sources are generated only in the fuel region and the magnitude of the source is set to be λΣ^f_p. The fixed source calculations are performed for preselected levels of optical lengths. In STREAM, default conditions of the fixed source calculation are 0.03 cm ray distance, 64 azimuthal angles and 3 polar angles. A flux converge criteria is that maximum point flux error <10⁻⁵. With the converged solution, a fuel-to-fuel collision probability is computed by Equation (9)(9) $$P_{ff}=\frac{{\Sigma }_t^f{\phi }_f}{q_f}.$$(9) which is the ratio of the total number of collided neutrons in the fuel region to the total number of neutron sources. Note that the magnitude of the source is not important to calculate the fuel-to-fuel collision probability. Any non-zero value can be used as source in the fuel. From the fixed source calculations, fuel-to-fuel collision probabilities for each level of optical length are tabularized with pairs of (X_l, P^MOC_{ff, l}) where l denotes the level of optical length and P^MOC_{ff, l} is the corresponding fuel-to-fuel collision probability from the MOC solution. The coefficients α_n and β_n in rational expressions of fuel-to-fuel collision probability are computed by the least square fitting to MOC solutions to minimize the difference between the MOC solution and the rational form of the escape probability as follows: (11) $${∥\Delta P_{ff}∥}^2=\sum _{l=1}^L{\left[P_{ff,l}^{\mathrm{MOC}}-\sum _{n=1}^N\frac{{\beta }_n{X}_l}{X_l+{\alpha }_n}\right]}^2,$$(11) where L is the total number of optical length levels.

Figure 3 shows the fuel-to-fuel collision probability from MOC and fitted rational expression.

Figure 3. Fuel-to-fuel collision probability.

Equations (10)(10) $${\Sigma }_t^K\approx \lambda {\Sigma }_p^K=\sum _{i\in K}{N}^i{\lambda }^i{\sigma }_p^i,$$(10) and (11(11) $${∥\Delta P_{ff}∥}^2=\sum _{l=1}^L{\left[P_{ff,l}^{\mathrm{MOC}}-\sum _{n=1}^N\frac{{\beta }_n{X}_l}{X_l+{\alpha }_n}\right]}^2,$$(11) ) are applicable to an infinite array of fuel pin or isolated fuel pin problems. Therefore, for practical fuel assembly problems with various types of fuel pins, the following procedures were developed:

1. The optimum rational approximation method is applied to each pin type.

2. The infinite array Dancoff factors (Γ_∞) for each pin type are computed by the enhanced neutron current method [22] with reflective boundary condition on unit pincells.

3. The individual Dancoff factors (Γ_i) of each fuel pin in the fuel assembly are calculated using the enhanced neutron current method.

4. The effective XSs from the optimum rational approximation method are adjusted using the two Dancoff factors as follows: (12) $${\sigma }_{x,i}^{\mathrm{opt}}={\sigma }_{x,\infty }^{\mathrm{opt}}\frac{{\sigma }_x\left(\lambda {\sigma }_p+{\Gamma }_i{\sigma }_e\right)}{{\sigma }_x\left(\lambda {\sigma }_p+{\Gamma }_{\infty }{\sigma }_e\right)},$$(12)

where σ^opt_{x, i} is the effective XS of pin i for reaction x using the optimum escape probability; σ^opt_∞ is the effective XS of a fuel pin in an infinite array using the optimum escape probability; σ_x is the XS; (λσ_p + Γ_iσ_e)is the background XS of resonance isotope in an individual fuel pin i; and (λσ_p + Γ_∞σ_e)is the background XS of resonance isotope in an infinite array of pins.

The effective XS is corrected by the ratio of XSs in individual pins and in an infinite array of pins to consider the interaction among fuel and non-fuel cells in the fuel assembly. This method makes it possible to perform fast resonance treatment for a fuel assembly problem because the fixed source calculation for assembly is needed only for computing individual Dancoff factors. The P_ff calculation is performed for each pin type with a reflective boundary condition.

Generally, the subgroup method codes [4,5] perform fixed source calculations for each category, energy group and subgroup level. Even though a relatively low number of iterations is not needed for each fixed source calculation, the total elapsed time for the fixed source calculation is 30%–50% of the total lattice transport computation time. Generally in deterministic neutron transport codes, most of the computing time is consumed by transport solver such as ray tracing and iteration calculations. The runtime can be briefly estimated by assuming that ray tracing is the most time consuming calculation: the number of transport sweeps for resonance calculations = (3 iterations per fixed source calculation) × (13 resonance energy groups) × (5 subgroup levels) × (2 categories) = 390. The number of transport sweeps for eigenvalue calculations = (10 iterations with acceleration) × (70 energy groups) = 700. The total elapsed time for the fixed source calculation = 390/(390 + 700) = 36%. There are a couple of acceleration techniques [23,24], and those enable to converge eigenvalue solutions within about 10 iterations.

The gray resonance treatment method [10] also requires the assembly-wise fixed source calculations for each energy group and level of total XS. On the contrary, the optimum rational approximation in this paper computes the individual Dancoff factor by assembly-wise fixed source calculation for each energy group only, and smaller size pin-wise fixed source calculations are performed for each energy group, level of optical length and pin type. It should be noted that 1–3 types of fuel pins are usually used in an assembly. The elapsed time of pin-wise fixed source calculation is negligible compared to the time for assembly fixed source calculation. The new method can reduce the computation time for effective XS and can also be applicable to general geometry problems.

4. Resonance treatment methodology for resonance isotopes in cladding region

Usually, the resonance treatments are applied only to the resonance isotopes in the fuel region since isotopes in the fuel region such as ²³⁸U are the dominant resonance absorbers in the reactor. However, there are also resonance isotopes in the cladding region like zirconium. Not only is the amount of Zr non-negligible in the reactor but also Zr has large resonances between 100 and 0.1 MeV. If the resonance scattering and the absorption of Zr is not treated properly, it can cause as much as 100–300 pcm errors in reactivity predictions. Traditionally, the resonance in Zr was either simply ignored or a typical background XS was used to generate the XS of Zr when running NJOY [1,6]. In this paper, a new approach for Zr resonance treatment is introduced to overcome the inaccuracy of the traditional methods.

One of the difficulties in modeling the Zr resonance effect is the fact that an analytic expression of rational approximation is not available for the cladding region due to its ring shape. In this paper, the annular shape of the cladding region is approximated as an infinite slab as shown in Figure 4.

Figure 4. Geometry approximation for cladding region.

The Roman's rational approximation [1] can now be used for isolated infinite slab geometry as follows: (13) $$p_{ff}=X\left(\frac{1.1}{X+1.4}-\frac{0.1}{X+5.4}\right).$$(13)

Similarly with Carlvik's rational approximation, the cladding-to-cladding collision probability in the lattice configuration can be derived from Equation (14)(14) $$\begin{array}{c}\hfill P_{cc}=p_{cc}+\frac{X{(1-p_{cc})}^2}{X(1-p_{cc})+A}=X\left(\frac{\beta }{X+{\alpha }_1}+\frac{1-\beta }{X+{\alpha }_2}\right).\end{array}$$(14) with the Dancoff factor for cladding. The subscript ff is replaced by cc in order to indicate that the target region is cladding. (14) $$\begin{array}{c}\hfill P_{cc}=p_{cc}+\frac{X{(1-p_{cc})}^2}{X(1-p_{cc})+A}=X\left(\frac{\beta }{X+{\alpha }_1}+\frac{1-\beta }{X+{\alpha }_2}\right).\end{array}$$(14)

The coefficients of P_cc are calculated as follows: (15) $${\alpha }_{1,2}=\frac{1.7A+1.89\mp \sqrt{A^2+4.536A+3.5721}}{\left(A+1\right)},$$(15) (16) $$\beta =\frac{\frac{5.8A+7.56}{A+1}-{\alpha }_1}{{\alpha }_2-{\alpha }_1},$$(16) (17) $$A=\frac{\Gamma }{1-\Gamma },$$(17) where Γ is the Dancoff factor.

The Dancoff factor of the cladding region is computed using the enhanced neutron current method [22]. Although the optimum rational approximation can be used for the resonance treatment in the cladding region too, it was observed that the Roman's rational approximation with the Dancoff factor could produce sufficiently accurate results.

5. Verifications

5.1. Accuracy improvement by extended resonance energy range

In this section, the impact of the extension of the resonance treatment energy range is evaluated. The extended resonance energy ranges for three isotopes are shown in Table 2.

Table 2. Extended resonance energy range.

Isotope	Extended energy range [eV]
^235U	1.1213–4
^242Pu	2.1–4
^241Am	0.3–4

The isotopic compositions of homogeneous problems are shown in Table 3. The compositions were chosen in such a way that the background XS of ²³⁸U became around 50–60 barn which is a typical range of background XS for PWR fuel pins. The lower boundary of the resonance treatment energy range of ²³⁸U is 4 eV because there is no resonance below 4 eV for ²³⁸U.

Table 3. Material composition.

Number density [#/barn-cm]
Case #	^1H	^238U	^235U	^242Pu	^241Am
1	2.5	1	0.04	–	–
2	2.5	1	–	0.0105	–
3	2.5	1	–	–	0.015

The calculated effective capture XSs of resonance isotopes are compared to the reference values calculated from a Monte Carlo simulation. The capture XS differences are shown in Figures 5–7. In the figures, “Extended LIB.” means XS library with resonance integral below 4 eV. It is easily noted that the extended library significantly improves the accuracy of computed effective XSs of all three resonance isotopes. From the results, it can be noted that the extended resonance integral can improve the accuracy in spite of the underlying assumptions which are ignoring upscattering, usage of same escape probability and intermediate resonance parameters with those of the last energy group 27 (4–9.877 eV).

Figure 5. Case 1: ²³⁵U capture XS error.

Figure 6. Case 2: ²⁴²Pu capture XS error.

Figure 7. Case 3: ²⁴¹Am capture XS error.

5.2. Accuracy improvement by optimum rational approximation

In this section, the accuracy of the new optimum rational approximation for fuel-to-fuel collision probability is demonstrated by comparing it with two conventional methods: (1) Wigner's one-term rational approximation and (2) Carlvik's two-term approximation. A model pincell problem is developed for the verification calculations. A UO₂ pin problem with 3.9 wt% uranium is designed for the verification of the cladding resonance treatment method. Figure 8 shows the geometry of the model problem and the isotopic compositions are shown in Table 4.

Figure 8. Geometry of pincell problem.

Table 4. Material composition of pincell problem.

Region	Isotope	Number density [#/barn-cm]
Moderator	^1H	0.0442326
	^16O	0.0221163
Cladding	^90Zr	0.0205952
	^91Zr	0.0044913
	^92Zr	0.0068651
	^94Zr	0.0069571
Fuel	^235U	0.0009035
	^238U	0.0219827
	^16O	0.0457725

Dancoff factors are used to compute coefficients in the rational approximations of the two conventional methods, and the enhanced neutron current method is used to compute the Dancoff factors. In the optimum rational approximation, a three-term approximation as well as two-term approximation is introduced with the same principle of minimizing the least square fitting error compared with the exact values of the probabilities. The coefficients from each method are presented in Table 5. The calculation results are presented for group 13 (1.503 × 10⁴ to 2.478 × 10⁵ eV) which is the first resonance group. The reference values of probabilities calculated by the MOC module are shown in Figure 9 for 150 points of optical lengths.

Table 5. Coefficients of rational approximations (group 13).

Coefficients	Wigner's one-term	Carlvik's two-term	Optimum two-term	Optimum three-term
β1	1.00000	1.06781	1.27600	1.19595
β2	–	−0.06781	−0.27600	−0.45000
β3	–	–	–	0.25405
α1	0.66902	0.90744	1.05000	1.05263
α2	–	4.42354	2.44000	4.24660
α3	–	–	–	5.20000

Figure 9. Fuel-to-fuel collision probability calculated by MOC.

Figure 10 shows the difference of each rational approximaiton compared against the reference values from MOC calculations. All the methods show very good agreement in both extremes of the optical thicknesses, i.e., at 2 × 10⁻³ and 3 × 10¹. In the intermediate range, however, the Wigner and Carlvik's approximations show relatively large differences. In that range, the difference of Wigner's approximation is larger than 1% and the maximum difference is about 7%. The Carlvik's approximation shows a much smaller difference compared to the Wigner's approximation, but it still shows a difference larger than 1%. The optimum two-term rational approximation always shows less than 1% difference. The optimum three-term rational approximation can further reduce the difference to less than 0.3%.

Figure 10. Comparison of fuel-to-fuel collision probabilities.

The comparison of effective multipication factors is summarized in Table 6. The comparison of the effective multiplication factor may not proper aspect to judge the accuracy of each method because many approximations in equivalence theory can cause errors and it should be totally included in the final product of the effective multiplication factor. However, the difference of multiplication factors can be assumed to be caused by differences in fuel-to-fuel approximation. The eigenvalue, using Wigner's rational approximation, shows 622 pcm difference compared to the MCNP result and about 700 pcm difference compared to the other multi-term rational approximation. Wigner's one-term approximation showed large errors in fuel-to-fuel collision probability (the maximum difference of which is 7%) and it can cause a large error in the multiplication factor. The eigenvalue using Calvik's two-term rational approximation shows 127 pcm difference compared to the MCNP results and about 60 pcm difference compared to the optimum two-term and three-term results. There is no significant difference between optimum two-term and three-term results. The difference compared to the MCNP results is about 60 pcm and difference between them is 15 pcm. The multiplication factors show converging behavior as using more accurate rational approximation in view of the fuel-to-fuel collision probability.

Table 6. k_eff comparison for pincell problem.

	Reference MCNP6	Wigner's one-term	Carlvik's two-term	Optimum two-term	Optimum three-term
keff	1.39020 (0.00007)	1.39642	1.38893	1.38967	1.38953
Diff. [pcm]	–	622	−127	−53	−67

In the optimum rational approximations, the least square fitting process is necessary. The accuracy of fuel-to-fuel collision probablilty can be increased as the number of terms increases. However, it should be noted that many rational terms require finer grids of background XSs in the XS library, and an elapsed time in the least square fitting process will be increased. The elapsed time also depends on the least square fitting algorithm. Therefore, it is important to decide the proper number of rational terms. There is about 15 pcm difference between eigenvalues of the optimum two-term and three-term. It can be assumed that the accuracy of fuel-to-fuel collision probability of optimum two-term is enough. In the STREAM code, the optimum two-term rational approximation is used as a default option for fast execution time as well as high accuracy.

5.3. Application of resonance treatment to isotopes in cladding region

The identical problem described in Figure 8 and Table 4 is solved in order to verify the resonance treatment method in the cladding region. The absorption reaction rates and XSs in the cladding region are compared with three different approaches. The first approach assumes the infinite dilution condition of the resonance self-shielding effect on the background cross section of the cladding material. The second approach is to use a representative background XS for Zr isotopes in a typical PWR pin. The third approach is to explicitly model the resonances in the cladding region as described in Section 4. The solutions of the three approaches are compared to that of the Monte Carlo calculation by MCNP6 in Table 7.

Table 7. Reactivity difference contributions of Zr absorption reaction rates [Unit : pcm].

	w/o Zr treatment	σb = 300barn	w/ Zr treatment
^90Zr	−13	0	2
^91Zr	−156	89	2
^92Zr	−14	11	1
^94Zr	−7	5	0
Total	−190	105	5

With no resonance treatment in the cladding region, the reactivity difference caused by Zr isotopes is in the order of 190 pcm. It is noted that most of the difference comes from ⁹¹Zr. When a typical value of background XS, 300 barn, is used for Zr isotopes, the total difference is reduced to 105 pcm which is still a non-negligible magnitude of difference. It should be noted that, with the explicit resonance treatment in the cladding region, most of the difference disappears and the total difference decreases to 5 pcm. Figures 11 and 12 show additional comparisons of reaction rates and XSs of ⁹¹Zr. The explicit resonance treatment of Zr isotopes shows the smallest difference in reaction rates and effective cross sections.

Figure 11. Reaction rate comparison for ⁹¹Zr.

Figure 12. Capture XS comparison for ⁹¹Zr.

5.4.Code-to-code comparison

In this section, more practical problems are designed to further investigate the performance of the combined resonance treatment methods presented in this paper. The first model is a set of UO₂ pincell problems. The geometry information and material composition are presented in Tables 8 and 9. There are variations in pin radii and pin pitches to cover a range of pincell geometry. The composition variations of 0.7–5 wt% enrichments are selected to cover the practical range of fuel design for commercial PWRs. The natural zirconium in the cladding region requires the resonance treatment. Light water is used as a moderator. Temperature is constant at 293.6K in every region of the cell.

Table 8. Geometry description of UO₂ pincell problem.

	Geometry type
	A	B	C
RFuel [cm]	0.49	0.42	0.38
RCladding [cm]	0.55	0.47	0.44
Pin Pitch [cm]	1.47	1.28	1.20

Table 9. UO₂ pincell problem description.

		Fuel : UO2		Cladding: Nat. Zr	Moderator : H2O
Case #	Geometry type	Enrichment [wt%]	Density [g/cm^3]	Density [g/cm^3]	Density [g/cm^3]
1	A	0.7
2		1.8
3		2.9
4		3.9
5		5.0
6	B	0.7
7		1.8
8		2.9
9		3.9	10.96	6.52406	0.661
10		5.0
11	C	0.7
12		1.8
13		2.9
14		3.9
15		5.0

The solutions of the model problems are compared against the reference solutions from MCNP6. The optimum two-term rational approximation is used for the isotopes in the fuel region and the resonance isotopes in the cladding region are treated explicitly by the cladding resonance treatment method. In addition to MCNP6 and STREAM, the same problems are solved by KENO and NEWT neutron transport codes in the SCALE6.1 code package [3] for a demonstration of the level of accuracy of STREAM. KENO has MG and continuous energy libraries. In this comparison, MG library is used to compare the performance of resonance treatment methods. For clarification, KENO with MG library is named as KENO-MG in this paper. The information of each computational code is summarized in Table 10. The main difference between SCALE codes and STREAM is the resonance treatment methods. BONAMI/CENTRM modules in SCALE codes are used to compute the effective XSs. CENTRM performs PW flux calculations with 1 D geometry, and condenses the PW XS to MG XS. In the calculation, a rectangular pincell is approximated as a cylindrical model by Wigner–Seitz approximation. KENO uses the Monte Carlo method and NEWT uses the extended step characteristics (S_n) method for the eigenvalue transport calculation.

Table 10. Comparison of KENO-MG, NEWT and STREAM.

	KENO-MG	NEWT	STREAM
Nuclear data	ENDF/B-VII.0
XS processing code	AMPX		NJOY
Energy group	238-group		70-group
Main resonance treatment method	Equivalence theory and PW flux calculation		Equivalence theory
Method to treat spatial self-shielding effect	Wigner one-term approximation (BONAMI) and direct 1-D PW calculation (CENTRM)		Optimum two-term approximation
Geometry in resonance treatment	1D cylindrical cell		Exact 2D cell
Transport solver	Monte Carlo	Extended step characteristics	Method of characteristics
Transport solver condition	20,000 histories/	24 in quadrature set/	0.01 cm ray distance/
	25 inactive cycles/	6 order of Sn levels	3 polar angles/
	500 active cycles		128 azimuthal angles
Scattering order	P5	P2	Transport corrected P0

The calculation results are summarized in Table 11. The root mean square (RMS) difference of the STREAM code is 45 pcm while the RMS difference of the KENO-MG and NEWT codes are 582 and 648 pcm, respectively. The standard deviation of the differences of STREAM is 44 pcm which is much smaller than 81 pcm and 90 pcm of KENO-MG and NEWT. KENO-MG and NEWT show similar results with each other because they use an identical resonance treatment method. However, the difference of SCALE codes and STREAM is quite large. There may be many sources of errors. First of all, the geometric assumption in the resonance treatment process can cause errors. The error will be included in the calculation of the Dancoff factor, background XS and PW flux. Second, the difference XS processing program AMPX can cause the difference. MCNP and STREAM use NJOY for XS processing. For the test problems in Table 11, the solutions of the STREAM code are more accurate in an average sense and also the spread of difference is smaller than in the other two codes.

Table 11. k_eff results comparison for UO₂ pincell problem 1–15.

Case #	MCNP6 (SD.)	KENO-MG* (Diff.^**)	NEWT (Diff.^**)	STREAM (Diff.^**)
1	0.91509	0.91074	0.91047	0.91576
	(0.00002)	(−435)	(−462)	(67)
2	1.21181	1.20595	1.20567	1.21205
	(0.00003)	(−586)	(−614)	(24)
3	1.32194	1.31581	1.31535	1.32180
	(0.00003)	(−613)	(−659)	(−14)
4	1.38006	1.37372	1.37336	1.37962
	(0.00003)	(−634)	(−670)	(−44)
5	1.41667	1.41050	1.40988	1.41594
	(0.00003)	(−617)	(−679)	(−73)
6	0.91349	0.90931	0.90868	0.91415
	(0.00002)	(−418)	(−481)	(66)
7	1.21384	1.20791	1.20737	1.21428
	(0.00003)	(−593)	(−647)	(44)
8	1.32541	1.31873	1.31842	1.32552
	(0.00003)	(−668)	(−699)	(11)
9	1.38420	1.37782	1.37699	1.38394
	(0.00003)	(−638)	(−721)	(−26)
10	1.42092	1.41460	1.41367	1.42045
	(0.00003)	(−632)	(−725)	(−47)
11	0.90844	0.90430	0.90369	0.90911
	(0.00002)	(−414)	(−475)	(67)
12	1.21493	1.20931	1.20849	1.21546
	(0.00003)	(−562)	(−644)	(53)
13	1.32937	1.32329	1.32234	1.32957
	(0.00003)	(−608)	(−703)	(20)
14	1.38952	1.38336	1.38229	1.38943
	(0.00003)	(−616)	(−723)	(−9)
15	1.42695	1.42079	1.41970	1.42671
	(0.00003)	(−616)	(−725)	(−24)
RMS. Diff. [pcm]		582	648	45
SD.[pcm]		81	90	44

*Standard deviation of MG Monte Carlo code KENO-MG results are 14–21 pcm.

^**Difference is calculated by (k_eff − k^MCNP6_eff) × 10⁵ (unit: pcm).

5.5. Mosteller benchmark

In addition to the 15 pincell problems, the Mosteller benchmark problem [25] is solved for further verification. The benchmark consists of not only UO₂ fuel but also reactor-recycle mixed oxide fuel (MOX) fuel and weapons-grade MOX fuel. Each type of fuel has a wide range of enrichments in hot zero power (HZP) and hot full power (HFP) conditions. Therefore, the benchmark is suitable for verifying the resonance treatment methods implemented in STREAM. In the verification, MCNP6 is used as a reference. However, the official version of MCNP6 does not have the capability to model the resonance upscattering event. The DBRC method has been implemented in MCNP5 by modifying the code in order to compare the results of STREAM. In the verification, three options of STREAM are used:

1. STREAM: Original STREAM with MG library generated by NJOY (asymptotic scattering kernel).

2. STREAM-SVT: STREAM with SVT factor.

3. STREAM-DBRC: STREAM with DBRC factor.

The eigenvalue results from STREAM and STREAM-SVT are compared to those from MCNP6, and the eigenvalue from STREAM-DBTC is compared to MCNP5-DBRC. The calculation results of the UO₂ pincell at HZP and HFP conditions are summarized in Tables 12 and 13.

Table 12. k_eff for UO₂ fuel (Mosteller benchmark).

	HZP				HFP
Enrichment	MCNP6	STREAM	STREAM - SVT		MCNP6	STREAM	STREAM - SVT
[wt%]	(SD.)	(Diff.*)	(Diff.*)		(SD.)	(Diff.*)	(Diff.*)
0.711	0.66563	0.66648	0.66662		0.65982	0.66038	0.66074
	(0.00001)	(85)	(99)		(0.00002)	(56)	(92)
1.6	0.96085	0.96181	0.96201		0.95268	0.95323	0.95376
	(0.00002)	(96)	(116)		(0.00002)	(55)	(108)
2.4	1.09903	1.09988	1.10011		1.08996	1.09033	1.09092
	(0.00002)	(85)	(108)		(0.00002)	(37)	(96)
3.1	1.17704	1.17773	1.17797		1.16745	1.16768	1.16830
	(0.00002)	(69)	(93)		(0.00002)	(23)	(85)
3.9	1.23972	1.24021	1.24046		1.22983	1.22980	1.23044
	(0.00002)	(49)	(74)		(0.00002)	(−3)	(61)
4.5	1.27514	1.27555	1.27579		1.26504	1.26494	1.26560
	(0.00002)	(41)	(65)		(0.00002)	(−10)	(56)
5	1.29950	1.29983	1.30007		1.28929	1.28912	1.28978
	(0.00002)	(33)	(57)		(0.00003)	(−17)	(49)
RMS. Diff. [pcm]		69	90		–	35	81

*Difference is calculated by (k_eff − k^MCNP6_eff) × 10⁵ (unit: pcm).

Table 13. k_eff for UO₂ fuel with the resonance upscattering correction (Mosteller benchmark).

	HZP			HFP
Enrichment	MCNP5 - DBRC	STREAM - DBRC		MCNP5 - DBRC	STREAM - DBRC
[wt%]	(SD.)	(Diff.*)		(SD.)	(Diff.*)
0.711	0.66531	0.66603		0.65878	0.65946
	(0.00001)	(72)		(0.00001)	(68)
1.6	0.96033	0.96115		0.95121	0.95187
	(0.00002)	(82)		(0.00002)	(66)
2.4	1.09847	1.09914		1.08827	1.08878
	(0.00002)	(67)		(0.00002)	(51)
3.1	1.17644	1.17695		1.16569	1.16605
	(0.00002)	(51)		(0.00002)	(36)
3.9	1.23909	1.23940		1.22799	1.22811
	(0.00002)	(31)		(0.00002)	(12)
4.5	1.27450	1.27472		1.26323	1.26324
	(0.00002)	(22)		(0.00002)	(1)
5	1.29878	1.29900		1.28741	1.28740
	(0.00002)	(22)		(0.00002)	(−1)
RMS. Diff. [pcm]		55		-	43

*Difference is calculated by (k_eff − k^{MCNP5 − DBRC}_eff) × 10⁵ (unit: pcm).

The RMS differences of STREAM are 69 and 35 pcm for HZP and HFP, respectively. STREAM-SVT show slightly large RMS differences which are 90 and 81 pcm for each temperature, respectively. When the resonance upscattering is considered, STREAM-DBRC has 55 and 43 pcm differences compared to results from MCNP5-DBRC. All three of STREAM methods show good agreement in eigenvalue calculations. The maximum difference compared to corresponding MCNP results is 108 pcm.

Table 14 and Figure 13 show the FTCs for UO₂ Fuel, which were calculated by Equation (18)(18) $$\mathrm{FTC}=\left(\frac{1}{k_{\mathrm{HZP}}}-\frac{1}{k_{\mathrm{HFP}}}\right)\frac{{10}^5}{300}[\mathrm{pcm}/\mathrm{K}].$$(18) . (18) $$\mathrm{FTC}=\left(\frac{1}{k_{\mathrm{HZP}}}-\frac{1}{k_{\mathrm{HFP}}}\right)\frac{{10}^5}{300}[\mathrm{pcm}/\mathrm{K}].$$(18)

Table 14. FTC for UO₂ fuel (Mosteller benchmark).

Enrichment	MCNP6	STREAM	STREAM-SVT	MCNP5-DBRC	STREAM-DBRC
[wt%]	(SD.)	(Diff.*)	(Diff.*)	(SD.)	(Diff.^**)
0.711	−4.41	−4.62	−4.45	−4.97	−4.99
	(0.02)	(−4.77%)	(−0.91%)	(0.01)	(−0.40%)
1.6	−2.98	−3.12	−3.00	−3.33	−3.38
	(0.01)	(−4.85%)	(−0.74%)	(0.01)	(−1.60%)
2.4	−2.52	−2.65	−2.55	−2.84	−2.89
	(0.01)	(−5.18%)	(−1.13%)	(0.01)	(−1.46%)
3.1	−2.33	−2.44	−2.34	−2.61	−2.65
	(0.01)	(−4.71%)	(−0.68%)	(0.01)	(−1.32%)
3.9	−2.16	−2.28	−2.19	−2.43	−2.47
	(0.01)	(−5.22%)	(−1.20%)	(0.01)	(−1.68%)
4.5	−2.09	−2.19	−2.10	−2.33	−2.38
	(0.01)	(−5.02%)	(−0.80%)	(0.01)	(−1.84%)
5	−2.03	−2.13	−2.05	−2.27	−2.31
	(0.01)	(−4.88%)	(−0.70%)	(0.01)	(−2.01%)

*Difference is calculated by (FTC- FTC^MCNP6)/ | FTC^MCNP6| × 10² (unit: %).

^**Difference is calculated by (FTC - FTC^MCNP5-DBRC)/ | FTC^MCNP5-DBRC| × 10² (unit: %).

Figure 13. Doppler coefficients for UO₂ Fuel (Mosteller benchmark).

STREAM has differences of 4.77%–5.22% in FTC compared to MCNP6. It was reported that the deterministic transport codes compute more negative FTC compared to the Monte Carlo codes [13,25,26]. The difference comes from using different scattering kernels between two methods, which can be noted in results in STREAM-SVT. The differences in FTC are reduced in the results of STREAM-SVT. STREAM-SVT shows differences of 0.70%–1.20%. STREAM-DBRC also shows good agreement in FTC compared to MCNP5-DBTC. The differences are 0.40%–2.01%.

The eigenvalue and FTC results of reactor-recycle MOX pincells are summarized in Tables 15–17 and Figure 14.

Table 15. k_eff for reactor-recycle MOX fuel (Mosteller benchmark).

	HZP				HFP
PuO2 Content	MCNP6	STREAM	STREAM - SVT		MCNP6	STREAM	STREAM - SVT
(wt%)	(SD.)	(Diff.*)	(Diff.*)		(SD.)	(Diff.*)	(Diff.*)
1	0.94483	0.94601	0.94620		0.93520	0.93610	0.93662
	(0.00002)	(118)	(137)		(0.00002)	(90)	(142)
2	1.02072	1.02226	1.02247		1.00973	1.01099	1.01154
	(0.00002)	(154)	(175)		(0.00003)	(126)	(181)
4	1.07586	1.07770	1.07790		1.06415	1.06566	1.06620
	(0.00002)	(184)	(204)		(0.00003)	(151)	(205)
6	1.10453	1.10639	1.10659		1.09265	1.09422	1.09475
	(0.00002)	(186)	(206)		(0.00003)	(157)	(210)
8	1.12796	1.12965	1.12984		1.11611	1.11745	1.11796
	(0.00003)	(169)	(188)		(0.00002)	(134)	(185)
RMS. Diff. [pcm]		164	186		–	134	184

*Difference is calculated by (k_eff − k^MCNP6_eff) × 10⁵ (unit: pcm).

Table 16. k_eff for reactor-recycle MOX fuel with resonance upscattering correction (Mosteller benchmark).

	HZP			HFP
PuO2 content	MCNP5 - DBRC	STREAM - DBRC		MCNP5 - DBRC	STREAM - DBRC
(wt%)	(SD.)	(Diff.*)		(SD.)	(Diff.*)
1	0.94436	0.94536		0.93380	0.93478
	(0.00002)	(100)		(0.00002)	(98)
2	1.02020	1.02158		1.00829	1.00960
	(0.00002)	(138)		(0.00002)	(131)
4	1.07537	1.07702		1.06272	1.06426
	(0.00003)	(165)		(0.00003)	(154)
6	1.10408	1.10573		1.09139	1.09286
	(0.00003)	(165)		(0.00003)	(147)
8	1.12756	1.12901		1.11489	1.11613
	(0.00002)	(145)		(0.00003)	(124)
RMS. Diff. [pcm]		145		–	132

*Difference is calculated by (k_eff − k^{MCNP5 − DBRC}_eff) × 10⁵ (unit: pcm).

Table 17. FTC for reactor-recycle MOX fuel (Mosteller benchmark).

PuO2 content	MCNP6	STREAM	STREAM - SVT	MCNP5 - DBRC	STREAM - DBRC
(wt%)	(SD.)	(Diff.*)	(Diff.*)	(SD.)	(Diff.^**)
1	−3.63	−3.73	−3.60	−3.99	−3.99
	(0.01)	(−2.68%)	(0.81%)	(0.01)	(0.02%)
2	−3.55	−3.63	−3.52	−3.86	−3.87
	(0.01)	(−2.27%)	(0.89%)	(0.01)	(−0.32%)
4	−3.41	−3.49	−3.39	−3.69	−3.71
	(0.01)	(−2.50%)	(0.47%)	(0.01)	(−0.57%)
6	−3.28	−3.35	−3.26	−3.51	−3.55
	(0.01)	(−2.12%)	(0.71%)	(0.01)	(−1.13%)
8	−3.14	−3.22	−3.14	−3.36	−3.41
	(0.01)	(−2.68%)	(0.08%)	(0.01)	(−1.41%)

*Difference is calculated by (FTC - FTC^MCNP6)/ FTC^MCNP6×10² (unit: %).

^**Difference is calculated by (FTC - FTC^MCNP5-DBRC)/ FTC^MCNP5-DBRC×10² (unit: %).

Figure 14. Doppler coefficients for reactor-recycle MOX Fuel (Mosteller benchmark).

For the reactor-recycle MOX fuel problems, STREAM, STREAM-SVT and STREAM-DBRC show the order of 100–200 pcms differences. It is noted that there exists k_eff bias for reactor-recycle MOX fuels. STREAM has differences of about 100 pcm larger in the MOX fuels than in UO₂ fuels. The MOX fuels contain large amounts of ²⁴⁰Pu, and STREAM overpredicts the eigenvalues. The bias can be attributed to the resonance interference effects of ²³⁸U and ²⁴⁰Pu isotopes. The accurate treatment of the resonance interference effect is still an open problem.

Similar behavior is also noted in that STREAM-SVT shows more accurate FTC than STREAM compared to MCNP6. STREAM-DBRC and MCNP5-DBRC show comparable results in FTC. The differences are 0.02%–1.41%.

The k_eff and FTC results of weapons-grade MOX pincells are presented in Tables 18–20 and Figure 15.

Table 18. k_eff for weapons-grade MOX fuel (Mosteller benchmark).

	HZP				HFP
Enrichment	MCNP6	STREAM	STREAM -SVT		MCNP6	STREAM	STREAM - SVT
[wt%]	(SD.)	(Diff.*)	(Diff.*)		(SD.)	(Diff.*)	(Diff.*)
1	1.09027	1.09087	1.09110		1.08116	1.08133	1.08194
	(0.00002)	(60)	(59)		(0.00002)	(17)	(47)
2	1.18087	1.18119	1.18142		1.16980	1.16972	1.17035
	(0.00002)	(32)	(29)		(0.00002)	(−8)	(28)
4	1.24963	1.24963	1.24987		1.23728	1.23693	1.23756
	(0.00002)	(0)	(5)		(0.00002)	(−35)	(4)
6	1.28626	1.28608	1.28631		1.27369	1.27316	1.27377
	(0.00003)	(−18)	(−16)		(0.00002)	(−53)	(−15)
RMS. Diff. [pcm]		28	34		–	52	28

*Difference is calculated by (k_eff − k^MCNP6_eff) × 10⁵ (unit: pcm).

Table 19. k_eff for weapons-grade MOX fuel with resonance upscattering correction (Mosteller benchmark).

	HZP			HFP
Enrichment	MCNP5 - DBRC	STREAM - DBRC		MCNP5 - DBRC	STREAM - DBRC
[wt%]	(SD.)	(Diff.*)		(SD.)	(Diff.*)
1	1.09001	1.09012		1.07981	1.07978
	(0.00002)	(11)		(0.00002)	(−3)
2	1.18057	1.18039		1.16838	1.16808
	(0.00002)	(−18)		(0.00002)	(−30)
4	1.24928	1.24885		1.23584	1.23530
	(0.00002)	(−43)		(0.00003)	(−54)
6	1.28603	1.28532		1.27240	1.27158
	(0.00003)	(−71)		(0.00002)	(−82)
RMS. Diff. [pcm]		43		–	51

*Difference is calculated by (k_eff − k^MCNP5-DBRC_eff) × 10⁵ (unit: pcm).

Table 20. FTC for weapons-grade MOX fuel (Mosteller benchmark).

PuO2 content	MCNP6	STREAM	STREAM - SVT	MCNP5 - DBRC	STREAM - DBRC
(wt%)	(SD.)	(Diff.*)	(Diff.*)	(SD.)	(Diff.^**)
1	−2.56	−2.70	−2.59	−2.89	−2.93
	(0.01)	(−5.51%)	(−1.23%)	(0.01)	(−1.37%)
2	−2.67	−2.77	−2.67	−2.95	−2.98
	(0.01)	(−3.73%)	(−0.04%)	(0.01)	(−1.03%)
4	−2.65	−2.74	−2.65	−2.90	−2.93
	(0.01)	(−3.32%)	(−0.07%)	(0.01)	(−0.90%)
6	−2.55	−2.63	−2.55	−2.78	−2.80
	(0.01)	(−3.04%)	(0.06%)	(0.01)	(−0.93%)

*Difference is calculated by (FTC - FTC^MCNP6)/ FTC^MCNP6×10² (unit: %).

^**Difference is calculated by (FTC - FTC^MCNP5-DBRC)/ FTC^MCNP5-DBRC×10² (unit: %).

Figure 15. Doppler coefficients for weapons-grade MOX Fuel (Mosteller benchmark).

In the weapons-grade MOX fuel problems, all of STREAM has the order of 50 pcm difference on average. It should be noted that most of the Pu in the weapons-grade MOX fuel is ²³⁹Pu, and is different to the reactor-recycle MOX fuel. The behavior of FTC for this MOX fuel is in the same manner as UO₂ and reactor-recycle MOX fuels. STREAM-SVT and STREAM-DBRC compute comparable results with MCNP6 and MCNP5-DBRC, respectively. The difference is about 1% on average.

5.6. Fuel assembly problem

As the last verification problem, a fuel assembly (FA) problem is selected from the FA problem-2B of VERA benchmark [11]. In both STREAM and MCNP6 inputs, an assembly inter-gap is ignored in this paper for simplification. Figure 16 shows geometry of 17 × 17 fuel assembly and Figures 17–19 show descriptions of each pincell. Table 21 shows material information.

Figure 16. 17 × 17 fuel assembly.

Figure 17. Description of fuel rod.

Figure 18. Description of Instrument tube.

Figure 19. Description of guide tube.

Table 21. Description of material.

Region	Material	Temperature
Fuel	3.1 wt% UO2	600 K
Cladding	6.56 g/cc of Zircaloy-4	600 K
Moderator	0.661 g/cc of H2OBoron: 1300 ppm	600 K

Table 22 and Figure 20 show the calculation results of eigenvalues and fission source distributions of MCNP6 and STREAM codes. The STREAM code shows -5 pcm of eigenvalue difference, and a fission source distribution difference of 0.09% and 0.22% for RMS and maximum differences, respectively. The average statistical error of the fission source distribution of MCNP6 results is 0.06%. The difference in the distribution is quite small, and less than 1% difference is acceptable for computing the fission source distribution generally. It can be noted that STREAM has a 50-pcm difference on average for UO₂ pincell problems in previous sections, and STREAM also shows a similar order of difference for FA problem.

Table 22. Results for fuel assembly problem.

keff			Fission Source Distribution
MCNP6	STREAM	Diff. [pcm]	RMS. Diff. [%]	Max. Diff. [%]
1.18297 (0.00003)	1.18292	−5	0.09%	0.22%

Figure 20. Comparison of fission source distributions from MCNP6 and STREAM. (IT: instrument tube; GT: guide tube).

6. Conclusions

The resonance treatment methods implemented in the newly developed lattice physics code STREAM are presented. First, the extended resonance integral method improves the accuracy of the effective XSs below the 4 eV where many resonances of heavy isotopes exist. In particular, the extended resonance integral method reduces hundreds or thousands barns of difference in the effective cross section of ²⁴¹Am and ²⁴²Pu. Second, the optimum rational approximation makes it possible to analyze general geometry with equivalence theory. In addition, the optimum rational approximation improves the accuracy in the eigenvalues by about 500 pcm compared to Wigner's one-term approximation and 70 pcm compared to Carlvik's two-term rational approximation. Third, it is possible to treat resonances in cladding material explicitly by the cladding resonance treatment method. The resonances in cladding have an impact in eigenvalue greater than 190 pcm, and the difference is reduced by the cladding resonance treatment method. The equivalent correction methods are also discussed to treat resonance upscattering events, which are an important feature of high temperature fuels. Each method can improve the accuracy of effective XSs and it is verified through detailed comparisons of XS, reaction rate and fuel escape probability. The accuracy of the new neutron transport code STREAM with newly developed resonance treatment methods is verified with many verification problems. The verification problems include various geometries and enrichments. Well-known benchmark problems such as the Mosteller benchmark and the VERA benchmark are solved. For the UO₂ pincell problems, STREAM can compute eigenvalues with a 50-pcm difference, on average. For the Mosteller benchmark, STREAM shows ∼100 pcm of difference in UO₂ and weapons-grade MOX fuels, and 100–200 pcms in reactor-recycle MOX fuels. STREAM shows about 1% difference in FTC. In the fuel assembly benchmark, it is confirmed that the accuracy of eigenvalue in the fuel assembly geometry is comparable to that in the pincell geometry for the STREAM code. Through various verification problems, STREAM shows around 100 pcm of reactivity difference.

Acknowledgements

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP).

Additional information

Funding

This work was also partially supported by KETEP, which is funded by the Korean government Ministry of Trade, Industry and Energy [grant number 20131610101850].