Low-Noise Synthetic Turbulence Tailored to Lateral Periodic Boundary Conditions

: The present work is dedicated to turbulence synthesis tailored to lateral periodic boundary conditions for direct noise computations through compressible large eddy simulations. Synthetic turbulence can be essential for aeroacoustic applications when computing airfoil turbulent inﬂow noise or for accurately capturing the behavior of boundary layers. This behavior determines both trailing edge noise and complex ﬂow structures such as laminar separation bubbles. For airfoil simulation purposes, spanwise periodic boundary conditions are usually considered. If synthetic perturbations are injected without observing the periodicity rule, strong spurious pressure waves are emitted and pollute the entire computational domain. In this work, the random Fourier modes method for turbulence generation is adapted in order to respect the spanwise periodicity constraint right at the computational domain inlet. This approach does not affect the turbulence properties such as the spectral shape and the turbulent kinetic energy decay. Since the emphasis is put on the generation and convection of the turbulence, only the turbulence convection region between the inlet and the airfoil is considered in this paper, without the airfoil. Two geometrical conﬁgurations are tested: the ﬁrst one is a simple box with a constant mesh size, and the second one concentrates the ﬁne cells on the area in front of the airfoil. In the second conﬁguration, the computational cost is reduced by up to 25%, but more spurious noise is present because of interpolation areas between different grids using the Chimera method. Finally, the results’ reproducibility is assessed using different turbulence realizations. T.R.; T.R.; T.R.; T.R.; and and P.L.; visualization, T.R.; supervision, and P.L.; administration, P.L. and T.R.;


