Verification and Validation of the SP_L Module of the Deterministic Code AZNHEX through the Neutronics Benchmark of the CEFR Start-Up Tests

Abstract

A new module for the AZtlan Nodal HEXagonal (AZNHEX) code, which is part of the AZTLAN Platform, was recently developed based on the Simplified Spherical Harmonics ($S{P}_L$) scheme to deal with the challenges presented in small fast reactor cores, such as the China Experimental Fast Reactor (CEFR), with high leakage and significant scattering effects. For the verification and validation process, we generated nodal homogenized macroscopic cross-sections (XS) through a full heterogeneous core model using the stochastic code SERPENT and subsequently, these XS were employed in AZNHEX. To verify the $S{P}_L$ implementation, several mesh sensitivity exercises were performed demonstrating that the $S{P}_L$ module was implemented successfully. Furthermore, to validate the code with this new implementation, we modeled some exercises contained in the CEFR benchmark with AZNHEX and compared the results with the experimental data available. The final results show a great improvement compared with the original diffusion solver reducing the deviations significantly from experimental data. In conclusion, it is shown and discussed the relevance of improved numerical models (transport approximations instead of diffusion) for the deterministic calculations of small fast reactors.

1. Introduction

The AZTLAN Platform [1] emerged in Mexico as a national project to develop a software platform for the analysis and design of nuclear reactors. This project is an initiative led by the National Institute of Nuclear Research (ININ) and its objective is to take significant steps towards positioning Mexico, in the midterm, at a competitive international level in the development of software for nuclear reactor analysis and modeling. The AZTLAN Platform project consists of a thermal-hydraulic code (namely AZtlan THErmohydraulics Core Analysis (AZTHECA)), a neutron transport code (namely AZtlan TRANsport (AZTRAN) [2]), two neutron diffusion codes, namely AZtlan KInetics in Neutron Diffusion (AZKIND) [3] and AZtlan Nodal HEXagonal (AZNHEX) [4], and a code for sensitivity and uncertainty analysis (namely AZtlan Tool for Uncertainty and SensItivity Analysis (AZTUSIA) [5]). Most of these codes have already passed through a test phase by means of different benchmarks, comparing the produced results with those obtained by codes globally used in the nuclear community and, when available, with experimental data.

The AZNHEX code [4] was originally developed as a 3D neutron diffusion code for nuclear reactor core analysis with hexagonal-Z geometry. Two new solvers are under development for AZNHEX: a Discrete Ordinates ($S_N$) solver [6] and a Simplified Spherical Harmonics ($S{P}_L$) solver [7]. The present work shows the results of a verification and validation (V&V) exercise of the AZNHEX new $S{P}_L$ solver by simulating the exercises contained in the Neutronics Benchmark of the China Experimental Fast Reactor (CEFR) Start-Up Tests [8] organized by the IAEA and proposed by the China Institute of Atomic Energy (CIAE). The experimental data of some of the CEFR start-up tests were used for this objective. Several institutes and universities around the world have participated in the reactor core modeling for these exercises, and many other researchers have validated other deterministic codes [9,10], which indicates that a proper verification and validation test with this benchmark can be done.

2. Description of China Experimental Fast Reactor

The China Experimental Fast Reactor (CEFR) is a 65 MWth pool-type sodium-cooled fast reactor constructed and operated by the China Institute of Atomic Energy (CIAE). It is the first fast reactor in China, which reached criticality for the first time in 2010. Figure 1 shows a cross-section of the reactor core, including the color identifiers of the in-core materials.

When fully loaded, a total of 79 fuel subassemblies (SAs) are present in the core, but it can reach criticality with only 72 and one control rod inserted at a certain height. The fuel SAs positions are originally occupied by Mock-Up fuel SAs, which are specially made SAs with similar geometry as regular fuel SAs but with no fuel present, and which are replaced by fuel SAs in stages. Other CEFR core components include 8 control rods SAs (3 safety rods, 3 shim rods, and 2 regulatory rods), 338 stainless steel SAs, 230 shielding SAs, and 1 neutron source SA [8].

3. Description of the Codes and Models

3.1. AZNHEX Code Background

The original AZNHEX code solves the neutron diffusion equation with elements of hexagonal-Z geometry. To do this, the code takes a hexagonal element and uses the Gordon-Hall transformation [11,12] to turn it into 4 squares in which the Raviart-Thomas-Nédélec order 0 (RTN-0) nodal method [13,14] can be applied for the spatial dependency of the equation. The interested reader can find a detailed description of the implementation in [2,15,16,17]. For the energy and time dependency, the code uses a multi-group approach and the $\theta$ method, respectively.

A previous effort on the code Verification and Validation process was done by means of the Benchmark for Neutronic Analysis of Sodium-cooled Fast Reactor Cores with Various Fuel Types and Core Sizes [18] organized by the Nuclear Energy Agency (NEA/OECD). The results showed good concordance with other participants in the Benchmark and were published in [2].

This previous effort showed that the results obtained by the diffusion code are good enough in a large-size core, nevertheless when it was used in a smaller core the deviation from the reference values was much higher. This can be since the diffusion model does not properly represent what happens in areas with high neutron leakage, which are much more present in a small core compared with a larger one. This proved the need for an improved solver for a wider range of accurate problem solving and, thus, justify the implementation of improved numerical models in the AZNHEX code.

$S{P}_L$ Implementation in AZNHEX

