Neutronics and burnup analysis of VVER-1000 LEU and MOX assembly computational benchmark using OpenMC Code

A handful of computational benchmarks that incorporate VVER-1000 assemblies having low enriched uranium (LEU) and the mixed oxide (MOX) fuel have been put forward by many experts across the world from the Nuclear Energy Agency. To study & scrutinize the characteristics of one of the VVER-1000 LEU & MOX assembly benchmarks in different states were considered. In this work, the VVER-1000 LEU and MOX Assembly computational-benchmark exercises are performed using the OpenMC software. The work was intended to test the preciseness of the OpenMC Monte Carlo code using nuclear data library ENDF/B-VII.1, against a handful of previously obtained solutions with other computer codes. The kinf value obtained was compared with the SERPENT and MCNP result, which presented a very good similarity with very few deviations. The kinf variation with respect to burnup upto 40 MWd/kgHM was obtained for State-5 by using OpenMC code for both the LEU and MOX fuel assembly. The depletion curves of isotope concentrations against burnup upto 40 MWd/kg/HM were also generated for both the LEU and MOX fuel assembly. The OpenMC results are comparable with those of benchmark mean values. The neutron energy vs flux spectrum was also generated by using OpenMC code. Based on the OpenMC results such as kinf, burnup, isotope concentrations and neutron energy spectrum, it is concluded that the OPenMC code with ENDF/B-VII.1 nuclear data library was successfully implemented. It is planned to use OpenMC code for calculation of neutronics and burnup of the VVER-1200 reactor to be commissioned in Bangladesh by 2023/2024.


Introduction
The main purpose of nuclear reactor theory is calculation of distribution of neutrons in the reactor core. From knowledge of it, we can determine the rate of fission reaction occurring in a nuclear reactor and hence reactor power and operating point (sub critical, critical or supercritical) hence, stability of fission chain reaction can be inferred. Generally, neutronic analysis is performed based on ''Deterministic" and ''Stochastic" methods. In deterministic methods the transport equation is solved as a differential equation. In stochastic methods such as Monte Carlo, discrete particle histories are tracked and averaged in a random walk directed by interaction probabilities. In Bangladesh, two Russian design VVER-1200 (2400 MWth) type nuclear reactors are under construction and to be commissioned by 2023/2024. A program has been undertaken at the Department of Nuclear Science and Engineering of Military Institute of Science and Technology, Dhaka, Bangladesh to introduce some Monte Carlo computer codes such as MCN-PX (Louis and Amin 2016), SERPENT (Mercatali et al. 2021), SuperMC (Wu et al. 2015), and OpenMC (Romano et al. 2015) Romano and Forget 2013) for neutronics calculations of different types of nuclear reactors. As a part of this program, the OpenMC code has been adopted to perform neutronics analysis of benchmark problems for two types of fuel assembly (LEU and MOX) that are typically of the advanced Russian designs. Benchmarks are specified problems used to test nuclear data libraries & validate computer codes (Alhassan et al. 2014). These act as a learning tool for engineers that enable them to perform neutronic calculations for various types of nuclear reactors. To ease the sharing of prevailing information & intelligence about the forthcoming designs of VVER-1000 LEU & MOX fuel, a VVER-1000 LEU & MOX fuel assembly computational benchmark was put forward by a group of professionals at OECD/NEA (OECD NEA 2002). The benchmark was analyzed using six codes obtained from five institutions, of which two were developed using continuous energy Monte Carlo method, i.e. MCNP-4B & MCU & for the rest, the collision probability method was used. The average values obtained from these codes are denominated as "Benchmark Mean" which is mentioned below in the result section. Additionally, the same problem was analyzed using various codes & data libraries over the past several years. Some of them are EXCEL (Thilagam 2009), APOLLO2 & TRIPOLI4 (Petrov 2013), MCNPX (Louis and Amin 2016), VISWAM (Khan et al. 2016) & GETERA (Abuqudaira and Stogov 2018). Having been launched in December 2012, OpenMC is a fairly new Monte Carloparticle transport code (Romano and Forget 2013) which gives the user an advantage of obtaining results based on the average result of three different methods namely track length, collision probability & absorption. One can also use the results obtained from any individual method instead of the average one. In this study, ENDF/B-VII.1 was used. ENDF/B-VII.1 is an evaluated nuclear data library that contains all necessary cross-section data to perform a neutronic analysis. This library has nuclear data for 423 nuclides (ENDF/B-VII.1 2012). The work is intended to perform neutronic, burnup and isotope concentrations analysis of the VVER-1000 computational benchmark using OpenMC code.

