Turbulence-induced vibrations prediction through use of an anisotropic pressure fluctuation model

. In nuclear fuel rod bundles, turbulence-induced pressure ﬂuctuations caused by an axial ﬂow can create small but signiﬁcant vibrations in the fuel rods, which in turn can cause structural eﬀects such as material fatigue and fretting wear. Fluid-structure interaction simulations can be used to model these vibrations, but for aﬀordable simulations based on the URANS approach, a model for the pressure ﬂuctuations must be utilised. Driven by the goal to improve the current state-of-the-art pressure ﬂuc-tuation model, AniPFM (Anisotropic Pressure Fluctuation Model) was developed. AniPFM can model velocity ﬂuctuations based on anisotropic Reynolds stress tensors, with temporal correlation through the convection and decorrelation of turbulence. From these velocity ﬂuctuations and the mean ﬂow properties, the pressure ﬂuctuations are calculated. The model was applied to several test cases and shows promising results in terms of reproducing qualitatively similar ﬂow structures, as well as predicting the root-mean-squared pressure ﬂuctuations. While further validation is being performed, the AniPFM has already demonstrated its potential for aﬀordable simulations of turbulence-induced vibrations in industrial nuclear applications.


Introduction
A particularly important topic in the field of nuclear safety is the behaviour of fuel rods.Fuel rods are submerged in a coolant, such as water in current-day reactors or liquid metal in advanced fast reactors, which flows axially along the fuel rods.While the axial flow leads to efficient cooling of the fuel rods, it can also cause Turbulence-Induced Vibrations (TIV) of the fuel rods, which originate from pressure fluctuations in the flow.This phenomenon plays a critical role in terms of nuclear safety, as it can cause structural effects such as fatigue problems, fretting wear, and stress corrosion cracking [1,2].The phenomenon has been studied since the start of nuclear reactor development in the 1950s, and it has been the root of several incidents [1,3].In the previous century, an emphasis was put on experiments to increase the understanding and knowledge of TIV, as well as using semi-empirical laws to establish a relation between the amplitude of vibration and relevant flow parameters [4,5].More recently, due to the larger availability of computational resources, the use of Fluid-Structure Interaction (FSI) simulation applied to fuel rods has received more interest.In one approach the turbulence is (partially) scale resolved (LES/DNS), but * e-mail: noutvdbos@gmail.comdue to the excessive computational requirements, these simulations typically rely on a simplified approach for the fluid-structure interaction, at moderate Reynolds numbers [6][7][8].Due to the large Reynolds numbers associated with the axial flow along fuel rods, LES/DNS is too expensive for industrial use of simulations of fuel rods or fuel assemblies.For reference, the LES calculations of De Ridder [9] took 2200 days in equivalent computational time to simulate 0.1 s of 1/10th of the length of a fuel rod.Extending this to a fuel assembly is not realistic.Furthermore, these simulations are typically one-way coupled.However, due to the large density ratios that are typically found in nuclear reactors, two-way FSI coupling is preferred, in order to capture the effect of the changing flow field on the fuel rods.A computationally much cheaper method is to use URANS for the fluid domain [10,11], which would allow simulations of fuel assemblies.However, URANS causes severe underprediction in the vibration amplitude.This is because URANS only calculates forces resulting from the mean flow, as it does not resolve the instantaneous turbulent pressure fluctuations.To remedy this, the Pressure Fluctuation Model (PFM) was proposed by Kottapalli et al. [12], which assumes that the turbulence statistics from URANS are sufficiently accurate to construct synthetic turbulent pressure fluctuations to be superimposed on the Reynolds-averaged pressure.This would artificially reintroduce the turbulence-induced forcing to the structure, at a lower computational cost than scale-resolving approaches.From comparison to experimental data, it was found that the PFM shows vibration amplitudes in the same order of magnitude, however, it did not give the required accuracy yet [12].In the current paper, an improved pressure fluctuation model, called AniPFM, is proposed.The new model tackles several limiting assumptions of the PFM, such as the assumption of isotropic turbulence, as well as the method of time correlation.With this, the AniPFM the accuracy in the prediction of vibration amplitudes through FSI-simulations of TIV. Figure 1 summarizes how the AniPFM interacts with the other components in FSI simulations.