Proposed by Gelbard in the decade of 1960s [19,20,21], the $S{P}_L$ approximation is derived from the Spherical Harmonics Theory (the $P_L$ approximation) applied in the one-dimensional transport equation. This way, the original proposal consists of the simple replacement of the spatial partial derivatives in the one-dimensional form of the $P_L$ approximation by multivariable operators [19,22]. The resulting system of equations is made-up of a total of $L + 1$ expressions per energy group, which is considerably easier to solve (even in time-dependent cases) than the ${(L+1)}^2$ system (per energy group) [23] of the three-dimensional, full spherical harmonics method.

The so-called canonical form of the $S{P}_L$ approximation, which is derived from the previously mentioned $L + 1$ system, has been implemented in the AZNHEX code, re-shaping it to take the form of the neutron diffusion equation. The process that generates this diffusion-like equation also involves the calculation of “modified” cross-sections (XS) and diffusion coefficients, and the amount of these new values will depend on the order L selected by the user. Starting from the canonical form of the $S{P}_L$ approximation (Equation (1)):

$$-\frac{{\mu }_n^2}{{\Sigma }_{t_g}}{\nabla }^2{\mathsf{\Phi }}_n^g(\overrightarrow{r})+{\Sigma }_{t_g}{\mathsf{\Phi }}_n^g(\overrightarrow{r}) = \sum _{g^{\prime }=1}^G\left[{\Sigma }_{s_0}^{g^{\prime }\to g}{\varphi }_0^{g^{\prime }}(\overrightarrow{r})+S_0^{g^{\prime }}(\overrightarrow{r})\right]\mathrm{where}n = 1,\cdots ,\frac{L+1}{2},$$

(1)

and with the scalar neutron flux defined as follows (Equation (2)):

$${\varphi }_0^g(\overrightarrow{r}) = \sum _{n=1}^{(L+1)/2}{\omega }_n{\mathsf{\Phi }}_n^g(\overrightarrow{r}),$$

(2)

where the value of $\mathsf{\Phi }$ depends on the $S{P}_L$ index. If the sum is truncated after the first term when $L = 1$, $\mathsf{\Phi }$ results equal to the scalar neutron flux ${\varphi }_0^g$ (diffusion).

By replacing Equation (2) in the general form of $S{P}_L$ aproximation (Equation (1)) and expanding the neutron source term (without considering external sources), Equation (3) is obtained:

$$-\frac{{\mu }_n^2}{{\Sigma }_{t_g}}{\nabla }^2{\mathsf{\Phi }}_n^g(\overrightarrow{r})+{\Sigma }_{t_g}{\mathsf{\Phi }}_n^g(\overrightarrow{r}) = \sum _{g^{\prime }=1}^G{\Sigma }_{s_0}^{g^{\prime }\to g}\sum _{m=1}^{(L+1)/2}{\omega }_m{\mathsf{\Phi }}_m^{g^{\prime }}(\overrightarrow{r})+\frac{1}{k}{\chi }_g\sum _{g^{\prime }=1}^G{\nu }_{g^{\prime }}{\Sigma }_{f_{g^{\prime }}}\sum _{m=1}^{(L+1)/2}{\omega }_m{\mathsf{\Phi }}_m^{g^{\prime }}(\overrightarrow{r});$$

(3)

it is possible to manipulate Equation (3) to resemble a neutron diffusion-like equation for heterogeneous (Equation (4)) and homogeneous (Equation (5)) media since the goal of this is to be able to keep using the diffusion solver with only minimal modifications. Thus:

$$\begin{array}{cc}\hfill & -\nabla ·{\mathtt{D}}_i(\overrightarrow{r})\nabla {\mathsf{\Phi }}_i(\overrightarrow{r})+{\Sigma }_{\mathtt{R}i}(\overrightarrow{r}){\mathsf{\Phi }}_i(\overrightarrow{r}) = \sum _{j\ne i}^{((L+1)G)/2}{\Sigma }_{\mathtt{s}i\to j}(\overrightarrow{r}){\mathsf{\Phi }}_j(\overrightarrow{r})+{\chi }_i\sum _{j=1}^{((L+1)G)/2}\nu {\Sigma }_{\mathtt{f}j}(\overrightarrow{r}){\mathsf{\Phi }}_j(\overrightarrow{r})\hfill \end{array}$$

(4)

and

$$\begin{array}{cc}\hfill & -{\mathtt{D}}_i{\nabla }^2{\mathsf{\Phi }}_i(\overrightarrow{r})+{\Sigma }_{\mathtt{R}i}{\mathsf{\Phi }}_i(\overrightarrow{r}) = \sum _{j\ne i}^{((L+1)G)/2}{\Sigma }_{\mathtt{s}i\to j}{\mathsf{\Phi }}_j(\overrightarrow{r})+{\chi }_i\sum _{j=1}^{((L+1)G)/2}\nu {\Sigma }_{\mathtt{f}j}{\mathsf{\Phi }}_j(\overrightarrow{r})\hfill \end{array}$$

(5)

where: $i = 1,\cdots ,\frac{(L+1)G}{2}$,

$$\begin{array}{cc}\hfill & {\mathtt{D}}_{i(n,g)} = \frac{{\mu }_n^2}{{\Sigma }_{t_g}} = 3{\mu }_n^2{D}_g,\hfill \\ \hfill & {\Sigma }_{\mathtt{R}i(n,g)} = {\Sigma }_{t_g}-{\omega }_n{\Sigma }_{s_0}^{g\to g},\hfill \\ \hfill & \nu {\Sigma }_{\mathtt{f}i(n,g)} = {\omega }_n{\nu }_g{\Sigma }_{f_g},\hfill \\ \hfill & {\chi }_{i(n,g)} = {\chi }_g,\hfill \\ \hfill & {\Sigma }_{\mathtt{s}j(m,g^{\prime })\to i(n,g)} = {\omega }_n{\Sigma }_{s_0}^{g^{\prime }\to g},\hfill \end{array}$$