Benchmark description
The benchmark model has two fuel assemblies (LEU & MOX) of the VVER-1000 reactor. Each of the assemblies consists of 331 elementary cells of four types for LEU assembly & six types for MOX assembly inside a hexagonal lattice. Twelve Gd 2 O 3 pins are located inside each assembly as a burnable absorber at completely different positions. The layout of both the assembly types obtained using OpenMC is shown in Figs 1 and 2. The assembly has a uniform enrichment of 3.7% in 235 U except in Gd, bearing pins that contain 3.6% wt. of 235 U and 4% wt. of Gd 2 O 3 . The assembly lattice pitch & the hexagonal cell pitch are respectively 23.6 cm & 1.275 cm. All three types of unit cells along with their dimensions for each assembly are given in Figs 1 and 2. The cladding and structural materials are of a Zircaloy. Both material & geometry specifications used in the assemblies are taken from the benchmark report (OECD NEA 2002).

Geometry description
The cell type geometry specifications are given in Table 1. The fuel and non-fuel cells are shown in Fig. 3.

Methodology
The models were represented in the OpenMC using python (python 3.7) code in Jupyter notebook. Initially, one of each of the different types of rods was created. To obtain the assemblies, they were placed in a hexagonal lattice with a lattice cell pitch of 1.275 cm and an assembly pitch of 23.6 cm. Modeling was done utilizing Boolean operations to define different zones within the cells. The hexagonal lattice and two different planes in the z-axis with reflecting boundary conditions bounded the geometry, which is equivalent to the geometry being infinite in the z-axis. For thermal scattering at low energies, S(α,β) table was used. The benchmark demands a solution for a variety of states, encompassing both hot and cold conditions, as shown in Table 2. Assumptions taken during the OpenMC process are: • Pin-by-pin Model • Reflective boundary condition in all directions (x, y & z) • Finite boundary at z-axis with reflective boundary conditions. • Cross-section data library: ENDF/B-VII.1 • 1000 batches with 100 inactive batches & 11000 particles in each batch.
In OpenMC code, there are three ways to calculate k-eigenvalue, including track-length estimator, collision estimator, and absorption estimator, respectively. They are expressed in the following equations: where W is the total weight starting each generation (or batch), w j is the pre-collision weight of the particle as it enters event j, d j is the length of the j th trajectory, and υΣ f and Σ a are macroscopic neutron production cross section and absorption cross section.
The OpenMC Monte Carlo code was used to calculate the average k inf based on the combined collision estimator, track length estimator, & absorption estimator (Wu et al. 2015).

Infinite multiplication factor
The infinite multiplication factor was calculated in eigenvalue mode of the OpenMC (Version 0.12.2) monte carlo code. To obtain accurate results based on the initial guess value for the fission source distribution, analysis of the iteration method source convergence is necessary. A study into convergence of Monte Carlo criticality analysis has proved that the Shannon entropy of the fission source distribution, H src , is an effective parameter for identifying the convergence of the fission source distribution (Brown 2006). The OpenMC monte carlo code includes new tools for determination and plotting of Shannon entropy of the fission source distribution to assess problem convergence.
To ensure that the result converges, the Shannon entropy was measured at zero burnups, and an approximation of the inactive batches was determined. The Shannon entropy versus generation curve as shown in Fig. 4, aids in determining the number of inactive batches. The number of inactive batches is taken to be 100, before which the entropy curve starts showing constant behavior.
The OpenMC results of the present study were compared to the results of SERPENT code as well as the MCNP results (Mercatali et al. 2015). The k inf values obtained by using OpenMC code are given in Tables 3 and 4 at zero burnup for both LEU and MOX fuel assembly, respectively. Tables 3 and 4 reveal that OpenMC and other  codes (such as SERPENT and MCNP) have a high degree of similarity in terms of nature and value. Table 3 compares the infinite multiplication factor computed by OpenMC with ENDF/B-VII.1 nuclear data library with those of the literature values obtained from SERPENT code using ENDF/B-VII nuclear data library at zero burnup for LEU fuel assembly, and shows that OpenMC and SERPENT agree quite well in case of LEU fuel assembly. As demonstrated in Table  3, there is a little variation due to the use of different nuclear data libraries. The comparison of these two nuclear data libraries can help explain the causes. In the case of MOX fuel assembly (Table 4), the OpenMC results are comparable with those of the results obtained from SERPENT code. Using PREPRO code (Diop et al. 2007), a detailed comparison was done between the two libraries. In ENDF/B-VII.1, there were around 170 variations, with 120 of them differing by 1% or more of total cross-sections from the same isotope in ENDF/B-VII.0. Because these methods utilized various different data libraries such as ENDF/B-VI.0, JEFF2.2, JEFF 3.1, JENDL 3.2, ENDF/B-VII.0, etc., the result differs even more from those shown in benchmark problem & benchmark mean results. Another explanation is that neutron data for fuel material, non-fuel material, and the S(α, β) table was used for 579 K instead of 575 K for S2-S4 stats due to the unavailability of data at these two temperatures. Table 5 shows a comparison of different results achieved by different codes for the cold state (State-5) for LEU and MOX fuel assembly. The JEFF 3.3 nuclear data library was also used in the OpenMC code to calculate the k inf for the State-5 only ( Table 5). The difference between the maximum and the minimum k inf value for the same state is roughly 680 pcm for LEU and about 1695 pcm for MOX fuel, according to the table. All of the values for the various codes were collected from various sources that are cited in the references. We can conclude from this comparison that the OpenMC results are acceptable.

