Verification of the code DYN3D for calculations of neutron flux fluctuations

Insufﬁciently explained magnitudes and patterns of ﬂux ﬂuctuation observed mainly in KWU PWRs are recently investigated by various European institutions. Among the numerical tools used to investigate the neutron ﬂux ﬂuctuations is the time-domain reactor dynamics code DYN3D . As DYN3D and comparable codes have not been developed with the primary intention to simulate low-amplitude neutron ﬂux ﬂuc-tuations, their applicability in this ﬁeld has to be veriﬁed. In order to contribute to the veriﬁcation of DYN3D for the simulation of neutron ﬂux ﬂuctuations, two special cases of perturbations of the neutron ﬂux (a localized absorber of variable/oscillatory strength and a travelling oscillatory perturbation) are considered with DYN3D on the one hand and with the frequency-domain neutron noise tool CORE SIM as well as analytical frequency-domain approaches, respectively, on the other hand. The obtained results are compared with respect to the distributions of the amplitude and the phase of the induced neutron ﬂux ﬂuctuations. The comparisons are repeated with varied amplitudes and frequencies of the perturbation. The results agree well both qualitatively and quantitatively for each of the conducted calculations. The remaining deviations between the DYN3D results and the reference results exhibit a dependence on the perturbation magnitude, which is attributed to the neglect of higher-order terms (linear theory) of the perturbed quantities in the calculation of the reference solutions.


Introduction
Investigation of neutron flux fluctuations has attracted new attention since numerous pressurized water reactors have exhibited insufficiently interpreted cycle-by-cycle increases or decreases of the magnitudes of these fluctuations.The phenomenon has become relevant especially in KWU (Kraftwerk Union) built reactors (Seidl et al., 2015), as their power limitation system performs automatic power reduction measures if given criteria with regard to high neutron flux fluctuations are met.Besides the changes of the amplitudes, neutron flux fluctuations of KWU PWRs exhibit spatial correlations that are still not entirely understood.
In the process of enhancing the methods of neutron noise analysis and as one measure to improve the understanding of the mentioned phenomena, neutron flux fluctuations have been simulated with dedicated deterministic tools such as CORE SIM (Demazière, 2011a) (e. g.Hoang et al., 2016), CORE SIM+ (Mylonakis et al., 2021b) (e. g.Mylonakis et al., 2021a), and the neutron noise solvers in APOLLO3 (Rouchon et al., 2017a;Rouchon et al., 2020), which operate in the frequency domain, and as FEMFFUSION (Vidal-Ferràndiz et al., 2020), which operates in the time domain and frequency domain, and also with reactor dynamics codes such as DYN3D (Rohde et al., 2016;Kliem et al., 2016) (e. g.Rohde et al., 2018;Viebach et al., 2020) or S3K (Grandi, 2011) (e. g.Verma et al., 2021), which also operate in the time-domain.As the latter have not been developed for the analysis of neutron noise, their applicability in this field is also question of recent research.From the physical point of view, those codes are applicable.However, numerical issues might be present, i. e. the setting of the numerical parameters of the codes, probably chosen for the transient simulation of incident scenarios, may require a readjustment.
One approach for verification of a time-domain code for neutron noise calculations is to compare its results with those of suitable frequency-domain calculations for certain benchmark cases.
Such analyses have been presented, e. g. by Olmo-Juan et al. (2019) and also by Vidal-Ferràndiz et al. (2020).In this context, the matter of the paper at hand is a comparison of results of the code DYN3D with those of analytical considerations and those of the tool CORE SIM (see the work by Viebach et al. (2019b) for preliminary work of this study).
As DYN3D results are compared with results of analytical calculations and with results of CORE SIM, for which verification cases exist, the comparison addresses a validation of DYN3D for simulations of neutron-flux fluctuations.The foundation of each of the compared approaches is diffusion theory, such that the comparison can focus on probing the applicability of time-domain codes in the given context.For a comparison of diffusion theory with transport theory in this context, it is referred to the works of Bahrami and Vosoughi (2018) and Gong et al. (2021).For Monte-Carlo based calculations of neutron noise, it is referred to the contributions of Yamamoto (2013), Rouchon et al. (2017b), and Demazière et al. (2020).It also has to be mentioned that thermal hydraulics and its feedback on the neutron kinetics is beyond the scope of the paper at hand.
For the comparisons, three different reactor models are used: a homogeneous one-dimensional fine-mesh reactor, the same with a coarse mesh, and a heterogeneous three-dimensional reactor.The models are prepared both for DYN3D calculations and for reference calculations.The reference calculations are performed with an analytical approach for the one-dimensional models and with CORE SIM for the three dimensional model.Sinusoidal perturbations are applied at fixed frequencies and fixed magnitudes as a localized absorber of variable strength or a travelling perturbation, respectively, to each of the models.Each perturbation case is calculated twice: first within DYN3D and second with the analytical approach or within CORE SIM, respectively.The induced neutron flux fluctuations are considered in terms of the spatial distributions of their amplitude and their phase.
The paper is structured as follows.After a brief description of the used codes, the considered reactor models are introduced.Afterwards, it is shown how the systems are perturbed and how the specific perturbation cases are realized.After some remarks on the preparation of the analytical solutions and the postprocessing of the DYN3D results, the comparisons are shown each for five different sets of perturbation frequency and magnitude.A discussion closes the main part of the paper.The paper finishes by drawing conclusions.

