Research Article
Print
Research Article
Neutronic calculations for the VVER-1000 MOX core computational benchmark using the OpenMC code
expand article infoMd Imtiaj Hossain‡, Abdus Sattar Mollah‡, Yasmin Akter‡, Mehraz Zaman Fardin‡
‡ Military Institute of Science and Technology, Dhaka, Bangladesh

Corresponding author: Md Imtiaj Hossain ( imtiajhossain854@gmail.com )

Academic editor: Maria Shchurovskaya
© 2023 Md Imtiaj Hossain, Abdus Sattar Mollah, Yasmin Akter, Mehraz Zaman Fardin.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation: Hossain MI, Mollah AS, Akter Y, Fardin MZ (2023) Neutronic calculations for the VVER-1000 MOX core computational benchmark using the OpenMC code. Nuclear Energy and Technology 9(4): 215-225. https://doi.org/10.3897/nucet.9.91090
Open Access

Abstract

The goal of this study is to perform neutronic calculations of the VVER-1000 MOX core computational benchmarks with an OpenMC code along with ENDF/B-VII.1 nuclear data library. The results of neutronic analysis using the OpenMC Monte Carlo code for the VVER-1000 MOX core, containing 30% mixed oxide fuel with low enriched uranium fuel, are presented in this study. As per the benchmark report, all six states are considered in the present study. The keff values, assembly average fission reaction rates, and pin-by-pin fission rates were calculated as per benchmark criteria. In addition, 2D thermal and fast neutron-flux distribution were also generated. The reactivity results and neutron flux distribution were compared with other results in which benchmark analysis was performed using the same core geometry and it showed great similarity with slight deviation. This shows that the modeling of the VVER-1000 MOX core was done successfully using OpenMC. Because OpenMC was successfully used for neutronics calculation of the VVER-1000 whole core, it may be mentioned here that OpenMC code can also be utilized for neutronics and other reactor core physics analyses of the VVER-1200 reactor which is to be commissioned in Bangladesh in the upcoming year.

Keywords

VVER-1000 MOX core, OpenMC, Low Enriched Uranium (LEU) Assembly, Mixed Oxide (MOX) fuel Assembly, Neutron Flux, Burnup

Introduction

The nuclear reactor, which is the center of a nuclear power plant, generates thermal power that is then converted to electric power for use in the economy by a variety of means. To avoid any unfortunate situation occurring, it is necessary to execute several core parameter calculations continuously. Calculations of multiplication factors, reactivity coefficients, fuel temperature (Doppler) and poison effect on reactivity, burnup, reactivity and isotopic concentration changes with burnup, fast and thermal neutron flux density, axial and radial power peaking of the core, fission rates distribution, power distribution of the core, etc. are among the crucial calculations (Lamarash 1988). To guarantee the integrity of the nuclear reactor core during operation, these calculations are carried out and evaluated regularly. The neutronic behavior of fuel assemblies and the core of a nuclear reactor with various combinations of fuel with different enrichments, moderator materials, and non-fuel structural components has been studied by using a suitable neutronic simulation code.

An OECD-NEA paper contains a comprehensive list of benchmarks that can be used to carry out this type of verification (Gomin et al. 2005). Some of the well-known Monte Carlo neutron transport programs, including MCNP (X-5 Monte Carlo Team 2008), SERPENT (Leppanen 2013), MONK (Richards et al. 2015), KENO (Petrie and Landers 1998), SuperMC (Wu et al. 2015) and TRIPOLI (Nimal and Vergnaud 1990), are currently gaining popularity as the greatest sources of information for computations involving reactor core physics. Unfortunately, a lot of these codes, which are frequently utilized as sources for neutronic calculations, are not easily available, and their dissemination is frequently restricted. However, some codes linked to reactor physics neutronic analysis, such as the OpenMC Monte Carlo code and the deterministic code DRAGON, are freely available and are increasingly used in code-to-code comparisons (Islam et al. 2022). An OpenMC code was used in our earlier research (Imtiaz et al. 2022; Nasim et al. 2022; Khan et al. 2022) to investigate “A VVER-1000 LEU and MOX Assembly Computational Benchmark” and predict the neutronic and burnup behavior at the lattice level.