Introduction
Low-Reynolds number (i.e., Re ≤ 5 × 10 5 [1]) aerodynamic applications such as small wind turbines or drones appear to be more and more present for civil or military use. The flow on blades can be subject to complex phenomena such as laminar boundary layer separation, transition towards turbulence, boundary layer reattachment or even stall [1]. Boundary layer separation and reattachment form the well-known laminar separation bubble (LSB) that can lead to wing stall when the LSB bursts at high angles of attack. In this situation, the underlying mechanisms still need to be understood [2,3]. Furthermore, the noise emitted by drones at low Reynolds numbers is a subject of growing interest [4,5]. It is also the case for higher-Reynolds number applications such as full-scale wind turbines [6].
For both low-and high-Reynolds number applications, inflow turbulence is of significant importance. In the case of a full-scale wind turbine, inflow turbulence is responsible for broadband low-frequency noise [7,8]. For low-Reynolds number applications, inflow turbulence not only affects emitted noise but also strongly impacts the aerodynamics. Indeed, it is experimentally suggested that the transition towards turbulence of the LSB on an airfoil, its length and its burst are highly affected by inflow perturbations [2,3]. The effect of inflow turbulence on the LSB separation, transition and reattachment points has also been numerically studied using different turbulent intensities or integral length scales. For example, Hosseinverdi et al. [9] generated an LSB on a flat plate and injected turbulence built with Fourier modes, while Schmidt et al. [10] injected turbulence based on a modified digital filter approach to disturb an LSB present on an SD7003 airfoil.
For acoustic predictions or LSB insights, numerical simulations such as large eddy simulations (LES) are a good compromise between resolution and computational cost. Previous studies showed the importance of generating realistic inflow turbulence in simulations. There are still two major challenges concerning synthetic turbulence. The first one is that most methods are unable to accurately reproduce experimental turbulence high-order statistics at the computational domain inlet [11]. The second challenge concerns direct noise computation (DNC) which requires that synthetic turbulence methods generate low levels of spurious noise [12].
The different categories and methods of synthetic inflow turbulence were reviewed in Dhamankar et al. [11]. While library-based turbulence methods are efficient in terms of statistics' reproducibility, they can only be used for specific applications and can be very expensive. The recycling methods avoid a long turbulence adaptation distance but present the drawback of generating spurious modes in the flow [11]. Furthermore, these methods exploit turbulent boundary layer similarity laws for rescaling. The last category of methods is based on synthetic turbulence generators and includes several approaches. These approaches are still missing correct third-order turbulence statistics but have the advantage of being directly embedded in solvers, without running any other simulation in parallel. Kraichnan [13] first proposed a spectral approach with a Gaussian distribution for the mode amplitudes. Smirnov et al. [14] and Batten et al. [15] extended this work by adding anisotropy and inhomogeneity properties. Among others, Karweit et al. [16] or Bechara et al. [17] generated isotropic turbulent fields with a deterministic energy distribution using a theoretical spectrum. This method, called random Fourier modes (RFM), consists of building turbulence with a realistic prescribed energy spectrum and has the advantage of naturally producing an incompressible velocity field, in order to synthesize low-noise turbulence. Sescu et al. [18] later showed that a synchronicity between the turbulence perturbations and the mean flow is needed to avoid the generation of spurious vorticity waves. The RFM method was used with two-dimensional structures to compute the leading edge noise of wavy airfoils with the Euler equations in Clair et al. [19].
The synthetic eddy method (SEM) [20] is another technique to produce turbulence and was initially used for channel flow incompressible simulations. This method has been modified to introduce realistic scale inhomogeneities in a flat-plate turbulent boundary layer in the wall-normal direction [21], and to reduce spurious noise emission [18] by making the resulting fluctuating velocity field incompressible. Later, Kim and Haeri [22] extended this work to retrieve a von Kármán spectrum by optimizing the method with constraint parameters. They applied it to turbulence interaction noise on a wavy airfoil with the Euler equations.
Another group of methods for the generation of inflow perturbations was introduced by Klein et al. [23]. With an incompressible solver, the authors used white noise digital filtering to create turbulent fields. In [23], these perturbations were employed to study the influence of the integral length scale, the mean flow velocity profile and the turbulent intensity on plane jets and on the interface of a water film ejected into air. Later, Ewert [24] introduced the random particle mesh method (RPM) to obtain a solenoidal velocity field using digital filtering on a fluctuating vector stream function. He successfully applied this method to compute slat noise. At the same time, Xie et al. [25] improved the computational efficiency of the digital filtering method. They applied it to a plane channel flow and to an urban-type configuration, using an incompressible solver. Kim et al. [26] then modified this last method by adding a mass flux correction at the inlet and by making the velocity field solenoidal through the Pressure-Implicit with Splitting of Operators (PISO) algorithm corrector steps. This allowed them to delete the spurious pressure waves by making the turbulence generation method consistent with the incompressible solver even if the adaptation distance increased. Kim and Xie [27] used this turbulence generation method to study the effect of freestream turbulence on a NACA0012 airfoil dynamic stall. They worked in the context of a small or medium-size wind turbine by choosing a chord-based Reynolds number equal to 1.35 × 10 5 . Recently, Gea-Aguilera et al. [28] extended the RPM method by replacing the white noise field by a scalar randomly taking the value ±1 for each synthetic eddy, which resulted in a computational cost improvement. This method is called advanced digital filtering and was used to compute airfoil leading edge noise by solving the linearized Euler equations.
For DNC, silent synthetic turbulence is sought but is still difficult to obtain. Indeed, in compressible codes, the velocity is often prescribed at the inlet without the thermodynamic variables, that leads the code to adapt the flow field and to generate undesired pressure fluctuations, as explained in Gloerfelt and Robinet [12]. Another source of spurious noise is found in airfoil aeroacoustic simulations. Indeed, spanwise periodic boundary conditions are commonly used to limit the computational cost, but the inflow turbulence does not necessarily respect this condition, that leads the solver to generate undesired pressure waves. Kim and Haeri [22] solved this problem for their modified SEM by duplicating the generated turbulence virtual boxes on each side of the span limits. Clair et al. [19] used in their Euler-based solver a two-dimensional RFM method built with a wavenumber component in the periodic direction equal to zero. Indeed, at a low Mach number, the wavenumbers with a non-zero component in the spanwise direction have a negligible impact on the far-field noise radiated in the mid-span plane, if the airfoil span is at least three times greater than its chord [19,29]. However, when the Mach number increases, the modes whose wavenumber component in the spanwise direction is non-zero become important for the far-field noise calculation.
The aim of the present work is to eliminate the spurious noise triggered by the spanwise non-periodicity of three-dimensional synthetic turbulence. In this paper, synthetic turbulence is generated using the RFM method inside Code_Safari [30], an EDF R&D code. It enables one to carry out relaxation filtering LES [31,32] using high-order finite-difference schemes and selective filters. In the first step, the RFM method is used inside Code_Safari to generate divergence-free turbulence at the inlet. This method was especially chosen because of its straightforward capability to impose a desired turbulent spectrum, its simplicity of implementation and its ability to synchronize the turbulent fluctuations with the mean flow [18]. In the second step, an extension of the RFM method called RFM-P is proposed. The RFM-P approach forces the random Fourier modes to respect the periodicity condition in the lateral direction, meaning that less spurious noise is generated. In the context of airfoil simulations, the lateral direction is chosen along the span. This paper is organized as follows. First, the flow solver and the synthetic turbulence methods are described in Section 2. Then, spatially decaying homogeneous isotropic turbulence (HIT) simulations using both the RFM method and the RFM-P extension are detailed in Section 3. Two geometrical configurations are proposed. The first one is used for a reference calculation, and the second one is built for the purpose of computational cost reduction. Since this study focuses on the turbulence generation and convection, the airfoil is not included in the present computations. Only the region between the inlet and a virtual airfoil is considered. The turbulence properties are analyzed to determine the RFM-P performance compared to the RFM one.