(6)

with: $dn,m=1,\cdots ,\frac{L+1}{2}$; where $g,g^{\prime }=1,\cdots ,G$;

the terms ${\mu }_n$ and ${\omega }_n$, which were introduced in Equations (1) and (2) and used to calculate the values of Equation (6), are constants that are calculated by solving the equation systems which come from the original $L + 1$ system and that led to the canonical form of the $S{P}_L$ approximation. In the particular case of the implementation in AZNHEX, only $S{P}_1$, $S{P}_3$, $S{P}_5$, and $S{P}_7$ were used, and their proper ${\mu }_n$ and ${\omega }_n$ values are given in

Table 1. Values of constants ${\mu }_n$ y ${\omega }_n$ for given $S{P}_L$ aproximation order.

n = 1,⋯,(L+1)/2	μn	ωn
$L=1$
1	$1/\sqrt{3}$	1
$L=3$
1	0.339981043584856	0.652145154862546
2	0.861136311594053	0.347854845137454
$L=5$
1	0.238619186083197	0.467913934572691
2	0.661209386466265	0.360761573048139
3	0.932469514203152	0.171324492379170
$L=7$
1	0.183434642495650	0.362683783378360
2	0.525532409916329	0.313706645877890
3	0.796666477413627	0.222381034453374
4	0.960289856497536	0.101228536290376

.

As we mentioned before, the application of the $S{P}_L$ method in AZNHEX is aimed to odd values of L up to 7 (i.e., $S{P}_1$, $S{P}_3$, $S{P}_5$, and $S{P}_7$). This way, considering that the number of XS and diffusion coefficient values will increase by two times for the $S{P}_3$ approximation, three times for the $S{P}_5$ approximation, and four times for the $S{P}_7$ approximation, the resulting system of diffusion-like equations for each odd order of L will increase the same way. To solve one of these systems, the AZNHEX code will implement an artificial energy mesh structure in which matrix size will vary according to the number of XS and diffusion coefficient new values (see [7] for full details on this implementation).

This method allows solving any three-dimensional energy-dependent problem without the need of using relatively high computational power and just by adding some minor modifications to the input file to be solved. Even if the user wants to implement the $S{P}_7$ approximation in the resolution of a given problem, it is expected that the required number of computational resources will be significantly less in comparison with a transport code. As a matter of example, if we consider a simple comparison of the number of unknowns to be found in a transport methodology, for instance in a Discrete Ordinates ($S_N$) solver, with G energy groups and $NC$ unknowns per discrete ordinate, if N is the order of discrete ordinate, the total number of unknowns ($TNU$) to be solved would be: $TNU = NC\times (N+2)\times N\times G$. In the case of $S{P}_L$ approximation, $TNU = NC\times G\times (L+1)/2$, being L, the order of $S{P}_L$ approximation. In this case, assuming that the $S_N$ and $S{P}_L$ approximations are equivalent when $N = L + 1$, for an approximation $S_8$ in discrete ordinates $TNU = 80\times NC\times G$ meanwhile, for $S{P}_7$, $TNU = 4\times NC\times G$, i.e., 20 times less unknowns like in $S_8$ transport approach (8 times in the case of $S_4$ vs. $S{P}_3$, and 16 times for the case of $S_6$ vs. $S{P}_5$). Of course, this analysis is done just in the number of unknowns of the problem independent from solvers and parallel algorithms that could be implemented in both methodologies. However, just having 20 times more unknowns as is the case for the transport $S_8$ implies that there is surely a big difference in computational time from what AZNHEX ($S{P}_L$) would require.

Since the $S{P}_L$ approximation has its basis on a transport scheme, the results that are obtained through the implementation of this method are in general more accurate than the ones that are generated by using the conventional diffusion scheme. Nevertheless, it is important to mention that the design of this $S{P}_L$ solver has been made neglecting the higher-order transport effects of the expansion of the neutron scattering kernel. Additionally, the required boundary and interface conditions, which are specific for every odd order of the $S{P}_L$ approximation presented here, have not been included in the formulation of the $S{P}_L$ solver.

3.2. CEFR Core Modelling

3.2.1. Model for Serpent

Since AZNHEX is not able to generate its own XS, it is necessary to generate them by using another code. To perform this task, the Serpent code version 2.1.30 [24] was used due to its capability to generate macroscopic XS that can be employed to feed deterministic codes, which gives the opportunity to verify and validate the AZNHEX code by simulating more realistic scenarios and demonstrate its calculation capabilities. Additionally, the Serpent code is widely used in reactor calculations due to its high precision in 3D simulations and great capabilities and advantages over other Monte Carlo codes, which allow obtaining accurate reference solutions.

The 3D model of the CEFR reactor core was built considering the material specifications contained in the benchmark technical documentation such as the type of geometry of the core elements, their dimensions, the isotopic composition of the materials, etc. Since the experiments are carried out at different temperatures, thermal expansion effects have an impact on the materials, and they should be considered in the simulations. To simplify these calculations, we assumed that all the materials of the reactor core have the same linear expansion coefficient of $1.8\times {10}^{-5}{/}^{\circ }$C, which is the value for the stainless steel established in the benchmark specifications [8] and it is present in all the SAs. It also corresponds to the highest coefficient of the materials present in the core. To keep the same masses, the densities of all the materials were smeared except for the sodium coolant, whose density was always obtained using Equation (7) in agreement with technical specifications given in [8]. It is defined as:

