On the simulation of neutron noise using a discrete ordinates method

The method of discrete ordinates is investigated for neutron noise simulations in the frequency domain. For this purpose, the solver NOISE-SN is developed and used to simulate two neutron noise problems that are respectively derived from the two-dimensional systems described in the neutron transport simulation benchmarks C4V and C5G7. In the ﬁrst problem based on the C4V system, NOISE-SN is compared to the diffusion-based simulator CORE SIM+. These results show that NOISE-SN and CORE SIM+ calculate similar spatial distributions of neutron noise, although signiﬁcant differences can be found at the location of the perturbation and at locations with strong variations of material properties, where the discrete ordi-nates method is expected to be more accurate than diffusion theory. Then NOISE-SN calculations are performed to test different S N approximations, and the ﬁctitious source method that may be applied to mitigate possible numerical artifacts, known as the ray effect. In the second problem based on the C5G7 system, the choice of a low order of discrete ordinates in NOISE-SN leads to unphysical values of the neutron noise because of the ray effect. The increase of the order of discrete ordinates or introducing a ﬁctitious source in the equations to be solved alleviates the issue. The second option is shown to remove the ray effect without a high order of discrete ordinates and thus without too expensive calculations, even though the strength of the ﬁctitious source needs to be tuned carefully to avoid very slow convergence rates.


