Monte Carlo-based dynamic calculations of stationary perturbations

Capitalizing on some earlier work, this paper presents a novel Monte Carlo-based approach that allows estimating the neutron noise induced by stationary perturbations of macroscopic cross-sections in the frequency domain. This method relies on the prior computation using Monte Carlo of modified Green’s functions associated to the real part of the dynamic macroscopic cross-sections, mimicking equivalent subcritical problems driven by external neutron sources. Once such modified Green’s functions are estimated, the neutron noise induced by any type of perturbations can be recovered, by solving a linear algebra problem accounting for the interdependence between the real and imaginary parts of the governing balance equations. The newly derived method was demonstrated on a large homogeneous test system and on a small heterogeneous test system to provide results comparable to a diffusion-based solver specifically developed for neutron noise applications. The new method requires the specification by the user of the real part of the Fourier transform of the macroscopic cross-sections. This is accomplished using ACE-formatted cross-section files defined by the user. Beyond this input data preparation, no change to the Monte Carlo source code is necessary. This represents the main advantage of the proposed method as compared to similar efforts requiring extensive modifications to the Monte Carlo source code.


INTRODUCTION
Nuclear reactor core monitoring using non-intrusive techniques is becoming increasingly important, due to the ageing fleet of reactors and, correspondingly, operational problems being more frequently encountered.In this respect, noise diagnostics, which relies on the analysis of the deviations of the local neutron flux from its time-averaged value, represents one of the most powerful techniques.Nevertheless, most of the diagnostic tasks are based on unfolding from the detector readings the nature and possibly location of the driving perturbation or anomaly.This unfolding requires the prior determination of the reactor transfer function, i.e. the function giving the response of the system to localized perturbations.
The determination of the reactor transfer function can be carried out using either time-domain or frequency-domain approaches, diffusion or transport-based methods, and deterministic or probabilistic strategies.The focus of this paper is on frequency-domain Monte Carlo methods.Most of the past work in this respect relies on extensive modifications of the Monte Carlo source code [1][2][3][4].Some earlier work performed by one of the present authors demonstrated the possibility to perform Monte Carlo-based noise calculations without modifying the source code [5].The strategy developed in this earlier work was based EPJ Web of Conferences 247, 21003 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124721003 on mimicking the noise calculations with two static subcritical calculations: one for the real part of the balance equations and one for the imaginary part of the balance equations.
In this paper, the limitations encountered during the preliminary work are addressed, namely, the coupling between the real and imaginary parts of the balance equations and the handling of possible negative reaction rates in such balance equations.More specifically, the neutron noise equations in the frequency domain are, as before, re-arranged into two sets of coupled balance equations, one for the real part of the balance equation and one for the imaginary part of the balance equation.Nevertheless, the coupling between those balance equations is treated as external neutron noise sources, for which a modified set of Green's functions is correspondingly derived.By taking spatial integrals of the actual noise source and those coupling terms multiplied by the modified Green's functions, the neutron noise can be determined.The main advantage of the new framework lies with the fact that those modified Green's functions only involve real quantities, are equivalent to balance equations obtained for static subcritical systems driven by external neutron sources and can be, again, readily obtained from Monte Carlo without source code modifications.Furthermore, possible negative real/imaginary quantities only appear in the actual noise source and coupling terms and can thus be easily handled.Various test cases demonstrate the viability of this new procedure and are reported in this paper using the probabilistic code MCNP6.2.This paper is structured as follows.The derivation of the new method is first reported in Section 2. Section 3 presents two representative test cases: a large homogeneous system and a small heterogeneous system.In both cases, the Monte Carlo solutions are compared to the ones obtained using the CORE SIM tool [6], which is a diffusion-based deterministic solver specifically developed and verified for noise calculations [7].Some conclusions and recommendations for future work are then presented.

THEORY AND METHOD
For the sake of simplicity, the following derivation of the method is exemplified in the framework of diffusion theory and for two energy groups.The notations used are standard and are defined in, e.g., [5].

Governing balance equations for the neutron noise
The balance equations describing the effect of stationary fluctuations of the macroscopic cross-sections on the neutron flux in the frequency-domain are obtained in the following manner.The time-dependent conservation equations are first considered, and all time-dependent terms are written as sums between their steady-state values and their deviations from steady-state.Neglecting second-order terms, removing from the balance equations the ones corresponding to the steady-state conditions of the system and performing a temporal Fourier transfer lead, in a two-group formulation, to balance equations valid for describing the first-order effect of the stationary fluctuations of macroscopic cross-sections.Such equations are given as: , 1 The solution to the noise problem can also be found using the Green's function technique, according to which the induced neutron noise is given by integrals between the noise source as: The Green's function is obtained when the right-hand side of Eqs. ( 1) or ( 2) is replaced by thus gives the neutron noise induced in the spatial point r and energy group g due to a Dirac-like perturbation in the spatial point a r and energy group g a .

