# Application of the SPH equivalence technique to the CABRI full core in non-fundamental mode

Alain Hébert, Julien Taforeau, Jean-Jacques Ingremeau

► To cite this version:

Alain Hébert, Julien Taforeau, Jean-Jacques Ingremeau. Application of the SPH equivalence technique to the CABRI full core in non-fundamental mode. PHYSOR 2020 – International Conference on Physics of Reactors: Transition to a Scalable Nuclear, American Nuclear Society, Mar 2020, Cambridge, United Kingdom. pp.03019, 10.1051/epjconf/202124703019. irsn-04106635

HAL Id: irsn-04106635

https://irsn.hal.science/irsn-04106635

Submitted on 25 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
APPLICATION OF THE SPH EQUIVALENCE
TECHNIQUE TO THE CABRI FULL CORE
IN NON-FUNDAMENTAL MODE

Alain Hébert¹, Julien Taforeau², and Jean-Jacques Ingremeau²

¹Polytechnique Montréal
P.O. Box 6079, Station Centre-Ville, Montréal QC

²Institut de Radioprotection et de Sûreté Nucléaire (IRSN)
31 avenue de la Division Leclerc, 92260, Fontenay-aux-Roses, France

alain.hebert@polymtl.ca, julien.taforeau@irsn.fr, jean-jacques.ingremeau@irsn.fr

ABSTRACT

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 quasi-Newtonian gradient search based on the Limited memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method. Numerical results are presented with Hafnium bars withdrawn or inserted.

KEYWORDS: SPH equivalence, CABRI reactor, Limited memory Broyden-Fletcher-Goldfarb-Shanno method, unconstrained optimization

# 1. INTRODUCTION

The objective of the CABRI International Programe (CIP) conducted by the Institute for Radio-protection 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.[1]

<figure><img src="image_2.png" /><figcaption>Figure 1: Full-core SPH equivalence over the CABRI core.</figcaption></figure>

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 fictitious absorption rate $L_g = d_g B^2 \phi_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 \text{ 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:

* Fuel temperature $T_f$: (293.15 K, 553.15 K, 923.15 K, 1573.15 K, 2000.0 K)

* $^3$He pressure in transient rods $P_{He3}$: (0.5 bar, 1.0 bar, 2.0 bar, 8.0 bar, 15.0 bar)

* Hafnium bars insertion $\alpha_{bar}$: (0.0, 1.0).

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 macro-calculation. 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].
## 2. 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 $\Lambda(\beta_g)$ corresponding to a non-conservative boundary is defined as

$$ \Lambda(\beta_g) = \frac{1(1 - \beta_g)}{2(1 + \beta_g)} \qquad (1) $$

where $\beta_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

$$ -J_g(r) \cdot N(r) + \Lambda(\beta_g) \phi_g(r) = 0 \quad \text{if } r \in \partial V_{\beta} \qquad (2) $$

where $\partial V_{\beta}$ 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

$$ \Phi_{\ell,g}(r) \cdot N(r) - 2\Lambda(\beta_g) \sum_{\substack{\ell'=0 \\ \ell' \text{ even}}}^{N-1} \mathcal{M}_{\ell,\ell'} \phi_{\ell',g}(r) = 0 \quad \text{if } r \in \partial V_{\beta} \qquad (3) $$

where $\ell$ is an odd integer. The classical Marshak coefficients in Eq. (3) are defined in term of Legendre polynomials $P_\ell(\mu)$ as

$$ \mathcal{M}_{\ell,\ell'} = (2\ell' + 1) \int_{0}^{1} d\mu P_{\ell}(\mu) P_{\ell'}(\mu) . \qquad (4) $$

The dual formulation of the Marshak albedo boundary condition, consistent with Raviart-Thomas finite element methods, is written

$$ 2\Lambda(\beta_g)\phi_{\ell,g}(r) - \sum_{\substack{\ell'=1 \\ \ell' \text{ odd}}}^{N} \mathcal{M}_{\ell,\ell'}[\Phi_{\ell',g}(r) \cdot N(r)] = 0 \quad \text{if } r \in \partial V_{\beta} \qquad (5) $$

where $\ell$ 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 written

$$ \tilde{\Lambda}_g = \mu_{M+1,g} \Lambda_g^* \qquad (6) $$

where M is the total number of macro-regions, $\Lambda_g^*$ is the albedo function of the reference calculation in macro-group g and $\mu_{M+1,g}$ is the additional SPH factor assigned to the albedo function
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 flux-volume normalization condition which consists to preserve the averaged flux in the lattice. This normalization condition, satisfied in each macro-group g, is written

$$ \sum_{m=1}^{M} \int_{V_m} d^3 r \, \tilde{\phi}_g(\mathbf{r}) \equiv \sum_{m=1}^{M} \tilde{F}_{m,g} = \sum_{m=1}^{M} F_{m,g}^*, \quad g \le G \quad (7) $$