$$\rho = 950.0483-0.2298T-14.604\times {10}^{-6}{T}^2+5.6377\times {10}^{-9}{T}^3,$$

(7)

where T is the temperature expressed in °C and $\rho$ is the density in kg/m${}^3$.

3.2.2. Cross-Section (XS) Generation

To simplify the XS generation, the subassemblies (either fuel, control, reflector, or shielding) are grouped into rings of SAs of the same type. This methodology reduces considerably the amount of data that is handled. Additionally, all the subassemblies are divided into 14 axial zones; hence, for every ring of SAs, 14 sets of XS are generated. These XS are then extracted from Serpent through a Python script and are then used in AZNHEX. The core discretization can be seen in Figure 2.

For the XS energy multigroup structure, we used six energy groups when generating the XS sets, in agreement with the optimization of multi-group energy structures study presented in [25]. This number of energy groups is appropriate to reduce the computational time in nodal calculations with AZNHEX. The energy bins structure is shown in

Table 2. Neutron energy groups employed.

Group	Upper Limit (MeV)
1	2.000000 × 10${}^1$
2	1.353400 × 10${}^0$
3	5.23400 × 10${}^{-1}$
4	6.73790 × 10${}^{-2}$
5	3.35460 × 10${}^{-3}$
6	7.48520 × 10${}^{-4}$

.

As an example, the group constants generated with Serpent in the axial zone 7 (from bottom to top as it can be seen in Figure 2 above) for the Fuel subassemblies grouped in Ring 1 (closest to the center of the core), are shown in

Table 3. Group constants generated with Serpent.

Material	Energy Group			Group Constants
		${\mathit{D}}_{\mathit{g}}$	${\mathbf{\Sigma }}_{{\mathit{r}}_{\mathit{g}}}$	$\mathit{\nu }$∗${\mathbf{\Sigma }}_{{\mathit{f}}_{\mathit{g}}}$	$\mathit{\kappa }$∗${\mathbf{\Sigma }}_{{\mathit{f}}_{\mathit{g}}}$	${\mathit{\chi }}_{\mathit{g}}$	${\mathbf{\Sigma }}_{\mathit{s}{\mathit{g}}^{\prime }\to \mathit{g}}$
Fuel SA’s (Ring 1)	1	2.968	4.290 × 10${}^{-2}$	2.205 × 10${}^{-2}$	2.031 × 10${}^2$	5.794 × 10${}^{-1}$	1.220 × 10${}^{-1}$	2.474 × 10${}^{-2}$	8.885 × 10${}^{-3}$
							3.425 × 10${}^{-4}$	2.323 × 10${}^{-6}$	4.666 × 10${}^{-7}$
	2	2.150	2.978 × 10${}^{-2}$	1.529 × 10${}^{-2}$	2.023 × 10${}^2$	2.825 × 10${}^{-1}$	0.00	1.739 × 10${}^{-1}$	2.223 × 10${}^{-2}$
							4.088 × 10${}^{-4}$	3.151 × 10${}^{-6}$	1.542 × 10${}^{-7}$
	3	1.527	1.407 × 10${}^{-2}$	1.747 × 10${}^{-2}$	2.022 × 10${}^2$	1.303 × 10${}^{-1}$	0.00	0.00	2.352 × 10${}^{-1}$
							4.892 × 10${}^{-3}$	1.206 × 10${}^{-6}$	6.873 × 10${}^{-8}$
	4	1.066	1.811 × 10${}^{-2}$	2.767 × 10${}^{-2}$	2.022 × 10${}^2$	7.495 × 10${}^{-3}$	0.00	0.00	0.00
							3.248 × 10${}^{-1}$	1.253 × 10${}^{-3}$	8.235 × 10${}^{-7}$
	5	6.931 × 10${}^{-1}$	5.512 × 10${}^{-2}$	7.978 × 10${}^{-2}$	2.022 × 10${}^2$	8.053 × 10${}^{-5}$	0.00	0.00	0.00
							0.00	5.400 × 10${}^{-1}$	3.866 × 10${}^{-3}$
	6	7.796 × 10${}^{-1}$	9.172 × 10${}^{-2}$	1.488 × 10${}^{-1}$	2.022 × 10${}^2$	9.843 × 10${}^{-6}$	0.00	0.00	0.00
							0.00	0.00	3.634 × 10${}^{-1}$

.

4. Results

4.1. Verification and Validation Exercise

Since the fuel loading and criticality experiment is the one (among all Start-up experiments) that has direct reactivity values (and that can also be directly calculated with nodal solvers), we chose it for the verification and validation process of AZNHEX and to perform a sensitivity mesh refinement study. In this section, the fuel loading criticality exercise of the benchmark is described in detail, and the results obtained, first with Serpent and then with $S{P}_L$ implementation in AZNHEX and their comparison with the experimental data are presented.

Description of Fuel Loading and Criticality Experiment

Before the start-up of the reactor, the core is preliminarily loaded with mock-up fuel SAs in the active fuel positions. The reactor moves to first criticality by replacing, step by step, the mock-up SAs with fuel SAs. In the sub-critical extrapolation process, the number of fuel SAs to be loaded is determined by extrapolation of reciprocal of count rate and following safety requirements; the process is described in [26] and shown in Figure 3.

