A time and frequency domain analysis of the effect of vibrating fuel assemblies on the neutron noise

The mechanical vibrations of fuel assemblies have been shown to give rise to high levels of neutron noise, triggering in some circumstances the necessity to operate nuclear reactors at a reduced power level. This work analyses the effect in the neutron field of the oscillation of one single fuel assembly. Results show two different effects in the neutron field caused by the fuel assembly vibration. First, a global slow variation of the total reactor power due to a change in the criticality of the system. Second, an oscillation in the neutron flux in-phase with the assembly vibration. This second effect has a strong spatial dependence that can be used to localize the oscillating assembly. This paper shows a comparison between a time-domain and a frequency-domain analysis of the phenomena to calculate the spatial response of the neutron noise. Numerical results show a really close agreement between these two approaches.


Introduction
Being able to monitor the state of nuclear reactors while they are running at nominal conditions is a safety requirement. The early detection of anomalies gives the possibility to take proper actions before such problems lead to safety concerns or impact plant availability. The CORTEX project (Demazière et al., 2017b), funded by the European Commission in the Euratom 2016-2017 work program, aims at developing an innovative core monitoring technique that allows detecting anomalies in nuclear reactors, such as excessive vibrations of core internals, flow blockage, coolant inlet perturbations, etc. The technique is based on using the inherent fluctuations in neutron flux recorded by in-core and ex-core instrumentation, referred to as neutron noise, from which the anomalies will be detected and differentiated depending on their type, location and characteristics. The method is non-intrusive and does not require any external perturbation of the system.
In recent years, operational problems in some pre-Konvoi Pressurized Water Reactors (PWR) were encountered. More specifically, an abnormal increase of the neutron noise level was noticed. The vibration of fuel assemblies is suspected to be a possible cause of increased neutron noise due to the correlation existing between the stiffness of the installed fuel assemblies (FA) and the noise amplitudes (Seidl et al., 2015). The objective of this paper is twofold. First, the assessment of the capability of time and frequency domain codes to simulate the propagation in the neutron field of FA vibrations. Second, to evaluate the order of magnitude and spatial shape of the neutron noise behaviour when a single FA is vibrating. This contribution investigates the influence of lateral deflection of one FA by moving the whole homogeneous material back and forth following a sinusoidal oscillation.
Although the study of the noise induced by vibrating absorbers began in 1948 (Weinberg and Schweinler, 1948), applications to the diagnostics of internals vibration started in the 1970s in the framework of one-dimensional (Antonopoulos-Domis, 1976) and two-dimensional simplified models (Pázsit and Analytis, 1980;Pázsit and Glöckler, 1983). More recently, the localization of unseated fuel assemblies (Demazière, 2006) was also undertaken. Moreover, the properties of the neutron noise seen by ex-core detectors and induced by an FA vibration were investigated (Tran et al., 2015). Other studies considered FA vibrations by randomly modifying the size of the water gap which surrounds the FA of interest (Chionis et al., 2017). Also, the influence of coherent lateral deflections for all FAs on the neutron flux was considered in (Seidl et al., 2015) and (Viebach et al., 2018). In the last years, some research aimed at comparing time domain and frequency domain calculations was reported (Chionis et al., 2017;Viebach et al., 2019;Olmo-Juan et al., 2019).
Mechanical vibrations of FA have been studied from a structural point of view in (Park et al., 2003) and (Fry et al., 1984). Natural frequencies reported in literature range from 0.8 Hz to 24.5 Hz depending on the mode and the idealized form of bearing. The amplitude of the vibrations ranges can reach up to 1 mm. Fuel elements can be mechanically described in a first approach by a damped spring with hysteresis. There are possibilities for amplitudes larger than 1 mm in case of hypothetical synchronous motions and forced excitations. The usual FA vibration model cover motions in the parameter field from 0.5 to 25 Hz and amplitudes up to 5 mm in a sinusoidal movement. This work is based on FA vibration of 1 Hz and 1 mm of amplitude as an illustrative example.
This paper is organized as follows. First, the time domain analysis of the effect of the FA vibration is explained in Section 2. Then, Section 3 describes briefly the frequency domain analysis and how the vibrating fuel assembly is included in this context. After that, some numerical results are given in Section 4 to verify the proposed methodologies. Finally, Section 5 summarizes the main conclusions of the study.