Methodology
The AniPFM is based on the formulation given of the PFM of Kottapalli et al. [12] and Senthooran et al. [13], though major modifications have been made to incorporate anisotropy and time correlation, using approaches from Billson et al. [14], and Shur et al. [15].Details of the new AniPFM are given in this section.

Pressure fluctuations
Similarly, to the method of Kottapalli et al. [12], the divergence operator is applied to the Reynolds averaged and instantaneous Navier-Stokes equations, to obtain the following Poisson equation: Equation (1) shows that the pressure fluctuation p solely depends on the mean velocity u i , the Reynolds stresses u i u j , and the velocity fluctuations u i .Only the velocity fluctuations are unknown.The governing equation for computing the pressure fluctuations is given by equation (1) where the velocity fluctuations are modelled.
Note that equation (1) assumes that the full spectrum of pressure fluctuations is computed, but in reality, only the fluctuations that can be resolved by the mesh are taken into account.Thus, the effect of unresolved turbulence on pressure fluctuations is not accounted for in this method.

Velocity fluctuations
Constructing the velocity fluctuations consists of three main parts, namely the construction of dimensionless velocity fluctuations, correlating the fluctuations in time, and finally scaling the fluctuations to match the input Reynolds stresses.The dimensionless fluctuations w t (x) are created by a Fourier decomposition, as shown in equation (2).Here, q n is the mode amplitude, σ n is the direction vector, κ n is the wavenumber vector, and φ n is a random phase shift with a uniform distribution, the subscript n denotes that it is for the nth mode.

Turbulent kinetic energy spectrum
The amplitude for each mode is determined by the nondimensionalized Von-Karman spectrum, shown in equation (3) [15].Here, k e is the eddy wavenumber, which is the wavenumber which contains the highest energy density.The variable k η refers to the Kolmogorov wavenumber, and f cut is a cut-off frequency based on the maximum resolvable wavenumber.As shown in equation ( 4), the amplitude is based on the relative weight of the energy of mode n.Note that even though the AniPFM can represent anisotropic turbulence, the energy spectrum is based on isotropic turbulence, thus it is still expected that the results near walls are not as accurate as in fully isotropic conditions.
The eddy wavenumber is based on the RANS length scale as shown in equation ( 5), where C l is a calibration coefficient set to 3.0, based on isotropic turbulence.
The cut-off function is shown in equation ( 6) [15].For the cut-off length, several definitions were investigated.However, it was found that the cut-off length as defined by Shur et al. [15] (shown in Eq. ( 7)) gives the most accurate results.

Wavenumber and direction vector
As shown in equation ( 2), the velocity fluctuations are determined by a sum of wavenumbers.As the sum is finite, it requires the specification of a start and final wavenumber, as well as the distribution of discrete wavenumbers in the resolved wavenumber space.The starting wavenumber is based on a conservative dimensional analysis, shown in equation ( 8).In the case of a zero-mean velocity, the RANS length scale is used to determine the starting wavenumber.Additionally, the starting wavenumber must also be lower than the wavenumber with the maximum energy density.Finally, the geometrical considerations must also be taken into account.This is done by introducing a user length scale, which denotes the maximum length that can be captured given the geometry of the problem.These requirements are fulfilled by the expression for the starting wavenumber, as shown in equation (9).
The highest wavenumber is based on the cut-off wavenumber (see Eq. ( 10)), as the energy spectrum quickly goes to zero after this point.In between the lowest and highest wavenumber, the wavenumbers are logarithmically distributed.
The wavenumber vectors are defined similarly to [12,14], as shown in Figure 2. The mathematical description is displayed in equation (11).Here, θ n and ψ n are uniformly distributed random variables, described by equation (12).
The wavenumber direction vector is determined from the wavenumber vector.From applying the divergence operator to equation ( 1), it can be found that k n • σ n = 0 must hold in order to adhere to the continuity equation.This is the only requirement.To achieve this, σ n is defined as the normalised cross-product between a random vector ζ n , and the wavenumber vector, as shown in equation ( 13).
Since the dot-product of a vector with the cross product of the same vector is always equal to zero, the continuity condition is met.