Various benchmark problems may be used to extensively assess the core of a VVER reactor. For this investigation, a VVER-1000 full core containing 30% MOX fuel was used as a benchmark problem which was obtained from a benchmark analysis (Gomin et al. 2005) performed by a group of reactor physics experts at the Nuclear Energy Agency with MCNP-4c, MCU and RADAR codes. The benchmark problem specifies the different parameters to be calculated. Several other studies have been conducted by researchers utilizing different codes to accomplish the same computations, such as (Thilagam et al. 2009) who performed the VVER-1000 MOX core computational benchmark using indigenous codes EXCEL, TRIHEX-FA, and HEXPIN. OpenMC is a relatively new and freely accessible Monte Carlo particle transport code (Romano and Forget 2013) that allows users to find the criticality (keff) based on the average of three separate approaches such as track length, collision probability, and absorption. The ENDF/B-VII.1 data library, which contains all of the required cross-section data to perform a neutronic analysis, was employed in our investigation. Nuclear data for 423 nuclides are available in this collection (ENDF/B-VII.1 2012).

Model description

The designed model is a VVER-1000 reactor full core which contains 30% mixed oxide fuel alongside low-enriched uranium fuels. The modeling was done in OpenMC in a jupyter notebook with Python 3.9. The core includes both fresh and burned fuel from various burnups, which are arranged in a periphery-to-center pattern inside the core. Because fresh fuel can achieve a higher burnup compared to once or twice-burned fuel and produces a lot of power compared to other burned fuel, the neutron flux associated with this assembly is likewise a lot higher. The core has seven different types of fuel assemblies as mentioned in the benchmark problem, which are as follows:

  • Fresh UOX fuel assembly
  • 15 MWD/KgHM burned UOX fuel assembly
  • 32 MWD/KgHM burned UOX fuel assembly
  • 40 MWD/KgHM burned UOX fuel assembly
  • Fresh MOX fuel assembly
  • 17 MWD/KgHM burned MOX fuel assembly
  • 33 MWD/KgHM burned MOX fuel assembly

Each assembly contains 331 elementary cells of various types such as different enriched fuel, gadolinium pins, guide tubes, and central tubes and for state-6, some control rods are inserted in some specific assemblies, as mentioned later.

Various Assumptions were taken during the modeling process, they are:

  • Reflective boundary condition in the z-axis, transmission boundary condition between assembly boundaries, and vacuum boundary condition at the outermost surface of the core.
  • 8,000 batches with 150 inactive batches and 80,000 particles per batch were observed.

Several steps must be followed to model a full core. At the very beginning, each sort of elementary cell was designed. These elementary cells include fuel cells, fuel cells with gadolinium absorbers, guide tube cells, central tube cells, and absorber rod cells. Fig. 1 and Table 1 show a description of the cell geometry.

Table 1.
Download as
CSV
XLSX

Cell type geometry specification

Cells Name Cell Radius (cm)
Fuel cell R1 = 0.386
R2 = 0.455
Central tube cell R1 = 0.55
R2 = 0.63
Guide tube cell R1 = 0.55
R2 = 0.63
Guide tube with absorber rod R1 = 0.35
R2 = 0.41
R3 = 0.55
R4 = 0.63
Figure 1. 

Fuel and non-fuel cells.

Following the design of the elementary cells, seven different types of fuel assemblies were designed, and for state six, an additional five fuel assemblies containing absorber rod cells were created. Figs 2, 3 depict the two basic types of fuel assemblies (LEU and MOX) designed using OpenMC.

Figure 2. 

LEU Assembly.

Figure 3. 

MOX assembly.

Figure 4. 

Core geometry description.

After the successful modeling of the assemblies necessary for modeling the whole core, the core description from the benchmark report was followed and the full VVER-1000 core consisting of 163 fuel assemblies was designed. The 1/6th portion of the geometry description and the full core which was modeled using OpenMC is shown in Figs 4–6.

Figure 5. 

Geometry description for the reflector region.

Figure 6. 

VVER-1000 full core.

After successful modeling of the whole core, various parameters were calculated for analysis purposes using OpenMC. There are six states described in the benchmark report in which the calculation was performed. The operational states’ description is given in Table 2.

Table 2.
Download as
CSV
XLSX

The reactor states for both assemblies