Time domain analysis
For a given transient, the balance of neutrons inside a nuclear reactor core can be approximated using the time dependent neutron diffusion equation in the two energy groups approximation without up-scattering. This model is of the form of where K is the number of delayed neutron precursors groups considered and the matrices are defined as where φ 1 and φ 2 are the spatial and time dependent fast and thermal neutron fluxes, respectively. The diffusion coefficients and cross-sections, D g , Σ 12 , Σ ag , νΣ f g , g = 1, 2, appearing in the equations depend on the reactor materials and as a result they are position and time dependent functions. β k is the yield of delayed neutrons in the k-th precursors group and λ k is the corresponding decay constant. Both coefficients are related to the delayed neutron precursors. Effective delayed precursors coefficient are calculated as β eff = K k=1 β k and λ eff = β eff K k=1 The spatial discretization used is a high-order finite element (FEM) as the one presented in Vidal-Ferràndiz et al. (2014). Once the spatial discretization has been selected, a semi-discrete version of the time dependent neutron diffusion equation is solved. Since the system of ordinary differential equations resulting from the discretization of the neutron diffusion equations is, in general, stiff, implicit methods are necessary. Particularly, a first order backward method is used as the one in (Ginestar et al., 1998) and (Vidal-Ferràndiz et al., 2016). The implemented finite element code to integrate the time dependent neutron diffusion equations is named FEMFFUSION.
The semi-discrete equations are of the form, where L and M are the matrices obtained from the spatial discretization of operators L and M. Φ and C k are the vector of coefficients that represent the neutron flux and the precursors concentration, respectively, in terms of the polynomials used in the FEM. The matrices X and [ṽ −1 ] are defined as where the matrix P is the mass matrix of the spatial discretization, which appears due to the fact that the polynomial basis used in the finite element method is not orthogonal. The neutron precursors equation (4) is integrated with an explicit scheme as In the same way, Euler's backward method is used in equation (3) and taking into account equation (5), we get a system of linear equations, where the matrices are defined as, , and the coefficientâ isâ This system of equations is large and sparse and has to be solved for each new time step with a high accuracy due to the small size of the perturbation compared to the reactor dimensions. The code uses a fixed time step of ∆t = 10 −4 s and the tolerance for the relative error to the solve the linear systems is set to 10 −12 . The preconditioned GMRES method (Saad, 2003) was chosen as the solver for the linear system of equations.
To solve the problem given by equation (6) and to update the neutron precursors concentration with equation (5), a matrix-free methodology was implemented (Kronbichler and Kormann, 2012). Only the diagonal blocks of matrix T are stored using a sparse matrix format in order to be able to use an Incomplete LU factorization as preconditioner. The same strategy was used to solve the steady-state equation in (Vidal-Ferràndiz et al., 2019).
Finally, to compare the time domain computations with frequency domain analysis and to evaluate the neutron noise, we must define the flux noise as Then, we can apply a Fourier transform to the temporal flux noise, This transform is applied numerically using the Fast Fourier Transform (FFT) algorithm over the time dependent results. In the case of FA vibrations, the frequency space can be cut for the frequency of the assembly oscillation because the noise spectrum is composed basically of only this frequency.

Principles of frequency domain simulations
The CORE SIM tool  solves the neutron diffusion equation in the frequency domain and was already successfully used to study other neutron noise sources (Demazière and Andhill, 2005;Demazière and Pázsit, 2009;Pázsit and Demazière, 2010). This tool uses a finite differences scheme to solve the two energy groups diffusion equation in the frequency domain for the first-order neutron noise approximation. The tool, specifically made for noise calculations, was also verified against many (semi-)analytical solutions . It must be taken into account that FEMFFUSION allows any number of delayed neutron precursor groups while CORE SIM only uses 1 group of precursors, so effective values, β eff and λ eff , are used. However, as shown in this work, this approximation does not have a remarkable effect on the accuracy of the results for usual natural FA vibration frequencies because these frequencies fall in the middle of the so-called plateau region (Demazière et al., 2017a).
CORE SIM uses a spatial discretization scheme based on central finite differences and assumes that the scalar neutron flux in the middle of the nodes is equal to the node-averaged scalar neutron flux (box-scheme). Version 1.2 of the code only allows a uniformly refined mesh. This constraint enforces to use a fine mesh of less than 2 cm wide in order to obtain accurate results for light water reactor applications.
The first-order neutron noise theory is based on splitting every time dependent term, expressed as U(r, t), into their mean value, U 0 , which is considered as the steady-state solution, and their fluctuation around the mean value, δU as U(r, t) = U 0 (r) + δU(r, t).
The fluctuations are assumed to be small compared to the mean values. This allows to neglect second-order terms (δU(r, t) × δU(r, t)) = 0. Also, the fluctuations of the diffusion coefficients are neglected and δD g = 0 is assumed. This approximation was demonstrated to be valid for light water reactor applications (Larsson and Demazière, 2009). Thus, the first-order neutron noise equation can be written as ) where The neutron noise equation is an in-homogeneous equation with complex quantities that has to be solved after the steady-state solution is obtained, as φ 1 and φ 2 represent the steady state fast and thermal neutron fluxes, respectively. The related static eigenvalue problem must be solved with the same spatial discretization as the frequency domain neutron noise equation to get coherent results.