Time correlation and scaling
The procedure for calculating the non-dimensional fluctuations w t (x) has been fully defined in the previous subsections.The next step is to create a space-timedependent velocity signal v t (x, t), from the purely spacedependent signal.Two main phenomena have to be taken into account when constructing a time-correlated velocity field.Namely, the convection of turbulent eddies and the decorrelation due to the production and dissipation terms.
A method similar to Billson et al. [14] is used to introduce this time dependency.First, the fluctuations are convected by solving for equation (14).Here, v m−1 t are the non-dimensional velocity fluctuations generated at timestep m − 1, and U j is the Reynolds-averaged velocity as calculated by the accompanied URANS simulation.Then in the second part, a new solution v m t is calculated from a combination of the (convected) previous solution v m−1 t , and a newly generated field w m t , as shown in equation (15).
The coefficients a and b are defined in equation (16), where τ is the dissipation timescale determined from the URANS simulation, and f τ is a modification factor for fine-tuning the correlation.For the latter, a value of f τ = 17 is used, in line with what was used in [14] for the simulation of a 3D jet.Tests showed that this value also gives a satisfactory correlation in the simulations of interest.The coefficients a and b are defined such that the squared mean properties of v t are still respected, i.e.
Finally, the last step is to scale the dimensionless fluctuations, such that they replicate the Reynolds stresses.Since v t 2 is equal to the Kronecker delta, the fluctuations need to be scaled with a tensor a ij such that a ij a ji = R ij , as shown in equation ( 17).This can be realised with a Cholesky decomposition of the Reynolds stress tensor, which is shown in equation (18).
With this method, the AniPFM can reconstruct anisotropic Reynolds stresses.For flows with a constant pressure gradient, such as channel flows, linear Eddy viscosity models show isotropic Reynolds stresses.In order to improve the accuracy of these models, a correction is used to transform the isotropic tensor into an anisotropic tensor, based on the nonlinear Eddy viscosity model of Wilcox [16].The Wilcox correction is shown in equation (19).

Model overview
This concludes the full model of the AniPFM.In Figure 3, the full model is summarised.First, the non-dimensional velocity fluctuations are calculated, based on the energy spectrum.Then the time correlation is performed.After this, the velocity fluctuations are computed by scaling with the Cholesky tensor.Finally, p can be computed.

Simulation set-up
Two cases were simulated, in order to validate the AniPFM.First, a Homogeneous Isotropic Turbulent (HIT) box was simulated, which was used to verify the AniPFM in isotropic conditions, since an isotropic energy spectrum was used as input.Second, a turbulent channel flow was simulated, with the purpose of testing the accuracy of the AniPFM versus the PFM in anisotropic turbulence.The setup of these two cases is discussed in this section.The AniPFM is implemented in OpenFOAM 8, and thus all subsequent simulations have been performed with the same software.

Homogeneous isotropic turbulent box
A box of L•L•L is created, which is discretized by a cartesian mesh of N • N • N cells.All boundaries are periodic, and the domain has a zero-mean velocity.The experiment of Comte-Bellot and Corrsin [17] has been replicated, as well as the DNS of Gotoh and Fukayama [18].In Table 1,  Re λ = 287 of Gotoh and Fukayama is replicated; the setup is also shown in Table 1.