Burnup effects
The variation of k inf with respect to burnup (MWd/ kgHM) is shown in Fig. 5 for LEU and MOX fuel assembly, respectively for State-1. The time burup steps were considered upto 10 MWd/kgHM to model Gd depletion accurately. The S1 state was followed upto 40 MWd/ kgHM using the power density 108 MW/m 3 given in the benchmark problem. The benchmark values are also included in Fig. 5 for comparison purposes. At burnup of 8 MWd/kgHM the peak value of k inf is seen in Fig. 5 for LEU fuel assembly. In the case of MOX fuel assembly, the peak value of k inf is observed at burnup of ~10 MWd/ kgHM. The OpenMC results for the both LEU and MOX fuel assembly are comparable with those of benchmark mean values.

Reactivity effects
The reactivity effect was computed using the k inf values obtained from various reactor operational states at zero burnup for LEU and MOX fuel assembly is given in

Figure 5. Variation of k inf vs burnup
Tables 6 and 7. The OpenMC reactivity values are comparable with those of benchmark mean value including SERPENT value (Mercatali et al. 2015).

Isotopic composition changes with burnup
Figs 6-22 exhibits the results of various radionuclide concentrations for LEU and MOX fuel assembly with respect to burnup for the operating poisoned condition (S1). Predictor Integrator, a first-order predictor algorithm, was employed for depletion analysis. The concentration of neutron poisons Xe-135 and Sm-149 must be in equilibrium for state-1, which is an operating poisoned state.
The main difference between State-1 and State-2 is that the former requires the xenon concentration and neutron flux to remain in equilibrium, whereas the latter does not.
The purpose of building this forced equilibrium before running the simulation is to account for the xenon oscillation that occurs during a depletion analysis which has a significant impact on the reactivity (Isotalo et al. 2013).
Burnup calculation is carried out for this state only. Isotopic composition changes of nuclides 235 U, 236 U, 238 U, 239 Pu, 240 Pu, 241 Pu, 242 Pu, and 149 Sm in cell-1 and cell-24 (as shown in Figs 2 and 3) of the LEU and MOX fuel assembly as a function of burnup is depicted in Figs 6-22. These changes were obtained from burnup calculation in State-1 which is compared with the Benchmark mean value. In addition, for both the LEU and MOX fuel assemblies, two additional isotopes, 155 Gd and 157 Gd, are compared in cell-24. Because of its large thermal fission cross-section of approx. 583 barns, U-235 is used as the primary fuel in an LEU fuel assembly. As the fission process progresses and burnup increases, the concentration of U-235 in LEU assembly falls. Because Pu is the major fuel in MOX assembly, 239 PU shows similar behavior which can be seen from Fig. 8 below. In an LEU assembly, all U-235 atoms do not cause fission. Radiative capture yields U-236 in a small fraction of U-235, hence the concentration of U-236 rises with burnup. The same goes for the MOX assembly since it also contains U-235 along with Pu. Despite the fact that U-238 makes up the majority of the fuel element material in an LEU assembly, its overall fission cross-section is too low compared to that of U-235. However, due to its fast fission & thermal fission (very small amount) and radiative capture, the concentration of U-238 drops as burnup increases. By radiative capture of neutrons, the fertile atom U-238 generates U-239 in an LEU assembly. U-239 quickly emits a beta particle to transform into Np-239 and Np-239, in turn, releases a beta particle to transform into Pu-239, which is relatively stable. Thus, the concentration of Pu-239 increases with burnup in LEU assembly. And in MOX assembly, since the main fuel is Pu-239, as burnup increases the fuel gets used up and its concentration decreases. In an LEU assembly, initially, there's no existence of Pu-240, but as the burnup increases, its concentration also builds due to the neutron capture of previously generated Pu-239. In MOX assembly, the concentration of Pu-240 rises by a large margin with the increase of burnup because the radiative capture of neutrons by U-238 produces Pu-240 after subsequent beta decays. By neutron capture, some Pu-240 nuclei may in turn form Pu-241. In both the LEU and MOX assemblies, the concentration of Pu-241 rises in a similar manner to that of Pu-240. But in LEU assembly, Pu-241 just takes a little longer to grow its concentration as Pu-239 does not exist in this assembly. Pu-242 comes into being when due to the neutron capture of Pu-241. Since neither of the two assemblies had an initial concentration of Pu-242, variation of Pu-242 concentration with burnup in both the LEU and MOX assemblies indicates a similar scenario. However, because of the presence of Pu-241 in the MOX core, the concentration of Pu-242 in MOX assembly rises earlier and is greater than that of LEU assembly. The thermal neutron absorption cross-sections of Gd-155 and Gd-157 are very large. As a result, in LEU assembly, they burn extremely quickly. However, due to the harder neutron spectrum in MOX assembly, Gd isotopes burn slowly. Gd-157 depletes faster in both assemblies than Gd-155, because of its larger absorption cross-section.
The isotopic concentrations for cell-1 are shown in Figs 6-13.