Description of the used codes
The used version of the code DYN3D solves the two-group diffusion equation in the time domain using a nodal expansion method with transverse integration for the spatial dependence and taking into account six groups of delayed neutron precursors.The used version of the code CORE SIM solves the linearized two-group diffusion equation in the frequency domain using finite-differences methods (mesh centered) for the spatial dependence and taking into account one group of delayed neutron precursors.The considered analytical solutions of the one-dimensional problems are evaluated using Python scripts.Thus, DYN3D is compared with a code or scripts, respectively, with methodically different approaches for the treatment of the time and space dependency.This enables a more generic probing of the DYN3D results than comparing them against a time-domain nodal methods code.The comparison of a non-linearized method (DYN3D) and a linearized method (CORE SIM; analytical approach, see Section 2.2.6) is only possible as the applied magnitudes of the perturbations are chosen sufficiently small (see also Section 2.2.4) such that higher-order terms (i.e., in this paper, dR a;2 Á d/ 2 ) are small compared to the linear terms.

Description of the models
Two one-dimensional models and one three-dimensional model are considered, all based on the OECD/NEA and U.S. NRC PWR MOX/UO2 Core Transient Benchmark (Kozlowski and Downar, 2003).The spatial layouts of the models are shown in Fig. 1.In order to comply with the restrictions of the considered analytical solution and the CORE SIM solution, only one group of delayed neutron precursors is considered (using b ¼ P 6 i¼1 b i and b=k ¼ P 6 i¼1 b i =k i for the condensation of the fractions b i ; i ¼ 1; 2; . . .; 6 of delayed neutron precursors to one fraction b and of the decay constants k i ; i ¼ 1; 2; . . .; 6 to one constant k).For the same reason, upscattering, selfscattering, and assembly discontinuity factors are disregarded.For the transient part of the DYN3D calculation, the thermal-hydraulics feedback is switched off.In the mentioned preliminary study (Viebach et al., 2019b), the physical models in DYN3D were not subject to such simplifications.

One-dimensional models
The one-dimensional models have homogeneous two-group constants.These are derived from the heterogeneous threedimensional model (see Section 2.2.2) by homogenization and numeral simplification.The spatial setup corresponds to the axial dimension of the three-dimensional case.For the coarse-mesh model (see Fig. 1a), the partition into 19 nodes of the threedimensional model is kept.For the fine-mesh model (see Fig. 1b), the number of nodes was increased to the round number of 500, increasing the number of nodes by one order of magnitude, maintaining that, except for the system boundaries, no node boundaries coincide with those of the coarse-mesh model, and keeping a finite number of digits for the node length.The two models are identical except for nodal partitioning.Zero-flux boundary conditions are applied.The models are set to critical by adjusting the thermal fission cross section mR f;2 .The used set of cross sections is given in Table 1. 1

