Coupling the GUARDYAN code to subchanﬂow

. GUARDYAN is a GPU-based dynamic Monte Carlo code developed at the Budapest University of Technology and Economics, Hungary. Dynamic Monte Carlo 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 with coupling to Thermal-Hydraulic 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 describes some convergence studies regarding reaching the initial equilibrium state. 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 ( L 2 ) of power release, fuel, and coolant temperatures for both static and dynamic convergence. Dynamic mode low-sample simulations provided surprisingly stable global L 2 measures, correct fuel temperatures, and power release, while coolant temperatures were oﬀ, without any indication of the incorrectness of the result. Static convergence showed an alternating, ﬂuctuating L 2 behavior that did not aﬀect the ﬁnal stable state.


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 [1][2][3]. 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. * e-mail: legrady@reak.bme.hu 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.
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 SALOME 1 platform. The two codes are compiled together, thus data is directly exchanged between SCF 1 https://www.salome-platform.org/ and GUARDYAN without the need for generating data files for data swapping.

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 (T f ), coolant temperature (T c ), and coolant density (ρ) data of the previous time step (i − 1) are used in the MC simulation for the calculation (F N ) of the pin-wise power release: The thermal-hydraulics module returns (F TH ) the new temperature and density data based on the power release: 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, sourceconvergence 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.
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 esti- mated. Thereafter, during the active cycles (GK A) nodewise 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 k eff .
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×3assembly 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.
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.

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 (2 20 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 (2 17 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 sub-sequent time steps: 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 L 2 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.
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.
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 L 2 evolution as seen in Figure 9.
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.

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 2 20 , 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.
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%.
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.

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.
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.

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 (L 2 ) of power release, fuel, and coolant temperatures for both  static and dynamic convergence. Dynamic mode lowsample simulations provided surprisingly stable global L 2 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 L 2 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.