States State name Fuel temperature (K) Non-fuel temperature (K) Reflector temperature (K) Moderator in the fuel assembly Water hole, water gap, and downcomer material Absorber rod
State 1 Working state 1027 575 560 M575B1.3 M60B1.3 -
State 2 State with constant temperature 575 575 560 M575B1.3 M560B1.3 -
State 3 Cold state with high boron content 300 300 300 M300B2.8 M300B2.8 -
State 4 Working state without boron 1027 575 560 M600B0 M560B0 -
State 5 State with constant temperature without boron 575 565 560 M560B0 M560B0 -
State 6 State with control rods inserted 565 565 560 M553B0 M553B0 Inserted

Here, in Table 2 columns 4 and 5, MxBy represents the Moderator at temperature x with y*1,000 ppm of boron contents.

Methodology

The models were represented in the OpenMC using python (python 3.7) code in jupyter notebook in the latest version of OpenMC (OpenMC 0.13.0). OpenMC has the feature to model hexagonal geometry which was used to design each type of assembly separately. Different materials in different regions inside the assemblies were defined using Boolean operation for modeling which is also known as constructive solid geometry. A hexagonal prism with an assembly pitch of 23.6 cm was used to bind the geometry giving it a hexagonal shape. For each pin cell, 1.275 cm of cell pitch was used. Two separate planes on the z-axis with reflecting boundary conditions, which is equivalent to the geometry being infinite on the z-axis, were defined. After completing the design of the seven types of assemblies, they were placed inside another hexagonal prism to produce the core. The core consists of a total of 163 fuel assemblies. To account for the thermal scattering at lower energies, S(α,β) table was provided. A total of six states were considered for the calculation of different parameters which are given in Table 2. State 6 is a special state where all of the control rods were inserted in their respective positions.

Results and discussion

Convergence test

Computing a value known as the Shannon entropy of the fission source distribution, Hsrc, has been done in research work to evaluate the convergence of the fission source distribution for the Monte Carlo method (Brown 2006; Ueki and Brown 2002). The behavior of the Shannon entropy curve in a Monte Carlo simulation is very important as the constant behavior of the entropy curve indicates the convergence of the simulation and also the number of inactive batches that should be ignored at the very start of each simulation process. The Shannon entropy of the discretized fission source distribution for a batch is given by (Brown 2006):

Hsrc=-∑J=1NPJ·ln2(PJ)

where Ns is the number of grid boxes in the superimposed mesh, and PJ = (number of source sites in J-th grid box)/(total number of source sites). Hsrc varies between 0 for a point distribution to ln2(Ns) for a uniform distribution.

The Shannon entropy curve is shown in Fig. 7a, according to which 100 inactive batches were decided for our simulation. Fig. 7b is a plot of effective multiplication factor vs generation or batches. This plot also showed a near-constant Shannon entropy per generation after a few inactive batches at the beginning. By observing this curve, the set of the number of batches and particles to be simulated in each batch was determined, which greatly increased the acceptability of our result.

Figure 7. 

Shannon Entropy and keff vs generation for state 1.

Figure 8. 

Shannon Entropy and keff vs generation for state 6.

Effective multiplication factor

The VVER-1000 whole core benchmark was first introduced and the effective multiplication factor was calculated for six different states from state-1 to state-6. The result obtained from OpenMC was compared with other results from benchmark reports such as MCNP-4c, MCU, which used the MCUDAT-2.1 data library as the basic data library, and MCNP5 (Lüle et al. 2015), which used ENDFB66 data library and our results showed very good similarity with the benchmark results through three codes.

Table 3 suggests that the obtained keff values agree well with other Monte Carlo codes, MCNP5, MCNP4C, MCU, and Benchmark Mean. The percent deviation in keff values between computed and benchmark Mean values for states S1–S5 vary from -0.375 percent to -0.507 and for S6 the variance is +0.535. The variations observed between the obtained values and other results from different codes are most likely due to the usage of different cross-section libraries.

Table 3.
Download as
CSV
XLSX

keff for states S1–S6