Modelling a vibrating fuel assembly in the frequency domain
Now, we describe how to model the vibration of an FA as a perturbation of cross sections in the frequency domain, following the derivation originally introduced for modelling the vibrations of control rods (see (Pázsit, 1988)), thereafter the modelling of core barrel vibrations (see e.g. (Sunde et al., 2006)), and more recently, the modelling of vibrating fuel assemblies (see (Jonsson et al., 2012)). This model is often referred to as the " /D model". One assumes that the cross section, Σ α , at the interface x = b between two material regions, as the one shown in Figure 1, is described as: where H is the unit step function, Σ I α and Σ II α are the cross sections at region I and II, respectively. Therefore, a vibrating assembly can be described as two in-phase moving interfaces (one moving interface on one side of the assembly and another moving interface on the other side of the moving assembly). For the sake of simplicity, only one moving interface is first considered in the derivations presented hereafter. An interface moving as b(t) = b 0 + A sin(ω p t), results in: Using a first order Taylor expansion, the cross section perturbation can be expressed as and, in the frequency domain, the perturbation is written as follows Since in numerical tools like CORE SIM the perturbation is introduced node-wise, one could assume that the perturbed region is If the mesh used does not match the perturbed region, the perturbation must be rescaled.
To ensure the validity of the first order Taylor approximation, a numerical fast Fourier transform of the time dependent cross sections is calculated. The numerical FFT was obtained for a perturbation of 1 Hz. Figure 2 shows the frequency spectrum of the perturbation at x = b 0 using the FFT. The amplitude of the cross section perturbation is maximum at 1 Hz, ω p , and the other frequencies are at least three times smaller while the first order approximation only uses the perturbation at 1 Hz. Figure 3 shows the spatial dependence of the perturbation |δΣ α | at 1 Hz for the FFT and the first order approximation. It can be seen that the first order approximation has the same perturbation integral as the numerical FFT. This statement is based on the comparison of numerically integrated values of the area under the perturbation curve. This shows that the first order approximation can be introduced using only one node, thus being the most convenient way to introduce the vibrating FA.

