Issue
EPJ Nuclear Sci. Technol.
Volume 9, 2023
Euratom Research and Training in 2022: the Awards collection
Article Number 11
Number of page(s) 9
Section Part 1: Safety research and training of reactor systems
DOI https://doi.org/10.1051/epjn/2023002
Published online 10 February 2023

© D. Legrady et al., Published by EDP Sciences, 2023

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

GUARDYAN (GpU Assisted Reactor Dynamic ANalysis) is a GPU-based dynamic Monte Carlo code developed at the Budapest University of Technology and Economics, Hungary [13]. As an emerging calculation technology, Dynamic Monte Carlo (DMC) [4] computes the neutron population evolution by calculating the direct time dependence of the neutron histories in multiplying systems. Some well-established Monte Carlo codes have DMC versions [5, 6] with coupling to Thermal-Hydraulic (TH) solvers. GUARDYAN has the computational advantage of applying GPUs, thus calculation burden can be carried by commonly available hardware and is capable of handling power plant size systems for kinetics problems. GUARDYAN has been recently coupled to the subchannel thermal-hydraulics code SUBCHANFLOW in order to carry out dynamic calculations with TH feedback. This paper is the first to report on this coupling, further, we describe some convergence studies regarding reaching the initial equilibrium state. Comparison is given of reaching the steady state using dynamic or static (k-mode) calculations.

2. The main features of GUARDYAN

GUARDYAN is specifically written to perform dynamic Monte Carlo calculations using GPUs. Geometry setup uses constructive solid geometry, where cells with homogenous material distributions are bounded by second-order surfaces. Within a cell, temperatures may continuously change thanks to the Woodcock (delta-) tracking algorithm for free path sampling. The code is written for the Linux environment in C++ and CUDA. Code verification involved rigorous comparison to well-established Monte Carlo and deterministic codes. Dynamic validation involved the reproduction of rapid reactivity insertions and withdrawals at the BME Training Reactor, while the current highlight being the reproduction of ex-core detector signals of a rod-drop experiment of a VVER-440/V213 nuclear power plant of unit 4 at the Paks NPP at the beginning of campaign number 32 with realistic burnup data considered [2], see Figure 1

thumbnail Fig. 1.

Power distribution in the core of Paks NPP Unit 4 campaign 32-rod drop experiment beginning state (left), halfway (middle), and end state (right).

The dynamic versions of the CPU-based codes TRIPOLI-4 [5] and Serpent-2 [6] are also coupled to SUBCHANFLOW and transient simulations for realistic small reactors with comparison to measurements or benchmark scenarios for larger, NPP-sized cores have been reported. Performance comparison is challenging due to different computing architectures.

3. Coupling to thermal-hydraulics

GUARDYAN has been coupled to the thermal-hydraulic subchannel code SUBCHANFLOW (SCF) [7]. Based on the COBRA-TF set of codes [8, 9] and COBRA-EN, the code SCF was developed at the Karlsruhe Institute of Technology. The development included a modernized software implementation in Fortran 95, substantial improvement in user-friendliness, and more robust error handling but foremost a refined two-phase flow handling. Empirical correlations of pressure drop, void generation, heat transfer, etc. are incorporated into the basic conservation equations of energy, momentum, and mass and implemented in SCF making it capable of accurate modeling of light water reactors. Water properties are updated to comply with IAPWS-97 (International Association for the Properties of Water and Steam). Other liquid and gas coolants can also be specified such as liquid sodium or lead, or gases like Helium or air. Hexagonal and rectangular lattices can be specified for regular plant geometries.

In the course of the McSafe [10] project SCF was coupled to the DMC version of TRIPOLI and Serpent-2. The SCF features developed to ease coupling with these external codes provided an almost readily applicable way for coupling to GUARDYAN. For a flexible solution only subroutine-level modifications (like ISO_C_BINDING) were used and we did not utilize the SALOME1 platform. The two codes are compiled together, thus data is directly exchanged between SCF and GUARDYAN without the need for generating data files for data swapping.

3.1. Coupling scheme

The thermal-hydraulics solver expects to receive input data on the power release per pin axial node resulting from the Monte Carlo simulation. GUARDYAN considers interactions for fission power release other sources of energy releases were not considered. The probability of fission interaction probability was favorably biased towards more frequent sampling with proper weight compensation to keep the expected value intact.

The coupling of thermal-hydraulics and neutronics follows the Operator Splitting Semi-Implicit (OSSI) method. Essentially, fuel temperature (Tf), coolant temperature (Tc), and coolant density (ρ) data of the previous time step (i − 1) are used in the MC simulation for the calculation (FN) of the pin-wise power release:

P k = F N ( T f i 1 , T c i 1 , ρ i 1 ) . $$ \begin{aligned} P^{k}=F_{N} \left({T_{f}^{i-1}, T_{c}^{i-1} ,\rho ^{i-1}}\right). \end{aligned} $$(1)

The thermal-hydraulics module returns (F TH) the new temperature and density data based on the power release:

[ T f i , T c i , ρ i ] = F TH ( P i ) . $$ \begin{aligned} \left[{T_{f}^{i},T_{c}^{i},\rho ^{i}}\right]=F_{\rm TH} \left({P^{i}}\right). \end{aligned} $$(2)

The relaxation parameter was not introduced to average previous time step results to obtain the new one, though this may have improved the convergence rate, especially with the GK mode. The thermal-hydraulic calculations are executed separately for each subchannel. Fuel and coolant temperatures are used for calculating temperature-dependent cross-sections.

GUARDYAN coupled dynamic calculations are executed in three distinct phases, see Figure 2, GD stands for GUARDYAN Dynamic mode. In the first, source-convergence phase the equilibrium dynamic flux shape is calculated. Starting from uniform distribution dynamic transport is pursued with instant decay of precursors for a millisecond (in thermal systems). At the end of the time step the node-wise power distribution is normalized to the nominal target value and the thermal-hydraulics module is called. A fully converging static TH calculation is carried out until the TH equilibrium is reached. After these temperature and density values are returned, the loop of these two steps repeats until convergence is reached. Criticality is set by the negative feedback of the system controlled by the moderator Boron concentration.

thumbnail Fig. 2.

Coupling scheme of GUARDYAN Dynamic (GD) mode to thermal-hydraulics (TH) solvers.

In the second step, the delayed neutron precursors are generated in a pure neutronics run, during which neutron histories are followed and at fission sites precursor samples are banked until the designated memory space is filled. Temperatures and densities are kept constant, node-wise power release is not estimated and the thermal-hydraulic module is not called.

In the third step, the transient calculation starts with precursor decaying with properly sampled delayed neutron time dependence, and at each coupling time step the TH module is called and temperature and density parameters are exchanged. Node-wise power release is estimated but it is not normalized anymore to the nominal power. The TH module now only calculates the length of the coupling time step in contrast to the first step where static TH results are obtained.

General purpose MC codes usually operate in k-eigenvalue mode. For aiding comparison with these codes GUARDYAN can also calculate static k eff. In this GUARDYAN k-mode (GK) calculation neutrons are followed from their birth at a fission site until the next induced fission or a termination (escape, failed Russian Roulette) event. For the coupling scheme, see Figure 3 inactive cycles are run first (GK I) to reach neutron source convergence, now the power release is not estimated. Thereafter, during the active cycles (GK A) node-wise power release is estimated. After normalization to nominal power, data is given to the TH solver. With the new temperature and density data the next inactive cycle starts. This loop repeats until stopped, the last step offers the estimation for keff.

thumbnail Fig. 3.

Coupling scheme of static calculation (GUARDYAN, GK) and thermal-hydraulics (TH).

4. Convergence studies on the TMI Minicore benchmark

4.1. The TMI Minicore model

Based on the Three Mile Island TMI-1 reactor a 3 × 3-assembly simplified core model is presented in the NEA benchmark set [11] intended for benchmarking kinetic calculations. The 9 fuel assemblies contain 15 × 15 fuel pins made of 4.12% enriched UOX fuel. In each assembly, four Gd-containing burnable poison pins are located. The fuel length is 353 cm, and the assembly is 21.64 cm wide. At the center of the assemblies, void tubes can be found. The core is surrounded by a 21.64 cm thick layer of a mixture of steel and borated water. The central assembly has the specialty of housing 16 control rods made of Ag–In–Cd. At Hot Zero Power control rods are fully inserted, see Figure 4.

thumbnail Fig. 4.

TMI Minicore cross-section material models.

The SCF model of the TMI Minicore places a subchannel between each set of four fuel pins, filled with water. The inlet temperature is kept constant at 565 K, outlet pressure is 15.5 MPa, and mass flow rate is set to 773.64 kg/s. The power is defined as 140.94 MW. The core axial length was divided into 30 axial nodes, within a pin further radial divisions were not created.

4.2. Convergence to steady state in dynamic mode