State OpenMC (OP) MCNP5 MCNP4C MCU BM* ∆K
(ENDF/B-VII.1) (ENDF/B-VI.6) (JEF2.2) (MCUDAT 2.1) (OP−BM)/(OP) × 100%
1 1.0337 ±0.006 1.03614 ±0.007 1.03770 ±0.007 1.03341 ±0.013 1.03769 -0.386
2 1.0465 ±0.006 1.04339 ±0.010 1.05132 ±0.010 1.04719 ±0.012 1.04989 -0.315
3 0.9294 ±0.009 0.93397 ±0.011 0.93416 ±0.011 0.93237 ±0.01 0.93286 -0.367
4 1.1310 ±0.004 1.13511 ±0.010 1.13871 ±0.010 1.1339 ± 0.012 1.13781 -0.432
5 1.1472 ±0.004 1.14333 ±0.010 1.15400 ±0.010 1.14932 ±0.012 1.15302 -0.507
6 1.0506 ±0.002 1.03914 ±0.010 1.04729 ±0.011 1.04267 ±0.009 1.04498 +0.535

*Benchmark mean value was obtained from MCNP-4C, RADAR, and MCU codes as per the benchmark report.

Assembly average fission reaction rates

The thermal output of the VVER-1000 core is roughly 3000 MW. The overall power is distributed across the 163 assemblies that constitute the core. Each assembly or pin within an assembly does not produce the same amount of power, and the power that it produces also changes depending on its enrichment and composition. The reactor power is proportional to fission reaction rates. Assembly average fission reaction rates for assemblies 1 through 28 were determined, along with their standard deviation, and compared to the findings from data from the literature review’s MCNP4C, MCU, and Radar (Gomin et al. 2005), HEXPIN (Thilagam et al. 2009), CNUREAS (Lüle et al. 2015), etc. Based on these data, it is clear that the result achieved using OpenMC is readily acceptable, as the maximum and minimum deviation value ranges from +7.9% to -9.7% for the six states below, given in Figs 9–14. It should be noted that not all of the findings shown here were obtained using the same data library. Also, the results presented here were generated by multiplying each data by a thousand. Each program made use of different data libraries, each with a different number of nuclides data. As a result, a little deviation is unavoidable. Due to the working principle and modeling approximation of different types of codes, even employing the exact same data library can result in a slight variation. As per the benchmark report, the pin-to-pin fission rates distribution of selected fuel assemblies 3, 21, and 27 for state-1 is shown in Figs 15–17. The deviation (%) between OpenMC and MCU is also shown in Figs 15–17 for comparison purposes. The OpenMC results are comparable with those of results from MCU data.

Figure 9. 

Assembly average reaction rates (×1000) and deviation (from MCNP & MCU) for state 1.

Figure 10. 

Assembly average reaction rates (x1000) and deviation (from MCNP & MCU) for state 2.

Figure 11. 

Assembly average reaction rates (x1000) and deviation (from MCNP & MCU) for state 3.

Figure 12. 

Assembly average reaction rates (x1000) and deviation (from MCNP & MCU) for state 4.

Figure 13. 

Assembly average reaction rates (x1000) and deviation (from MCNP & MCU) for state 5.

Figure 14. 

Assembly average reaction rates (x1000) and deviation (from MCNP & MCU) for state 6.

Figure 15. 

Pin-by-pin fission reaction rates for fuel assembly no 3 at state 1.

Figure 16. 

Pin-by-pin fission reaction rates for fuel assembly no 21 at state 1.

Figure 17. 

Pin-by-pin fission reaction rates for fuel assembly no 27 at state 1.

Neutron flux density spectrum

The neutron flux density spectrum was obtained from the flux tally via OpenMC. The flux spectrum is a 2D slice plot. Since the current version of OpenMC can’t generate an isometric plot of the neutron flux density, a 2D plot was generated for State-1, and State-6 only. The four slice plots of thermal fast-flux spectrum plots are shown in Figs 18–19. The VVER-1000 reactor is a PWR with a variety of fuel assemblies that have different multiplication characteristics owing to changes in enrichment and burnup. In PWRs, the out-in loading pattern is used, with the fresh fuel batch on the periphery of the core and the intermediate and high burnup batches in the center. During refueling, the highest burnt fuel assemblies are discharged, with fresh fuel loaded at the periphery and other batches inserted within. The right figure from Fig. 18 represents that the thermal flux is comparatively more towards the periphery of the core in comparison to other core positions. Consequently, the power production is more towards the periphery which aids in preventing power from peaking at the center of the core. The thermal absorption of neutrons by fissile nuclides increases as the thermal flux increases towards the periphery, increasing the fission reaction rate and hence enhancing the fast neutron flux, as illustrated in Fig. 18. The 2D plot of the neutron flux density spectrum was shown.