When the core approaches to criticality, a super-critical extrapolation (that uses the control rods to reach criticality by the period method) is then used. In this extrapolation process, 72 fuel SAs were loaded. Control rods, with exception of one regulating rod that is fully inserted, are totally withdrawn to the out-of-core position. The regulating control rod is then withdrawn step by step to three different positions to reach super-criticality. At each position, a positive period is obtained. Based on that, the critical position of the control rod is predicted by extrapolation (based on the calculated control rod worth curve). Finally, the control rod is put to the predicted critical position, and the reactor’s clean-core criticality state is reached. For the CEFR, the final clean-core criticality state was reached with 72 fuel SAs and the regulating control rod at the position of 70 mm with a measured sodium temperature of 245 ${}^{\circ }$C. Three start-up detectors located near the active core (temporarily installed) were used to get the count rate throughout the criticality approaching process [26].

4.2. Verification and Validation of Serpent Model

Since calculations with codes directly provide the value of $k_{eff}$, we used the following equation to report the reactivity in each case:

$$\rho = (k_{eff}-1.0)/k_{eff};$$

(8)

The reference calculation was the critical case. A study on the selection of the evaluated data library was done to ensure that the numerical results calculated with both Serpent and AZNHEX, would have minimum differences in comparison with those obtained experimentally. Thus, eight different libraries were tested with the same Serpent model (for the critical case). The evaluated data files were: BROND-3.1 [27], CENDL-3.1 [28], ENDF/B-VII.1 [29], ENDF/B-VIII.0 [30], JEFF-3.1 [31], JEFF-3.1.2 [32], JEFF-3.2 [33] and JENDL-4.0 [34].

As mentioned in Section 4.1, the core was modeled with 72 fuel SAs inserted and all control rods (CRs) out (position 500 mm) except for the Regulating Rod 2 (RE2) which is inserted in position 70 mm. This state was selected because the core is expected to reach criticality under those conditions, as experiments showed [8].

The exercise showed that the ENDF/B-VIII.0 library is the best suited for this case, resulting in a reactivity of only 4 pcm above the expected reactivity value (${\rho }_{exp} = 0.0$ for the critical case), as shown in Figure 4.

For this reason, we selected the ENDF/B-VIII.0 library for all the numerical calculations done in this work.

The results of the Serpent model for the expected calculations in terms of reactivity, with an estimated value of 0.0 in the final critical position, are presented in

Table 4. Reactivity values and deviation in pcm from Serpent calculations and experimental data.

RE2 Position	Exp. Measurement	Serpent	Absolute Dev.
	${\mathit{\rho }}_{\mathit{exp}}$ (pcm)	${\mathit{\rho }}_{\mathit{Serpent}}$ (pcm)	${\mathit{\rho }}_{\mathit{exp}}-{\mathit{\rho }}_{\mathit{Serpent}}$
190 mm	40	48.0	−8.0
170 mm	34	41.0	−7.0
151 mm	25	28.0	−3.0
70 mm	0	4.0	−4.0

.

4.3. Verification and Validation of $S{P}_L$ Implementation in AZNHEX

Table 5. Reactivity values and deviation in pcm from AZNHEX calculations and experimental data.

RE2 Position	Exp. Measurement	AZNHEX (pcm)
	${\mathit{\rho }}_{\mathit{exp}}$	${\mathit{\rho }}_{{\mathit{SP}}_1}$	Dev	${\mathit{\rho }}_{{\mathit{SP}}_3}$	Dev	${\mathit{\rho }}_{{\mathit{SP}}_5}$	Dev	${\mathit{\rho }}_{{\mathit{SP}}_7}$	Dev
190 mm	40	−2031.34	2071.34	356.52	−316.52	−143.08	183.08	−123.57	163.57
170 mm	34	−2032.81	2066.81	355.13	−321.13	−144.24	178.24	−124.73	158.73
151 mm	25	−2038.18	2063.18	349.97	−324.97	−149.21	174.21	−129.70	154.70
70 mm	0	−2065.01	2065.01	324.15	−324.15	−174.29	174.29	−154.75	154.75

shows the reactivity values obtained from the AZNHEX $k_{eff}$ results employing all the $S{P}_L$ approximations and the deviation of these results with respect to the experimental data for the supercritical and critical steps.

The original version of AZNHEX only included the diffusion approximation ($S{P}_1$) and, as expected, the results obtained with the $S{P}_1$ approximation have the highest discrepancies compared with the experimental values and with other $S{P}_L$ approximations since the $S{P}_1$ method corresponds to diffusion and is not able to properly represent the real conditions in this kind of small and high leaking fast reactors. It is clear that as the index of the $S{P}_L$ approximation increases, the results are more accurate when compared with the reference results, with differences below 200 pcm in the $S{P}_5$ and $S{P}_7$ cases.

In Figure 5, the same results of Table 5 are presented but in $k_{eff}$ values, where it is possible to appreciate the difference in employing different $S{P}_L$ approximations.

Regarding the computational times, in

Table 6. Comparison of computational times for different $S{P}_L$ approximations.

	${\mathit{SP}}_{\mathbf{1}}$	${\mathit{SP}}_{\mathbf{3}}$	${\mathit{SP}}_{\mathbf{5}}$	${\mathit{SP}}_{\mathbf{7}}$
Time [s]	237.7	2260.4	4244.9	6546.5
Time factor vs. $S{P}_1$	1	9.5	17.8	27.5
Time factor vs. $S{P}_3$	–	1	1.8	2.9