Three-dimensional model
The three-dimensional model in DYN3D was prepared in an earlier work by Beckert and Grundmann (2008).The spatial layout is shown in Fig. 1c.The wide-range two-group macroscopic crosssection library was generated using the lattice code HELIOS (Wemple et al., 2008).Vacuum boundary conditions are applied.
The model in CORE SIM is derived from the steady-state solution of the DYN3D model.The converged nodal distributions of the group constants are directly used for CORE SIM.The kinetic parameters, which are global entities in CORE SIM, were homogenized before the DYN3D steady-state calculation.In order to increase the spatial accuracy of the CORE SIM calculations, the mesh given by the homogeneous regions (i.e. the DYN3D nodes) is refined by a factor of 7 in each dimension, giving 7 Â 7 Â 7 subregions per node region. 2oth DYN3D and CORE SIM perform steady-state calculations before calculating the neutron fluctuations in the transient calcula-tion or the noise calculation, respectively.The relations between the DYN3D and the CORE SIM runs and the corresponding flow of information are illustrated in Fig. 2.

Perturbation of the system
Neutron flux fluctuations are calculated in the time domain and in the frequency domain node indices and t and f denoting time and frequency, respectively.
In the one-dimensional case, the z-direction is considered (only i z exists).In the following, the equations are only written with respect to the three-dimensional case for the sake of brevity.The systems are perturbed by modification of the nodal cross sections.For the sake of simplicity, in the work at hand, only the thermal absorption cross section R abs;2 is modified: with R abs;2;0 being the steady-state value of the cross section, dR abs;2 being the time dependent perturbation, A being the amplitude, f 0 being the frequency, and a being the phase.The transformation of the perturbation to the frequency domain gives with i denoting the imaginary unit and d . . .ð Þ denoting the Diracdelta function.The amplitude is expressed in terms of a relative amplitude A rel ,

Considered perturbation cases
Two different perturbation cases f dR abs;2 are considered: A localized absorber of variable strength (see Section 2.2.4.1) in the central node and a vertically travelling perturbation (see Section 2.2.4.2) in the central channel.Both cases are simulated for each of the three reactor models.And for both cases, the perturbation magnitude and the perturbation frequency are varied giving the five perturbation sets that are shown in Table 2.The entry values, especially of the second column and the fourth column, are motivated by the relevant ranges of measured neutron flux fluctuations in KWU PWRs recently investigated (see Section 1), which are around 1 Hz and some percents, respectively.
For the travelling perturbation, the cross sections of multiple nodes are perturbed, giving an integral perturbation magnitude that is large compared to that of a single-node perturbation.Therefore, as large-amplitude perturbations conflict with the linearization of the neutron kinetics that is assumed in the analytical solution (see Section 2.2.6) and in the CORE SIM methodology (see Section 2.1), the respective magnitude A rel of the travelling perturbation is chosen smaller by two orders compared to the ones Fig. 1.Layout of the considered reactor models.The one-dimensional models (Fig. 1a and Fig. 1b) are used for DYN3D and the analytical calculations.The three-dimensional model (Fig. 1c) is used for DYN3D and for CORE SIM calculations.Reflector areas are represented in dark gray (only Fig. 1c).Regions considered for the perturbations are represented in shades of light gray.

Table 1
Group constants used in the one-dimensional simulations.
In the following, the construction of the two perturbation cases by the distributions A i x ; i y ; i z À Á and a i x ; i y ; i z À Á (see Eq. ( 2)) is described in detail.
2.2.4.1.Absorber of variable strength at fixed location.The perturbation that represents a localized absorber of variable strength is given by with the triple i x;p ; i y;p ; i z;p À Á denoting the node where the perturbation is applied.A perturbation of the nodes in the center or next to the center of the core, respectively, is considered, i. e.
Assuming a perturbation travelling velocity v and an equal distance Dz between the center of axially adjacent nodes, the time shift at layer i z;p relative to the lowermost position i With specifying the perturbation frequency f 0 and choosing the phase of the lowermost position as a i x;p ; i y;p ; 1 , the time shifts of the axial layers can be expressed by the perturbation phases as Thus, a perturbation travelling vertically up the channel defined by i x;p ; i y;p À Á is described here by   For the three-dimensional calculations, the nodes of the bottom and top reflectors are excluded from the explicit perturbation.