Neutron flux distribution
Figs 33 and 34 present the neutron energy spectrum with respect to fluxes at three burnup cycles i.e. Beginning of Cycle (BOC), Middle of Cycle (MOC) and End of Cycle (EOC) for State-1 of LEU and MOX fuel assembly, respectively. The typical neutron energy spectrum is comprised of three kinds of spectra; thermal Maxwellian spectrum in thermal energy range, 1/E spectrum in epi-thermal energy range and fission spectrum in fast energy range. The shape of the typical spectrum has been given as results from neutron transport phenomena which are supposed to be observed in the LWR. Emitting neutrons by nuclear fission and moderation by light water are being considered in the phenomena for generation of such neutron energy spectrum (Figs 33 and 34). The neutron energy flux spectrum for each energy group has been calculated by dividing the neutron flux by the width of every energy group, φ(cm −2 . s −1 .eV −1 ) = Φ(cm −2 .s −1 )/ΔE(eV) by using OpenMC code.
As may be seen from the Figs 33 and 34, good agreements are noticed between the neutron energy flux spectrum for BOC, MOC and EOC with slight variation among them especially MOX fuel assembly. No such neutron energy spectrum is available in the benchmark report to compare with OpenMC results. For comparison purpose, a typical neutron energy spectrum of a Light Water Reactor (Tanaka et al. 2015) is shown in Fig. 35. The neutron energy spectrum (Figs 33a and 34a) of VVER-1000 for State-1 (LEU and MOX fuel assembly) is almost identical shape of light water reactor (Fig. 3). Hence, the OpenMC code was successfully used for obtaining a neutron energy spectrum.

Conclusion
The k inf values were calculated for VVER-1000 LEU & MOX assemblies that are typically of the advanced Russian designs in different reactor operating states using OpenMC code with nuclear data library ENDF/B-VII.1. The k inf values were also calculated against fuel burnup upto 40 MWd/kgHM. In addition, the isotope composition was also calculated for burnup upto 40 MWd/kgHM as per benchmark requirements. The calculated results were compared with the benchmark mean values along with the literature data. The OpenMC results showed very good agreement with the benchmark mean values alongwith other literature values. The neutron energy spectrum was successfully generated by using OpenMC code for both LEU and MOX fuel assembly for State-1. It is concluded that the OpenMC code along with the nuclear data library ENDF/B-VII.1 was successfully implemented at the Department of Nuclear Science and Engineering Department, MIST. In Bangladesh, two Russian design VVER-1200 (2400 MW th ) type nuclear reactors are under construction and to be commissioned by 2023/2024. Based on the experience achieved for implementation of OpenMC code in the field of neutronics and burnup calculations, it is planned to calculate k inf or to perform burnup calculations for VVER-1200 and other PWR and BWR by using OpenMC along with SuperMC.

Authorship contribution statement
Md.