, a comparison among the different times of $S{P}_L$ approximations is done. It can be seen that, there is a significative difference between $S{P}_1$ and $S{P}_3$ (and, of course, with the other $S{P}_5$ and $S{P}_7$ approximations). The reason is due to the fact that, in the original version of AZNHEX (pure diffusion implementation), the scattering matrix is numerically adjusted in order to avoid up-scattering terms, for this reason, the inner iteration process finalized in significant less steps. However, in the methodology used in this $S{P}_L$ implementation, and described with detail in [7], the introduction of higher orders of $S{P}_L$ implies the creation of pseudo energy-groups, i.e., the original matrix (diffusion) is expanded with the new $S{P}_L$ coefficients. The order of the matrix is multiplied by 2, 3, and 4, respectively for $S{P}_3$, $S{P}_5$, and $S{P}_7$ including “up-scattering” terms between pseudo groups. The consequence of this artificial up-scattering effect, results in an increase of the inner iterations, and thus of the computational times almost for 10 times between $S{P}_1$ and $S{P}_3$. It could be significant amount of time, however, such implementation allows us to move from more than 2000 pcms of difference to 350 pcms. In the final line of Table 6, it is possible to see the time factor between the $S{P}_5$ and $S{P}_7$, comparing with $S{P}_3$, in this case, the result is very consistent with the size of the matrix to be solved in each case.

In Figure 6, the normalized radial power distribution obtained with Serpent and AZNHEX $S{P}_1$ to $S{P}_7$, with a $2 \times 2$ mesh refinement, is shown. The mesh refinement applied in AZNHEX is described in the following subsection. The values were normalized to 1, this is, the average power value is equal to 1. The figure presents also, the path line to track the power. As it can be seen, the maximum value is located in the first assembly, which is near to the center of the core, for all the calculations. The AZNHEX results present higher discrepancies in the border of the core (point 5) which can be related to the high leakage effect in this interphase.

Mesh Sensitivity Analysis with AZNHEX $S{P}_L$

An extra analysis was performed for this exercise in which different mesh refinements were implemented to approximate even better the AZNHEX results to the experimental data. These refinements were done radially and axially along the reactor core. For a deeper analysis, besides the combination of radial and axial refinements, both refinements were also done separately. Figure 7 and Figure 8 show the way these radial and axial refinements are done in AZNHEX.

The

Table 7. Reactivity values and deviation in pcm from AZNHEX calculations with radial refinement and experimental data.

RE2 Position		190 mm	Dev	170 mm	Dev	151 mm	Dev	70 mm	Dev
Exp.
Measurement		40 (pcm)		34 (pcm)		25 (pcm)		0 (pcm)
	$S{P}_1$_$1\times 1$	−2031.34	2071.34	−2032.81	2066.81	−2038.18	2063.18	−2065.01	2065.01
	$S{P}_1$_$2\times 1$	−1649.18	1689.18	−1650.81	1684.81	−1656.31	1681.31	−1683.73	1683.73
	$S{P}_1$_$4\times 1$	−1508.35	1548.35	−1510.34	1544.34	−1515.89	1540.89	−1543.48	1543.48
	$S{P}_3$_$1\times 1$	356.72	−316.72	354.74	−320.74	349.77	−324.77	323.95	−323.95
	$S{P}_3$_$2\times 1$	702.06	−662.06	700.42	−666.42	695.17	−670.17	668.81	−668.81
AZNHEX	$S{P}_3$_$4\times 1$	827.93	−787.93	826.24	−792.24	821.08	−796.08	794.40	−794.40
	$S{P}_5$_$1\times 1$	−143.20	183.20	−144.21	178.21	−149.22	174.22	−174.30	174.30
	$S{P}_5$_$2\times 1$	202.30	−162.30	201.03	−167.03	195.93	−170.93	170.28	−170.28
	$S{P}_5$_$4\times 1$	327.31	−287.31	326.01	−292.01	320.88	−295.88	295.01	−295.01
	$S{P}_7$_$1\times 1$	−124.15	164.15	−125.16	159.16	−130.17	155.17	−155.24	155.24
	$S{P}_7$_$2\times 1$	221.92	−181.92	220.65	−186.65	215.56	−190.56	189.99	−189.99
	$S{P}_7$_$4\times 1$	346.73	−306.73	345.43	−311.43	340.30	−315.30	314.57	−314.57

shows the reactivity values obtained from the AZNHEX $k_{eff}$ results employing the $S{P}_L$ approximations with radial refinement and the deviation of these results from the experimental data.

It can be seen in Table 7 that larger differences are obtained with the radial refinements in comparison with the no refinement cases ($1\times 1$) for all $S{P}_L$ approximations in the four cases of approach to criticality process, getting the best results with the $S{P}_7$_$1\times 1$ refinement. Additionally,

Table 8. Reactivity values and deviation in pcm from AZNHEX calculations with axial refinement and experimental data.

RE2 Position		190 mm	Dev	170 mm	Dev	151 mm	Dev	70 mm	Dev
Exp.
Measurement		40 (pcm)		34 (pcm)		25 (pcm)		0 (pcm)
	$S{P}_1$_$1\times 1$	−2031.34	2071.34	−2032.81	2066.81	−2038.18	2063.18	−2065.01	2065.01
	$S{P}_1$_$1\times 2$	−2302.14	2342.14	−2303.70	2337.70	−2309.08	2334.08	−2336.42	2336.42
	$S{P}_1$_$1\times 4$	−2374.52	2414.52	−2376.13	2410.13	−2381.54	2406.54	−2408.99	2408.99
	$S{P}_3$_$1\times 1$	356.72	−316.72	354.74	−320.74	349.77	−324.77	323.95	−323.95
	$S{P}_3$_$1\times 2$	99.90	−59.90	98.90	−64.90	92.91	−67.91	66.96	−66.96