Previously developed Monte Carlo-based method
Eqs. (1) or ( 2) cannot be directly solved in existing codes that do not handle complex arithmetic as formulated, since they involve complex quantities as the result of the Fourier transform.This is particularly true for Monte Carlo codes.Nevertheless, splitting each complex quantity into a real part (denoted with the superscript r ) and an imaginary part (denoted with the superscript i ), four coupled balance equations are obtained: EPJ Web of Conferences 247, 21003 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124721003 Noticing that each of the above equations are equivalent to subcritical systems driven by external neutron sources, the four equations above were proposed in [5] to be solved in an iterative manner, as described hereafter.Starting with an assumed imaginary part identically equal to zero for both the fast induced neutron noise (12) and ( 13), which thus in addition to the driving noise sources also contain an artificial noise source spatially distributed over the entire system geometry.Eqs. ( 12) and ( 13) still correspond to inhomogeneous equations where the right-hand sides only contain the sum of the actual noise sources and of the artificial ones, and can thus be solved in a source-driven type of calculations.Once the imaginary parts of the fast induced neutron noise EG and of the thermal induced neutron noise EG are determined, they can be inserted back into the right-hand sides of Eqs. ( 10) and ( 11), and the process can be repeated until convergence.
In order to be able to perform such source-driven calculations in Monte Carlo as explained above, full control of the macroscopic cross-sections, as appearing in Eqs. ( 10)-(13), is required.Some Monte Carlo tools offer such a possibility via so-called ACE-formatted cross-section files (with ACE standing for A Compact Evaluated Nuclear Data File).In the reported work, MCNP6.2 was used [8], together with the method to generate group-wise diffusion-like cross-sections for MCNP6.2, as originally proposed in [9].The generated cross-sections were utilized to obtain the solution to the static problem, which was then used to generate new cross-section ACE-files in the dynamic problem and to define the driving neutron source.

New Monte Carlo-based method
As already pointed out in [5], solving subcritical-equivalent problems using Monte Carlo as highlighted above requires the reaction rates in such subcritical-equivalent problems to be positive.In the most general case, there is no guarantee that the real and imaginary parts of either the dynamic macroscopic crosssections or the induced neutron noise are positive, thus possibly leading to negative reaction rates.This problem was alleviated in [5] by only considering frequencies of 1 Hz, for which the imaginary part of the various terms and of the balance equations can be neglected.The novelty of the work reported hereafter resides in rewriting Eqs. ( 10)-(13) in a manner that leads to positive reaction rates when solving subcriticalequivalent systems.This is achieved by recasting those equations into: , , Eq. ( 14) can be solved using the Green's function technique by first replacing the right-hand side of Eq. ( 14) by a Dirac-like perturbation E a r r in either the thermal and the fast energy group, respectively.
The solution to Eq. ( 14), i.e. , 1,2 g g EG , corresponds in such a case to a Green's function i , , . The solution to the original noise problem, i.e. the solution to Eq. ( 14) with the actual right-hand side of this equation, is then given by integrals between the right-hand side of this equation, , and the Green's function i , , It is important to emphasize that the Green's function i  9) to Eq. ( 17) lies with the fact that the reaction rates appearing in the matrix A are all real and it could be demonstrated that those reaction rates all have the same sign as when solving a subcritical system.Determining the modified Green's function i , , can thus be carried out without difficulty using Monte Carlo techniques, as long as the macroscopic cross-sections can be specified using ACE-formatted cross-section files.
Because of the spatial variables r and a r appearing in the modified Green's function i , , , a spatial discretization of the system being modelled is required prior to computing such modified Green's functions.This means that the system has to be replaced by a set of homogenized regions using properly averaged macroscopic cross-sections.In essence, the method proposed in this paper is similar to the fission matrix method used in Monte Carlo simulations, where the production of neutrons in a given homogenized region due to neutrons born in any homogenized region is determined.After spatial discretization, Eqs. ( 17) and ( 18) can be recast into: where all quantities appearing in the equation above are generically written as vectors and matrices taking both the spatial and energy discretization of all balance equations into account.The implementation details can be found in [10].Eq. ( 19) can finally be rewritten as: Eq. ( 20) represents a linear system of equations that can easily be solved once the discretized modified Green's functions i X G have been determined from Monte Carlo runs.It also means that the coupling between the real and imaginary parts of the induced neutron noise is resolved externally to the Monte Carlo simulations.Likewise, the actual effect of the considered noise source is computed outside of the Monte Carlo runs.

VERIFICATION OF THE DEVELOPED METHOD
In the following, the proposed method is tested on two configurations: a large homogeneous system and a small heterogeneous system.In both cases, the solutions computed with the new method are compared to the ones determined using the deterministic code CORE SIM.

Homogeneous test case
A one-dimensional slab of thickness 201 cm is considered hereafter, with vacuum boundary conditions applied, as represented in Fig. 1.The static macroscopic cross-sections, point-kinetic data, and noise source definition are given in Table 1.1.52 10 cm.s v q 5 1 2 4.13 10 cm.s