To scout, the convergence properties of the dynamic coupling scheme, calculations for the TMI Minicore benchmark were carried out with two different parameter settings. Run #1 used more samples (220 neutrons, about 1 000 000) but with less frequent coupling with an interval size of 5 ms. Run #1 was based on the settings found in référence [12]. Run #2 employed a significantly smaller sample size (217 neutrons, equaling about 100 000) with frequent coupling stepping carried out at each 100 μs. With Run #2 the intention was to create a statistically unstable, noisy solution that is coupled very frequently to the TH solver inflating the chance that the coupled calculation becomes statistically unstable and observably fails.

For assessing the convergence state a measure has been defined as the L 2 distance of the distributions in two subsequent time steps:

L 2 ( X l X l 1 ) = i [ X i l X i l 1 ] 2 N , $$ \begin{aligned} L_{2} \left( {X^{l}-X^{l-1} } \right)=\frac{\sqrt{\mathop \sum \nolimits _{i} \left[ {X_{i}^{l} -X_{i}^{l-1} } \right]^{2}} }{\sqrt{N} }, \end{aligned} $$(3)

where l is the time step number and X i is a certain quantity in the space segment i, with i having the maximum of N. X i may stand for node-wise power release, temperature, or density. When L 2 becomes stationary as a function of l, the convergence to a steady state may have been reached. As a result of a stochastic calculation, a fluctuation around the stationary value appears, reaching a steady state is meant within the stochastic error.

The equilibrium state of the neutron population in GD mode can be attributed to an equilibrium state of the neutron flux resulting from an initial source convergence process similar to power iteration for equilibrium fission source distribution. This would reflect also in a steady power evolution as a function of time steps.

Figure 5 shows the initial source convergence. For Run #2, power stabilizes after a few coupling time steps without a visible long-term trend. Fuel temperature and similarly coolant temperature achieve a stable value after 20–30 time steps. Larger sample numbers (and consequently lower relative variances in the power estimates) in Run #1 allow for faster convergence to the initial ground state in terms of the number of time steps, see also Figure 6, showing the data of Figure 5 in an L2 sense. Average coolant temperature and power converge to the same value, however, converged fuel temperatures differ. This indicates that stochastic perturbation of the released power has a nonlinear effect on the equilibrium state, triggering an ongoing detailed analysis of TH-MC coupling, see référence [13]. Lower sample numbers of Run #2 result in stable-looking estimates, devoid of large fluctuations of the averages, with little indication that the convergence may be false.

thumbnail Fig. 5.

Source convergence in dynamic mode, power is normalized to a single starting neutron.

thumbnail Fig. 6.

L2 norm evolution of the GD source convergence for Run #1 (red) and Run #2 (blue).

A map of the pin-wise fuel temperature distribution shows a visible difference in the estimation stability of higher and lower sample numbers, see Figure 7.

thumbnail Fig. 7.

Pin-wise fuel temperature distribution at the converged state, for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step.

The apparently noisy node-wise temperature distribution in Figure 7 hardly suggests that a converged equilibrium state may be reached, or a stable longer-term calculation is feasible. Low sample numbers and short time steps warrant very noisy node-wise power release see Figure 8 that could average in longer time spans, and through that, may provide relaxed, smoother input for the TH calculation with larger coupling time steps. That not being the case, surprisingly, the TH solver tames the noisy input into a globally smooth L2 evolution as seen in Figure 9

thumbnail Fig. 8.

Node-wise power estimates for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step.

thumbnail Fig. 9.

Dynamic simulation average power, fuel, and coolant temperatures, for Run #1 (red) and Run #2 (blue).

Looking at the time evolution of the integral quantities (see Fig. 9) of the time-dependent system after the source convergence and delayed neutron sampling have already been done, power seems to be stable within the statistical fluctuations: smoother for Run #1 (red curve), more undulating for Run #2.

Temperature evolution is very smooth. Coolant temperature seems to be closely matching, but fuel temperatures systematically differ by about 20 degrees. It seems that the statistical poorness of the sampling does not directly translate into the statistical ill-being of the TH integral quantities but may result in an unconverged, false solution.

4.3. Convergence to steady state in k-mode

An equilibrium critical state may also be reached by GK calculation. Starting sample number was chosen to be 220, at each coupling step the neutronics source convergence was run for 100 inactive and 400 active cycles. During the inactive cycles, the estimates of the power release were disregarded but the updated temperature values took effect. Therefore, the new temperature and neutron fields did not go in entwined progress like in the GD mode. Figure 10 shows the alternating convergence to an equilibrium end state. Different random number generator seeds were chosen to show that the large stepwise fluctuations are not caused by randomness and that the converged state is common. The inactive cycles disconnect the previous neutronic state with the next coupling state, hence the alternating deviation from the end state appears. The converged temperature values agree with those found in the GD mode.