Extraction of the results
The results were prepared such that amplitudes and phases can be compared.For the frequency domain, the fluctuations f d/ 1 ; f d/ 2 are available through the CORE SIM output, i. e. absolute value and phase of complex numbers given on the fine-mesh domain.In order to provide a spatial resolution compatible with that of DYN3D, the CORE SIM output is coarsened to the original nodal domain by applying the arithmetic mean over the respective sub-regions.The coarsened distributions of absolute values are normalized to the coarsened steady-state distributions / 1;0 ; / 2;0 À Á to provide the distribution of relative fluctuation amplitudes.
For the time domain, evaluation of amplitudes and phases was performed via the following procedure.After the steady-state calculation, the transient calculations were performed with a temporal length of a bit more than 2 Á 1=f 0 such that each simulation covers two full periods.By comparison of the simulated values of the neutron flux of the first period and those of the second period, the drift of the local mean values was quantified.For each position (node), the drift is approximated as a linear function and then used to detrend the neutron flux values.Afterwards the detrended timeseries of the local neutron flux were interpolated by cubic splines.
The minimum value min / 1=2 and the maximum value max / 1=2 of the fast and the thermal flux, respectively, evaluated for the first of the two simulated periods are used to determine the relative fluctuation amplitude: The phase a / 1=2 of the fluctuation is determined by finding the instant t 0 when the fluctuating neutron signal passes the average of the minimum and the maximum, i. e.
for the second time.It has to be emphasized that two different relative units are defined: % amp is the unit for the relative fluctuation amplitudes d/ 1=2 =/ 1=2;0 , and % dev is the unit for the relative deviation of the amplitudes of the DYN3D solutions relative to the corresponding analytical or CORE SIM solution.
It should also be emphasized that, different from the frequencydomain approach, the approach applied for calculating and extracting the time-domain results only approximates stationary neutron oscillations.Sufficiently long time series have to be simulated with DYN3D, if also the oscillation of the concentrations of the delayed neutron precursors is required to be settled.However, in this work, only a few oscillation periods are modelled in DYN3D.This might result in non-stationary conditions with respect to most notably the delayed neutrons.For future studies of stationary oscillations with DYN3D, the calculation of long time series (in the range of 100 s) is planned.

Preparation of the analytical solution of the one-dimensional case
For the perturbations in the one-dimensional models (see Section 2.2.1), solutions to the problems can be found using the Green's function approach as presented by Pázsit and Demazière (2010).According to the given context, the mathematical problems are the one-dimensional time-dependent two-group diffusion equation with one group of delayed neutron precursors, homogeneous group constants, and added perturbations as given in Section 2.2.4.After setting higher-order terms to zero (linear theory), i. e. dR a;2 Á d/ 2 ¼ 0, and performing a Fourier transform of the remaining time-dependencies, the equations to be solved read: and all other quantities carrying their usual meaning.Here, only the perturbation of the cross section of the thermal absorption is considered.In order to enable solving the problem by the Green's function approach, in Eq. ( 16), the inhomogeneity is replaced by the Dirac-delta function d z À z 0 ð Þ and the neutron flux fluctuation is given by the Green's function The components G 1 and G 2 of the Green's function are composed of sine and sinh functions (see the works of Pázsit and Demazière (2010) and Demazière (2011b)).The fluctuation of the neutron flux is given by the convolution of the inhomogeneity of Eq. ( 16) with the Green's function: In order to comply with the nodal perturbations defined in Eqs. ( 1)-( 6), the perturbation of the cross section of the thermal absorption dR a;2 z 0 ; x ð Þ is given by means of a rectangle function, i. e. in case of the localized absorber of variable strength: with z l i z;p À Á and z u i z;p À Á denoting the lower and upper boundary of the perturbed node i z;p , respectively.For the case of the travelling perturbation (see Eq. ( 11) and Eq. ( 12)), the expression given in Eq. ( 21) has to be formulated with respect to the amplitude and the phase of the perturbation for all nodes of the perturbed channel.
The neutron flux fluctuations (Eq.( 20)) are calculated on a spatial mesh that was refined with respect to the node mesh by a factor of 100 for fine mesh (see Fig. 1b) and factor of 500 for the coarse mesh (see Fig. 1a).The nontrivial integration domain z 0 2 z l i z;p À Á ; z u i z;p À Á Â Ã (see Eq. ( 21)) is refined by a factor 181 for the fine mesh and 581 for the coarse mesh.The integration in Eq. ( 20) is carried out by simple rectangle rule.Finally, the neutron flux fluctuations given on the respective refined meshes are averaged over the respective node regions to provide quantities that can be compared with the corresponding DYN3D results.