Figure 18. 

Thermal and fast neutron flux density (state-1).

Figure 19. 

Thermal and fast neutron flux density (state-6).

Here, the symmetric behavior of the core is seen. State-6 is a very special state, where control rods are inserted in some specific places inside the guide tubes in some selected assemblies. Hence, the thermal and the fast neutron flux density are less in the middle of the core due to control rod insertion, as illustrated by Fig. 19.

The benchmark report lacks a neutron energy spectrum for comparison. In their 2009 article, Thilagam et al. (2009) published 2D thermal and fast neutron spectrum by using HEXPIN code. The 2D neutron flux distribution calculated using OpenMC and those acquired using the HEXPIN code are equivalent. It is clear from this comparison that the OpenMC algorithm is appropriate for collecting the neutron energy spectrum for the whole VVER-1000 core.

Conclusions

The OpenMC code was used in this investigation to calculate the effective multiplication factor for states one through six, assembly average fission reaction rates, and pin-by-pin fission reaction rates. In addition, 2D thermal and fast neutron flux density distributions were calculated. Following that, the obtained results were contrasted with those from MCU and MCNP as well as other findings from the literature values. It was evident from the comparisons of keff values that OpenMC had been successfully implemented for the model mentioned in the OECD benchmark problem. The assembly average fission reaction rates also showed slight deviation from other assemblies, as shown in the result sections. The absence of the three mm water layer right outside the core could be one of the causes. A very substantial discrepancy was seen at interior assemblies compared to periphery assemblies for state six as well. As can be seen from the obtained neutron flux spectrum in various states as well as assembly average fission reaction rates, OpenMC demonstrated a very good capability in performing neutronic calculations for VVER geometry-based nuclear reactors, with the exception of the minor deviations caused by modeling errors and the use of a new data library.

Authorship contribution statement

Md. Imtiaj Hossain: Methodology, Data collection, Formal analysis, Writing – original draft. A. S. Mollah: Supervision, Conceptualization, Results interpretation, Writing – review & editing. Yasmin Akter: Resources, Data analysis, Writing. Mehraz Zaman Fardin: Resources, Literature review, Writing. All authors reviewed the results and approved the final version of the manuscript.

Acknowledgments

The authors really acknowledge the efforts of the Department of Nuclear Science and Engineering, Military Institute of Science and Technology, Dhaka, Bangladesh for their academic support. The authors thank the anonymous referees for their critical reading of the paper and for the improvements they suggested.