Flow Solver
In this work, an EDF R&D finite-difference code, called Code_Safari [30], is used to carry out DNC through LES. Considering a Newtonian fluid and a Cartesian mesh x = (x, y, z) transformed into a curvilinear one ξ = (ξ, η, ζ), the compressible equations to be solved are ∂ ∂tÛ where t is the time,Û = U/J with U = (ρ, ρu, ρe) T the conservative variables vector and J is the curvilinear transformation matrix Jacobian. The density is ρ, u = (u, v, w) T denotes the velocity vector in the Cartesian coordinates, and e is the total specific energy. Here, the equation of state is assumed to be the ideal gas law. Non-viscous fluxes are F ξ , F η , F ζ , while the viscous ones are F ν ξ , F ν η , F ν ζ . Their expressions can be found in Daude et al. [30]. To avoid the dissipation and dispersion of acoustic waves, the spatial discretization is carried out by an optimized 4th-order 11-point scheme [33], and the time advancement is realized by an optimized 2nd-order 6-step Runge-Kutta scheme [33]. From the first to the third points in the orthogonal direction with respect to the boundaries, the Tam and Webb dispersion-relation-preserving (DRP) [34] backward scheme is used for the spatial discretization. For the fourth and the fifth points, the DRP centered scheme is chosen. The dissipation of the non-resolved scales and of the grid-to-grid oscillations is ensured by a 6th-order explicit selective filter optimized on 11 points instead of subgrid models [35]. In the orthogonal direction with respect to the boundaries, filters based on centered decreasing-order schemes are used for numerical robustness from the first to the fifth points. The DRP centered scheme is applied for the fourth and fifth points, and standard schemes are chosen from the first to the third points. In this paper, outflow boundary conditions designate the Bogey and Bailly boundary conditions [36].
The Tam and Dong inflow boundary conditions [37] are implemented on 5 points in order to minimize the spurious noise when an inlet disturbance is introduced [12]. These conditions are given by where p is the fluctuating pressure assumed equal to the acoustic pressure, ρ and u are the fluctuations of density and velocity, r is the radius in the spherical coordinates, V g is the phase velocity and (ρ i , u i , p i ) are the desired fluctuating quantities at the inlet. In Equation (2), ρ i and p i are set to zero since their expressions are unknown, although this is a source of spurious noise, as noticed by Gloerfelt and Robinet [12]. As it will be seen in Section 2.2.1, u i is analytically known so that the right-hand side term of Equation (2) can be semi-analytically computed before being introduced in the equation, the Jacobian matrix between the Cartesian and the curvilinear grid being numerically computed. Finally, Equation (2) is computed at each substep of the Runge-Kutta algorithm to reduce the emitted spurious noise [12,37]. Code_Safari is parallelized using the OpenMPI Message Passing Interface (MPI) library. Periodic boundary conditions are treated as simple MPI communications using ghost points. Overset meshes with 8th-order interpolation are built with Overture [38], developed in the Lawrence Livermore National Laboratory. Details about the computational platform and libraries are given in Section 3.1.4.

Inflow Synthetic Turbulence
This section describes the synthetic turbulence algorithms implemented in Code_Safari. Firstly, the classical random Fourier modes (RFM) method is described in Section 2.2.1. Secondly, the RFM-P extension is described in Section 2.2.2; it takes into account the lateral periodic boundary conditions when computing the turbulence wavenumbers.

Classical Random Fourier Modes Method
The RFM method, based on [39], builds an incompressible turbulent velocity field by adding up N Fourier modes. Figure 1 schematically shows one mode and gives notations for the local spherical base (e r n , e θ n , e φ n ).  The turbulent velocity field u can be written as where • denotes the scalar product, U x is the mean flow field,ũ n is the nth mode amplitude, ψ n is the nth mode random phase, σ n is the nth mode direction and k n = k nx , k ny , k nz is the nth mode wavenumber. To limit spurious noise, the incompressibility condition ∇ • u = 0 is used, that leads to k n • σ n = 0.
Thus, the vector σ n must be in the (e θ n , e φ n ) plane. Furthermore, Kim and Haeri [22] and Sescu and Hixon [18] showed that the mean flow and the disturbances must be synchronized to generate minimum spurious noise. This condition means that the angular frequency of the disturbance must be written as k n • U, as it is written in Equation (3). Other details about the RFM method can be found in Appendix A. The RFM method generates a frozen divergence-free isotropic homogeneous turbulence with imposed turbulent kinetic energy k t and a longitudinal integral length scale Λ f . The chosen spectrum is the modified von Kármán one [17], given by Equation (A6). The wavenumber amplitudes are logarithmically distributed, and the wavenumber components are randomly computed using the spherical coordinate system angles, as seen in Appendix A.