Results
In the DYN3D simulations, it turned out that neutron-kinetics time step of Dt NK ¼ 10 À5 s suited the demands of the considered cases for each of the three models.In each simulation, the perturbation was updated with a time step of Dt TH ¼ 10 À4 s.In order to optimize the initial conditions of the transient calculations, the precision for the fission source iteration in steady-state calculation was set to f;steady ¼ 10 À14 (default: 2 Á 10 À6 ).In the transient calculation, tightening the precision for the fission source iteration to f;transient ¼ 10 À8 (default: 10 À5 ) was sufficient.
Each of the DYN3D runs was performed on one core of an Intel (R) Core(TM) i5-6440HQ CPU @ 2:20 GHz.For two seconds of simulated neutron fluctuations, computing times ranged between about one minute for the 1D coarse mesh cases, about 15 min for the 1D fine mesh cases, and about two and a half hours for the 3D cases.Each of the CORE SIM runs was performed on two cores of an Intel(R) Xeon(R) CPU E7-4850 v3 @ 2:20 GHz.For two seconds of simulated neutron fluctuations, computing times amounted to around two and a half days.
In the following, the results are shown for all considered cases.The solutions are presented only for perturbation set A, while the deviations are shown for all perturbation sets in order to assess the change of the deviations due to parameter variation.

One-dimensional models
Figs. 3 and 4 show the comparison of the DYN3D calculation and the analytical solution of the problems for the fine onedimensional mesh.Figs. 5 and 6 show the corresponding results for the coarse one-dimensional mesh.
With respect to the amplitude distributions, the results of parameter set A obtained with DYN3D and the analytical calculation agree with each other (see the left parts of Figs.3a, 4a, 5a,  and 6a).A closer inspection of the deviations reveals maximum Fig. 3. Neutron flux fluctuations for the fast ( f ) and thermal ( t ) energy group simulated with DYN3D ( d3d ) and given by an analytical solution ( an ).Localized absorber of variable strength in 1D on a fine mesh, shown for perturbation set A (f ¼ 1 Hz; A rel ¼ 1%).The deviations between both solutions are given for all considered perturbation sets (A-E).Except for the region around the perturbation in Fig. 3a, the f-and the t-curves of each set overlap everywhere.In Fig. 3b, the f-and the t-curves of each set overlap everywhere and the curves of set A are partly covered by those of set D (an arrow points to the hidden curves of set A). Fig. 4. Neutron flux fluctuations for the fast ( f ) and thermal ( t ) energy group simulated with DYN3D ( d3d ) and given by an analytical solution ( an ).Travelling perturbation in 1D on a fine mesh, shown for perturbation set A (f ¼ 1 Hz; A rel ¼ 0:01%).The deviations between both solutions are given for all considered perturbation sets (A-E).The set E* is equal to the set E but with f;transient ¼ 10 À10 .Except for the set E*, the f-and the t-curves of the deviations of each set overlap everywhere.deviations less than one percent.For each considered case, the deviations have their minimum values in the center of the reactor.The deviations exhibit a quantitative dependence on the perturbation amplitude (sets A vs. B vs. C) and on the perturbation frequency (sets A vs. D vs. E).However, only in case of the travelling perturbation, set E develops a qualitative change of the distribution of the amplitude deviation.The effect was mended by a readjustment of precision of the fission source iteration to f;transient ¼ 10 À10 as shown by the lines of parameter set E*. Systematic effects are hard to find for the considered variations.Nevertheless, for the local absorber (see Figs. 3a and 5a), a systematic effect is an increase of the deviations for increasing the perturbation magnitude.
Also with respect to the phase distributions, the results of parameter set A obtained with DYN3D and the analytical calculation agree with each other (see the left parts of Figs.3a, 4b, and 6b), except for the localized absorber simulated for the coarse mesh, where differences are visible (see the left part of Fig. 5b).With regard to the parameter variations, the phase exhibits clearer dependences.Smaller perturbation amplitudes and higher frequencies lead to better agreements of the phases.