AZNHEX	$S{P}_3$_$1\times 4$	28.99	11.01	27.99	6.01	22.00	3.00	−4.00	4.00
	$S{P}_5$_$1\times 1$	−143.20	183.20	−144.21	178.21	−149.22	174.22	−174.30	174.30
	$S{P}_5$_$1\times 2$	−398.58	438.58	−399.59	433.59	−404.63	429.63	−430.85	430.85
	$S{P}_5$_$1\times 4$	−470.20	510.20	−471.21	505.21	−476.26	501.26	−502.51	502.51
	$S{P}_7$_$1\times 1$	−124.15	164.15	−125.16	159.16	−130.17	155.17	−155.24	155.24
	$S{P}_7$_$1\times 2$	−378.43	418.43	−379.43	413.43	−384.47	409.47	−410.68	410.68
	$S{P}_7$_$1\times 4$	−449.01	489.01	−451.03	485.03	−456.07	481.07	−481.31	481.31

shows the reactivity values obtained from the AZNHEX $k_{eff}$ results employing the $S{P}_L$ approximations with only axial refinement and the deviation of these results concerning the experimental data. It is possible to observe that, with exception of $S{P}_3$ approximation, the axial refinement has significant reactivity differences. The main reason, could be related with the increase in the dimensions ratio between axial and radial nodes. The aspect ratio is defined as the ratio of the longest dimension to the shortest dimension of a quadrilateral element. In many cases, as the aspect ratio increases, the inaccuracy of the solution increases, and in general, an element yields best results if its shape is compact and regular, a more detailed description of this phenomenon is presented in [35].

In order to decrease the aspect ratio and to make shape compact and more regular, a combination of radial and axial refinements was used. In

Table 9. Reactivity values and deviation in pcm from AZNHEX calculations with radial and axial refinements and experimental data.

RE2 Position		190 mm	Dev	170 mm	Dev	151 mm	Dev	70 mm	Dev
Exp.
Measurement		40 (pcm)		34 (pcm)		25 (pcm)		0 (pcm)
	$S{P}_1$_$1\times 1$	−2031.34	2071.34	−2032.81	2066.81	−2038.18	2063.18	−2065.01	2065.01
	$S{P}_1$_$2\times 2$	−1916.62	1956.62	−1918.37	1952.37	−1923.88	1948.88	−1951.81	1951.81
	$S{P}_1$_$4\times 4$	−1846.55	1886.55	−1848.37	1882.37	−1853.94	1878.94	−1882.17	1882.17
	$S{P}_3$_$1\times 1$	356.72	−316.72	354.74	−320.74	349.77	−324.77	323.95	−323.95
	$S{P}_3$_$2\times 2$	450.05	−410.05	448.43	−414.43	443.18	−418.18	416.24	−416.24
AZNHEX	$S{P}_3$_$4\times 4$	509.14	−469.14	507.44	−473.44	502.13	−477.13	475.00	−475.00
	$S{P}_5$_$1\times 1$	−143.20	183.20	−144.21	178.21	−149.22	174.22	−174.30	174.30
	$S{P}_5$_$2\times 2$	−49.11	89.11	−50.67	84.67	−55.78	80.78	−81.86	81.86
	$S{P}_5$_$4\times 4$	9.11	30.89	7.63	26.37	2.46	22.54	−23.91	23.91
	$S{P}_7$_$1\times 1$	−124.15	164.15	−125.16	159.16	−130.17	155.17	−155.24	155.24
	$S{P}_7$_$2\times 2$	−28.41	68.41	−29.79	63.79	−34.89	59.89	−60.94	60.94
	$S{P}_7$_$4\times 4$	29.99	10.01	28.54	5.46	23.34	1.66	−3.01	3.01

the reactivity values obtained from the AZNHEX $k_{eff}$ results employing the $S{P}_L$ approximations by using both radial and axial refinements and the deviation of these results from the experimental data are shown.

It is possible to see that a better consistency in the results is obtained when combined refinements (axial and radial) are used, mainly in the highest order $S{P}_L$ approximations ($S{P}_5$ and $S{P}_7$) obtaining the best results with the $S{P}_7$_$4\times 4$ refinement approximation, having differences below 11 pcm.

The computational cost of such refinement is presented in

Table 10. Comparison of computational times for different refinements.

Refinement	${\mathit{SP}}_{\mathbf{1}}$		${\mathit{SP}}_{\mathbf{3}}$		${\mathit{SP}}_{\mathbf{5}}$		${\mathit{SP}}_{\mathbf{7}}$
	Time [s]	Factor	Time [s]	Factor	Time [s]	Factor	Time [s]	Factor
$1\times 1$	237.7	1	2260.4	1	4244.9	1	6546.5	1
$2\times 2$	1764.95	7.4	18,009.6	7.9	33,045.0	7.8	49,684.02	7.6
$4\times 4$	13,454.95	56.6	138,142.54	61.1	245,042.81	57.7	355,767.57	54.3

