MULTI-PHYSICS TRANSIENT SIMULATIONS WITH TRIPOLI-4

In this work, we present the first dynamic calculations performed with the Monte Carlo neutron transport code TRIPOLI-4 R © with thermal-hydraulics feedback. For this purpose, the Monte Carlo code was extended for multi-physics capabilities and coupled to the thermal-hydraulics subchannel code SUBCHANFLOW. As a test case for the verification of transient simulation capabilities, a 3x3-assembly mini-core benchmark based on the TMI-1 reactor is considered with a pin-by-pin description. Two reactivity excursion scenarios initiated by control-rod movement are simulated starting from a critical state and compared to analogous simulations performed using the Serpent 2 Monte-Carlo code. The time evolution of the neutron power, fuel temperature, coolant temperature and coolant density are analysed to assess the multi-physics capabilities of TRIPOLI-4. The stabilizing effects of thermal-hydraulics on the neutron power appear to be well taken into account. The computational requirements for massively parallel calculations are also discussed.


INTRODUCTION
Multi-physics simulations of nuclear reactor cores involve the coupling of neutron transport to other disciplines, such as thermal-hydraulics or thermo-mechanics. So far, such problems have mostly been analyzed using deterministic neutron-transport codes [1][2][3]. In recent years, intensive research efforts have focused on multi-physics calculations with Monte Carlo neutron transport codes in view of constructing high-fidelity simulation tools for the simulation of reactivity-induced transients in PWRs [4][5][6]. The goal is to provide reference solutions for neutron transport with multi-physics feedbacks.
In this context, that we have implemented an external coupling scheme between the Monte Carlo code TRIPOLI-4 R [7] and the thermal-hydraulics subchannel code SUBCHANFLOW [8]. The coupling was realized by developing a multi-physics interface for TRIPOLI-4, which makes it possible to perform criticality calculations with feedback [9] as well as time-dependent coupled calculations via the use of the kinetic transport methods that have been recently conceived in TRIPOLI-4 [10]. In this work, we present the TRIPOLI-4 multi-physics capabilities for the simulation of transients and we illustrate them with a few configurations that were benchmarked against similar Serpent 2 calculations [11].

TRIPOLI-4 MULTI-PHYSICS CAPABILITIES
In order to be able to perform multi-physics calculations with TRIPOLI-4 and thermal-hydraulics, we have developed an external program, called "supervisor", that drives a TRIPOLI-4 calculation step by step. The main role of the supervisor is to orchestrate data exchange between the coupled codes. The ICoCo API [12] from the SALOME platform serves as an interface to facilitate the coupling. The data exchanged between the two codes are represented in the form of fields on discretized meshes and are transferred in memory using the MEDCoupling library [13,14].
TRIPOLI-4 scores the power distribution on a fuel-centered mesh. For SUBCHANFLOW, two meshes are defined: a fuel-centered mesh receives the power distribution from TRIPOLI-4, and a coolant-centered mesh is used for the resolution of the conservation equations. Note that the TRIPOLI-4 and SUBCHANFLOW meshes must have the same bounding box, but they do not need to superimpose exactly or to be numbered in the same way, as long as spatial interpolation is performed between them. The power distribution scored by TRIPOLI-4 is transferred to SUB-CHANFLOW using the interpolation tools provided by the MEDCoupling library. The subchannel code computes the updated properties of the fuel (temperatures) and the moderator (temperatures and densities), which are transferred back to TRIPOLI-4. Again, the MEDCoupling interpolation tools are used to transfer the SUBCHANFLOW results to the corresponding cells of the TRIPOLI-4 mesh. The TRIPOLI-4 material temperatures and densities are then updated; note that this requires the TRIPOLI-4 geometry to be discretized in such a way that each mesh cell corresponds to an individual volume. The temperature dependence of the cross sections is taken into account with stochastic interpolation in the Monte-Carlo code.
Prior to a dynamic calculation, a criticality calculation with feedback must be performed to generate the initial reactor state. During this phase, relaxation is imposed on the power distribution in order to achieve a smooth convergence. When the convergence has been reached, it is possible to store the fission sources in a file: the position, energy, direction and weight of each particle are written and can be used as the initial state for another calculation; in the case of dynamic calculations, fission sources are converted into neutron and precursor sources for the dynamic phase by an additional power iteration [10]. For the time-dependent calculation, we formally apply the explicit Euler discretized scheme for the system of coupled neutron transport and thermal-hydraulics equations.
In general, the total CPU time is largely dominated by TRIPOLI-4 and parallelism is unavoidable. In a parallel calculation, we have chosen to have only one instance of SUBCHANFLOW for all the processors, so that the SUBCHANFLOW input is represented by the average power distribution calculated by all the parallel units. The rationale for this choice is that the thermal-hydraulics equations are non-linear and the propagation of the statistical fluctuations through the coupling scheme is not straightforward. Averaging the power distribution over all the parallel units minimizes the fluctuations on the input to the thermal-hydraulics solver and the resulting bias.
The TRIPOLI-4 parallelism scheme is the following. One parallel unit, the monitor, is in charge of orchestrating the other TRIPOLI-4 parallel units, of running the supervisor and SUBCHANFLOW. Another parallel unit (the scorer) is in charge of collecting the scores. The simulators are in charge of the TRIPOLI-4 simulation. When all simulators have completed a time step, the scorer collects the results and sends the averages to the monitor. The monitor sends the power distribution to SUBCHANFLOW, before finally launching the SUBCHANFLOW calculation. We emphasize the fact that the SUBCHANFLOW run can only start once all the simulators have completed the simulation of each time step.