References

  • Abuqudaira TM, Stogov YV (2018) Neutronic calculations for the VVER-1000 LEU and MOX assembly computational benchmark using the GETERA code. Journal of Physics: Conference Series 1133: 012018. https://doi.org/10.1088/1742-6596/1133/1/012018
  • Brown D (2012) Release of the ENDF/B-VII.1 Evaluated Nuclear Data File. ANS Winter Meeting Town and Country Hotel and Resort, San Diego, CA November 11–15. https://doi.org/10.2172/1043398
  • Brown FB (2006) On the use of shannon entropy of the fission distribution for assessing convergence of Monte Carlo criticality calculations. PHYSOR-2006, ANS Topical Meeting on Reactor Physics, Organized and hosted by the Canadian Nuclear Society. Vancouver, BC, Canada, September 10–14.
  • Chadwick MB, Obložinský P, Herman M, Greene NM, McKnight RD, Smith DL, Young PG, MacFarlane RE, Hale GM, Frankle SC, Kahler AC, Kawano T, Little RC, Madland DG, Moller P, Mosteller RD, Page PR, Talou P, Trellue H, White MC, Wilson WB, Arcilla R, Dunford CL, Mughabghab SF, Pritychenko B, Rochman D, Sonzogni AA, Lubitz CR, Trumbull TH, Weinman JP, Brown DA, Cullen DE, Heinrichs DP, McNabb DP, Derrien H, Dunn ME, Larson NM, Leal LC, Carlson AD, Block RC, Briggs JB, Cheng ET, Huria HC, Zerkle ML, Kozier KS, Courcelle A, Pronyaev V, van der Marck SC (2006) ENDF/B –VII.0: Next Generation Evaluated Nuclear Data Library for Nuclear Science and Technology. Nuclear Data Sheets 107(12): 2931–3060. https://doi.org/10.1016/j.nds.2006.11.001
  • Hossain MdI, Akter Y, Fardin MZ, Mollah AS (2022) Neutronics and burnup analysis of VVER-1000 LEU and MOX assembly computational benchmark using OpenMC Code. Nuclear Energy and Technology 8(1): 1–11. https://doi.org/10.3897/nucet.8.78447.figure11
  • Islam A, Rahim TA, Mollah AS (2022) Inter-code comparison of computational VERA depletion benchmark using OpenMC. OpenMC-ONIX and DRAGON, Atom Indonesia (in press). https://doi.org/10.17146/aij.2022.1191
  • Khan GR, Nasim AS, Erfan KA, Mollah AS (2022) Verification of Monte Carlo Code OpenMC using VVER-1000 MOX Fuel Assembly against Criticality Benchmark Data. International Journal of Integrated Sciences & Technology 4S(2022): 7–13.
  • Khan SA, Jagannathan V, Kannan U, Mathur A (2016) Study of VVER1000 OECD LEU and MOX computational benchmark with VISWAM Code system. Nuclear Energy and Technology 2(4): 312–334. https://doi.org/10.1016/j.nucet.2016.11.008
  • Lamarsh J, Baratta A (2013) Introduction to Nuclear Engineering, 3rd Edition, Pearson, USA.
  • Leppanen J (2013) Serpent- A Continuous-energy Monte Carlo Reactor Physics Burnup Calculation Code. VTT Technical Research Centre of Finland.
  • Louis HK, Amin E (2016) The Effect of Burnup on Reactivity for VVER-1000 with MOXGD and UGD Fuel Assemblies Using MCNPX Code. Journal of Nuclear and Particle Physics 6(3): 61–71.
  • Mercatali L, Venturini A, Daeubler M, Sanchez VH (2015) SCALE and SERPENT solutions of the OECD VVER-1000 LEU and MOX burnup computational benchmark. Annals of Nuclear Energy 83: 328–341. https://doi.org/10.1016/j.anucene.2015.03.036
  • Mercatali L, Beydogan N, Sanchez-Espinoza VH (2021) Simulation of low-enriched uranium burnup in Russian VVER-1000 reactors with the Serpent Monte-Carlo code. Nuclear Engineering and Technology 53: 2830–2838. https://doi.org/10.1016/j.net.2021.03.014
  • Nasim AS, Khan GR, Erfan KA, Mollah AS (2022) Study on processing and validation of ENDF/B-VIII nuclear data library by criticality benchmark of PWR pin cells using NJOY21 and OpenMC. International Journal of Integrated Sciences & Technology 4S(2022): 1–6.
  • Petrie LM, Landers NF (1998) KENO V.a: An Improved Monte Carlo Criticality Program With Supergrouping, NUREG/CR-0200, Revision 6, Volume 2, Section F11, ORNL/NUREG/CSD-2/R6, September.
  • Thilagam L, Jagannathan V, Sunil sunny C, Subbaiah KV (2009) VVER-1000 MOX Core Computational Benchmark analysis using indigenous codes EXCEL, TRIHEX-FA and HEXPIN. Annals of Nuclear Energy 36(10): 1502–1515. https://doi.org/10.1016/j.anucene.2009.08.001
  • Thilagam L, Karthikeyan R, Jagannathan V, Subbaiah KV, Lee SM (2010) Inter-comparison of JEF-2.2 and JEFF-3.1 evaluated nuclear data through Monte Carlo analysis of VVER-1000 MOX Core Computational Benchmark. Annals of Nuclear Energy 37(2): 144–165. https://doi.org/10.1016/j.anucene.2009.11.009
  • Ueki T, Brown FB (2002) Stationarity Diagnostics Using Shannon Entropy in Monte Carlo Criticality Calculation I: F test. Transactions of the American Nuclear Society 87: 156.
  • X-5 Monte Carlo Team (2008) MCNP – A General MC N-Particle Transport Code, Version 5, Vol. 836.
login to comment