. It is possible to see that, in general, moving from no refinement ($1\times 1$) to full refinement ($4\times 4$) increases the computational time with a factor of approximately 60 whereas for an intermediate refinement ($2\times 2$) the computational time just increases 7.5 times. From Table 9, for $S{P}_5$ and $S{P}_7$, $2 \times 2$ results are very acceptable and the computational time is also suitable, thus, a $2 \times 2$ refinement could be seen as an optimal election.

5. Discussion

Constant improvement is part of the software development. Even though the original version of the AZNHEX code was successfully tested in a previous set of exercises, it failed when used in a smaller core; therefore, an improvement on the mathematical solver was needed. In this work, a recently developed $S{P}_L$ solver (based on the Spherical Harmonics Theory), along with its validation process, were presented. As part of the objectives of the AZTLAN Platform, the designing and incorporation into AZNHEX of this $S{P}_L$ solver are aimed to increase the numerical precision of the code itself, making AZNHEX a reliable computational tool for designing and analyzing nuclear reactors made up of subassemblies with hexagonal-z geometry such as the Russian Lead-Cooled Fast Reactor, BREST [36] under development, the Korean Sodium-Cooled Fast Reactor, KALIMER [37] and the Advanced Lead Fast Reactor European Demonstrator, ALFRED [38]. This is especially relevant since the newest generations of nuclear reactors are characterized by this hexagonal-z geometry of its subassemblies.

The structure of the $S{P}_L$ solver incorporated into AZNHEX allows its relatively easy implementation in other diffusion codes of the AZTLAN Platform, since there is only the need of modifying the values of the XS and the diffusion coefficients into an artificial array of energy groups.

On the other hand, the methodology employed to generate the XS data with Serpent was very useful to simplify the amount of data to handle for the AZNHEX models in conjunction with the use of 6 neutron energy groups. This resulted in a considerable reduction of the computing time in both codes and in to getting very accurate results with AZNHEX, as can be seen in this paper.

About the results themselves, a great improvement was seen when using the new $S{P}_L$ solver in AZNHEX, getting very similar values to those obtained with Serpent and with the experimental ones. Referring to the data shown in Table 5, it can be seen that the major improvement in the precision of the results occurs between the $S{P}_1$ and the $S{P}_3$ approximations. Although the $S{P}_5$ and the $S{P}_7$ schemes offer results that are closer to the reference data, the implied computational cost of solving a given problem through such approximations is, in most cases, unnecessarily expensive for the desired precision. Furthermore, the data presented in Table 7, Table 8 and Table 9 allow considering in detail the impact of the radial and axial refinements. Although a reduction of the deviation with increased $S{P}_L$ refinements was expected, many values within these results present a different behavior. For example, in Table 7, it can be seen that, for the diffusion approximation ($S{P}_1$), the deviation is smaller as the radial refinement increases. Nevertheless, for the $S{P}_3$, $S{P}_5$, and $S{P}_7$ schemes, the absolute deviation shows a trend to increase along with the radial refinement. Besides, the absolute deviations for the $S{P}_7$ approximation, considering superior radial refinements, are bigger than the equivalent values for the $S{P}_5$ approximation. Additionally, for the reactivity data shown in Table 8, the impact of the axial refinement behaves somewhat differently. While the absolute deviation increases alongside the axial refinement for every order of the $S{P}_L$ approximation but the $S{P}_3$ approximation, such deviations are smaller as the odd order of L increases.

In order to tackle such inconsistencies (maybe related to the big differences on radial and axial sizes), a combined radial and axial refinement was used. In this case, a better consistency in the results is obtained mainly in the highest order $S{P}_L$ approximations ($S{P}_5$ and $S{P}_7$) where the best results were obtained, as expected, with the $S{P}_7$_$4\times 4$ refinement approximation.

However, it is important to clarify that some of the previously mentioned unexpected variations in the reactivity deviations could be attributed to the absence of the aforementioned superior transport effects of the scattering kernel and, to a lesser extent, to the lack of the required interface and boundary conditions. Since one of the main advantages of the $S{P}_L$ approximation is its higher precision in the description of the scalar neutron flux in highly heterogeneous media, the absence of the anisotropic behavior of the neutron scattering has a significant effect in small cells where the non-consistent change in the spatial dimensions could affect the flux gradient. We planned further investigations in this direction [39].

6. Conclusions

Despite the unexpected changes in the reactivity deviations, it can be concluded that the implementation of the $S{P}_L$ module in the current version of AZNHEX fulfilled the objective to improve the estimation of the results for a small fast reactor such as the CEFR. More accurate estimations were obtained using $S{P}_5$ and $S{P}_7$ approximations in most of the exercises, where the differences concerning the experimental values were around 100 pcm. This indicates that higher L values lead to relatively closer values to the reference data having also lower computational times in comparison with a stochastic code.

Further development of the $S{P}_L$ solver described in this work involves the restructuring of its mathematical background to appropriate include the higher order elements of the expansion of the scattering kernel with the required interface and boundary conditions, to describe with more precision the neutron scalar flux and every involved anisotropy term. It is expected that such a restructuring process allows to correct the variations of the reactivity deviations, despite the implemented order of the $S{P}_L$ approximation.

Additionally, the restructuring of the $S{P}_L$ solver includes the addition of a time-dependent module, which will allow analyzing phenomena such as the variations of the scalar neutron flux and the reactivity deviations as the rods within the core move from one position to another.

Finally, a Discrete Ordinates $S_N$ module is also being implemented in AZNHEX. Further work will be performed to compare these new solvers against stochastic and experimental data.