Three-dimensional model
In order to give an overview about the solution of threedimensional model, Fig. 7 shows the steady-state neutron flux calculated with DYN3D and the deviation against the steady-state neutron flux calculated with CORE SIM for the mid axial layer (i z ¼ 10).The solutions qualitatively agree.For the fuel region, quantitative differences are in the range of some percents.However, the reflector region exhibits deviations up to about 17%.The respective multiplication factors are given in Table 3 accounting for a difference of about 53:2 pcm.Fig. 5. Neutron flux fluctuations for the fast ( f ) and thermal ( t ) energy group simulated with DYN3D ( d3d ) and given by an analytical solution ( an ).Localized absorber of variable strength in 1D on a coarse mesh, shown for perturbation set A (f ¼ 1 Hz; A rel ¼ 1%).The deviations between both solutions are given for all considered perturbation sets (A-E).In Fig. 5a, except for the region around the perturbation, the f-and the t-curves of each set overlap everywhere.In Fig. 5b, the f-and the t-curves of each set overlap everywhere.Both in Fig. 5a and in Fig. 5b, the curves of the deviations of set A and set D are partly covered by those of set E (arrows point to the hidden curves of set A and set D). Fig. 6.Neutron flux fluctuations for the fast ( f ) and thermal ( t ) energy group simulated with DYN3D ( d3d ) and given by an analytical solution ( an ).Travelling perturbation in 1D on a coarse mesh, shown for perturbation set A (f ¼ 1 Hz; A rel ¼ 0:01%).The deviations between both solutions are given for all considered perturbation sets (A-E).The set E* is equal to the set E but with f;transient ¼ 10 À10 .In Fig. 6a, except for set E and set E*, the f-and the t-curves of the deviations overlap everywhere.In Fig. 6b, the f-and tcurves of all sets overlap everywhere and the curves of set B and set E are partly covered by those of set E* (arrows point to the hidden curves of set B and set E).
Fig. 8 presents the comparison between DYN3D and CORE SIM for the localized absorber of variable strength, and Fig. 9 presents it for the travelling perturbation.The results are compared along the central channel, where the perturbations are applied, and along the mid fuel assembly row at mid height of the reactor.They qualitatively agree.
On a quantitative base, the amplitudes differ in the range of percents with maximum deviations of about 6% dev , which is found at the perturbed node in the case of the localized absorber (see Fig. 9 and in the axial reflector in the case of the travelling perturbation (Fig. 9a).
The deviations of the amplitudes show little dependence on the varied parameters.The corresponding lines of the localized absorber (Fig. 8a and Fig. 8c) coincide.For the travelling perturbation, the deviations are slightly dependent on the perturbation frequency.An increase of the frequency leads to slightly increased deviations.
The deviations of the phases are in the range of hundredth of rad.They exhibit a marginal dependence on the perturbation fre-quency and a stronger dependence on the perturbation amplitude.A decrease of the perturbation amplitude leads to a decrease of the deviation in the solutions obtained with DYN3D and CORE SIM.

Discussion
For each considered reactor model, for each perturbation case, and for each perturbation set, the DYN3D results showed qualitative agreement with the respective reference results and quantitative agreement with maximum deviations of the induced fluctuation magnitudes ranging between about 0:005% dev and about 6% dev .As indicated by the perturbation set E* for the onedimensional models (Fig. 4 and Fig. 6), the numerics parameters play a crucial role for the correctness of the results.It also has to be emphasized that the postprocessing necessary for the timedependent simulations is a root for certain ambiguity: Since DYN3D takes into account the arising nonlinear terms (dR abs;2 Á d/ 2 ), the calculated neutron flux does not strictly oscillate as a sine. 3The distortion, i. e. deviation from the sine curve, is the more pronounced, the higher the perturbation amplitude was chosen.For a sufficiently strong perturbation, the system diverges.This  was observed, e. g., with A rel ¼ 2% for the travelling perturbation with the coarse-mesh one-dimensional reactor during preliminary calculations.The chosen means to extract the amplitude and the phase from the DYN3D results (see Section 2.2.5) gave the most consistent for the considered simulations.
The sensitivity on the perturbation amplitude can be seen in Fig. 3a and Fig. 5a, i. e. for the amplitudes of the neutron flux for the one-dimensional simulations of the localized absorber.In contrast, for the corresponding three-dimensional simulations, the amplitudes of the neutron flux are insensitive to variations of the perturbation amplitude and also of the perturbation frequency.However, the dependence is hidden by the comparatively large deviations between the DYN3D and CORE SIM results due to the differing methodology for the treatment of the spatial dependence in the neutron diffusion equation (i.e. nodal vs. finite differences).As indicated in Section 2.2.2, a mesh refinement of the CORE SIM model results to a reduction of the deviations by several percents.The results obtained here may be compared with those of Viebach et al. (2019b).There, a refinement factor 2 is used, and the maximum deviations amount to about 15% dev .
Elaborating on the effects of variations of the perturbation parameters on the phases, a low perturbation amplitude results to small deviations, and, except for the three-dimensional simulation of the travelling perturbation, an increase of the perturbation frequency also leads to a reduction of the deviations.As described above, the higher order terms of the diffusion equation taken into account by DYN3D distort the sine-shape expected for the timedependence of the simulated neutron flux fluctuations.Therefore, the deviations between the DYN3D solutions and the respective reference solutions from the analytical considerations and CORE SIM with respect to the phase are the larger, the stronger the applied perturbation is.In contrast to the deviations of the amplitudes, the deviation of the phase is sensitive also to variations of the perturbation frequency.This effect may occur because the determination of the phase is more prone to influences of gradually non-sinusoid time dependencies.With a low perturbation frequency, the time of the perturbation to act with unchanged sign on the neutron flux is longer than for a high frequency.Therefore, the higher-order terms have more time to develop their effects on the time-dependence of the neutron flux.
The investigations shown here, demonstrate that the effects of large-scale as well as finely localized perturbations of the material can be well simulated with DYN3D.The considered range of the perturbation amplitude can be extended to even smaller ampli-Fig.8. Neutron flux fluctuations for the fast ( f ) and thermal ( t ) energy group simulated with DYN3D ( d3d ) and simulated with CORE SIM ( cs ).Localized absorber of variable strength in 3D, shown for perturbation set A (f ¼ 1 Hz; A rel ¼ 1%).The deviations between both solutions are given for all considered perturbation sets (A-E).Except for the region around the perturbation, the f-and the t-curves of all sets overlap.Both in Fig. 8a and Fig. 8c, the curves of all sets overlap.Both in Fig. 8b and Fig. 8d, the curves of set A are partly covered by those of set D (arrows point to the hidden curves of set A). tudes, promising even better agreement with the linearized tools.For extending the considered range of the perturbation frequency, the perturbation magnitude applied in the study has to be paid attention for.As described above, for low perturbation frequencies, the perturbation has more time to act, leading to more pronounced deviations of the time dependence of the neutron flux from the ideal sine-shape form.For increasing perturbation frequencies, the perturbation amplitude is less relevant.However, the timestep in DYN3D needs to be refined.In this context, the stop criteria for the iterations in DYN3D (e. g. f;transient )) have to be mentioned, since they have to be more strict for small-amplitude perturbations.In general, when applying a time-domain reactor dynamics code to neutron noise, preparatory analyses as shown here should be conducted that include the intended ranges of the frequency and amplitude of the perturbation.
As already mentioned, reactor dynamics codes like DYN3D have not been primarily developed for the analysis of small-amplitude neutron flux-fluctuations.However, the analysis at hand shows that, with some readjustments of the numerical parameters, these codes may also be applied to questions related to neutron noise.In fact, the application of such codes brings some advantages.As discussed, operation in time-domain does not require the neglect of higher-order terms, which is commonly utilized for the transfor-mation into frequency space.Thus, if the considered perturbations are not of sufficiently small amplitude, time-domain codes may still correctly simulate the response of the neutron flux, while linearized tools fail.In particular, for the time-domain codes, the question whether higher order contributions are negligible or not has only minor relevance, as these codes take them into account anyway.Furthermore, time-domain codes can handle more than one frequency at once.Nevertheless, it should be emphasized that frequency-domain approaches are faster, as they directly give the fluctuations at the perturbation frequency.The calculation of the neutron noise at the perturbation frequency is equivalent to a source calculation that does not require any time discretization.Another advantage of the reactor dynamics codes is their sophisticated treatment of multiphysics aspects of a nuclear reactor, mainly thermal hydraulics and heat conduction, as feedback may play a major role for neutron flux fluctuations.However, feedback is not included in the investigations shown here.A verification of relevant feedback mechanisms in the context of small-amplitude neutron-flux fluctuations give rise to another paper and is matter of future work.It should be mentioned that not only the various relevant physics fields as neutron kinetics or thermal hydraulics but also their mutual coupling has to be subject of the verification process.
Fig. 9. Neutron fluctuations for the fast ( f ) and thermal ( t ) energy group simulated with DYN3D ( d3d ) and simulated with CORE SIM ( cs ).Travelling perturbation in 3D, shown for perturbation set A (f ¼ 1 Hz; A rel ¼ 0:01%).The deviations between both solutions are given for all considered perturbation sets (A-E).The curves of the deviations of set A and set B overlap with those of set C. In Fig. 9c and in Fig. 9d, except for the region around the travelling perturbation, the f-and the t-curves overlap.