TEST CASE
The selected configuration, illustrated in Fig. 1, is a 3x3 mini-core based on the TMI-1 reactor [15]. Each fuel assembly consists of 15x15 rods, made of 4.12% enriched UOX, and also contains four (Gd 2 O 3 +UO 2 ) burnable poison pins. An instrumentation tube is located in the center of each assembly. Assemblies are surrounded by a reflector made of borated water and stainless steel (not shown in Fig. 1). The central assembly contains 16 extra control rods composed of a Ag-In-Cd core and Inconel cladding. The active length of 353.06 cm is divided into 30 axial slices. The coolant inlet temperature is set to 565 K, the outlet pressure is 15.51 MPa and the inlet mass flow rate is 773.64 kg s −1 . The power of the system in stationary conditions is normalized to 140.94 MW. All the fuel and coolant compositions are individualized in the TRIPOLI-4 geometry, so as to be able to independently update their temperatures and densities. For the sake of simplifying the comparison with the published Serpent 2 calculations, we assume in this work that all the power is deposited in the fuel. In order to compute the fuel rod temperature, each axial slice of the fuel is divided into ten radial rings and the heat diffusion equation is solved with a finite-volume method. The fuel temperature returned by SCF is the the volume average over the radial nodes.
The calculations have been run on the CEA Cobalt supercomputer at the TGCC (Très Grand Centre de Calcul, Bruyères-le-Châtel, France) with 1000 parallel units during 24 hours.
Preliminary criticality calculations with feedback were performed to prepare the system on a critical state, with the power being normalized to 140.94 MW. Control rods are full inserted in the core.
With a boron concentration of 1305.5 ppm, the multiplication factor is k eff = 1.00018 ± 8 × 10 −5 , showing that the system is close to a critical state. The resulting thermal-hydraulics fields were stored at the end of the calculation: fuel temperatures, coolant temperatures and coolant densities, as well as the fission sources (position, energy, direction and weight of the fission neutrons for each one of the 1000 parallel units), to be used by the parallel units in charge of the dynamic calculations.