Random Fourier Modes Method Tailored to Periodic Boundary Conditions
In this section, the RFM method is adapted to take into account spanwise periodicity boundary conditions for airfoil simulations. This extension will be called RFM-P in the following, for making a distinction from the RFM method. As the wavenumber components are randomly chosen to respect isotropy in the classical RFM method, spurious noise can be generated when using periodic boundary conditions, as it will be seen in Section 3.2.1. Indeed, a discontinuity of the turbulent velocity is present between the two opposite periodic boundaries. To enforce periodicity, the inflow turbulence is modified by the solver, leading to significant spurious pressure, vorticity and entropy fluctuations. To correct this issue, the spanwise component of the wavenumbers is first chosen in the RFM-P extension as a multiple of the smallest positive wavenumber component respecting the periodicity. Thus, no modification of the inflow turbulence is needed, meaning that the spurious fluctuations disappear. If z denotes the periodic direction, then the only possible values for the z-component of the wavenumbers are where L z is the computation domain length in the z direction and N z is the number of admissible positive wavenumber components along z. It is defined as where E(•) is the integer part and k Max z is the maximum wavenumber that the mesh can accept, based on the smallest mesh size and on the numerical scheme employed. For example, when using the optimized finite-difference 11-point scheme of Bogey and Bailly [33], this wavenumber must be set to k Max z = 2π 4.6∆ z to accurately resolve every injected wave, where ∆ z is the mesh size in the z-direction. The same notation is adopted for the maximum wavenumbers k Max x and k Max y along x and y. Then, the RFM-P algorithm computes a set of N targeted wavenumber amplitudes k Search n with a logarithmic distribution following Equation (A3). In the second step, randomly chosen wavenumber components k nx and k ny are combined with the admissible values of k nz given by Equation (5) to retrieve the targeted modes. The flow chart for the RFM-P approach is described in Figure 2. Note that the index p is used in this figure, such that the search base corresponding to the wanted wavenumber is written as k Search p . The pseudo-random wavenumbers are generated using the intrinsic function RANDOM_NUMBER of Fortran Intel compilers. The seed is initialized using the system clock.
The first step of the algorithm consists in generating a maximum number of modes with a non-zero z-component N k z =0 , using the values computed by Equation (5). Most of the time, it is impossible to obtain N k z =0 = N. Thus, in the second step, a certain number of modes N k z =0 with k nz equal to zero are generated to reach N modes. Eventually, there will be N k z =0 three-dimensional modes and N k z =0 two-dimensional modes. Furthermore, to ensure a maximum of three-dimensional modes, at least one mode must be found for each possible non-zero z-component from Equation (5) or the algorithm will stop. Several iterations (Iteration in Figure 2) can be needed, but only one is usually sufficient so that N iterMax in Figure 2 can be set to ten. Note also that it is not possible to impose a strict number of modes using the RFM-P extension. Indeed, it becomes difficult to obtain the low-wavenumber amplitudes when the z-component is enforced because it means that the two other wavenumber components must be very small. Thus, an effective number of modes N Eff exists. To be sure that N Eff is close enough to N, several iterations (denoted 2DIterationMax in Figure 2) are needed when adding two-dimensional modes. This number of iterations depends on N and on the mesh size. Typically, for the present simulations and N = 200 modes, 2DIterationMax = 500 is sufficient to obtain | N Eff −N N | ≤ 1 %. While the RFM method uses random θ n and φ n values for computing the wavenumber components, RFM-P starts from the components to compute the angles in a second step using θ n = arccos k nz k n and,  The other quantities such as α n , σ n , ψ n andũ n are computed the same way as in the RFM method described in Appendix A, and the final turbulent fluctuations are calculated using Equation (3). Contrary to the RFM method that is built to respect isotropy, it is clear that the RFM-P approach breaks this isotropy condition by selecting a discrete set of wavenumber values in the periodic direction. This lack of isotropy is illustrated in Figure 3 where the wavenumber distribution is represented for both RFM and RFM-P approaches using a domain length in the periodic direction L z = 0.25 m, a search base of N = 200 modes distributed between 25 and 500 m −1 and k Max The number of modes N Eff = 198 is obtained, with N k z =0 = 114 and N k z =0 = 86. While the RFM method distributes the modes randomly, the RFM-P extension concentrates the modes along the axis k x = k y = 0. It can also be noticed that many two-dimensional structures (k z = 0) are obtained. This constraint is, however, thought to be a good approximation when taking a small slice of an airfoil, since most of the structures along z will be aerodynamically seen as two-dimensional structures. One can notice that even with the RFM method, the largest structures will be two-dimensional due to the small span length (2Λ f ), meaning that the isotropy is also broken. When computing leading edge noise in the acoustic far-field, Clair et al. [19] and Atassi [29] explained that only the modes following the relation k nz ≤ k nx M √ 1−M 2 will be cut on for an observer in the mid-span plane if the span is at least three times greater than the chord. The mean flow Mach number is denoted M = U x /c, where c is the speed of sound and U x is the mean flow velocity along x. At a low Mach number, this means that only the modes with a low-wavenumber component in the periodic direction will play a role in the leading edge noise. Thus, taking only two-dimensional structures as Clair et al. [19] did could be a good approximation. However, when M increases, more and more modes will be cut on, meaning that it can be necessary to inject three-dimensional structures to correctly compute the leading edge noise in the mid-span plane, as conducted in the RFM-P approach.

Simulations of Spatially Decaying Homogeneous Isotropic Turbulence
In this section, the RFM-P extension is compared to the RFM method, using a case of spatially decaying homogeneous isotropic turbulence (HIT). First, the numerical configurations and the different simulations are described, along with the processing methods. In the second step, the results of both approaches are compared. The aim is to study the turbulence properties between the inlet plane and the close proximity of the airfoil leading edge. Therefore, only the domain between the inlet and the airfoil leading edge is considered in the following. Two meshing configurations are studied. The first one, called Single Grid, serves as a reference since it is a simple box with a constant mesh size. The second configuration, called Multi Grids, is a modified version of Single Grid in order to reduce the computational cost. Details are given in Section 3.1.1. A comparison of the turbulence properties between both configurations and both synthetic turbulence generation approaches will be performed. Since the turbulence is stochastic, the present results may depend on the particular turbulence realization. Thus, the results of two different realizations are presented for the Single Grid configuration only. Some data and results of this section are available in the supplementary materials at the end of the conclusion.