where $\tilde{F}_{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

$$ \sum_{m=1}^{M} \frac{F_{m,g}^*}{\mu_{m,g}} = \sum_{m=1}^{M} F_{m,g}^*, \quad g \le G \quad (8) $$

where $\mu_{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

$$ \mathcal{F}\{\phi(\mathbf{r})\} = \sum_{m=1}^{M+2} \sum_{g=1}^{G} (f_{m,g}\{\phi(\mathbf{r})\})^2 \quad (9) $$

where the components $f_{m,g}\{\phi(\mathbf{r})\}$ are the $M+2$ conditions to satisfy in each macro-group. They are defined as

$$ f_{m,g}\{\phi(\mathbf{r})\} = \begin{cases} \frac{\langle\Sigma_a, \phi\rangle_{m,g} P_{a,tot}^* - P_{a,m,g}^*}{\langle\Sigma_a, \phi\rangle \Delta_{a,m,g}} & \text{if } m \le M \\ \sqrt{M} \left( \frac{\langle\Lambda + \Sigma_a, \phi\rangle_g P_{f,tot}^* - \frac{L_g^* + P_{a,g}^*}{\Delta_{L,g}}}{\langle\nu\Sigma_f, \phi\rangle \Delta_{L,g}} \right) & \text{if } m = M+1 \\ \frac{1}{F_g^*} \sum_{j=1}^{M} \frac{F_{j,g}^*}{\mu_{j,g}} - 1 & \text{if } m = M+2 \end{cases} \quad (10) $$

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) $\nu$-fission rates obtained from the full-core reference transport calculation;
$\Delta_{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;

$\Lambda(r)$ = albedo function defined on the non-conservative boundaries $\partial V$ of the domain;

$\Delta_{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}^* = \sum_m P_{a,m,g}^*$, $P_{a,tot}^* = \sum_g P_{a,g}^*$, $P_{f,tot}^* = \sum_m \sum_g P_{f,m,g}^*$ and $F_g^* = \sum_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.

### 3. NUMERICAL RESULTS

We investigated the benefit of a quasi-Newtonian SPH technique on two calculation points of the 3-parameter grid for the 2D CABRI full core. We selected parameters $\alpha_{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 simpied 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:

$$P_i^{\text{fiss}} = \frac{\sum_i V_i}{V_i \sum_i P_i}$$

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:

$$\epsilon_{\max} = \max_i \frac{|P_i^{\text{fiss}} - P_i^{\text{fiss}*}|}{P_i^{\text{fiss}*}} \quad \text{and} \quad \bar{\epsilon} = \frac{1}{V_{\text{core}}} \sum_i \left[ \frac{|P_i^{\text{fiss}} - P_i^{\text{fiss}*}|}{P_i^{\text{fiss}*}} \right] V_i$$

where $P_i^{\text{fiss}*}$ is computed using the reference fission rates (as computed by DRAGON5) and $V_{\text{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.
<table><caption>Table 1: Summary of SPH convergence for the CABRI core</caption><thead><tr><th rowspan="2"></th><th colspan="2">Hafnium bars withdrawn</th><th colspan="2">Hafnium bars inserted</th></tr><tr><th>no SPH</th><th>after SPH</th><th>no SPH</th><th>after SPH</th></tr></thead><tbody><tr><td>ε<sub>max</sub></td><td>21.4</td><td>0.177</td><td>18.3</td><td>0.239</td></tr><tr><td>ε̄</td><td>7.99</td><td>0.068</td><td>10.6</td><td>0.071</td></tr><tr><td>K<sup>*</sup><sub>eff</sub> - K<sup>sp3</sup><sub>eff</sub></td><td>-310.4 pcm</td><td>11.9 pcm</td><td>-43.5 pcm</td><td>-3.1 pcm</td></tr><tr><td>F{φ(r)}</td><td>3.78</td><td>6.57×10<sup>-4</sup></td><td>6.08</td><td>8.12×10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Percent error on the fission rates – Hafnium bars withdrawn.</figcaption></figure>

4. 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 non-fundamental 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.
<figure><img src="image_4.png" /><figcaption>Figure 3: Percent error on the fission rates – Hafnium bars inserted.</figcaption></figure>

ACKNOWLEDGEMENTS

This work was supported by a grant from the Natural Science and Engineering Research Council of Canada.

REFERENCES

[1] A. Hébert, “DRAGON5 and DONJON5, the contribution of École Polytechnique de Montréal to the SALOME platform,” *Ann. nucl. Energy*, **87**, 12–20 (2016).

[2] V. Laboure, Y. Wang, J. Ortensi, S. Schunert, F. Gleicher, M. DeHart and R. Martineau, “Hybrid Super Homogenization and Discontinuity Factor Method for Continuous Finite Element Diffusion,” *Ann. nucl. Energy*, **128**, 443–454 (2019).

[3] A. Hébert, “A general presentation of the SPH equivalence technique in non-fundamental mode cases,” *Ann. nucl. Energy*, **141**, 1–10 (2020).

[4] L. Armijo, “Minimization of functions having Lipschitz continuous first partial derivatives,” *Pacific journal of mathematics*, Vol. 16, **1**, 1–3 (1966).

[5] J. E. Dennis and J. Moré, “Quasi-Newton methods, motivation and theory,” *SIAM Review*, Vol. 19, **1**, 46–89 (1977).

[6] J. Nocedal, “Updating Quasi-Newton Matrices With Limited Storage,” *Math. of computation*, Vol. 35, **151**, 773–782 (1980).