One-dimensional benchmark
In order to test the numerical tools developed for the FA vibration analysis, a simple onedimensional benchmark is defined. The benchmark is composed of 11 assemblies of 25 cm wide where the vibrating assembly is placed in the middle of the reactor as Figure 4 shows. The cross sections are defined in Table 1 and zero flux boundary conditions are imposed on the left and right frontiers of the reactor. Kinetic data are shown in Table 2. The problem is made critical before starting the time dependent calculation by diving νΣ f g by the previously calculated multiplicative factor of the problem, k eff = 0.9792500.
The movement of the central assembly is defined as where x i (t) represents the position of the vibrating assembly along time, originally placed in x i0 . Figure 6a shows the relative power evolution for an oscillation of A = 1 mm of amplitude and a frequency of ω p 2π = 1 Hz along 10 periods calculated with the time-domain code FEMFFU-SION. A sinusoidal evolution in the relative power can be seen with a small amplitude, around 7.87e-8 and with a constant increment along time. This increment is caused because the reactor becomes slightly supercritical when the central assembly moves from its starting position due to this assembly is in the middle of the reactor and this assembly is less reactive than the surrounding assemblies. Figure 6b displays the static multiplicative factor through the positions travelled during one period. It can be seen that the change in the k eff is less that 1.2e-9. These Figures also show that the sinusoidal global results present twice the frequency of the mechanical FA oscillation, because the reactor is symmetric and the moving assembly is in the centre of the reactor. In other words, the global reactor results of the assembly moving to the right are equal to the results when the assembly is moving to the left provoking a 2 Hz oscillation in the power and the multiplicative factor. However, the space dependent neutron noise has the same frequency as the mechanical FA vibration. The behaviour of the total power is analogous with the one studied analytically in a point kinetic reactor model in (Akcasu, 1958) and (Ravetto, 1997) from a sinusoidal change in reactivity where the linear increase of the power is caused by the delayed neutrons. However, thermal-hydraulic feedback will eliminate the slow increment in the total power because, in an operating nuclear reactor, the temperature coefficients of reactivity are negative (Stacey, 2007) and the increment in the power is very small.
In Figures 6a and 6b, two non-equidistant meshes are compared to each other. In both cases, a locally refined mesh around the moving FA is used. The coarsest mesh uses 47 cells with the following sizes [4 * 25.00, 24.5, 4 * 0.1, 10 * 0.02, 4 * 0.1, 24.0, 4 * 0.1, 10 * 0.02, 4 * 0.1, 24.5, 4 * 25.00] cm. A detail of the refined mesh near the vibrating assembly is shown in Figure 5. This model of refinement is based on the initial assembly configuration adding 20 cells where the movement of the assembly is located and 4 transition cells at each side of this area to make a smooth change of the refinement level. The other locally refined mesh uses 94 cells where each of the previous cells is split into two cells. Also, a uniform mesh with 17600 cells is utilized for (1/cm) (1/cm) (1/cm) (1/cm) Fuel 1 1.40343 1.17659e-2 5.62285e-3 2.20503e-3 1.60795e-2 2 0.32886 1.07186e-1 1.45865e-1 5.90546e-2 Vibrating Assembly 1 1.40343 1.17659e-2 5.60285e-3 2.19720e-3 1.60795e-2 2 0.32886 1.07186e-1 1.45403e-1 5.88676e-2 Reflector 1 0.93344 2.81676e-3 0.00000e+0 0.00000e+0 1.08805e-2 2 0.25793 8.87200e-2 0.00000e+0 0.00000e+0  The results with 47 cells mesh can thus be considered as spatially converged. The refinement model will be still valid when the oscillation amplitude is increased and the size of the cells is changed accordingly. Figure 7 shows the evolution of the reactor power and the k eff using a smaller number of uniform spatial mesh cells. When using too coarse meshes, the computed solution deviates significantly from the converged solution. These Figures indicate the necessity of using local refinements around the oscillating assembly and accurate numerical solvers to correctly integrate the time dependent neutron diffusion equation.   Figure 8 shows the power evolution for different oscillation amplitudes from 0.3 mm to 3 mm while the frequency is fixed to 1 Hz for a mesh with 47 cells. The size of the mesh cells has been adapted to the oscillation amplitude ensuring that the obtained results are spatially converged. Obviously, as the oscillation amplitude increases, its effect in the power increases. Figure 9 represents the spatial dependence of the neutron flux disturbance around the steady-state value at four different simulation times during the first period of oscillation. It can also be verified that the neutron noise response to the vibration disturbance is basically a sinusoidal function at the FA oscillation frequency.  Figure 9 because the principal response of the neutron noise is a sinusoidal function at the FA vibration frequency. Figure 11 displays the results for the phase. A close agreement is observed between the time domain and the frequency domain methodologies.