Turbulent channel flow
The results of the TCF are compared to the DNS results of Abe et al. [19].The simulation with Re τ = 640 is used for comparison.In this case, the flow between two infinitely long and wide stationary plates is simulated.The simulation domain used is equal to 6δ • 2δ • 3δ, where δ is the mid-channel height.The mesh has a size of N x • N y • N z cells, which are kept as variables.The mesh distribution is uniform in the stream-and spanwise-direction, and it is geometrically expanding to the mid-channel plane in the wall-normal direction.The mesh in the wall-normal direction is configured such that y + ≈ 1 for the first grid cell from the wall.The setup is summarised in Table 2, and   an example mesh is shown in Figure 4. Simulations were performed with the k − ω SST turbulence model, along with the Wilcox correction.However, for several simulations, no URANS calculation was done, but rather the mean flow properties of the DNS of Abe et al. were used as an input to the AniPFM.This was used to isolate any errors that originate due to the AniPFM, and not due to the input.

Results discussion
The results of both cases are elaborated upon in this section.In Table 3, the V&V done for all variables the test cases is shown.In this section, the most important results from this verification are shown and discussed.

Homogeneous isotropic turbulent box
In Figure 5, the simulated velocity fluctuations of the HIT box of Comte-Bellot Corrsin are shown, with a mesh size of N = 128.The velocity fluctuations on the left are from the AniPFM, whereas the results on the right are simulated through means of a Large-Eddy Simulation [20].
From Figure 5, it was concluded that the AniPFM can qualitatively reconstruct similar flow structures as highresolution methods.
In Figure 6, the energy spectra generated by both the AniPFM and the PFM of Kottapalli et al. are compared to the experimental spectrum of Comte-Bellot & Corrsin.The AniPFM results show a very good resemblance with respect to the experimental results, right up to the cutoff wavenumber.Compared to the results of the PFM, the peak of the energy spectrum is better predicted, and there is no unphysical accumulation of energy near the cut-off wavenumber.The latter difference is due to the fact that a cut-off filter is used in the AniPFM, which makes sure that the unresolved energy is not redistributed over the resolved part of the spectrum.
The effect of the cut-off filter can also be found in the root-mean-squared pressure fluctuations.With a mesh of N = 64, the AniPFM predicted p rms within a 2.4% error with respect to the results of Goth et al., whereas the PFM predicted a 3.7% error.The AniPFM showed an error of 1.1% when the mesh was refined N = 128.

Turbulent channel flow
For the turbulent channel flow, there are several sources of errors that cause a discrepancy between the AniPFM results and experimental/DNS data.Therefore, it is important that each source is carefully evaluated.
There exists an uncertainty range in the results of the AniPFM, which is due to the random numbers that are used throughout the simulation.For the chosen temporal correlation scheme, the random numbers gave an uncertainty range of ±0.25% for p rms , for a confidence interval of 95.4%.The number of modes did not seem to have a large effect on p rms , however, with a lower number of modes, the uncertainty went up, as fewer random numbers were used per iteration.For the presented simulations, 1024 wave number modes were used.
The amount of energy that is resolved by the AniPFM depends on the fineness of the mesh.There is no sub-grid model that models the effect of unresolved velocity fluctuations on the pressure fluctuations.Thus, the statistics of the pressure fluctuations only converge if a large part  of the velocity fluctuations is resolved.Since p rms is the quantity of interest, this quantity is evaluated for several meshes.This is shown in Figure 6, in which the results were obtained using the k − ω SST turbulence model, together with the Wilcox correction.It can be seen that the meshes converge to a final solution when using finer meshes.At the second finest mesh, the results near the wall are within the given uncertainty range, thus deeming the solution converged.
The modelling error was found by performing a simulation of the turbulent channel flow on the finest mesh from Figure 7, but with the RANS input variables (k, ε, u i and u i u j ) taken from DNS data.The replicated Reynolds stresses and p rms are shown in Figure 8.The Reynolds stresses are very closely approximated.The small underestimation is because the velocity fluctuations are not fully resolved.Nevertheless, it was found that these unresolved fluctuations had no effect on p rms .From Figure 8 (right), it can be seen that near the midchannel plane, the AniPFM very closely approximates p rms .This is because here the isotropic energy spectrum very closely approximates the actual energy spectrum.Near the wall, the energy spectrum is not approximated as accurately, due to the larger anisotropy in the flow.Thus, a larger discrepancy in p rms was found near the wall, with a maximum error of 4.4%.
When the k−ω SST turbulence model, together with the Wilcox correction was used as URANS input for the turbulent channel flow, an error of roughly 10% was found near the wall.This increase in error is due to the large difference in the Reynolds stress tensor between the selected RANS model and the DNS data.In comparison, the PFM showed an underestimation of roughly 47% at the wall.Thus, the AniPFM shows a sharp improvement in the prediction of p rms for anisotropic flows, with respect to the PFM.