Introduction
During normal, steady-state operations of a nuclear reactor, neutron flux measurements show small fluctuations around mean values.These fluctuations are referred to as neutron noise, and they may be driven by a variety of perturbations, such as mechanical vibrations of core internals, disturbances in the coolant flow, etc. From the analysis of the neutron noise, it is possible to identify anomalous patterns at an early stage so that appropriate actions can be taken to maintain plant availability.
Neutron noise analysis can consequently be used to develop monitoring techniques for enhanced nuclear reactor safety, see e.g.Demazière et al. (2018).For this purpose, the reactor transfer function, which describes the core response to a possible perturbation, is often required.The modeling of the reactor transfer function can be based on the neutron transport equation, while the possible perturbations are expressed in terms of changes in the macroscopic nuclear cross sections.As discussed in Pázsit (1992), a coarser estimation of the reactor transfer function may be sufficient for reactor diagnostics, and most of the work in the field of reactor neutron noise relies on neutron diffusion theory, e.g.Demazière (2011) and Mylonakis et al. (2020).However, recent efforts also focused on the development of stochastic methods (Yamamoto, 2013;Rouchon et al., 2017) and higher-order deterministic methods (Bahrami and Vosoughi, 2018) to provide more detailed analyses and assess the limitations of the diffusion approximation for neutron noise applications.The work presented in the paper investigates the method of discrete ordinates (or S N method) for neutron noise simulations in the frequency domain, using the solver NOISE-SN.Such a solver has already been tested on different numerical problems.For instance, as reported in (Vinai et al., 2021), NOISE-SN (indicated as ''the discrete ordinates solver developed by Chalmers") has been used to simulate the neutron noise induced by different types of neutron noise sources in a simplified nuclear fuel assembly and compared with Monte Carlo and deterministic transport solvers, such as TRIPOLI-4 Ò (Rouchon et al., 2017) and APOLLO3 Ò (Rouchon et al., 2021).The analysis of the calculated static neutron flux and neutron noise has shown a good agreement between NOISE-SN and the other transport solutions.
When applying the discrete ordinates method, the discretization of the angular variable cannot account for all possible directions of neutron travel.It can thus generate numerical artifacts, the so-call ray effect.Such an effect has been widely studied in sta-tic calculations and leads to anomalous oscillations of the scalar and angular flux, e.g., see Lathrop (1968) and Lathrop (1971).The problem may become severe for systems with isolated neutron sources and for strong absorbing media with low scattering ratios.In the paper, the issue is discussed for the case of frequencydomain neutron noise calculations.To mitigate these possible numerical distortions, the increase of the order of discrete ordinates and the introduction of a fictitious source in the equations to be solved are explored.The second option follows the method proposed by Miller and Reed (1977) and is attractive because it can lead to good results without using high orders of discrete ordinates and thus avoids large numbers of expensive transport sweeps and excessive requirements on computer memory.
The investigation is based on simulations of two neutron noise problems that are respectively derived from the 2-D system configurations described in the neutron transport simulation benchmarks C4V (Lefebvre et al., 1991) and C5G7 (Lewis et al., 2001).In the C4V neutron noise problem, NOISE-SN is compared with the diffusion-based noise simulator CORE SIM+, which was used to simulate various neutron noise numerical tests and experiments, e.g.Mylonakis et al. (2021).The objective is to study the possible differences between NOISE-SN and CORE SIM+.In addition, the effect of the order of discrete ordinates is analyzed and the algorithm based on the fictitious source method is tested.In the C5G7 neutron noise problem, the impact of the ray effect on the NOISE-SN simulations is shown and the two mitigation strategies (i.e., increasing the order of discrete ordinates and introducing a fictitious source in the discrete ordinates equation) are evaluated.
The paper is organized as follows.In Section 2, the transport neutron noise equation and the solver NOISE-SN are introduced.In Section 3, the NOISE-SN calculations of the C4V neutron noise problem are discussed.In Section 4, the C5G7 neutron noise problem is studied with NOISE-SN.In Section 5, conclusions are drawn.

The solver NOISE-SN
The solver NOISE-SN is based on the transport neutron noise equation in the frequency domain.The discretization of the equations is performed according to a finite difference, discrete ordinates, multi-energy group method.The algorithm for the numerical solution is accelerated using the CMFD method.The solver allows to simulate 2-D and 3-D systems with reflective and vacuum boundary conditions.For the mitigation of possible ray effects in 2-D problems, the equations are modified with a fictitious source term.

Neutron noise equation
In the case of a nuclear reactor, the multi-energy-group timedependent neutron transport equation with a generic number of families of delayed-neutron precursors is written as: where w is the angular neutron flux, / is the scalar neutron flux and C is the concentration of delayed-neutron precursors.The subscripts g and g 0 ¼ 1; Á Á Á ; G identify the energy groups and q ¼ 1; Á Á Á ; Q the families of precursors of delayed neutrons.In this work, neutron noise problems in stationary systems are considered, and hence the macroscopic fission cross sections are normalized with k eff to adjust the system to criticality.For steady-state conditions, the time derivatives in Eqs. ( 1) and ( 2) are zero and the static neutron flux w g;0 r; X ð Þ satisfies the equation: The neutron noise equation in the frequency domain used for NOISE-SN is derived in the case that linear theory is applicable.Small, stationary perturbations of the macroscopic nuclear cross sections are assumed.The other system parameters, including the effective multiplication factor k eff , are considered constant with respect to time.In the time-dependent equations, the neutron flux, the concentration of precursors and the nuclear cross sections are modelled as the sum of a static mean value and a fluctuating part.Eq. ( 3) is subtracted from Eqs. ( 1) and ( 2), the second-order perturbation terms are neglected and a temporal Fourier transform is performed.Then, the neutron noise equation in the frequency domain reads as: The quantities dw g and d/ g are the induced angular and scalar neutron noise, respectively, and take complex values.The variable x ¼ 2pf is the angular frequency of the driving perturbation.The imaginary unit is denoted as i.The term S g r; X; x ð Þin Eq. ( 4) is the neutron noise source, i.e.: S g r; X; x ð Þ¼À dR t;g r; The quantity dR is the perturbation associated with the nuclear cross section R.

Discretization of the equations
Eqs. ( 4) and ( 5) correspond to a fixed source problem whose solution requires the effective multiplication factor k eff and the neutron fluxes calculated under static conditions.Therefore, NOISE-SN consists of a static module and a module for the neutron noise simulations in the frequency domain.
The static module solves the criticality problem given by Eq. (3).The discretization of the equation is based on a standard discrete ordinates method (Lewis and Miller, 1984).The angular flux is evaluated along discrete directions, and the scalar flux is constructed from the angular flux via the Legendre-Chebyshev P N À T N quadrature set (Longoni, 2004).The spatial differencing relies on the diamond finite difference method.Possible negative values of the neutron flux are corrected using the set-to-zero fixup.The neutron scattering term is approximated by an L-order real spherical harmonics expansion (Handb. Nucl. Eng., 2010).
The frequency-domain neutron noise module calculates the complex-valued neutron noise using the static solution and the neutron noise source prescribed in terms of fluctuations of nuclear cross sections.The discretization of Eqs. ( 4) and ( 5) is similar to the static case.A set of discrete directions is selected, where the n-th discrete direction X n is defined by the direction cosines l n , g n and n n with respect to the orthogonal spatial coordinates x; y; z ð Þ.The equation is spatially discretized using the diamond finite difference method over a grid of rectangular cuboids.The centers of the grid cells are identified by triples of integer mesh indices I; J; K ð Þ associated with the coordinates x; y; z ð Þ.The faces of the cells are given by triples where one half is added to or subtracted from the integer mesh index related to the direction perpendicular to the cell face.The volumes of the cells are V I;J;K ¼ Dx I Â Dy J Â Dz K .Then, the discretized forms of Eqs.(4) and (5) read as: and S g;n;I;J;K x ð Þ ¼ ÀdR t;g;I;J;K x ð Þw g;n;I;J; In Eq. ( 6), the unknown is the angular neutron noise dw n .For the closure of Eqs. ( 6) and ( 7), the diamond finite-difference relationships between the values at the centers and at the faces of a cell are used, i.e.: In Eq. ( 6), the unknown is the angular neutron noise dw n .For the closure of Eqs. ( 6) and ( 7), the diamond finite-difference relationships between the values at the centers and at the faces of a cell are used, i.e.: In addition, the moments of the neutron noise d/ m l;g;I;J;K are obtained from the angular neutron noise with the following quadrature formula: The real spherical harmonics R m l (which are also needed for the neutron scattering terms in Eqs. ( 6) and ( 7)) are given in Handb.Nucl.Eng. ( 2010).The scalar neutron noise d/ g;I;J;K is the moment with l ¼ m ¼ 0. The directions and the relative weights in Eq. ( 9) are selected according to the Legendre-Chebyshev P N À T N quadrature set.The weights are normalized as: In Eqs. ( 6) ands ( 7), the fission spectrum v dyn g;I;J;K includes the contributions from both the prompt and delayed neutrons, i.e.: The expression of the noise source in Eq. ( 7) needs the static angular flux w g;n;I;J;K;0 , the flux moments / m l;g;I;J;K;0 and the effective multiplication k eff , which are provided via the static module.

CMFD method for acceleration
A typical inner-outer iterative procedure is used for NOISE-SN to solve the multi-energy-group static neutron transport equations and the multi-energy-group neutron noise equations.In an inner iteration, a transport sweep with respect to the discrete angular directions is performed to estimate the neutron fluxes for one energy group and construct the with-in group self-scattering term.After the inner iterations are completed from the highest to the lowest energy group, the neutron fluxes are used to evaluate the down-scattering term.Then, the fission source is updated in the outer iteration and a new set of inner iterations starts.In the outer iteration of the static calculation, the effective multiplication factor k eff is also updated.The procedure is repeated until a convergence criterion is met.
When solving the neutron transport equation with an innerouter iteration scheme, acceleration techniques are needed to improve the slow numerical convergence.Several studies showed that the Coarse Mesh Fine Difference -CMFD method allows an efficient acceleration of static and time-dependent neutron transport calculations, e.g., see Zhong et al. (2008) and Zhu et al. (2015).In the development of NOISE-SN, the CMFD method is found to also have a good performance for the frequency-domain neutron noise calculations, see details in Yi et al. (2021).The criticality calculations in the static module of NOISE-SN rely on a standard adCMFD method (Jarrett et al., 2016).A similar algorithm is followed for the CMFD acceleration of the neutron noise calculations in the frequency domain, and its details are discussed below.
After one inner-outer iteration, the estimated neutron noise is used to construct a low-order diffusion-like equation.This CMFD equation is spatially discretized over a grid that is coarser than the grid of the original problem and is solved to provide a correction of the neutron noise.If convergence is not reached, the scheme moves to the next inner-outer iteration, otherwise it stops.One full iterative loop then consists of one inner-outer iteration and the successive CMFD calculation.The notation is such that the index ITE þ 1=2 indicates the inner-outer iteration within the full ITE-th iterative loop.
For the CMFD problem, the mesh indices IC; JC; KC ð Þidentify the center of the cells of the coarser grid.The faces of the cells are given by triples in which one half is subtracted from or added to the mesh index associated with the direction perpendicular to the cell face.The CMFD equation is derived from Eq. ( 6) by applying to all the terms.The operation P n w n Á ð Þ represents an integration over the angular variable using the P N À T N quadrature set and the operation P The coarse mesh diffusion coefficients in Eq. ( 14) are calculated using the fine mesh static fluxes / m¼0 l¼0;g;I;J;K;0 : The coarse mesh parameters in Eq. ( 12) are computed on the fly using the scalar neutron noise d/ ITEþ1=2 g;I;J;K .The coarse mesh cross sections are given as: R dyn t;g;IC;JC;KC ¼ The fission spectrum term v dyn g;IC;JC;KC for the coarser mesh is also homogenized via the scalar neutron noise calculated in the innerouter iteration, i.e.: The neutron noise source term S g;IC;JC;KC is derived from a simple volumetric homogenization of Eq. ( 7) because it is defined by the problem and thus is fixed.
The CMFD equation is a fixed source problem.Then, considering the overall coarser grid, a system of linear equations like Eq. ( 12) is built and solved for the scalar noise dU If reflective boundaries are specified, incoming angular neutron noise quantities are also updated using the same procedure.
A convergence criterion is checked.The relative differences between the two last iterations for both the real and imaginary parts of the scalar neutron noise are calculated at each computational cell.If the maximum of the relative differences is below a specific value, the calculation is stopped.If it is not, the corrected neutron noise is used for the next ITE þ 1-th iterative loop.To avoid numerical instability of the CMFD algorithm in the neutron noise calculation, some inner-outer iterations without acceleration may be beneficial before the full iterative loops.
In the initial iteration, a diffusion-based guess is used for the transport sweeps.This guess is the solution of Eq. ( 12) using the static fluxes to homogenize the system parameters and with the correction factor b D equal to zero.

Fictitious source method for ray effect mitigation
In NOISE-SN, a fictitious source method can be used for the mitigation of the ray effects in 2-dimensional simulations.Accordingly, an additional, fictitious source term that preserves the rotational invariance characteristic of the transport equation is introduced to convert the discrete ordinates equation to a spherical-harmonics-like equation (Miller and Reed, 1977).Then the 2-dimensional version of Eq. ( 6) is rewritten as: The fictitious source is denoted as S Fict g;n;I;J and is multiplied by a factor c which can take values between 0 and 1.The choice of the factor c allows to tune the strength of the fictitious source.If c ¼ 0, the original discrete ordinate equation is solved.If c ¼ 1, the mitigation of ray effects is expected to be maximum, but the convergence rate is the worst.The theoretical model of the fictitious source developed by Lathrop (1971) is considered because of its relatively good convergence behavior in static calculations as compared to other types of sources.The model is coded following the procedure developed by Miller and Reed (1977), i.e.: The parameter N is the order of discrete ordinates and the spherical harmonics Y m l are related to the real spherical harmonics R m l , i.e.: The spatial derivatives in Eq. ( 24) are computed with a central finite difference scheme in each computational cell, i.e.: @ @x du 2mÀ1 The quantities du 2mÀ1 g;N;IAE1=2;J and du 2mÀ1 g;N;I;JAE1=2 are the moments at the cell surfaces and are estimated from the angular neutron noise at the cell surfaces, i.e.: The coefficients a m n;l are calculated with Y m l following the procedure given in (Miller and Reed, 1977).
Since the fictitious source depends on the angular neutron noise via Eqs.( 28) and ( 29), it is updated, similarly to the with-in group self-scattering term, after each inner iteration.In the construction of the CMFD equation, the fictitious source term vanishes when the integration with respect to the angular variable is performed using the P N À T N quadrature set.However, it is corrected before the start of a new full iterative loop as follows:

Simulation of a neutron noise problem in C4V
A neutron noise problem is defined using the system described in the C4V two-dimensional two-energy group neutron transport benchmark.The solution obtained from NOISE-SN is compared with the solution from the diffusion-based solver CORE SIM+ in order to identify possible differences between the discrete ordinates method and a less accurate approach such as the neutron diffusion approximation.The impact of the order of discrete ordinates on the NOISE-SN calculation and the correct use of a fictitious source for removing possible numerical oscillations are also tested.

Description of the problem and solution method
The C4V configuration consists of two UO 2 assemblies located in North-West and South-East positions and two MOX assemblies in north-east and south-west, see Fig. 1.Each assembly contains 17 Â 17 square cells.One square cell is 1:26 Â 1:26 cm and may be a fuel or a guide tube cell.The neutron cross sections for the fuel and the guide tube cells are given in two energy groups and for iso- In the neutron noise calculations, the fast velocity v 1 is equal to 1:9919 Â 10 7 cm=s and the thermal velocity v 2 is equal to 3:8184 Â 10 5 cm=s.One family of precursors of delayed neutrons is used, with the effective delayed neutron fraction b ¼ 0:0049 and the decay constant k ¼ 0:0797 s À1 .
The neutron noise source is a perturbation of the fuel cell indicated in red in Fig. 1.It is defined as a stationary fluctuation of the capture cross section in both energy groups.The amplitude of the source in each energy group is 5% of the respective nominal value, and the frequency is equal to 1 Hz.

Comparison between NOISE-SN and CORE SIM+
The static and neutron noise results obtained from NOISE-SN and CORE SIM+ are compared.The NOISE-SN and CORE SIM+ calculations are performed over a 170 Â 170 grid of square meshes equal in size (i.e., Dx ¼ Dy ¼ 0:252 cm).Thus, each fuel or guide tube cell of the C4V system is divided into 5 Â 5 identical, homogeneous squares.The order N of the discrete ordinates approximation is equal to 32.The CMFD acceleration of NOISE-SN makes use of a coarser grid whose computational cells correspond to the fuel and guide tube cells.Before each CMFD calculation, one inner iteration and one outer iteration are carried out.

Static results
The computed values of the effective multiplication factor k eff are very close, i.e., 0.91742 for NOISE-SN and 0.91755 for CORE SIM+.Figs. 2 and 3 show the static fluxes along the main diagonal across the MOX fuel assembly (bottom-left to the top-right corner of the system) and across the UO 2 fuel assembly (top-left to the bottom-right corner of the system), respectively.The mesh position index used for the x-axes of the plots counts the computational cells following the diagonal, from left to right.In general, the two solvers are in agreement.The higher relative differences are found for the thermal flux in the guide tubes of the MOX fuel assemblies, where they reach the value of 16.8%, see Fig. 2.These discrepancies are expected when comparing higher-order transport methods and the diffusion approximation in conditions of strong variations of the material properties, e.g., see (Vinai et al., 2021).The relative differences may also be significant at the boundaries, e.g., in Fig. 2, they are equal to 10.8% for the fast group and 6.3% for the thermal group.

Neutron noise results
As discussed in Section 2, the neutron noise is calculated in the frequency domain, so it is a complex quantity.In the current analysis, the results are presented in terms of amplitude and phase, taken along the main diagonal across the MOX fuel assemblies.This diagonal also goes through the perturbed fuel cell.Again, the mesh position index used for the x-axes of the plots counts the computational cells following the diagonal from left to right.
The neutron noise amplitude calculated with NOISE-SN and with CORE SIM+ are similar, see Fig. 4. The spatial distribution of the neutron noise resembles the spatial distribution of the static flux (compare Figs. 2 and 4) because the small dimension of the system favors a point-kinetic response.Yet, the local impact of the perturbation can be identified in the close surroundings of the noise source location (which corresponds to the mesh position index equal to 94 in Fig. 4).Consistently with the static case, the larger discrepancies between NOISE-SN and CORE SIM+ are found for the thermal noise amplitude in the guide tubes of the MOX fuel assemblies, being the maximum relative difference equal to 16%.At the neutron noise source location, the relative differences are also significant, i.e., 8.9% for the fast group and 14.3% for the thermal group.This outcome may be explained by the higher accuracy expected from a discrete ordinates method compared with a neutron diffusion method, e.g.see (Vinai et al., 2021).The values of the phase of the noise estimated with the two solvers are similar, as shown in Fig. 5.The phase is close to 180 because the neutron flux variation is opposite to the perturbation, i.e., the change of the capture cross section.

Order of discrete ordinates
Calculations are performed to verify the consistent behavior of the S N solver, using different orders N of discrete ordinates, namely equal to 8, 16 and 32.
In Fig. 6 the values are taken along the column of computational cells that are along the right boundary of the system and are counted from the top to the bottom (see Fig. 1).These cells are fuel cells, their material composition changes with the fuel assembly.The different orders of discrete ordinates give similar results for the amplitude (plots on the left) and the phase (plots on the right).When using the S8 approximation, numerical oscillations, even though very small, are found and they mainly affect the phase.These numerical oscillations can be eliminated by increasing the order of discrete ordinates, as the S16 and S32 calculations demonstrate.

Fictitious source method
As discussed in subsection 2.4, the introduction of a fictitious source into the neutron transport equation may suppress numeri-cal issues arising from a coarse discretization of the angular variable.This option is computationally cheaper for each iteration than a simple increase of the order of discrete ordinates.Although the numerical oscillations in the C4V neutron noise problem are very small, the fictitious source method is used for the order N ¼ 8 of discrete ordinates and tested.The S16 and S32 calculations are not tested with the fictitious source method since ray effects are negligible in these two cases.
In Fig. 7, the fast and the thermal neutron noise are obtained from a fictitious source of strengths c that are respectively equal to 0.0 (i.e., no fictitious source), 0.9 and 1.0.The noise is taken along the right boundary of the system, see subsection 3.3.The calculations with and without the fictitious source method provide similar results.The introduction of the fictitious source reduces the minimal numerical oscillations, e.g., see the phase of the noise in Fig. 7.In addition to smoother results, increasing the strength of the fictitious source causes a minor shift of the predictions.

Convergence
The number of full iterations required for the convergence of the calculations and the related wall clock time are reported in Table 1.This computational work is performed with a 1 Â 10core Intel Xeon E5-2630v3 processor.The convergence criterion is such that the maximum relative difference of the local scalar neutron noise between two consecutive iterations must be less than 1 Â 10 À6 .For simulations without the fictitious source (i.e., c ¼ 0) the convergence rate is insensitive to the order of discrete ordinates.On the other hand, the convergence of the S8 approximation deteriorates with stronger fictitious sources.

Simulation of a neutron noise problem in C5G7
A neutron noise problem is derived from the perturbation of the neutron capture cross section in one fuel pin of the twodimensional system defined in the C5G7 neutron transport benchmark (Lewis et al., 2001).The problem is simulated with NOISE-SN, using different orders of discrete ordinates.In addition, the fictitious source method is studied for the mitigation of numerical issues that may arise from the discretization of the angular variable.

Description of the problem and solution method
A detailed description of the C5G7 system is provided in (Lewis et al., 2001) and is summarized in the following.The configuration is 2-dimensional and includes four fuel assemblies arranged in a 2 Â 2 grid and a reflector/moderator water region on the east and south peripheries, see Fig. 8. Two UO 2 fuel assemblies are respectively in the north-west and south-east positions of the 2 Â 2 grid, and two MOX fuel assemblies are respectively in the north-east and south-west positions.Each fuel assembly consists of 17 Â 17 square cells that are equal in size and have a side length of 1.26 cm.These physical cells include a central circular region with a radius of 0.54 cm and a surrounding region.The circular region is either a fuel-clad mix, a guide tube or a fission chamber, while the surrounding region is moderator.In the modelling of the system, the physical cells are transformed according to the square cells of the computational spatial grid, see Section 4.2.A set of sevenenergy group cross sections is assigned to each kind of composition.The neutron kinetic parameters and the data for eight families of delayed neutron precursors are specified in (Hou et al., 2017).
The neutron noise source is defined as a perturbation in one of the MOX fuel pins, whose location is shown in Fig. 8.The perturbation is a fluctuation of the neutron capture cross section over all the energy groups.Its amplitude is 5% of the static value of the cross section and the frequency is equal to 1 Hz.

Computational spatial grids
The solution algorithm used in NOISE-SN is accelerated via a CMFD method, see subsection 2.3.Accordingly, a fine grid is used for the transport sweeps and a coarse grid for the CMFD calculation.
In the fine grid for the transport sweeps, a fuel, guide tube, or fission chamber cell is divided into 4 Â 4 computational cells, see Fig. 9.The central 4 computational cells are identical in size and they, as a whole, approximate the circular fuel pin, guide tube or fission chamber, preserving the actual area.The same spatial discretization is used for the reflector water region.Therefore, the fine grid for the entire system consists of 204 Â 204 computational cells.
For the CMFD calculation, the coarse grid is given by 51 Â 51 computational cells.The size of the computational cells is equal to the size of the fuel, guide tube and fission chamber cells.

Order of discrete ordinates
The neutron noise problem defined in the C5G7 system is simulated with different orders N of discrete ordinates, i.e., equal to 16, 32 and 64.
For the static calculations, one inner iteration and one outer iteration are performed before the adCMFD acceleration.For        the neutron noise calculations, two inner iterations and one outer iteration are performed before solving the CMFD equation.
According to the selected convergence criterion, the static calculations are stopped if the maximum relative differences of the local scalar neutron flux and k eff between two successive iterations are less than 1Â10 À6 .The effective multiplication factor k eff , the number of iterations required for a converged solution, and the wall clock time are reported in Table 2.The reference value for k eff is taken from (Smith et al., 2004).The calculations are performed with 1Â10-core Intel Xeon E5-2630v3 processor.The different approximations need the same number of iterations to converge.The largest discrepancy with respect to the reference k eff value is found for order N ¼ 64, i.e., 105 pcm.In Figs.10-12, the S16, S32 and S64 approximations are used to calculate the static flux, the neutron noise amplitude, and the neutron noise phase, respectively.The values are taken for the first and the seventh energy group, which are representatives of the fast and thermal energy groups, respectively.The plots are made along the diagonal that crosses the MOX fuel assemblies (from the bottomleft to the top-right corner of the core region) and that includes the perturbed fuel cell.The x-axis labels of these plots are such that the mesh index counts, along the diagonal, the computational cells starting from the bottom-left corner.The different orders of discrete ordinates give similar values.The neutron noise amplitude resembles the static flux because the system is small and tends to respond to small perturbations in a point-kinetic manner.Yet, a local effect is found at the location of the neutron noise source; compare Figs. 10 and 11 for the mesh indices equal to 74 and 75.
Although the spatial distributions of the noise calculated with different orders of discrete ordinates are consistent in most of the points, numerical issues are identified in the reflector region.
In some cells at the bottom-right corner of the system, the S16 approximation predicts the noise phase for the first energy group to be equal to 0 degrees, which is unphysical, see Fig. 13.In fact, the noise phase is expected to be 180 because a variation of the number of neutron captures induces an opposite variation in the neutron flux.Some numerical distortions, even if they are much smaller, are also found in the noise phase for the seventh energy group (and for the other energy groups).The reason for these incorrect values is the limited number of angular directions used in the calculation, i.e., the ray effects.The selection of a higher order of discrete ordinates such as N ¼ 64 improves the results, as the comparison between Figs. 13 and 14 shows.Nevertheless, the S64 approximation does not completely remove the numerical inconsistency.Therefore, discrete ordinates of larger orders, which lead to more expensive simulations, are needed.

Fictitious source
For the mitigation of the ray effects, an alternative to the increase of the order of discrete ordinates is the introduction of a fictitious source in the equations to be solved, as discussed in sub-  section 2.4.The fictitious source is included in the S16 neutron noise calculation.
To investigate the effect of a fictitious source on the calculations for the C5G7 neutron noise problem, a column of computational cells is taken along the lower part of the right boundary of the system, in the reflector region, where the S16 approximation predicts unphysical values of the noise phase (see Fig. 13).The fictitious source is selected to have strengths equal to 0.0 (no fictitious source), 0.8 or 0.9.In Fig. 15, the results are compared in terms of the real and imaginary parts of the neutron noise.The Ydirection mesh indices used in the x-axes correspond to the ones in Fig. 13, and the X-direction mesh equals to 204.The solution with unphysical phase values, i.e., calculated with S16 and the fictitious source of strength c ¼ 0, are due to relatively large oscillations of the real and imaginary parts of the neutron noise.Such oscillations are of a numerical type since they are reduced by increasing the order of ordinates or adding a fictitious source.Also, the use of a fictitious source combined with the S16 approximation leads to smoother results than the case of S64.The larger the fictitious source strength is, the smaller the oscillations are.
In Fig. 16 the spatial distributions of the neutron noise for the first and seventh energy group, over the entire system, are obtained from the S16 approximation combined with the fictitious source of strength c ¼ 0:9.Compared to Fig. 13, the noise phase is better predicted without taking unphysical values, although minor distortions in the noise phase, which are due to the small numerical oscillations of the real and imaginary part of the neutron noise, are still present in the reflector region.For complete elimination of the distortions in this problem, a higher order of discrete ordinates or a different fictitious source is necessary.

Convergence
Concerning the use of the fictitious source, tests are carried out to optimize the numerical performance of the solver for the simulation of the C5G7 neutron noise problem.First, the fictitious source of strength c ¼ 1:0 causes the convergence rate to be very slow, while choices such as c ¼ 0:8 or c ¼ 0:9 bring significant improvements.Then, a certain number of ordinary full iterative loops before introducing the fictitious source into the solution scheme is necessary to avoid numerical instabilities.In the current case, the initial 7 iterative loops are run without the fictitious source.In these loops, the CMFD equation is solved over the fine transport grid (instead of the coarser grid) to further enhance the stability of the algorithm.Calculations that use the fictitious source for all the energy groups or only for the most problematic ones such as the lowest and the highest, were also investigated and show similar convergence.This outcome suggests that the convergence rate is mainly governed by the energy groups that suffer from the most severe ray effects.
The convergence behavior of the algorithm for the real and imaginary parts of the scalar neutron noise is given in Fig. 17.The numerical performances of different types of calculations are compared, i.e., using the order of discrete ordinates N equal to 16, 32 and 64 and using the S16 approximation combined with a fictitious source of strength c equal to 0.8 and 0.9.The convergence criterion is such that the relative difference of the local scalar neutron noise between the two last iterations is less than 1 Â 10 À6 in any computational cell.
The number of iterations required for convergence and the wall clock time are shown in Table 3.These neutron noise calculations are performed with 2 Â 10-core Intel Xeon E5-2630v3 processors.The use of the fictitious source in the S16 simulation causes a higher number of iterations, but these iterations are computationally cheaper than the iterations for the S32 and S64 approximations.Therefore, mitigation of the ray effect can be obtained from the fictitious source method with a limited increase of the numerical effort.

Conclusion
The computer program NOISE-SN has been developed to solve the transport neutron noise equation in the frequency domain with the method of discrete ordinates.The solver has already been used to simulate numerical neutron noise problems and compared with other higher-fidelity solvers.In this work, two cases are studied to  further explore the performance of the method for neutron noise applications, i.e., the 2-dimensional systems defined in the C4V and C5G7 benchmarks are perturbed with a prescribed localized neutron noise source.
In the C4V problem, NOISE-SN is compared to the diffusionbased simulator CORE SIM+.These results show that NOISE-SN and CORE SIM+ calculate similar spatial distributions of neutron noise, although significant differences in values can be found in regions where diffusion theory is expected to be less reliable, such as at the location of the perturbation and at locations with strong variations of material properties.Then, NOISE-SN can provide more accurate solutions than the diffusion method.In addition, low orders of discrete ordinates are not affected by any severe ray effects and the NOISE-SN scheme with or without the fictitious source method gives very close results as expected.
In the C5G7 problem, a coarse discretization of the angular directions such as the S16 approximation is shown to estimate unphysical values of neutron noise.The choice of higher orders of discrete ordinates improves the results and avoids incorrect values, although ray effects are quite persistent and are still visible in the neutron noise calculated, e.g., with the S64 approximation.The fictitious source method allows to use low orders of discrete ordinates and reduce the impact of the ray effect.However, the strength of the fictitious source needs to be carefully adjusted to avoid very slow convergence rates.For instance, the S16 approximation with a fictitious source of strength 0.9 provided satisfactory results in terms of both reduction of ray effects and convergence.
As shown in the analysis of C4V and C5G7 problem, the application of the fictitious source affects the solution and the performance of the computational scheme.In problems where ray effects are significant, a parametric study is required to select an optimum fictitious source strength that allows an efficient balance between ray-effect mitigation and the convergence rate of the iterative process.
Future work is necessary for the evaluation of the discrete ordinates method when simulating 3-dimensional numerical neutron noise problems and neutron noise experiments.

Fig. 1 .
Fig. 1.System configuration of the C4V benchmark problem with the source location labeled in red.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Fast (left) and thermal (right) static flux along the diagonal crossing the MOX fuel assemblies of the C4V system.

Fig. 3 .
Fig. 3. Fast (left) and thermal (right) static flux along the diagonal crossing the UO2 fuel assemblies of the C4V system.

Fig. 4 .
Fig. 4. Amplitude results in the fast (left) and thermal (right) group for the neutron noise problem based on the C4V configuration with a perturbation frequency of 1 Hz.

Fig. 5 .
Fig. 5. Phase results in the fast (left) and thermal (right) group for the neutron noise problem based on the C4V configuration with a perturbation frequency of 1 Hz.

Fig. 6 .
Fig. 6.Simulation of C4V neutron noise problem with different orders of discrete ordinates; fast noise amplitude (top-left) and phase (top-right), thermal noise amplitude (bottom-left) and phase (bottom-right).

Fig. 7 .
Fig. 7. Simulation of C4V neutron noise problem with the S8 approximation and different strengths of the fictitious source; fast noise amplitude (top-left) and phase (topright), thermal noise amplitude (bottom-left) and phase (bottom-right).

Fig. 8 .
Fig.8.C5G7 system configuration with the location of the noise source.Fig.9.Discretization of the fuel, guide tube, or fission chamber cell; fine mesh for the transport sweeps.

Fig. 10 .
Fig. 10.Static flux with different orders of discrete ordinates along the diagonal that crosses the MOX fuel assemblies in the C5G7 system; first energy group (left) and seventh energy group (right).

Fig. 11 .
Fig. 11.Neutron noise amplitude with different orders of discrete ordinates along the diagonal that crosses the MOX fuel assemblies in the C5G7 system; first energy group (left) and seventh energy group (right).

Fig. 12 .
Fig. 12. Noise phase calculated with different orders of discrete ordinates along the diagonal that crosses the MOX fuel assemblies in the C5G7 system; first energy group (left) and seventh energy group (right).

Fig. 13 .
Fig. 13.Neutron noise calculated with the S16 approximation in the C5G7 system; first-energy group amplitude (top-left) and phase (top-right) and seventh-energy group amplitude (bottom-left) and phase (bottom-right).

Fig. 14 .
Fig. 14.Neutron noise calculated with the S64 approximation in the C5G7 system; first energy group amplitude (top-left) and phase (top-right) and seventh energy group amplitude (bottom-left) and phase (bottom-right).

Fig. 15 .
Fig. 15.Effect of the strength cof the fictitious source on the S16 calculations; real (left) and imaginary (right) part of the neutron noise for the first energy group close to the right boundary of the system.

Fig. 16 .
Fig. 16.Neutron noise calculated using the S16 approximation and a fictitious source of strength c ¼ 0:9; first-energy group amplitude (top-left) and phase (top-right) and seventh-energy group amplitude (bottom-left) and phase (bottom-right).

Fig. 17 .
Fig. 17.Convergence of the real (left) and imaginary (right) part of the scalar neutron noise for the C5G7 neutron noise calculations.

Table 1
Number of iterations and wall clock time required for convergence in the C4V neutron noise problem.

Table 2
Effective multiplication factor, number of iterations and wall clock time in the C5G7 static calculations.

Table 3
Number of iterations and wall clock time in the C5G7 neutron noise calculations.