Configurations
The inflow turbulence is based on Comte-Bellot and Corrsin's [40] results taken at their first measurement plane. They give the longitudinal integral length scale Λ f and the turbulent kinetic energy k t , whose values will be given in Section 3.1.2. These two quantities drive the configuration numerical setups. The Single Grid configuration is shown in Figure 4 and is first composed of a main grid, called the convection grid, with a uniform mesh size ∆. Then, at the end of the domain, a sponge zone is added to damp turbulent structures before they reach the exit boundary. The boundary conditions, described in Section 2.1, are also presented in this figure: IBC stands for inflow boundary condition, OBC for outflow boundary condition and PERBC for periodic boundary condition. The dimensions of the computational domain are based on similar configurations found in the literature. The domain length along x is set to L x = 9.5c, where c is the virtual airfoil chord. Indeed, Eljack [41] or Thomareis and Papadakis [42] showed by some convergence tests on a NACA0012 that this distance, between the inlet and the airfoil leading edge, is sufficient to have convergence on the pressure and friction coefficients, even when an LSB is present. The sponge zone length along x is equal to 1.5c, beginning at 9.5c. Similarly, the L y computational domain length, in the chord-normal direction, is set to 20c. Finally, the airfoil span is set to L z = 0.2c to limit the computational cost. The chord is set to c = 10Λ f , as it can been found in the experiments of Wang et al. [3] or Moreau et al. [43]. The data transfer between the convection grid and the sponge zone is ensured by an 8th-order interpolation [30]. The mesh size ∆ for turbulence convection up to the airfoil is set to ∆ = λ g /1.1 to have enough resolved structures, where λ g is the transverse Taylor micro-scale of the injected turbulence. Then, to decrease the computational cost, the Multi Grids configuration is proposed. The motivation behind this configuration is to reduce the size of the convection grid along the y axis, as shown in Figure 5a where an airfoil is drawn only for clarity purposes. The remaining computational domain is composed of a generation grid at the inlet with a mesh size equal to ∆ and of two acoustic grids with a mesh size equal to 2∆, since no turbulence is needed in these areas. The acoustic grids begin at 10Λ f downstream of the inlet plane to avoid an interaction with the generation grid. The turbulence is convected from the inlet plane towards the airfoil through the convection grid, while the smallest scales are dissipated when passing through the acoustic grids. As no airfoil is needed in the present study, Figure 5b displays the geometry considered in the present work. The configuration is basically the same but a sponge zone is added at the computational domain exit as in the Single Grid configuration. Boundary conditions and dimensions for acoustic and convection grids are given in the figure. The convection grid has a height of 40Λ f . Global dimensions are the same as those for the Single Grid configuration. In the acoustic grids, the small structures whose wavelengths respect λ < 4.6∆ are numerically damped by the filter during their convection by the mean flow. The advantage of this second configuration is that it halves the number of calculation points. Details such as the total number of points for both configurations are given in Table 1. Table 1. Geometrical and mesh data for both Single Grid and Multi Grids configurations. Λ f stands for the longitudinal integral length scale and λ g for the transverse Taylor micro-scale of the turbulence.

Number of Interpolation Points
Single

Simulation Parameters
For each configuration, three simulations are carried out. The first simulation is conducted with the RFM method described in Section 2.2.1, the second simulation with the RFM-P extension described in Section 2.2.2 and the third one with the RFM-P approach, but only with two-dimensional structures, meaning that k z = 0 for every mode. The simulations are respectively named RFM, RFM-P3D and RFM-P2D. The physical turbulence parameters are those measured by Comte-Bellot and Corrsin [40] in their wind tunnel first plane, that are k t = 0.0739 m 2 ·s −2 , Λ f = 0.024 m and λ g = 0.0048 m. This gives a Reynolds number of Re λ g = λ g √ 2/3k t ν = 72, where ν is the kinematic viscosity. This low Reynolds number could be problematic for subgrid models, particularly for the classical models which assume a large separation between small and large turbulent scales. The numerical filtering used in the present work has the advantage of overcoming this limitation, meaning that the mesh size only determines the spectral range that is resolved [32].
The modified von Kármán spectrum is discretized with N = 200 modes between 20 and 500 m −1 since the energy is concentrated in this wavenumber range. In addition, for the RFM-P extension, the maximum wavenumber in each direction is chosen according to the accurate limit of resolution of the 4th-order 11-point scheme such that k Max 4.6∆ = 313 m −1 . The accurate limit is a little bit exceeded so that more modes can be introduced, but it is important to note that all the modes are still resolved by the scheme and the mesh.
The Mach number is set to M = U x /c = 0.1, giving a turbulent intensity of , where t * L x is the normalized convection time based on L x and U x , ∆ t is the time step and N ite = 7.14 × 10 4 is the number of computation iterations. The numerical filtering strength is set to σ ν = 0.2, as in Bogey and Bailly [44]. As already seen in Section 2.2.2, an effective number of modes N Eff is obtained, whose values are given in Table 2, as well as the number of three-dimensional structures N 3D in the RFM-P3D simulations. The effective number of modes is close to the target value, with a maximum relative error of 2%.
Finally, the three simulations of the Single Grid configuration are conducted for another turbulence realization to study the influence of different mode draws, the second realization being independent from the first one. The corresponding simulation names end with _Realization2 in Table 2.

Probes and Processing Method
In this section, the sets of probes used to compute the resolved turbulent kinetic energy and the one-dimensional spectra are described. These quantities are calculated using data on the last 2t * L x . The resolved turbulent kinetic energy k r t along x is studied with a line of probes placed at y = z = 0 and with a 3∆ space between two successive probes. The longitudinal spectral power density E vv k y of the y-velocity v is computed on x = 0 and x = 90Λ f . On each x position at z = 0, 106 probes are placed between y = −10Λ f and y = 10Λ f with a ∆ space between each successive probe. The minimum and maximum y positions are chosen so that enough probes can be placed without being in the interpolation areas between convection and acoustic grids. The periodogram method is used to compute the spectra using all 106 points on each x position. Thus, these spectra have a spectral resolution of 13 m −1 and a sampling wavenumber of 1400 m −1 . They are computed at each iteration of the last 2t * L x and averaged on this time window.