Conclusion
In this paper, a new pressure fluctuation model was proposed, called AniPFM.This model improves the prediction of pressure fluctuations when using a URANS approach, which can be useful in particular for turbulenceinduced vibration prediction.Several aspects of the AniPFM were adjusted with respect to the previous PFM of Kottapalli et al., namely the energy spectrum cut-off filter, the replication of anisotropic Reynolds stresses, and the method for time correlation.
Two test cases were performed to assess the performance of the AniPFM: namely, a homogeneous isotropic turbulent box, and a turbulent channel flow.From the HIT case, it was found that the prediction for p rms was slightly improved due to the addition of the cut-off filter.From the turbulent channel flow, it was found that the introduction of anisotropy gave a slightly higher model error.This is because the energy spectrum is still based on isotropic turbulence.Nevertheless, the AniPFM showed an improvement in the prediction of p rms at the wall, compared to the PFM of Kottapalli et al.
While the initial results give an optimistic view, there are still several challenges that need to be tackled.Further validation must be done, both for fluid-only cases, as well as for relevant use cases.While there is no consensus if the AniPFM two-way coupled approach would be an improvement with respect to the one-way coupled LES approach for determining fuel pin vibration, the research does show to be of interest for strongly coupled cases, such as full fuel assembly FSI simulations.For this, both the computational costs and accuracy of the AniPFM must be compared to LES or Hybrid URANS/LES FSI simulations, to evaluate if the AniPFM is beneficial to use.Finally, there are still assumptions in the AniPFM which have to be evaluated, such as the use of an isotropic energy spectrum scaled to meet the Reynolds stress requirements.The construction of this spectrum is based on integral properties of the flow, with assumed isotropy, meaning that the energy distribution might not be well represented in areas of high anisotropy.
The next step is to apply the model to fluid flow cases that resemble nuclear fuel rod operating conditions more closely.While validation is still ongoing, the AniPFM shows potential for accurate, computationally cheap simulations of turbulence-induced vibrations in industrial nuclear applications.

Fig. 1 .
Fig. 1.A flow-chart of the global interactions between the several numerical solvers.

Fig. 3 .
Fig. 3. Flow chart of the different computational steps of the proposed AniPFM.

Table 1 .
Details of the simulation replicating Comte-Bellot and Corrsin, Gotoh and Fukayama at Re λ = 70.Turbulence Model k − ε k − ε Initial k [m 2 /s 2 ] the replicated set-up of Comte-Bellot and Corrsin are shown.The DNS of Gotoh and Fukayama was specified in dimensionless numbers.The simulation with

Fig. 5 .
Fig. 5.The instantaneous velocity fluctuations of the AniPFM, compared to LES results.

Fig. 7 .
Fig. 7.The root-mean-squared pressure fluctuations along the wall-normal coordinate for various meshes, versus the DNS results of Abe et al. [19].

Fig. 8 .
Fig. 8.The Reynolds stress profiles (left) and the RMS pressure fluctuations (right) along the wall-normal coordinate, versus the DNS results of Abe et al.[19].DNS data is used as input.

Table 2 .
Details of the channel flow simulation set-up.

Table 3 .
All variables that have been compared against experimental and/or DNS results.
CaseV&V variables HIT TKE, p spectrum, TKE spectrum, p rms u i u j time correlation TCF w/URANS input p rms , R ij TCF w/DNS input p rms , R ij , p spectrum, TKE spectrum, p frequency spectrum