TRANSIENT SIMULATIONS
The purpose of the following transient simulations is to verify the TRIPOLI-4 multi-physics capabilities. We therefore analyse the impact of thermal-hydraulics feedback on power excursions, and we also present the comparison between the TRIPOLI-4/SUBCHANFLOW results and the Serpent 2/SUBCHANFLOW results published in Ref. 11. Starting from the stored source (fission sources and thermal-hydraulics fields) mentioned above, the time evolution of the total power was followed over 5 s in 50 regularly spaced intervals by increments of ∆t = 0.1 s. As explained, at the end of each time interval, the neutron power is averaged over all the parallel units and transferred to SUBCHANFLOW, which then runs and solves the thermal-hydraulics equations for the time step. The next time step begins with the updated temperatures and densities. Population-control and variance-reduction methods specific to kinetic simulations are used, namely combing, forced decay and branchless collisions [16].
In order to reduce correlations between the batches of the dynamic calculations, we have chosen to perform ten additional power iterations at the beginning of each batch, starting from the fission sources and thermal-hydraulics fields of the previous batch. SUBCHANFLOW is called with the new power distribution, and the compositions are updated. Thus, each batch of the dynamic simulations begins with sensibly different fission sources, temperature and density fields, and correlations between the corresponding dynamic observables are considerably reduced.
For the first scenario, the eight control rods are progressively extracted by 40 cm between t = 0.3 s and t = 1.3 s, which makes the system prompt supercritical with ∆ρ ≈ 1.3 $. The rod extraction is discretized in 10 steps. A second similar simulation was performed, with a less severe scenario, extracting the rods by 30 cm, which corresponds to a reactivity insertion of about ∆ρ ≈ 0.5 $. Each simulator completed seven full batches for the 30 cm scenario and four full batches for the 40 cm scenario. The time evolution of the total power for both scenarios is shown in Fig. 2. The impact of the thermal-hydraulics feedback is clearly visible. First, the power is stable around 140.94 MW between t = 0 s and 0.3 s, as expected. Then, the power increases up to 220 MW for the 30 cm scenario and up to 1100 MW for the 40 cm scenario. Finally, the thermal-hydraulics feedback mechanisms absorb most of the reactivity and drive the system towards an equilibrium configuration around 150 MW. The power peak is much sharper for the 40 cm scenario. The results obtained with Serpent 2/SUBCHANFLOW are also presented (the standard error was not available and is not shown). We find a very good agreement between the TRIPOLI-4/SUBCHANFLOW and the Serpent 2/SUBCHANFLOW results.
The time evolution of the fuel temperature averaged over all the rods of the mini-core is also presented in Fig. 2. As the reactivity increases, the average fuel temperature increases, then it remains stable after about 2 s. At the same time, the coolant temperature increases and the coolant density decreases as shown by Fig. 3. The final equilibrium state is characterized by higher temperatures  for the 40 cm scenario than for the 30 cm. There are some differences between the two coupling schemes, but the agreement is good overall; the maximum difference for the average fuel rod temperature is 10 K, in the peak power; the difference is below 1 K for the average coolant temperature and below 1 kg/m 3 for the average coolant density. It is difficult to ascertain whether these discrepancies may be due to statistical fluctuations, since the temperature/density values at different times are strongly correlated.
The increase with time in the fuel temperature for the 40 cm scenario is shown in Fig. 4 at a pin-cell level for two slices. Slice 3 is located close to the bottom of the mini-core, where the neutron flux actually rises because of the rod extraction; slice 10 is located above the rod extraction. For this reason, the increase is much larger for slice 3: the temperature increases up to 900 K at the center of the mini-core. For slice 10, the temperature mainly increases in the surrounding assemblies, since at this height the center of the mini-core contains the absorbant part of the control rods. inactive for almost two hours because of one single simulator. The large fluctuations in calculation time probably reflect fluctuations in neutron population, which grow with the size of the neutron population itself and are amplified by the branching nature of the fission process. In order to reduce the total waiting time and increase the simulation efficiency, we have tried to enforce population control on a tight grid during this time step (every 0.01 s), but no significant improvement was observed. A tighter time grid may be required for this purpose.
Our coupling scheme exacerbates the negative impact of the fluctuations, because the thermalhydraulics solver is called only when all the simulators have completed their histories. One way to improve the efficiency of the simulation would be to split the calculation into smaller, truly independent replicas, with independent couplings to thermal-hydraulics. When all simulators of the same packet have completed the time step, SUBCHANFLOW could be run without waiting for the other packets, and at least these simulators could start the next time step. Therefore, we would still be waiting for some packets, but fewer simulators would stay inactive during this time.

CONCLUSIONS
We have performed the first dynamic simulations with the coupling scheme between the Monte Carlo neutron transport code TRIPOLI-4 and the thermal-hydraulics subchannel code SUBCHAN-FLOW. For this purpose, we have considered a 3x3 mini-core based on the TMI-1 reactor. The initial reactor state was calculated once with a criticality calculation including feedback. We reused it for the two dynamic calculations.
We have studied two reactivity-insertion scenarios, with the control rods being extracted by respectively 30 cm and 40 cm. For both scenarios, the feedback mechanisms make the power decrease and ultimately reach a new equilibrium state. We compared our results to the published Serpent 2/SUBCHANFLOW results and generally found good agreement.
The simulation of the peak power is very challenging because of the large variations on the population size. The 40 cm rod extraction especially induces a very large reactivity insertion (the system becomes prompt supercritical), which leads to large fluctuations on the population size and thus in CPU time among the parallel units. In particular, one unit was about 50 times slower than the median. In order to mitigate the impact of fluctuations in CPU time among parallel units, the simulation could be split up in independent replicas. However, this solution is in tension with the need to minimize the statistical fluctuations resulting from the Monte Carlo simulation.