Computation Platform and MPI Decomposition
Simulations are carried out on Occigen, a supercomputer belonging to the CINES (French National Computing Center for Higher Education). Only Haswell ES-2690V3@2.6Ghz processors are used with an Infiniband FDR 56 Gbit/s communication bus. The 2.0.1 OpenMPI library for Intel compilers is employed. The distribution of cores in simulations is conducted so that each core has roughly the same number of calculation points in each direction, both in Single Grid and Multi Grids configurations, except in the z direction where there are not enough calculation points to cut the domain into several processes. For this purpose, the Single Grid simulations are run with 180 cores spread on 8 nodes of 24 cores, while the Multi Grids simulations are carried out on 140 cores spread on 6 nodes of 24 cores. Thus, it is important to notice that the number of cores and their distribution on the grids are not the same between Single Grid and Multi Grids configurations for core layout optimization. For each configuration, the number of points per core in the three directions is given in Table 3. The numbers of interpolation points, due to the interface between acoustic and convection grids and the sponge zone, are given in Table 1 for both Single Grid and Multi Grids configurations. As the calculation point number decreases with the Multi Grids configuration, the number of interpolation points increases, which must be taken into account when estimating the computational cost. Table 3. Number of points per core for each configuration in every direction (see Figures 4 and 5).

Single Grid
Multi Grids   Similarly, the corresponding fluctuating pressure fields at z = 0 are shown in Figure 7. There is a high decrease in the spurious noise between RFM and RFM-P3D simulations. It is worth mentioning that the pressure is scaled between −50 and 50 Pa for the RFM simulations, while it is scaled between −0.5 and 0.5 Pa for the RFM-P3D simulations. This represents the main result of the paper, in that the RFM-P extension does not significantly affect the velocity field, while it strongly reduces the spurious noise. When considering only the RFM-P3D simulations in Figure 7, it can be seen that the pressure fields have higher levels in the Multi Grids configuration than in the Single Grid one. This issue can be attributed to the overset method where interpolations take place. As explained by Sharan et al. [45], numerical dispersive waves can be reflected between the different grids of the mesh. These waves can grow with time and pollute the computational domain. For the Single Grid configuration, noise is mainly radiated from the sponge zone. This will not be a problem for airfoil computations, since the sponge zone will be far behind the airfoil and will be extended to a longer distance to make it more efficient.  Figure 8 shows the decay of the resolved turbulent kinetic energy k r t using probes described in Section 3.1.3. In LES, only a part of the spectrum is calculated; thus, the resolved kinetic energy k r t is smaller than the total turbulent kinetic energy k t . In the present work, k r t is estimated as 0.055 m 2 ·s −2 at the inlet, with k t = 0.0739 m 2 ·s −2 . From Figure 8, it can be seen that RFM-P simulations are close to this value for both configurations, but the RFM simulations overestimate it. As the RFM method is not suitable for periodic boundary conditions, the solver must strongly adapt velocity and pressure fields to enforce the periodicity which causes the generation of high-amplitude undesired waves. It is important to note that the velocity fluctuations are not directly superimposed to the mean flow but computed through Equation (2). Thus, they can be modified right at the inlet if they do not respect the periodicity condition. This explain why the resolved turbulent kinetic energy takes different values on x = 0 depending on the simulation. Additionally, this issue is present because the span is very small (only 2Λ f , 11 computational points), that means that all the points of the domain are affected by the periodicity condition, contrary to the study [46], where the span was 5Λ f and the turbulent kinetic energy of every simulation started from the same value. However, the high-amplitude spurious noise problem was present even in this configuration.  The extreme values of the empirical decay exponent for grid-generated turbulence [47] are also plotted in Figure 8, showing that the spatial decay is really small in the simulations. For RFM-P simulations, 25% of the initial energy is lost when reaching the airfoil leading edge. This spatial loss should increase if the mean velocity decreases, giving more time for dissipation. Indeed, considering dk t dt = − [48], where is the turbulent kinetic energy dissipation rate, and assuming a frozen turbulence, the evolution of k t with respect to x is