thumbnail Fig. 10.

Convergence in static GK mode, evolution of the criticality (k eff), fuel, and coolant temperatures with the coupling steps.

Node-wise relative error figures per GK step stay below 2% in the middle axial section of the reactor, see Figure 11. Even for the topmost nodes, where a high leakage rate decreases estimation efficiency the nodewise relative standard deviation stays below 10%.

thumbnail Fig. 11.

Relative errors of the nodewise power estimates in percent for the lowest (left), central (middle), and topmost (right) nodes.

Relevance of these error figures may count if some local anomaly is sought, but stakes are higher if fitness for TH coupling is to be judged. In this respect, it may be more telling to investigate how the error in the power estimation translates into the error in temperatures. A local statistical power anomaly in a millisecond-scale time interval is averaged out by the much slower thermal-hydraulic processes. Moreover, the directed stream of the coolant flow may transport this statistical error to a different position in the core.

In Figure 12 the fuel temperature standard deviation is plotted against the standard deviation of the power estimates for each pin node. Correlation is not apparent, but fuel temperature errors are generally smaller than the causing statistical error in the power. As power estimation errors increase the amount of statistical inaccuracy transferred seems to decrease.

thumbnail Fig. 12.

Correlation of standard deviation of power and fuel temperature.

4.4. Comparison of k-code and dynamic convergence

Apart from the globally characteristic number as average fuel and coolant temperatures or average flux, node-wise data were compared, see Figure 13. Dynamic and static neutron fluxes should coincide for a critical reactor, and converged TH-states should not differ either.

thumbnail Fig. 13.

Comparison of dynamic and static convergence results: relative difference of power output (left), relative difference of coolant temperatures (middle), and relative difference of fuel temperatures (right).

Sub-percent pin-wise differences in all three regarded comparisons indicate that the converged end-states are the same. Apparently random deviations can be observed in the power and the fuel temperature, in the range of a few percent that corresponds to errors in the GK mode. The coolant temperature errors seem to show some weak patterns but the amplitude of the deviations stays way below a percent.

5. Conclusions

Coupling of the direct Monte Carlo code GUARDYAN to the subchannel code SUBCHANFLOW has been achieved and convergence to a steady state was investigated. A literature-suggested set of TH input settings and high sample numbers resulted in very low statistical errors of the power estimates and stable global measures (L2) of power release, fuel, and coolant temperatures for both static and dynamic convergence. Dynamic mode low-sample simulations provided surprisingly stable global L2 measures, correct fuel temperatures, and power release, while coolant temperatures were off, without any indication of the incorrectness of the result. Static convergence showed an alternating, fluctuating L2 behavior that did not affect the final stable state.

Further research should target the detailed analysis of the stability, robustness, and unbiasedness of coupled dynamic Monte Carlo and thermal hydraulics.


Conflict of interests

The authors declare that they have no competing interests to report.

Acknowledgments

The authors acknowledge the support and computational resources of the Wigner Scientific Computational Laboratory (WSCLAB) (the former Wigner GPU Laboratory).

Funding

The research reported in this paper is part of project no. BME-EGA-02, implemented with the support provided by the Ministry of Innovation and Technology of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021 funding scheme.

Data availability statement

Data associated with this article cannot be disclosed due to legal reasons.

Author contribution statement

David Legrady: manuscript preparation, research group leader, conceptualization.

Marton Pukler: algorithm testing, modeling.

Balint Panka: reactor geometry modeling.

Gabor Tolnai: simulations, modeling, programming, data analysis, graphs.