Conclusions
Two special cases of perturbations of neutron kinetics have been simulated with the time-domain reactor dynamics code DYN3D, i. e. a localized absorber of variable strength and a travelling perturbation.For both cases, the magnitude and the frequency of the perturbation were varied.Aiming on a verification of DYN3D for the simulation of neutron flux fluctuations, the results obtained for a one-dimensional fine-mesh homogeneous system, a onedimensional coarse mesh homogeneous system, and a threedimensional heterogeneous system were compared with respective reference solutions given by analytical calculations and calculations with the dedicated neutron noise tool CORE SIM.The deviations range between hundredth of percents for low perturbation amplitudes and percents for larger perturbation amplitudes.The observed sensitivity of the deviations to the magnitude of the perturbations is attributed to the neglect of the higher order terms of the neutron diffusion equation in the reference solutions.Thus, it is shown that DYN3D is well suited for the simulation of neutron flux fluctuations.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Sequence of calculations performed to obtain neutron flux fluctuations with DYN3D and with CORE SIM for the three-dimensional model.

Fig. 7 .
Fig. 7. Steady-state neutron flux normalized to the maximum flux value shown for the DYN3D calculation at mid axial layer (10 of 19) along with the deviation against the corresponding CORE SIM results.The upper and lower panel show the results for the fast group and the thermal group, respectively.

Table 2
Considered parameter sets for the applied perturbations by absorbers of variable strength and travelling perturbations.

Table 3
Comparison of k eff of the three-dimensional models.