Turbulent Kinetic Energy Evolution
If the turbulence temporal decay does not depend on the Mach number, its spatial decay does, as seen above. Here, M = 0.1 is sufficiently low to be able to capture some changes in the turbulence statistics during its convection along the domain length L x . Indeed, the convection time at a distance of L x = 95Λ f can be expressed as where τ η = ν is the Kolmogorov time scale and τ Λ f = Λ f √ 2/3k t is the decaying time for the largest scales. Thus, only the small scales have the time to evolve in this particular situation, meaning that the small spatial decay observed could only be a consequence of the relatively high Mach number of this particular situation in comparison with the Mach number of Comte-Bellot and Corrsin's experiments [40]. The ratio between the numerical and experimental turbulent kinetic energy spatial decays can be expressed as M Simu. = 0.37. Additionally, the same small energy decay was already observed in a previous study [46] using the same numerical tools but with a larger domain length L x = 250Λ f . As a consequence, the convection time was comparable to the large-scale decaying time. Thus, even when the large scales have the time to evolve during their convection, most of the dissipation is due to the small scales, meaning that the turbulent kinetic energy should remain nearly the same at a sufficiently long distance to reach the airfoil leading edge. The fact that the numerical decay is not in agreement with the empirical one comes from the synthetic turbulence generation method. Indeed, the third-order moment is not reproduced and this moment is responsible for the energy transfer from the large scales to the smallest, as explained by Dhamankar et al. [11] in their review. However, for airfoil simulations, this disagreement can be beneficial since the turbulent kinetic energy at the leading edge can be the same as or very close to the energy set at the inlet. Finally, Figure 8 shows that the energy using the Multi Grids configuration increases slightly compared to the one of the Single Grid configuration for RFM-P simulations. This increase is of 7% only. To explain it, Figure 9 compares the turbulent kinetic energy decay obtained for two different turbulence realizations using the Single Grid configuration. One can notice that the results between the two realizations are different and that the increase in the kinetic energy in the second realization is similar to the one obtained in Figure 8 with the Multi Grids configuration. The overestimation with the RFM method is still present compared to the other simulations. Thus, the difference in the kinetic energy between the Multi Grids and Single Grid configurations in Figure 8 is not necessarily the consequence of the different configurations but of the different realizations. The increase in the energy in Figure 8 for the Multi Grids-RFM-P2D simulation after 20Λ f is also due to the realization variability since the same behavior is visible in Figure 9 for Realization 2 in the Single Grid-RFM-P2D simulation. More realizations would be desirable to compute an uncertainty error in the method, but that would also be very computationally expensive.

One-Dimensional Spectra
One-dimensional spectra E vv along the y direction are obtained using probes described in Section 3.1.3. The top plots of Figure 10 show these spectra for the Single Grid and Multi Grids simulations at the box inlet, x = 0, while the bottom plots show the same spectra at x = 90Λ f . The resolution limit plotted as a vertical line corresponds to the smallest wavenumber that can be accurately calculated by the numerical scheme with the chosen mesh size. First, by comparing Figure 10a and Figure 10b, similar results are obtained between configurations Single Grid and Multi Grids at the inlet. The simulations are also compared to the measurements of Comte-Bellot and Corrsin [40], showing that both methods are capable of reproducing the desired experimental spectrum. However, the RFM method seems to overestimate the energy at the lowest wavenumber which can be explained by the same reason given for the kinetic energy overestimation in Section 3.2.2. The RFM-P extension tends to underestimate the energy at the lowest wavenumber. This could be due to the small domain length L z = 2Λ f , that may not be sufficient to introduce the largest scales. Indeed, using the same turbulence parameters, it was observed in [46] that the energy was correctly computed at the lowest wavenumber when the domain length along z was more than twice the present one. A different behavior is observed at intermediate wavenumbers (30 m −1 ≤ k y ≤ 60 m −1 ). On the one hand, the RFM and RFM-P3D approaches are close to each other. On the other hand, the RFM-P2D simulation shows lower energy levels compared to the two other simulations at x = 0. An explanation for this observation will be proposed later in this section.
Spectra right at the inlet are not the best criteria to assess the relevance of the RFM-P extension. Indeed, peaks in the spectra are present because of the discrete distribution of the energy between modes. It is necessary to look at the spectra near the end of the computational domain since this will be the location of the airfoil leading edge in future simulations. The adaptation distance between the inlet and the outlet allows the energy transfer between modes to be carried out by the Navier-Stokes equations. Spectra in Figure 10c,d are indeed smoother, and the three simulations give similar results. The high-wavenumber energy has been numerically dissipated between x = 0 and x = 90Λ f because of a lack of grid resolution, as the red dashed vertical line shows. Retrieving high-wavenumber energy is only possible by decreasing the mesh size. The RFM-P simulation with two-dimensional structures seems to be, again, less energetic than the other simulations at low wavenumbers, but only for the Single Grid configuration. Note that the Comte-Bellot and Corrsin [40] measurements correspond to the inlet plane (x = 0) for all plots of Figure 10. It is not possible to compare numerical spectra at x = 90Λ f with measurements at the same location because the mean flow velocity is different from the one used in the experiments.
Finally, Figure 11 shows the results between both planes for the Single Grid configuration and only for the RFM-P extension because the RFM method was shown to overestimate the low-wavenumber energy due to the periodicity problem. It can be seen that the energy at low wavenumbers is nearly the same between both positions, while the energy at high wavenumbers is dissipated. Spectra become smoother and non-resolved scales are attenuated, that explains the decay of the turbulent kinetic energy in Figure 8. It means that the large structures are purely convected and that the transfer towards the small scales does not take place, as already explained. Notice that smaller scales could be added to the simulation by reducing the mesh size. Thus, a full range of scales could be obtained, without dissipation between the inlet and the airfoil. As conducted for the kinetic energy, the results for a second turbulence realization are compared with those of the first realization for the Single Grid configuration in Figure 12. On the first plane, results between realizations seem to be quite similar for low wavenumbers only. The difference visible at high wavenumbers comes from the different distributions of modes between realizations. On the second plane, the difference between the two realizations can hardly be seen, and the RFM-P2D simulation still contains less energy at low wavenumbers. Looking at Figures 10-12 for x = 0, the RFM-P2D simulation is less energetic for k y ≤ 70 m −1 . This is a consequence of the algorithm of Figure 2 used to generate the wavenumbers. Indeed, setting k z = 0 allows greater values of k y to be reached in order to fill all the wavenumbers k = k 2 x + k 2 y + k 2 z of the desired spectrum. Thus, more energy is injected at higher values of k y compared to the RFM-P3D simulation where k z = 0. In conclusion, regarding spectra, the RFM-P3D simulation gives better results than the RFM one because there is no energy overestimation at the inlet for the lowest wavenumbers.

