APPLICATION OF THE SPH EQUIVALENCE TECHNIQUE TO THE CABRI FULL CORE IN NON-FUNDAMENTAL MODE

We developed a SPH equivalence technique in non-fundamental mode condition between a CABRI full-core model solved with the method of characteristics (MOC) in 2D and a simplified full-core model solved with the simplified P3 (SP3) method, linear anisotropic sources and discretized with Raviart-Thomas finite elements over a pure Cartesian mesh. The MOC and SP3 calculations are performed with DRAGON5 and DONJON5 codes, respectively. A three-parameter database is generated by DRAGON5 and is interpolated in DONJON5 as a function of the core condition. An objective function is set as the root mean square (RMS) error (MOC-SP3 discrepancy) on absorption distribution and leakage rates defined over the macro-geometry in DONJON5. Our algorithm is a quasiNewtonian gradient search based on the Limited memory Broyden-Fletcher-GoldfarbShanno (LBFGS) method. Numerical results are presented with Hafnium bars withdrawn or inserted.


INTRODUCTION
The objective of the CABRI International Programe (CIP) conducted by the Institute for Radioprotection and Nuclear Safety (IRSN) and operated by the Commissariatà l'Énergie Atomique et auxÉnergies Alternatives (CEA) is to study the behaviour of nuclear fuel rods and their cladding during a reactivity injection accident (RIA) in pressurized water reactors (PWR). Such an accident would result in a rapid, sudden and local increase in the neutron flux, which would induce an increase in nuclear power due to fission. The experimental program consists in exposing a section of an irradiated fuel rod under the thermohydraulic and neutronic power conditions encountered during a reactivity injection accident. Power transients are driven by depressurization of helium-3 (neutron-absorbing gas) within the CABRI reactor core, at the centre of which the rod to be tested is positioned.
The important neutron leakage due to the size and the complex geometry of the CABRI reactor (hodoscope instrumentation and water loop in the core) make the classical deterministic calculation scheme inadequate. Indeed, the usual 2-step procedure based on a transport calculation at the assembly level using fundamental mode assumption and then a full-core diffusion calculation fails to describe accurately the neutronic behaviour of the CABRI core.
To overcome this problem, we developed a modified 2-step procedure where a detailed 2D transport calculation of the CABRI full-core is performed in non-fundamental mode condition to produce homogenized and condensed cross sections to be used in a 3D full-core simplified P N (SPN) calculation. An intermediate macro-calculation is also defined to perform the SPH equivalence process. Figure 1 depicts the MOC geometry and the macro-geometry, respectively. The proposed procedure is implemented in the lattice DRAGON5 and full-core DONJON5 codes. [ Figure 1: Full-core SPH equivalence over the CABRI core.
The three types of flux calculations used in the modified 2-step procedure are: • The 2D calculation is performed in DRAGON5 using 26 energy groups and 18,000 nodes with flat source approximation and P 1 source anisotropy. A buckling approximation is used to represent axial leakage. DRAGON5 calculations are based on the method of characteristics (MOC) as implemented in modules SALT: and MCCGT:. Specific options are: -SHEM295 cross section library based on Jeff 3.3 used for the first level and condensed at 26 energy groups for the second level -Resonance self shielding based on subgroup method with CALENDF tables (PT option) as implemented in module USS: -The MOC calculation uses TISO 6 15.0 -The fuel pins are split in 50%/vol, 30%/vol, 15%/vol and 5%/vol subvolumes -Axial leakage for all loadings was simulated using a fictious absorption rate L g = d g B 2 φ g based on a multigroup leakage coefficient d g obtained using an homogeneous P 1 leakage model with an axial buckling fixed at B 2 = 0.001536 cm −2 , corresponding to the actual core height.
• The 2D SP3 macro-calculation is performed in DONJON5 using 8 energy groups and 361 Raviart-Thomas finite elements over a pure Cartesian mesh. A buckling approximation is used to represent axial leakage. DONJON5 calculations are based on SPN theory, as implemented in module TRIVAT:. Specific options are: -Bi-cubic Raviart-Thomas approximation with Gauss-Legendre integration -P 1 Scattering anisotropy.
• The 3D SP3 space-time calculation is intended to be used in multi-physics simulations. These calculations are beyond the scope of this paper.
The CABRI power transients are generated by the depressurization of pressurized 3 He contained in the transient rods located at the 4 corners of the core (see Fig. 1). The transient rods system is composed of 96 tubes spread around 4 fuel assemblies linked to a system of control and fast opening valves constituting a low and a high flow-rate channel. The initial transient state is obtained by the withdrawal of the Hafnium bars to the critical position. The Doppler coefficient varies with fuel temperature and core poisoning due to 3 He and to Hafnium bars. We have chosen to represent the complete core, including its graphite reflector, with a set of 117 homogenized cross section values, condensed over 8 energy groups. These cross sections are interpolated from a multi-parameter database (known as the DRAGON5/DONJON5 multicompo) as a function of three parameters: • Here, we are neglecting feedback effects due to coolant temperature and density, as these effects are delayed with respect to the increased of fuel temperature. A complete multi-physics simulation of a CABRI pulse is beyond the scope of this paper. A complete grid of parameters is used, so that 50 MOC calculations are required to construct the multicompo. This approach is expected to provide correct Doppler coefficients during the first 0.05 seconds of the transient simulation, provided that macro-calculation with homogenized/condensed cross sections are correctly reproducing MOC reaction rates over the 50 tabulated points of the multicompo.
The CABRI core is far to heterogeneous to avoid an equivalence step before performing the macrocalculation. A similar observation was made by Labouré et al in Ref. [2] for the high temperature test reactor (HTTR). Labouré proposed a SPH correction of the homogenized/condensed cross sections before performing the macro-calculation. Our approach is an extension of the Labouré proposal. We have implemented a SPH equivalence technique in non-fundamental mode condition between a CABRI full-core model solved with the method of characteristics (MOC) in 2D and a Cartesian SP3 macro-calculation, as depicted in Fig. 1. An objective function is set as the root mean square (RMS) error (MOC-SP3 discrepancy) on absorption distribution and leakage rates defined over the macro-geometry in DONJON5. Our algorithm is a quasi-Newtonian gradient search based on the memory limited Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method with backtracking-Armijo inexact line search, as described in Refs. [3] to [6]. The indirect contributions to the gradient vector are computed using generalized perturbation theory (GPT), using the approach described in Ref. [3].

THE SPH EQUIVALENCE TECHNIQUE IN NON-FUNDAMENTAL MODE CASES
Fundamental mode condition is the case where no neutron is leaking due to the boundary conditions. In the case where the macro-calculation over macro-group g is done in non-fundamental mode condition, it is proposed in Ref. [3] to apply a SPH correction on the albedo functions corresponding to boundaries with a non-conservative condition in the reference calculation. If the macro calculation is performed in diffusion approximation, the albedo function Λ(β g ) corresponding to a non-conservative boundary is defined as where β g is the albedo in macro-group g. The net current J g (r) escaping the domain at point r of the boundary is given by the albedo boundary condition as where ∂V β is the fraction of the domain where the non-conservative boundary condition is applied and N (r) is the outgoing normal unit vector.
The macro-calculation used with the CABRI core is based on simplified P N (SPN) approximation, so that Eq. (2) should be replaced by a Marshak boundary condition. The primal formulation of the Marshak albedo boundary condition is written where is an odd integer. The classical Marshak coefficients in Eq. (3) are defined in term of Legendre polynomials P (μ) as The dual formulation of the Marshak albedo boundary condition, consistent with Raviart-Thomas finite element methods, is written where is an even integer.
In order to preserve the neutron balance in macro-group g, cross section data and albedo functions must all be SPH corrected. The correction specific to albedo functions is writteñ where M is the total number of macro-regions, Λ * g is the albedo function of the reference calculation in macro-group g and μ M +1,g is the additional SPH factor assigned to the albedo function EPJ Web of Conferences 247, 03019 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124703019 in group g. This correction technique is proposed as an alternative to the discontinuity factor correction used in Ref. [2].
In fundamental mode condition and in cases where Eq. (4) is used, an infinity of SPH factor sets can satisfy the reference reaction rates in each macro-group g. A unique set is selected with the application of an arbitrary normalization condition. The simplest option is to use the fluxvolume normalization condition which consists to preserve the averaged flux in the lattice. This normalization condition, satisfied in each macro-group g, is written M m=1 Vm whereF m,g and F * m,g are the volume-integrated macro and reference fluxes in macro-region V m and macro-group g.

Equation (7) can be rewritten as
where μ m,g are the SPH factors assigned to cross sections in group g.
The proposed algorithm for obtaining the SPH factors is based on the minimization of a RMS functional where the gradient vector is obtained using generalized perturbation theory (GPT) and where the Hessian matrix is evaluated using the LBFGS recusion.
The RMS error on absorption distribution is an homogeneous functional of the flux defined as where the components f m,g {φ(r)} are the M + 2 conditions to satisfy in each macro-group. They are defined as with P * a,m,g = reference (or target) absorption rates obtained from the full-core reference transport calculation; P * f,m,g = reference (or target) ν-fission rates obtained from the full-core reference transport calculation; EPJ Web of Conferences 247, 03019 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124703019 Δ a,m,g = low limit absorption rates defined as max 10 −4 P * a,tot , P * a,m,g in order to avoid division by small numbers; L * g = reference leakage in macro-group g; Λ(r) = albedo function defined on the non-conservative boundaries ∂V of the domain; Δ L,g = low limit leakage defined as max 10 −4 P * f,tot , L * g + P * a,g in order to avoid division by small numbers; and where P * a,g = m P * a,m,g , P * a,tot = g P * a,g , P * f,tot = m g P * f,m,g and F * g = m F * m,g .
The condition m = M + 1 in Eq. (10) is based on the preservation of the effective multiplication factor of the core. The SPH normalization relations (6) are included in the RMS error in order to simplify the optimization process.

NUMERICAL RESULTS
We investigated the benefit of a quasi-Newtonian SPH technique on two calculation points of the 3parameter grid for the 2D CABRI full core. We selected parameters α bar = 0.0/1.0 (Hafnium bars withdrawn/inserted), T f = 318.2 K and P He3 = 11.3 bars. Due to the many heterogeneities of the CABRI core, no solution of the non-linear system exists and the optimization process converges towards a non-zero minimal value of the RMS functional. The fixed point iterative method diverges and the Newton method for unconstrained optimization is prohibitive in terms of CPU cost and memory consumption. Moreover, both BFGS and LBFGS methods return a non-descent direction if the recursive formula used to compute the Hessian is applied to many times. This issue was not observed on the simpled B&W Core XI investigated by Hébert in a previous study. [3] A simple fix to this issue is to restart the Hessian matrix recursion back from the steepest-descent case at each 10 iterations. In our previous publication, a complete error reduction was reported without having to restart the LBFGS algorithm.
Statistics are given on the effective multiplication factor and on the fission map accuracy. The fission densities are defined by the following relation: where P i is the total fission rate and V i is the volume of the region i. The maximum and averaged errors are respectively defined by: where P fiss * i is computed using the reference fission rates (as computed by DRAGON5) and V core is the total volume of the regions where the fission rate is not equal to zero.
After 400 iterations of the LBFGS method with restart, the RMS error defined by Eq. (9) is reduced below 10 −3 . A summary of the results is given in Table 1, showing the convergence of the SPH method. The relative error distribution on fission rates is shown in Figs. 2 and 3 for Hafnium bars withdrawn and inserted, respectively.

CONCLUSIONS
We have demonstrated the capability to reduce the condensation and homogenization error for a highly heterogeneous core using a new implementation of the SPH method consistent with nonfundamental mode cases. Many orders of magnitude error reduction was obtained for this core. In addition, this new approach paves the way to the modelling of power transient in the CABRI reactor and more generally to reactor cores with complex geometries and high leakage rates.