Macroscopic cross-section data Point-kinetic data
Arbitrary scaling S of the noise sources Noise source located at 0 0 x 0.0135 Hz f and 2 f X Q The results of the Monte-Carlo based simulations are represented in Fig. 2, on which the results of the deterministic tool CORE SIM are superimposed.Since the calculations are performed in the frequency  domain, the amplitude and phase of the induced neutron noise both for the fast and energy groups are plotted.The spatial discretization used in the Monte Carlo-based method was set to 1 cm.As can be seen in this figure, the induced neutron noise agrees very well between the two solutions.The local component of the neutron noise is particularly well reproduced.Some small deviation in the prediction of the global component can nevertheless be noticed, where the "buckling" of the amplitude of the global component is slightly different between the two solutions and where the slope of the phase of the global component also slightly differs in the thermal group.The maximum relative difference in amplitude is ca.10%, whereas the maximum difference in phase is ca. 2 deg.Additional simulations (not reported here for the sake of brevity) demonstrated that the observed deviation in amplitude decreases when decreasing the mesh size for the computation of the discretized modified Green's functions using Monte Carlo.Furthermore, when replacing the Monte Carlo-based modified Green's function with a modified Green's function estimated from CORE SIM, the agreement between the direct CORE SIM simulations and the ones using the modified Green's function estimated with CORE SIM as well is excellent.This homogeneous test case demonstrates the correctness of the Monte Carlo-based method.It should be mentioned that the frequency of the noise source (equal to 0.0135 as reported in Table 1) was chosen so that both the real and imaginary parts of the induced neutron noise need to be considered.

Heterogeneous test case
A one-dimensional system made of 11 fuel pins, each pin being of 1 cm thickness surrounded by water, is considered hereafter, with vacuum boundary conditions applied, as represented in Fig. 3.The spacing between each fuel pin is 0.25 cm, whereas the spacing between the outer periphery of the outermost fuel pins and the boundary of the system is 0.125 cm.The static macroscopic cross-sections, point-kinetic data, and noise source definition are given in Table 2.
The results of the Monte-Carlo based simulations are represented in Fig. 4, on which the results of the deterministic tool CORE SIM are superimposed.The spatial discretization used in the Monte Carlo-based method was set to 0.025 cm.Both the amplitude and the phase of the reference CORE SIM solution are reproduced by the Monte Carlo-based solution, even if the transport solution has larger spatial variation on small scales, as compared to the diffusion-based CORE SIM solution.As for the homogeneous case, a slightly different slope in the phase of the induced thermal neutron noise is observed.It should nevertheless be pointed out that the variation of the phase throughout the system is less than 0.6 deg, and thus the observed differences are in essence negligible.Although transport effects are expected and could explain the small differences noticed between the transport and diffusion-based solutions, a benchmark between the Monte Carlo-based approach and a deterministic transport solver would be necessary to properly assess the accuracy of the Monte Carlo solution strategy.

Macroscopic cross-section data Fuel region
Moderator region Arbitrary scaling S of the noise sources Noise source located at 0 0 x 0.159 Hz f and 2 f X Q

CONCLUSIONS
In this paper, a method allowing the estimation, using Monte Carlo techniques, of the neutron noise induced by stationary fluctuations of the macroscopic cross-sections was developed.This method is based on the prior estimation using Monte Carlo of a modified Green's function, which is thereafter used for recovering the neutron noise.The method was demonstrated to give qualitatively correct results, using a diffusion-based solution as the reference solution.Since no source code modification is required in the proposed method, this work thus represents a new alternative to recent advancements made in Monte Carlo-based software specifically modified for noise calculations.The main advantage of avoiding source code modification with the present method lies with the possibility to thus capitalize on the extensive verification and validation of the Monte Carlo software used for such simulations.It should nevertheless be emphasized that the proposed technique requires spatially discretizing the system being modelled and correspondingly estimating such modified Green's function.Depending on the size of the system and its geometrical complexity, the estimation of the modified Green's functions might represent a very large task. -

1 iEG and thermal induced neutron noise 2 i 1 r
EG , Eqs.(10) and (11) correspond to inhomogeneous equations where the right-hand sides only contain the noise sources, and can thus be solved in a Monte Carlo source-driven type of calculation.Thereafter, the computed real parts of the fast induced neutron noise EG and of the thermal induced neutron noise 2 r EG can be inserted into the right-hand sides of Eqs.
complex.The main advantage of reformulating the solution from Eq. (

Fig. 1
Fig. 1 Representation of the homogeneous test case.

Fig. 3
Fig. 3 Representation of the heterogeneous test case.

Fig. 4
Fig. 4 Results for the heterogeneous test case.All results were normalized to the same induced neutron noise in the fast energy group at the location of the noise source.