References

  1. B. Molnar, G. Tolnai, D. Legrady, A GPU based direct Monte Carlo simulation of time dependence in nuclear reactors, Ann. Nucl. Energy 132, 46–63 (2019) [CrossRef] [Google Scholar]
  2. D. Legrady, G. Tolnai, T. Hajas, E. Pazman, T. Parko, I. Pos, Full Core Pin-Level VVER-440 Simulation of a Rod Drop Experiment with the GPU-Based Monte Carlo Code GUARDYAN, Energies 15, 2712 (2022) [CrossRef] [Google Scholar]
  3. Official webpage of GUARDYAN, http://awing.reak.bme.hu/GUARDYAN [Google Scholar]
  4. B.L. Sjenitzer, J.E. Hoogenboom, Dynamic Monte Carlo method for nuclear reactor kinetics calculations, Nucl. Sci. Eng. 175, 94–107 (2013) [CrossRef] [Google Scholar]
  5. M. Faucher, D. Mancusi, A. Zoia, New kinetic simulation capabilities for Tripoli-4: Methods and applications, Ann. Nucl. Energy 120, 74–88 (2018) [CrossRef] [Google Scholar]
  6. J. Leppanen, Development of a dynamic simulation mode in Serpent 2 Monte Carlo code, in Proceedings of the International 426, Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2013), LaGrange Park, Illinois (2013), p. 427 [Google Scholar]
  7. U. Imke, V.H. Sanchez, Validation of the subchannel code SUBCHANFLOW using the NUPEC PWR Tests (PSBT), Sci. Technol. Nucl. Instal. 12, 465059 (2012) [Google Scholar]
  8. M.J. Thurgood et al., COBRA/TRAC – A thermal-hydraulic code for transient analysis of nuclear reactor vessel and primary coolant systems, in NUREG/CR-3046, 1983 [Google Scholar]
  9. D. Basile, R. Chierici, M. Beghi, E. Salina, E. Brega, COBRA-EN, an updated version of the COBRA-3C/MIT Code for thermal-hydraulic transient analysis of light water reactor fuel assemblies and cores, Report 1010/1 (revised 1.9.99), ENEL-CRTN Compartimento di Milano (2001) [Google Scholar]
  10. High-Performance Monte-Carlo Methods for SAFEty Demostration – From proof of concept to realistic Safety Analysis and Industry Applications, Grant Agreement ID 755097. https://cordis.europa.eu/project/rcn/211032/factsheet/en [Google Scholar]
  11. R.N. Bratton, M. Avramova, K. Ivanov, OECD/NEA benchmark for Uncertainty Analysis in Modeling (UAM) for LWRS – summary and discussion of neutronics cases (Phase I), Nucl. Eng. Technol. 46, 313–342 (2014) [CrossRef] [Google Scholar]
  12. M. Faucher, Coupling between Monte Carlo neutron transport and thermal-hydraulics for the simulation of transients due to reactivity insertions, PhD thesis, University of Paris-Saclay, 2019 [Google Scholar]
  13. T. Hajas, G. Tolnai, D. Légrády, Variance analysis of the coupling of thermal-hydraulics and point kinetics with stochastic noise term modeling dynamic Monte Carlo behavior for later use in GUARDYAN, in International Conference on Physics of Reactors 2022, 2022 May 15–20, Pittsburgh [Google Scholar]

Cite this article as: David Legrady, Marton Pukler, Balint Panka, Gabor Tolnai. Coupling the GUARDYAN code to subchanflow, EPJ Nuclear Sci. Technol. 9, 11 (2023)

All Figures

thumbnail Fig. 1.

Power distribution in the core of Paks NPP Unit 4 campaign 32-rod drop experiment beginning state (left), halfway (middle), and end state (right).

In the text
thumbnail Fig. 2.

Coupling scheme of GUARDYAN Dynamic (GD) mode to thermal-hydraulics (TH) solvers.

In the text
thumbnail Fig. 3.

Coupling scheme of static calculation (GUARDYAN, GK) and thermal-hydraulics (TH).

In the text
thumbnail Fig. 4.

TMI Minicore cross-section material models.

In the text
thumbnail Fig. 5.

Source convergence in dynamic mode, power is normalized to a single starting neutron.

In the text
thumbnail Fig. 6.

L2 norm evolution of the GD source convergence for Run #1 (red) and Run #2 (blue).

In the text
thumbnail Fig. 7.

Pin-wise fuel temperature distribution at the converged state, for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step.

In the text
thumbnail Fig. 8.

Node-wise power estimates for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step.

In the text
thumbnail Fig. 9.

Dynamic simulation average power, fuel, and coolant temperatures, for Run #1 (red) and Run #2 (blue).

In the text
thumbnail Fig. 10.

Convergence in static GK mode, evolution of the criticality (k eff), fuel, and coolant temperatures with the coupling steps.

In the text
thumbnail Fig. 11.

Relative errors of the nodewise power estimates in percent for the lowest (left), central (middle), and topmost (right) nodes.

In the text
thumbnail Fig. 12.

Correlation of standard deviation of power and fuel temperature.

In the text
thumbnail Fig. 13.

Comparison of dynamic and static convergence results: relative difference of power output (left), relative difference of coolant temperatures (middle), and relative difference of fuel temperatures (right).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.