Two-dimensional BIBLIS Reactor
The BIBLIS benchmark is selected to verify the frequency and time domain analysis for a vibrating FA in a two-dimensional reactor. It is a classical two energy group neutron diffusion problem taken as a benchmark for different numerical codes. The materials composition of the reactor is represented in Figure 12. The cross sections for each material are given in Table 3. Vacuum boundary conditions are applied at the boundary of the reactor. The full definition of the benchmark can be found in (Hébert, 1985).
The problem is made critical before starting the time dependent calculation by dividing νΣ f g by the previously calculated multiplicative factor of the reactor, k eff = 1.025482. In this work, the usual definition of 6 groups of precursors is used in the time domain calculation. For the CORE SIM calculation, 1 group effective values, β eff and λ eff , are used. The kinetic parameters used are given in Table 4.
The assembly placed in position (6, 6), marked with a cross, ×, in Figure 12, is selected to be oscillating along the x direction.
Due to the different scales of the problem, a fine mesh needs to be used to accurately solve the system. In the time domain analysis, a locally refined mesh in the x direction following the same scheme as in the one-dimensional benchmark with 869 cells is used. Cubic polynomials in the FEM are used for this problem. In the frequency domain analysis, a uniform mesh of 4624 cells is employed. If these fine meshes are not used, the effect of the FA vibration can be erroneously estimated.  Figure 13a shows the relative power evolution for an oscillation of 1 mm of amplitude and a frequency of 1 Hz along 3 periods. A sinusoidal change in the relative power with a really    small amplitude, about 3.8e-7, with a linear increment over time is observed. As it has been already mentioned, this increment is caused by the change of the criticality of the system when the assembly is moved from its starting position. Figure 13b displays the multiplicative factor of the reactor, k eff , through the positions travelled by the moving FA during one period. It can be seen that the change in the k eff is about 3.8e-8. In this case, as the vibrating assembly is no situated in the middle of the reactor, the k eff does not behave in a completely symmetric way. The behaviour of the power is analogous to the one previously studied in the one-dimensional benchmark. Figure 14 displays the spatial distribution of the amplitude of the neutron noise in the frequency domain calculated with the time dependent code FEMFFUSION and performing a FFT of the results and considering only the frequency of the perturbation. A wider influence of the perturbation on the fast flux can be seen. Otherwise, the thermal noise is localised in the surroundings of the oscillating FA. Two clear peaks can be observed in these Figures. They correspond to the perturbed region where the cross sections change along the FA movement. Figure 15 shows a comparison of the amplitude of neutron flux noise for the fast and thermal groups along y = 150.2969 cm between the FEMFFUSION code and the CORE SIM code, in other words, between the time domain and the frequency domain approaches. Figure 16 gives a comparison of the phase of the neutron noise between the two codes. A close agreement is observed for both the amplitude and the phase of the neutron noise, validating the time domain analysis against the frequency domain one. Figure 17 shows the fast and the thermal neutron noise magnitudes for three different oscillation amplitudes calculated with CORE SIM. Figure 18 displays the neutron noise magnitude associated with the fast and thermal fluxes for different vibrations amplitudes at the neutron detector, situated at (x = 104.0517, y = 150.2969) cm marked with a • in Figure 12. These results are calculated with FEMFFUSION and CORE SIM codes. In these Figures, a proportional dependency of the neutron noise magnitude with the amplitude of the FA vibration amplitude can be observed. Figure 19 shows the fast and the thermal neutron noise magnitudes for three different oscillation frequencies computed with CORE SIM. Even though the vibration frequencies are extreme as compared to realistic mechanical vibration of FA, the neutron noise results are similar. Figure 20 displays the neutron noise magnitude against the vibration frequency at the position of the neutron detector calculated with CORE SIM and FEMFFUSION with 1 group of delayed neutron precursors (1 gdnp) and 6 groups of delayed neutron precursors. This Figure is similar to the reactor transfer function shown in (Demazière et al., 2017a) where the usual mechanical FA vibration frequencies fall in the so-called plateau region, defined as the frequency range [λ eff , β eff /Λ 0 ]. This indicates that with respect to the amplitude of the neutron noise, the frequency-dependence of the reactor response tends to follow the one of the point-kinetics. As it can be seen for usual FA vibration frequencies, the neutron noise is basically independent of the vibration frequency. Since the effect of non-linearities is negligible, no resonance effect is predicted. The inclusion of thermalhydraulic feedback might nevertheless modify the frequency response of the system. Also, a close match between the time-domain and the frequency domain methodologies is obtained. Only at low frequencies, i.e. at frequencies smaller than 0.5 Hz, a discrepancy between calculations with different number of delayed neutron precursors groups taken into account can be observed.

Conclusions
The current research is an attempt to understand in more detail the coupling mechanism between the mechanical vibration of fuel assemblies and the generated neutron noise without thermal-hydraulic feedback. The problem combines different spatial scales, e.g. a 1 mm oscillation of a fuel assembly that measures around 20 cm in a reactor of some meters size. This implies that we need to work with a very high precision in the spatial discretization and in the tolerances allowed for the solvers.
Numerical results show two different effects in the neutron field caused by the fuel assembly vibration. First, a global slow variation of the power due to a change in the criticality of the system is observed. This effect is small and will be compensated by the thermal-hydraulic coupling because, in an operating nuclear reactor, temperature coefficients of reactivity are negative. Secondly, an in-phase oscillation in the neutron flux with the assembly vibration is demonstrated. The corresponding neutron noise is highly spatially dependent. For this second effect, this paper shows a comparison between a time domain analysis and a frequency domain analysis of the phenomena. Numerical results show a really close match between these two approaches. This paper also demonstrates that neglecting second order effects, as done in the frequency domain simulations, gives essentially the same results as not neglecting those, as done in the time domain simulations. Also, it must be noted that the frequency domain analysis takes much less computational time than the time domain analysis.
The simplified model for the vibration fuel assemblies, modelled as the movement of a homogeneous material is similar to the " /D model". This model is expected to be correct some diffusion lengths away of the vibrating assembly where the neutron detectors are located. In future works, experimental results in the framework of the CORTEX project are expected to be compared with this model in order to validate it.