Computational Cost
In this section, the computational costs of all simulations are presented. Section 3.1.4 should be kept in mind while reading the present results since the distribution and number of cores are not the same between configurations, meaning that the performance's results are approximate. Table 4 shows the computational cost in core hours for each simulation. The Multi Grids configuration seems to be less expensive which is due to the smaller mesh number in comparison with the Single Grid configuration. A gain of 25% in the computational cost can be achieved, which can be useful for expensive airfoil simulations.

Conclusions
To compute airfoil leading edge noise and to trigger boundary layer instabilities using compressible large eddy simulations, synthetic turbulence compatible with periodic boundary conditions is needed. When the classical random Fourier modes method is used to generate synthetic turbulence, spurious noise is present because this method does not respect the periodicity condition in the spanwise direction, since the wavenumber components are randomly chosen. In the present work, an adaptation of the random Fourier modes method for periodic boundary conditions, called the RFM-P extension, was developed and compared with the RFM method in a case of spatially decaying homogeneous isotropic turbulence, using two-dimensional or three-dimensional modes. The RFM-P algorithm computes modes whose wavenumber component in the lateral direction is compatible with the periodic boundary condition. Two geometrical configurations were considered to study the spatial decay of turbulence using both approaches. The first configuration, called Single Grid, is a simple box with a constant mesh size, and the second one, called Multi Grids, reduces the computational cost by up to 25% by decreasing the mesh size where detailed turbulence is not needed. Results show that the RFM-P extension permits a better control of the turbulent kinetic energy compared to the RFM method, while the spectral shape is not disturbed. Furthermore, the decay of the turbulent kinetic energy is small compared to the empirical decay exponents. This decrease is only due to the dissipation of the smallest scales, while larger scales are unchanged and purely convected towards the outlet. The RFM method produces undesired high-amplitude waves generated at the inlet to force the turbulent fluctuations to be periodic. The RFM-P extension, developed in this work, is effective in getting rid of these waves. The drawback of the second configuration is the interpolation using the Chimera method, which causes additional spurious pressure fluctuations at the inlet compared to the constant mesh size configuration. In the future, the RFM-P approach will be tested in airfoil simulations to study the effect of inflow turbulence on the broadband noise emitted by an airfoil at various angles of attack.
Finally, the authors of the present work are aware of the RFM-P extension's limitations, inherited from the RFM method, such as the spatial and temporal periodicities of the fluctuations or the inability to introduce a turbulent field directly into the computational domain. Gea-Aguilera et al. [28] showed that the advanced digital filtering and RFM methods resulted in similar airfoil leading edge noise predictions, in a two-dimensional configuration. Therefore, a comparison between the RFM-P extension, the incompressible SEM of Kim and Haeri [22] and the advanced digital filtering method of Gea-Aguilera et al. [28] in a three-dimensional configuration could be worthwhile.

Acknowledgments:
The authors would like to acknowledge the French supercomputing center CINES (Centre Informatique National de l'Enseignement Supérieur) for the A0052A10311 allocation in the GENCI (Grand Equipement National de Calcul Intensif) framework.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: This appendix provides some details about the RFM method. The formulation used in the present work is based on Gloerfelt et al. [39]. To meet the condition k n • σ n = 0 given in Equation (4), the vector σ n must be in the (e θ n ,e φ n ) plane, as seen in Figure 1. It is described in this plane by the angle α n , so that σ n = cos α n e θ n + sin α n e φ n . (A1) The e θ n , e φ n and e r n are expressed in the Cartesian coordinate system by The wavenumber amplitude k n is logarithmically distributed for a good discretization of the energetic low wavenumbers using k n = exp[ln k 1 + (n − 1)∆(ln k n )], n = 1, 2, ..., N, where k 1 is the lowest wavenumber and ∆(ln k n ) is given by where k N is the highest wavenumber. The parameters to set are the wavenumbers k 1 and k N , and the number of modes N. To respect isotropy, the probability densities for α n , θ n and φ n are p(α n ) = 1 2π , p(θ n ) = sin(θ n ) 2 and p(φ n ) = 1 2π [39]. The random phase is chosen so that p(ψ n ) = 1 2π . The mode amplitudeũ n is computed using to obtain the imposed turbulent kinetic energy k t , where E is the desired turbulent spectrum and ∆k n is the step between two successive wavenumbers. In this work, E is the modified von Kármán spectrum, that depends on the desired turbulent kinetic energy k t and the longitudinal integral length scale Λ f . This spectrum is written as [17] E(k) = α s 2k t 3 (k/k e ) 4 k e 1 + (k/k e ) 2 17/6 exp −2 where k η is the Kolmogorov wavenumber and α s and k e are constants. These constants take the following values: as computed in [49] to respect both the imposed turbulent kinetic energy k t and the integral length scale Λ f .