T Rigall

B Cotté

P Lafon

Airfoil Noise Numerical Simulations with Direct Noise Computation and Hybrid Methods Using Inflow Synthetic Turbulence

The interaction of inflow turbulence with an airfoil is significant for low Reynolds number applications, as it generates leading edge noise, but also strongly modifies the airfoil boundary layer properties and thus the produced noise. Inflow turbulence impact still demands to be better understood and is also a numerical challenge. Large Eddy Simulations are here carried out with both incompressible and compressible solvers. The Synthetic Eddy Method and the Random Fourier Modes method, used to generate synthetic turbulence, are studied in this work. In order to eventually perform airfoil simulations, the second method is modified to take into account spanwise periodic boundary conditions, to reduce spurious noise level in compressible simulations. First, properties of the methods are studied with a spatially decaying turbulence academic case. Even if the decay is not well represented, the methods show satisfactory properties in terms of isotropy, homogeneity and one-dimensional spectra. Furthermore, spurious noise levels are reduced thanks to the modified method. Preliminary airfoil simulations are then carried out with inflow turbulence for flow validation in presence of a laminar separation bubble. The drag coefficient is in very good agreement with experimental results. However, there are some discrepancies between simulations and experiments on the lift coefficient, which is extremely sensitive to perturbations. These discrepancies are attributed to a lack of experimental turbulence characteristic information, such as isotropy or integral length scale and to too short simulations at high angles of attack. * The Einstein convention is used : subscript repetition implies summation, subscript between brackets avoids it.

I. Nomenclature

= number of acceptable non-zero positive wavenumber components along z * PhD Student, ENSTA Paristech, Institute of Mechanical Sciences and Industrial Applications, tommy.rigall@ensta.fr † Assistant Professor, ENSTA Paristech, Institute of Mechanical Sciences and Industrial Applications, benjamin.cotte@ensta.fr ‡ Research Engineer, EDF R&D, Institute of Mechanical Sciences and Industrial Applications, philippe.lafon@edf.fr 1 R i j = Reynolds stress tensor (m 2 .s -2 ) R N or m u x u x = velocity longitudinal autocorrelation normalised function along x Re c = Reynolds number based on the airfoil chord and the upstream velocity Re λ g = Reynolds number based on the transverse Taylor micro-scale and the turbulent kinetic energy U = velocity field (m.s -1 ) U BC = Comte-Bellot and Corrsin's mean flow velocity before the wind tunnel contraction (m.s -1 ) U x = mean flow velocity along x (m.s -1 ) V B = virtual box volume around the inlet (m 3 ) a i j = lower triangular matrix of the Reynolds stress tensor Cholesky decomposition (m.s -1 ) (e r n , e θ n , e φ n ) = local spherical base for the Fourier modes f j = shape function in the j direction (m -3/2 ) f Λ f = shape function for present turbulence associated to Λ f (m -3/2 ) I t = turbulent intensity k e = wavenumber of energetic structures (m -1 ) k E xact n = wavenumber amplitude of the n th desired mode (m -1 ) k η = Kolmogorov wavenumber (m -1 ) k

x,y,z M ax = maximum wavenumber value acceptable in the x, y or z directions (m -1 ) k n = wavenumber of the n th mode (m -1 ) k t = mean turbulent kinetic energy (m 2 .s -2 ) k y , k z = wavenumber component in the y or z directions (m -1 ) k 1 , k N = minimum and maximum wavenumber amplitudes for the spectrum discretisation t = physical time (s) t * = kinematic viscosity (m 2 .s -1 ) ξ + , η+, ζ+ = non-dimensional wall mesh sizes

II. Introduction

Broadband noise from wind turbines benefits from a large range of studies [START_REF] Zhu | Modeling of Aerodynamically Generated Noise From Wind Turbines[END_REF][START_REF] Oerlemans | Prediction of wind turbine noise and validation against experiments[END_REF][START_REF] Tian | Wind turbine noise modeling based on Amiet's theory: Effects of wind shear and atmospheric turbulence[END_REF] but there are few models for the boundary layer separation noise or the stall noise prediction, which is stronger and can be amplitude modulated [START_REF] Oerlemans | An explanation for enhanced amplitude modulation of wind turbine noise[END_REF]. The appearance of domestic small wind turbines or the drone proliferation also rise the emitted noise question. These low Reynolds number systems become numerically affordable for wall resolved Large Eddy Simulations (LES).

For drones or small wind turbines working at low Reynolds numbers, typically Re ≤ 5 × 10 5 [START_REF] Alam | The ultra-low Reynolds number airfoil wake[END_REF], a Laminar Separation Bubble (LSB) can occur on wings, that is difficult to characterise both experimentally and numerically. The LSB structure is composed of a high recirculation area giving birth to large scale eddies on the suction side and in the wake of the wing. LSB have been largely studied experimentally and numerically. Some works [START_REF] Alam | The ultra-low Reynolds number airfoil wake[END_REF] deal with the lift increase followed by a stall, triggered by the LSB burst. Other authors address the acoustic retroaction loop between the LSB wake and its separation point [START_REF] Paterson | Vortex noise of isolated aircraft[END_REF] or the trailing edge impact on the LSB dynamics [START_REF] Thomareis | Effect of trailing edge shape on the separated flow characteristics around an airfoil at low Reynolds number: A numerical study[END_REF]. Finally, the considerable effects of the upstream turbulence on the separation, transition and reattachment points of the LSB are highlighted in some experiments [START_REF] Wang | Turbulent intensity and Reynolds number effects on an airfoil at low Reynolds numbers[END_REF]. All these studies bring to light that the LSB is very sensitive to perturbations and it is essential to better understand it, since such a structure can change the aerodynamic behavior of an airfoil. Numerical simulations can help a lot to understand the underlying mechanisms [START_REF] Thomareis | Effect of trailing edge shape on the separated flow characteristics around an airfoil at low Reynolds number: A numerical study[END_REF]. Furthermore, inflow turbulence generates leading edge noise and modifies the boundary layer properties which, in turn, changes the global noise produced. Inflow turbulence is thus essential for both acoustic and aerodynamic reasons.

In numerical simulations such as LES, inflow turbulence is also a major challenge if we want to reproduce experimental turbulence without getting too much spurious noise in Direct Noise Calculation (DNC). The differents categories and methods of inflow turbulence have been reviewed in Dhamankar et al. [START_REF] Dhamankar | Overview of Turbulent Inflow Boundary Conditions for Large-Eddy Simulations[END_REF]. While library-based turbulence methods are efficient in terms of statistics reproductivity, they can only be used for specific applications and can be very expensive. The recycling methods avoid to have a long distance for the adaptation of the turbulence to the solver but present the drawback of generating spurious modes in the flow [START_REF] Dhamankar | Overview of Turbulent Inflow Boundary Conditions for Large-Eddy Simulations[END_REF]. The last category is the synthetic turbulence generator which includes many different methods. They are still missing correct turbulence statistics but they have the advantage to be directly embedded in solvers, without running any other simulation in parallel. Recently, Gea-Aguilera et al. [START_REF] Gea-Aguilera | Synthetic turbulence methods for computational aeroacoustic simulations of leading edge noise[END_REF] used an advanced digital filtering approach with success to compute airfoil leading-edge noise with linearised Euler equations. Batten et al. [START_REF] Batten | Interfacing statistical turbulence closures with Large-Eddy Simulation[END_REF] used spectral approaches. They continued the work of Smirnov et al. [START_REF] Smirnov | Random flow generation technique for Large Eddy Simulations and particle-dynamics modeling[END_REF] by adding anisotropy and inhomogeneity properties to the first work of Kraichnan [START_REF] Kraichnan | Diffusion by a random velocity field[END_REF]. This method, called Random Fourier Modes (RFM), consists in building turbulence by summing Fourier components. The RFM method was used and simplified to compute leading edge noise of wavy airfoils with the Euler equations in Clair et al. [START_REF] Clair | Experimental and Numerical Investigation of Turbulence-Airfoil Noise Reduction Using Wavy Edges[END_REF]. Another approach has been proposed by Gloerfelt and Robinet [START_REF] Gloerfelt | Silent inflow condition for turbulent boundary layers[END_REF] who introduced unstable modes in a turbulent boundary layer, resulting in a very low spurious noise simulation. The SEM (Synthetic Eddy Method) [START_REF] Jarrin | A Synthetic-Eddy-Method for generating inflow conditions for large eddy simulations[END_REF] is another way to compute synthetic turbulence and has been used for channel flow simulations. This method has been several times modified [START_REF] Pamiès | Generation of synthetic turbulent inflow data for large eddy simulation of spatially evolving wall-bounded flows[END_REF] for turbulent channel flow, or for low noise condition [START_REF] Sescu | Toward low-noise synthetic turbulent inflow conditions for aeroacoustic calculations[END_REF] by making the resulting fluctuating velocity field incompressible. Later, Kim and Haeri [START_REF] Kim | An advanced synthetic eddy method for the computation of aerofoil-turbulence interaction noise[END_REF] extended this work to get a von Kármán spectrum by optimising the method with constraint parameters. They successfully applied it to the turbulence interaction noise on a wavy airfoil with the Euler equations. For DNC, silent synthetic turbulence is seeked but is still a tricky achievement. Indeed, in compressible codes, the velocity is often prescribed without modifying the thermodynamic variables, which forces the code to adapt the flow field with pressure fluctuations, as explained in Gloerfelt and Robinet [START_REF] Gloerfelt | Silent inflow condition for turbulent boundary layers[END_REF].

The present work aims to generate realistic and silent inflow synthetic turbulence in LES, in order to accurately capture LSB on airfoils as well as to compute airfoil noise. In this paper, two methods of synthetic turbulence are studied to check their properties. Preliminary LSB calculations on a NACA 0012 airfoil with inflow synthetic turbulence are then carried out. Synthetic turbulence is generated inside two EDF R&D codes. The first one, Code_Saturne, solves the incompressible Navier-Stokes equations with low order finite volumes while the second, Code_Safari, solves the compressible ones with high-order finite differences. The second code enables us to carry out DNC, i.e. to compute simultaneously flows and acoustics, but it is still highly computationally expensive for low Mach numbers and high Reynolds numbers. Incompressible simulations with Code_Saturne are a way to do hybrid calculations, in which flows are solved first and are then used as sources to compute the noise, thanks to acoustic analogies. This second approach has no limit at low Mach number flows which makes it complementary with DNC. The SEM is implemented in Code_Saturne and the RFM method is used in Code_Safari for its incompressibility property. A modified RFM method, called RFM-P, is proposed in this work, to be used with periodic boundary conditions in the context of airfoil simulations. It takes into account periodicity condition in the spanwise direction and forces the modes to respect this condition. Less spurious noise is generated but the isotropy is broken.

The paper is organised as follows. First, the flow solvers and the synthetic turbulence methods are described in Section III. Then, a spatially decaying Homogeneous Isotropic Turbulence (HIT) test case is carried out to assess the properties of the three methods in Section IV. Finally, preliminary results about airfoil simulations with LSB and inflow turbulence in Code_Saturne are shown in Section V. The lift and drag coefficients are compared to experimental results.

III. Numerical Methods : Flow Solvers and Inflow Synthetic Turbulence

A. Flow Solvers Two solvers are used in this work. Code_Saturne [START_REF] Archambeau | Code_Saturne : A Finite Volume Code for the computation of turbulent incompressible flows -Industrial Applications[END_REF] is a EDF-R&D collocated finite volume code which is here used to solve the incompressible Navier-Stokes equations by LES. The SIMPLEC algorithm ensures the pressure and velocity coupling. For LES, 2nd-order schemes are used both for time (Crank-Nicolson) and spatial discretisations. The role of non-resolved scales on resolved scales is played by functional sub-grid models.

The second code, Code_Safari [START_REF] Daude | A high-order finite-difference algorithm for direct computation of aerodynamic sound[END_REF], is a EDF-R&D finite difference code which is here used to solve the compressible Navier-Stokes equations by LES with high-order numerical schemes to be able to carry out DNC. To avoid dissipation and dispersion of acoustic waves, the time advancement is realised by an optimised 2 nd -order 6-step Runge-Kutta scheme [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF]. Likewise, the spatial discretisation is carried out by an optimised 4 th -order 11-point scheme [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF]. Unlike Code_Saturne, the dissipation of the non-resolved scales and of the grid to grid oscillations is ensured by a 6 th -order explicit selective filter optimised on 11 points [START_REF] Bogey | A shock-capturing methodology based on adaptative spatial filtering for high-order non-linear computations[END_REF]. The Bogey and Bailly outflow boundary conditions [START_REF] Bogey | Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases[END_REF] and Tam and Dong inflow boundary conditions [START_REF] Tam | Advances in numerical boundary conditions for computational aeroacoustics[END_REF] are implemented on 5 points in order to minimise the spurious noise when an inlet perturbation is introduced [START_REF] Gloerfelt | Silent inflow condition for turbulent boundary layers[END_REF] and to avoid any acoustic reflection at the boundaries. Close to the walls and the boundaries, centred decreasing order schemes are used both for spatial discretisation and numerical filtering for robustness. For Code_Safari, overset meshes with 8 th -order interpolation are built with Overture [START_REF] Henshaw | Mappings for Overture A Description of the Mapping Class and Documentation for Many Useful Mappings[END_REF], developped in the Lawrence Livermore National Laboratory.

B. Inflow Synthetic Turbulence

Synthetic turbulence is necessary to trigger boundary layer instabilities and to compute leading edge noise of an airfoil in a realistic way. Methods of turbulence generation used in each code are now described. In Code_Saturne, the Synthetic Eddy Method (SEM) developped by Jarrin et al. [START_REF] Jarrin | A Synthetic-Eddy-Method for generating inflow conditions for large eddy simulations[END_REF] is used. An illustration taken from the Code_Saturne user guide [START_REF]Code_Saturne version 5.0.8 practical user's guide[END_REF] of the SEM principle is shown on Fig. 1a). It consists in summing N coherent structures to form velocity fluctuations. First, an isotropic and normalised velocity fluctuation u N I at the spatial position x and at the physical time t is constructed as the sum of N structures following

SEM : Synthetic Eddy Method

u N I j (x, t) = V B N N i=1 i j f (j) (x -x i (t)) * , (1)
with i j taking randomly the values 1 or -1 and V B the virtual box volume around the inlet S of the computational domain, which are shown on Fig. 1a). More details about the virtual box can be found in Appendix VII.A. Finally, f (j) represents the shape function of the structures in each direction j. For the present work, this function should be the same in each direction to respect the isotropy condition and will be noted f Λ f from now on, as the function is associated to a turbulence with a longitudinal integral length scale Λ f . The centre of each structure is initially randomly generated inside the virtual box. Then, it is simply convected by the mean flow which is considered in the x direction in this work and is denoted by U c on Fig. 1a). Each time a centre exits the virtual box, a new structure is generated randomly along y and z on the virtual box inlet plane. Then, a Cholesky decomposition is used on the desired Reynolds stress tensor R i j = u i u j [START_REF] Lund | Generation of Turbulent Inflow Data for Spatially-Developing Boundary Layer Simulation[END_REF] to set the turbulent kinetic energy and a potential anisotropy to u N I , where u is the final fluctuation desired at the inlet. This fluctuation is computed by

u i x, t = a i j x, t u N I j x, t , (2)
where a i j is the lower triangular matrix of the R i j Cholesky decomposition. In the homogeneous and isotropic present framework, the matrix becomes a i j = 2 3 k t δ i j where k t is the known mean turbulent kinetic energy independent of time or space and δ i j is the Kronecker symbol.

For each calculation iteration, the velocity fluctuation u is computed on the inlet plane at x = x min using Eq. ( 2). In this work, we consider two shape functions to see their influence on the turbulence properties. The first one is a hat function calculated using

       f Λ f (x) = α 3 1 - | x | β 1 - | y | β 1 - | z | β if x, y, z ∈ [-β; β], = 0 otherwise, (3)
and the second one is a Gaussian function given by

         f Λ f (x) = α 3 exp - 1 2 x β 2 + y β 2 + z β 2 if x 2 + y 2 + z 2 ≤ 4β, = 0 otherwise, (4)
where, in both equations, the parameters α and β are determined as shown in Appendix VII.A in order to impose the values of k t and Λ f on the final generated turbulence. Thus, the SEM is designed to impose the homogeneity and isotropy properties as well as the turbulent kinetic energy and longitudinal integral length scale. A frozen turbulence is considered and the zero mean value is automatically ensured when the structures number N is sufficiently high to statistically cover the inlet plane. Furthermore, the correlation function, and consequently the turbulent spectrum, is given by the shape function. It means that the shape function should be specifically designed to give a desired spectrum.

RFM-P : Random Fourier Modes adapted to Periodic boundary conditions

In Code_Safari, two methods of synthetic turbulence are used. The first one is the classical RFM method and the second one is a modified RFM method, called RFM-P. It has been developped in the present work to cancel spurious noise due to periodicity conditions in numerical simulations. Both methods are now described.

Classical Random Fourier Modes method

The Random Fourier Modes (RFM) method has been initially introduced by Kraichnan [START_REF] Kraichnan | Diffusion by a random velocity field[END_REF]. It consists in summing N Fourier modes to form the velocity fluctuations. Following the formalism introduced in Bechara et al. [START_REF] Bechara | Stochastic Approach to Noise Modeling for Free turbulent flows[END_REF] and considering a mean flow U x along x, with u = U + u = U x e x + u the velocity field, the velocity fluctuation u is given by

u x, t = 2 N n=1 ũn cos k n • x + ψ n -k nx U x t σ n , (5)
where ũn , k n , ψ n and σ n are the amplitude, wavenumber, phase and direction of the n th mode and N is the number of modes.

It is important to note that the fluctuating velocity given by Eq. ( 5) is a solution of the linearised Euler equations; it is an advantage since Code_Safari uses the non-reflecting outflow conditions of Bogey and Bailly [START_REF] Bogey | Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases[END_REF] that are based on the linearised Euler equations. Then, to limit the noise generated by the turbulence, it is necessary to respect the incompressibility constraint ∇ • u = 0, as explained by Kim and Haeri [START_REF] Kim | An advanced synthetic eddy method for the computation of aerofoil-turbulence interaction noise[END_REF], which leads to the condition

k n • σ n = 0, (6)
giving a constraint to each mode direction σ n . Finally, the perturbation is convected by a mean flow and it has been shown by Sescu and Hixon [START_REF] Sescu | Toward low-noise synthetic turbulent inflow conditions for aeroacoustic calculations[END_REF] and by Kim and Haeri [START_REF] Kim | An advanced synthetic eddy method for the computation of aerofoil-turbulence interaction noise[END_REF] that the angular frequency of the perturbation must be equal to k n • U, as in Eq. ( 5). It means that the turbulent perturbation must be synchronised with the mean flow to generate minimum spurious noise. Figure 1(b) provides notations for the wavenumber k n and the direction σ n of the n th mode in spherical coordinates. Following Eq. ( 6), σ n must be in the (e θ n ,e φ n ) plane where the vectors e θ n and e φ n , completed by e r n , form the local spherical base. The angular position of σ n in the (e θ n ,e φ n ) plane is defined by the angle α n , so that σ n is given by

σ n = cos α n e θ n + sin α n e φ n . (7)
The e θ n , e φ n and e r n are expressed in the cartesian coordinate system by

e θ n = cos θ n cos φ n cos θ n sin φ n -sin θ n , e φ n = -sin φ n cos φ n 0 , e r n = sin θ n cos φ n sin θ n sin φ n cos θ n . (8)
The wavenumber k n can be written as k n = e r n k n , where k n is chosen with a logarithmic distribution in order to correctly discretise the energy in the low-wavenumber region. It follows the expression

k n = exp (ln k 1 + (n -1)∆ (ln k n )) , n ∈ {1, ..., N }, (9)
where k 1 is the lowest wavenumber, to be set as a parameter, and ∆ (ln k n ) is given by

∆ (ln k n ) = ln (k N ) -ln (k 1 ) N -1 , (10)
with k N the highest wavenumber, to be set as a parameter. Isotropy leads to a uniform probability density for α n and ψ n given by p(α n ) = 1 2π and p(ψ n ) = 1 2π , to avoid any favoured direction for σ n and to have a random phase. Likewise, the probability density for θ n and φ n are p(θ n ) = sin(θ n ) 2 and p(φ n ) = 1 2π [START_REF] Gloerfelt | Generation of inflow turbulence for aeroacoustic applications[END_REF], which is established by considering an elementary surface of a k n radius sphere. Finally, each mode amplitude is calculated by

ũn = E(k n )∆k n , (11)
where ∆k n is the step between two successive wavenumbers and E(k n ) is the energy associated to the n th mode. The way to obtain Eq. ( 11) is given in Appendix VII.B. The spectrum E chosen here is the von Kármán-Saffman [START_REF] Bailly | Turbulence[END_REF] one, which depends on the desired values of k t and Λ f , and is also given in Appendix VII.B. The link between the spectrum and both the turbulent kinetic energy and the integral length scale is also given in Appendix VII.B. Thus, the RFM is designed to impose the homogeneity and isotropy properties as well as the turbulent kinetic energy, the longitudinal integral length scale and the turbulent spectrum. A frozen turbulence is considered and the zero mean value is automatically ensured.

Random Fourier Modes for Periodic boundary conditions method

The second method, RFM-P, is now described. This method aims at reducing the spurious noise produced at the inlet when the turbulence is injected and when periodic conditions are set to opposite faces in a simulation. Indeed, the RFM generates random wavenumbers so that its components in the periodicity direction do not respect the periodicity condition in general. The code has to adapt the velocity field, which leads to large pressure fluctuations inside the computational domain. The RFM-P is thus designed to take into account this periodicity condition inside the turbulent velocity construction. The method lies on the same physical background as the RFM one but the algorithm is changed and is now explained. We consider in the following that the only periodic direction is along z. The algorithm can be straightforwardly extended to other cases.

First, a logarithmic distribution for the wavenumbers is still desired. We start with an exact set of N E xact modes amplitudes, k E xact n , which are calculated with Eq. ( 9). Then, the only non-zero acceptable values for the wavenumbers component in the periodic z direction, k nz , follow the equation

k nz = ± 2πn L z , (12)
and are generated, where n is an integer between 1 and N z Modes , the number of positive non-zero acceptable wavenumbers in the z direction and L z is the length of the computational domain at the inlet in the z direction. The inlet lies in the y-z plane for this example. The number N z M odes is computed by

N z Modes = E k z M ax 2π/L z , (13)
where E (•) represents the integer part and k z M ax is the maximum admissible value in the z direction due to the mesh size and to the numerical scheme employed. For example, using the finite difference 11-point scheme optimised in the phase space of Bogey and Bailly [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF], at least 4.6 mesh spaces by wavelength are needed to accurately resolve the wave which gives k z M ax = 2π 4.6∆ , with ∆ the mesh size. The next steps for one iteration are described by the pseudo-algorithm on Fig 2 . The number N k z =0 is the number of modes for which k z = 0 and N k z 0 the number of modes for which k z 0. The NiteMax number is the number maximal of iteration autorised to find modes, first, for each non zero k nz mode component, and then, for each k nz = 0 mode component. The global algorithm is repeated until a sufficient number of modes N k z 0 + N k z =0 has been found, with a certain relative error with N E xact Generally, only one iteration for the global algorithm, and 2 or 3 iterations to find at least one correct mode for each non-zero k nz are necessary when we accept that the total modes number is at 1% or less from N E xact . The final step is to compute the other variables α n , σ n , ψ n and ũtn necessary to build the velocity fluctuations and this is done the same way as in the RFM method. The fluctuating velocity is then contructed by summing the modes with Eq. ( 5). While the RFM method computes randomly φ n and θ n to have the wavenumber components, the RFM-P computes first the wavenumber components and then φ n and θ n from it.

Though the RFM-P should generate less spurious noise in case of periodicity condition, it is easy to see that the isotropy is broken. Indeed, we set the value of the wavenumbers components in the periodic direction and, to complete the wavenumber vector, we begin with the largest value for the periodic component in order to make sure it is used. Indeed, the largest value is the most difficult to get since it needs a combination of two small values in the two other directions to be in the range of the desired wavenumbers. The risk is to favour the z direction, thus breaking the isotropy. Figure 3 illustrates this issue by comparing the distribution of modes result between the RFM and the RFM-P methods in the y-z plane, both with 200 modes, L z = 0.25 m and k

x M ax = k y M ax = k z M ax = 466 m -1 .
It is clear that more wavenumbers are present on the high z values with the RFM-P than with the RFM, which means that the z direction appears to be favoured.

FOR each k nz , starting from the highest absolute value IterationNumber = 0

LOOP : WHILE (IterationNumber < NiteMax) a) N k z 0 = 0 b) Generate N E xact random values of k ix ∈ [-k x M ax ; k x M ax ] and k iy ∈ {-k y M ax , ..., k y M ax } c) IF at least one ln (k ni ) = ln k 2 ix + k 2 iy + k 2 nz is at 0.05%
N k z =0 = N E xact -N k z 0 modes, setting k nz = 0, IterationNumber = 0, GO TO LOOP with N k z =0 instead of N k z 0

IV. Synthetic Turbulence Validation in Homogeneous Isotropic Case

Before computing turbulence interaction with airfoils, it is necessary to make sure that the SEM, the RFM and the RFM-P methods meet the desired criteria in terms of realistic energy density spectrum, isotropy, homogeneity and a decrease in spurious noise for the RFM-P method. Moreover, the influence of numerical parameters on those properties must be assessed. For this purpose, a spatially decaying HIT test case is carried out for the three methods.

Experimental data of Comte-Bellot and Corrsin [START_REF] Comte-Bellot | Simple Eulerian time correlation fo full and narrow-band velocity signals in grid-generated, isotropic turbulence[END_REF] are used for comparison. They studied spatially decaying HIT in a wind tunnel and made measurements of turbulent kinetic energy k t , longitudinal integral length scale Λ f and transverse Taylor micro-scale λ g in three vertical planes at U BC t M = [42; 98; 171], with M the size of the mesh grid generating the turbulence in the experiments and U BC the mean flow velocity before the wind tunnel contraction. Only M = 5.08 cm is considered here. For the three methods, the synthetic turbulence in the inlet plane of the computational domain is set with the turbulence parameters measured by Comte-Bellot and Corrsin [START_REF] Comte-Bellot | Simple Eulerian time correlation fo full and narrow-band velocity signals in grid-generated, isotropic turbulence[END_REF] in the plane U BC t M = 42. The values are Λ f = 0.024 m = 0.47M and k t = 0.0739 m 2 .s -2 with a Reynolds number based on λ g equal to

Re λ g = √ k t λ g ν = 71.6
, where ν is the kinematic viscosity. Then, the mean flow convects the turbulence on the corresponding next plane U BC t M = 98 and the spatial spectrum is compared with the one measured by Comte-Bellot and Corrsin. It is important to note that the convection time between the two planes is significant since a spatially decaying turbulence is in reality a convected temporal decaying turbulence. Every simulations are run on 12 non-dimensional time

t * L x = t L x U x L x ,
where t L x is the necessary time to cross the entire computational domain of length L x at the simulation mean flow U x . All the statistics (turbulent kinetic energy and spectra) are calculated on the last 2 non-dimensional times.

A. Incompressible LES with the SEM -Turbulence Properties

Computational Set-up and Parameters

For the Code_Saturne incompressible LES using the SEM, the computational domain is a simple box. Its dimensions are [180M;5M;5M] = [381Λ f ; 10.6Λ f ; 10.6Λ f ] and the numerical mesh size is ∆ = λ g /1.5 = 3.23 × 10 -3 m. A small mesh size is chosen because of the low Reynolds number Re λ g . Indeed, the LES subgrid models used here are based on the assumption that an inertial slope is present in the spectrum but this one only appears clearly at higher Reynolds numbers [START_REF] Bailly | Turbulence[END_REF]. Because the separation between large and small scales is not obvious here, it is important to set a mesh grid size lower than the Taylor scale to activate the subgrid model approximately where the inertial range appears, and distant from the large scale spectral area. Indeed, the size of the Gaussian filter used in the Code_Saturne LES for the subgrid model is twice as big as the mesh size. To see the impact of the subgrid model on the turbulence resolution, the Smagorinsky model and the Germano and Lilly [START_REF] Bailly | Turbulence[END_REF] dynamic model are tested. For this simple test case, if the mesh size is well chosen, the differences between the results given by the models should be minor. However, for the airfoils simulations, the dynamic model will be needed to cope with the high differences between inlet turbulence and boundary layer turbulence. The influence of shape functions on the turbulence resolution will be tested as well as the effect of the time-step ∆ t through the Courant number defined as CF L = U x ∆ t ∆ . Periodic conditions are set on opposite lateral faces. A Dirichlet condition on the velocity and a zero Neumann condition on the pressure are imposed at the inlet. A Dirichlet condition on the pressure updated at each iteration is imposed at the outlet, in order to meet the condition ∂ 2 P ∂n∂τ = 0, with P the static pressure, n the unit normal and τ the unit tangent at each outlet face. The mean flow velocity is 12.7 m.s -1 as in the Comte-Bellot and Corrsin's wind tunnel contraction. Since the velocity before the contraction U BC is equal to 10 m/s in the experiments, the x = 0 plane in the simulation corresponds to the

U BC t M
= 42 plane of Comte-Bellot and Corrsin's experiments while the U BC t M = 98 experimental plane is reproduced at x = 71.12M = 150Λ f downstream of the inlet in the simulations to keep the same convection time between both planes. Three simulations, SEM1, SEM2 and SEM3, are run and their characteristics are gathered in Table 1. They are run with CFL = 1 which gives ∆ t = 2.54×10 -4 s and 3.398×10 4 iterations. Each simulation is characterised by a subgrid model and by a shape function. As a reminder, the shape functions are driven by Eqs. ( 3) and (4). Finally, the Smagorinsky constant C s is set to 0.18, as the theoretical value for HIT [START_REF] Bailly | Turbulence[END_REF] and the value of C s , notation extended to the dynamic model, is authorised to vary between 0, to prevent backscattering, and 0.18, to avoid unrealistic energy exchange. First, the resolved turbulent kinetic energy, averaged on each y-z plane is plotted for the three simulations on Fig. 4 and is compared with the empirical values for the decay slopes from Mohamed and Larue [START_REF] Mohamed | The decay power law in grid-generated turbulence[END_REF], whose values are between -0.95 and -1.43. The first obvious conclusion is that neither the subgrid model nor the shape function seems to have a high impact on the solution which tends to prove the robustness of the method and the correct mesh size choice. The dynamic model seems to have a better behaviour at the outlet than the Smagorinsky one, which must be linked to the boundary conditions. Whatever the simulation, the kinetic energy is almost constant up to 50Λ f which means that no dissipation takes place. Then, from 50Λ f , the energy starts to decay. These results are also found in the litterature and this distance is representative of an adaptation distance [START_REF] Dhamankar | Overview of Turbulent Inflow Boundary Conditions for Large-Eddy Simulations[END_REF]. It represents a problem for most of the synthetic turbulence methods. The problem lies in the initial spectrum choice and in the capacity of synthetic turbulence to accurately reproduce the 3 rd order moment (skewness) of the turbulence, which is responsible for the good quality of energy exchange between large scales and small scales. It means that the adaptation distance corresponds to a distance or to a time necessary for the energy cascade to appear and the dissipation to occur. The adaptation distance can be roughly estimated by the large eddy turn-over distance, as Jarrin et al. did [START_REF] Jarrin | A Synthetic-Eddy-Method for generating inflow conditions for large eddy simulations[END_REF], such that

Influence of the Subgrid Model and of the SEM Shape Function on the Turbulence Properties

L Λ f = U x Λ f 2/3k t (14)
which gives L Λ f = 57Λ f . This value is in agreement with the previous 50Λ f found on Fig. 4. However, the empirical laws for the decay slopes are really retrieved only later for x = 150Λ f which means that, to completely recover the correct statistics, the real adaptation distance is bigger than the one estimated by Eq.( 14). The surestimation and oscillations of the turbulent kinetic energy on Fig. 4 will be analysed in Section IV.A.3. A conclusion is that the turbulence seems to be simply convected with very little dissipation up to 50 Λ f . In other simulations, it has been seen that the results do not depend on the number of structures N if the inlet plane is statistically covered by those structures, i.e., if

N ≥ S i S Λ f
where S i is the inlet plane surface and S Λ f is the surface of a structure associated to Λ f . Figure 5 shows, for the three simulations and for both plane positions, the one-dimensional longitudinal autospectra in the y-direction. These spectra are calculated based on spatial data taken from probes distributed along the y axis at z = L y /2. The data are compared to the Comte-Bellot and Corrsin's ones. It is clear that the initial spectrum at the inlet, on Fig. 5(a), is far from the experimental one in terms of form. It is simply due to the choice of a Gaussian function for the shape function. Again, the subgrid model or the shape function are not important choices except at high wavenumbers where the Smagorinsky model is more dissipative than the dynamic one, which is the same conclusion as the one obtained by Kato et al. [START_REF] Kato | Numerical prediction of sound generated from flows with a low Mach number[END_REF] on their temporally decaying turbulence. The large structures should be correctly represented by the SEM with a Gaussian shape function. It is clear that the middle to high wavenumbers are badly represented,which can explain the low decay of k t on the Fig. 4. On Fig. 5(b), at x = 150Λ f , the energy seems to be better distributed accross the wavenumbers, even if the spectra are still far from the experimental ones. If the energy is still high, the redistribution takes place and the energy cascade begins to appear, which leads, a few Λ f later, to the theoretical turbulent energy decay.

Influence of CFL Number on the Turbulence Properties

The influence of the CF L is now studied through four new simulations, for which the corresponding CFL, time step and iteration number are gathered in Table 2. Simulations are still run on 12 non-dimensional times. To avoid expensive simulations, the mesh size is ∆ = λ g . The dynamic model and a Gaussian shape function are used since it has been seen that they give slightly better results on the kinetic energy. Figure 6 shows the turbulent kinetic energy decay along x/Λ f for these four simulations. It is clear that the CFL has a high impact on the results. Indeed, for CF L = 1, the energy is overestimated and it is only with CF L = 0.5 or lower that a convergence is reached. The adaptation distance is still present since it is linked to the initial spectrum but the overestimation and the oscillations of energy disappear. It is important to note that every curve starts at the same k t value at x = 0 but that the energy for CF L = 1 rises on the first cells. Figure 7 shows the one-dimensional longitudinal autospectra on the two planes under study for the four simulations. The simulation for CFL = 1 shows a higher energy at all wavenumbers at x = 0 which confirms the previous observations on the kinetic energy.

Finally, a more problematic effect of high CF L numbers can be seen on Fig. 8(a). It shows frequency spectra produced thanks to time-domain x velocity fluctuations data on a probe located at x = 25M, y = L y /2 and z = L z /2. For CFL = 1 and CFL = 0.5, high energy spurious peaks are present at respectively f = 300 Hz and f = 400 Hz. These could be an aliasing problem, since CFL = 1 corresponds to a time step equal to half the lowest time scale introduced in the simulation but Choi and Moin [START_REF] Choi | Effects of the computational time step on numerical solutions of turbulent flow[END_REF] show that a time step six times lower than the lowest time scale to be solved is necessary to be independent of the time discretisation, while simulating a plane channel flow. The authors showed that the problem was not linked to the time integration scheme. Later, Kato et al. [START_REF] Kato | Numerical prediction of sound generated from flows with a low Mach number[END_REF] carried out DNS simulations on temporally decaying HIT and they came to the same conclusion with a factor 2 between the smallest time scale and the time step, when seeking a convergence on the spatial spectrum, which is also found in the present work. Regarding temporal spectrum, the condition on the CFL seems to be more restrictive. The fact that a time step lower than half the lowest time scale is necessary could still be explained by aliasing due to the non-linearity of the Navier-Stokes equations. Indeed, no peak is present when a frequency spectrum is computed from the SEM fluctuations without being injected into the Navier-Stokes equations. The convective term of the Navier-Stokes equations causes the apparition of small structures which cannot be resolved by the mesh grid. This was studied by Chow and Moin [START_REF] Chow | A further study of numerical errors in large-eddy simulations[END_REF] which show that aliasing due to the convective term is the cause of high energy at high wavenumbers on a narrow spectral range. Kravchenko and Moin [START_REF] Kravchenko | On the Effect of Numerical Errors in Large Eddy Simulations of Turbulent Flows[END_REF] show that it has a bad impact on turbulence statistics such as fluctuating or mean velocity profiles. Thus, the energy overstimation on

Isotropy and Homogeneity of the SEM

The isotropy property of the SEM is now evaluated. The previous simulation with CFL = 0.25 in Table 2 is used. Figure 8(b) shows the comparison between the spatial spectra E vv y k y and E ww z (k z ) on the planes x = 0 and x = 150Λ f . It is easy to see that from the start to the end, the isotropy is perfectly respected. Finally, it has been verified that the homogeneity property is also perfectly respected thanks to time-domain spectra, at x = 10M and x = 150Λ f .

B. Compressible Simulations with the RFM and RFM-P Methods -Turbulence Properties and Spurious Noise

Computational Set-up and Parameters

For the compressible LES with Code_Safari using the RFM and RFM-P methods, the computational domain is composed of two 8 th -order interpolation overset boxes. The first is the study domain and has dimensions [250Λ f ; 5Λ f ; 5Λ f ]. The second is a sponge zone, extending from 245Λ f to 290Λ f with the same dimensions along y and z. The dimensions of the computational domain are different from those used in Section IV.A.1 because compressible calculations are more expensive. The mesh size is ∆ = λ g 1.1 in the box and it is progressively stretched with an exponential law in the sponge zone. With the explicit filtering, the previous problem of models working only in the inertial zone is not present; the mesh size choice has less importance. The CFL choice is based on numerical considerations, in order to limit diffusion and dispersion errors. With the schemes described in Section III.A, CF L = 0.95 is low enough. Bogey and Bailly outflow boundary conditions [START_REF] Bogey | Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases[END_REF] are used for the up, down and outflow faces, periodic conditions are set for the front and back faces (z direction). Tam and Dong inlet conditions [START_REF] Tam | Advances in numerical boundary conditions for computational aeroacoustics[END_REF] are used for the inlet face to prevent acoustic reflection and to introduce the synthetic turbulence, built with the RFM or the RFM-P methods. The centre of the Tam and Dong radiation conditions is placed at the middle of the domain. A mean flow Mach number M x = 0.1 has been taken to give sufficient time to the turbulence to decay but also to avoid numerical errors since Code_Safari is ill-conditionned at very low Mach numbers. We note τ η = ν the Kolmogorov time scale, with is the turbulent dissipation rate such that = (2/3k t ) 3/2 Λ f . We also note t C BC the convection time between the two first measurements planes in Comte-Bellot and Corrsin's experiments, and t 250Λ f the time to travel the entire study box. We have

τ η = 5.7 × 10 -3 s t 250Λ f = 0.18 s < t C BC = 0.2240 s. ( 15
)
Even if the simulation does not allow the turbulence to be in the same state as in the experiments at the end of the computational domain, the convection time should be high enough to let it decay and to see its properties evolve. With this CF L number and this Mach number, the time step is ∆ t = 1.13 × 10 -5 s which is ten times lower than the smallest Table 3 RFM method parameters for the three Code_Safari simulations. N is the modes number, k 1 and k N are the minimum and maximum wavenumber amplitudes for the spectrum discretisation, k

x,y,z M ax is the maximum wavenumber component value in each direction x,y and z (only for RFM-P) and N Modes z is the number of positive non zero wavenumber components along z (only for RFM-P) turbulent time scale to capture, so that no aliasing problem will occur. The aim here will be to see qualitative changes in the turbulence along its convection and not to do some quantitative comparisons with the experimental results, since the convection times are not the same. Three simulations have been run with different inlet turbulence conditions. For 12 non-dimensional simulation times, 1.85×10 5 iterations are needed. The filtering strength coefficient [START_REF] Daude | A high-order finite-difference algorithm for direct computation of aerodynamic sound[END_REF] is equal to 0.2.

N k 1 (m -1 ) k N (m -1 ) k x,
For the first simulation, the turbulence is generated using the RFM method, the second one uses the RFM-P method and the last one uses the RFM-P method with two-dimensional structures, i.e. k z = 0 m -1 only. This last simulation is a limit case to see the influence of the absence of the wavenumber along z on the isotropy property and is close to a typical airfoil simulation. Indeed, the airfoil spanwise is limited to a short distance (10% to 50% of the chord) which means that a lot of the turbulent structures will be seen as two-dimensionnal for the airfoil. Table 3 gathers information on modes number, minimum and maximum wavenumbers authorised for the turbulent spectra discretisation, the maximum value for the wavenumber components (only for RFM-P) and the number of positive non zero wavenumber along z (only for RFM-P simulations). We use 200 modes which is twice the number used in Gloerfelt and Le Garrec [START_REF] Gloerfelt | Generation of inflow turbulence for aeroacoustic applications[END_REF] but different numbers should be tested in further studies.

Turbulent Kinetic Energy Decay

Figure 9(a) shows the resolved turbulent kinetic energy averaged in each y-z plane comparison with theoretical slopes. Figure9(b) has same y axis than Fig. 6 for direct comparison with incompressible simulations. A simple conclusion is that the energy is very slowly decaying compared to theoretical slopes and to incompressible results. Indeed, only 20% of the energy is lost at 200 Λ f . Then, the maximum difference of energy between simulations is 8% and this is way less at short distance (x<100Λ f ). It means that the methods gives very similar results for the resolved kinetic energy decay. The behavior of the turbulence seems to change at 150Λ f . Indeed, the resolved energy seems to increase. This increase is more important with the RFM-P simulation with k z = 0 simulation than with the RFM-P one, the RFM simulation being an intermediate case. This is attributed to numerical errors due to the interaction between turbulent structures and outflow boundary conditions, up and down. Even if the convection time is high enough to see some changes, they are negligible, above all at short distances used for airfoils simulations (<50Λ f ).

Spatial Spectra

Figure 10 compares the one-dimensional spectra along y for the three simulations at x = 0 and x = 120Λ f . The results on x = 0 are compared with experimental data. It is not relevant to compare results on x = 120Λ f since the convection time is not the same as in the experiment. From this figure, it can be concluded that no major differences exist between the three methods on the spectra, except for the RFM-P -k z = 0 simulation where the energy at high wave numbers is lower than the two other simulations in the plane x = 120Λ f . In the three simulations, the -5/3 slope is retrieved.

Spurious Noise Produced by the Three Simulations

Figure 11 shows the sound pressure level calculated on the two last non-dimensional times at some probes located at L y /2, L z /2 and at different positions along x for the three simulations. It is clear that the RFM-P method creates 10 dB less than the RFM method at short distances (x < 50Λ f ). Likewise, the RFM-P -k z = 0 simulation creates 20 dB less than the RFM method. This figure shows that adapting the synthetic turbulence to the periodic conditions is a way to reduce the spurious noise produced by its introduction in the computational domain. Also, a way to reduce the noise is to place the airfoil sufficiently far from the inlet (at least 40Λ f ).

Turbulence Isotropy and Homogeneity

Figure 12 shows the one-dimensional longitudinal autospectra along y and z for the isotropy study on x = 0 and x = 120Λ f and for the three simulations. While the RFM is perfectly isotropic, the RFM-P shows a slight anisotropy at low wavenumbers, which tends to decrease with distance. The simulation RFM-P -k z = 0 is fully anisotropic with no energy along z which is fully coherent with the fluctuations contruction. These observations are coherent with the discussion in Section III.B.2 and with the Fig. 3. It has also been verified that the homogeneity property is perfectly verified with time-domain spectra in x = 0 and x = 120Λ f for every simulation.

C. Comparison between SEM and RFM-P

Figure 13 compares one-dimensional spectrum for incompressible using SEM and Code_Saturne, 200 structures and shape function with spectrum for the compressible simulations using the RFM-P method, Code_Safari and 200 modes. It is important to note that the spectra are computed at x = 0 to see only the inlet perturbation in each code and not the influence of the code on it. These spectra are compared to experimental data corresponding to the ideal spectrum. Furthermore, the theoretical spectrum for both methods are plotted on the figure. For the SEM, the theoretical spectrum is computed directly using the velocity expression of Eq (2) while, for the RFM-P, it is based on the von Kármán-Saffman spectrum E of Eq. ( 26) combined with the equation from Bailly and Comte-Bellot [START_REF] Bailly | Turbulence[END_REF] :

E vv k y = 1 2 ∫ ∞ k y E(k) k 1 - k 2 y k 2 dk, (16)
with k = | k | the wavenumber amplitude, k y the wavenumber along y and E vv the longitudinal autospectrum along y. The theoretical RFM-P spectrum is done with E discretised between 25 m -1 and 500 m -1 with 200 modes, same parameters as before, to give the spectrum injected in LES, without the non-resolved scales. That is why the theoretical spectrum on Fig. 13 draws away from the experimental one. The figure shows that the theoretical spectrum are respected by the simulations but the Gaussian spectrum for the SEM is far from the experimental one for high wavenumbers, while the spectrum used for the RFM is very close to it. The discrepancies between the compressible simulation using the RFM method, and the experiment, on the plane x = 0, are only due to the LES cut. It means that the RFM should have a better ability to reproduce quickly the energy cascade mechanism and so to reduce the adaptation distance, compared to the SEM. The Gaussian shape for the SEM explains why the adaptation distance is so long. The fact that the turbulent energy is not decaying in the compressible simulations with the RFM method is probably due to the high Mach number and to boundary conditions numerical errors rather than the turbulent spectrum prescription.

V. Preliminary Results About Turbulence Injection on a NACA 0012 Airfoil at Low Reynolds Number

In this part, preliminary results about flow validations on a NACA 0012 airfoil with inflow turbulence are shown. First, compressible simulations without turbulence are presented in Section V.A to emphasise the role of turbulence fluctuations on the LSB appearance. Then, preliminary incompressible simulations of LSB with inflow perturbations are carried out in Section V.B and the lift and drag coefficients are compared to experimental results.

A. Compressible Simulations without Inflow Turbulence

A first compressible simulation is performed with Code_Safari and without inflow turbulence, for a NACA 0012 with an angle of attack α = 15 • . With a chord-based Reynolds number Re c = 10 4 and an upstream Mach number M x = 0.2, the LSB burst which should normally appear on the suction side, according to Alam et al. [START_REF] Alam | The ultra-low Reynolds number airfoil wake[END_REF], is not present in the simulation. Figure 14 illustrates that while visualising the vorticity magnitude at a non dimensional time t * c = 336, based on the chord and on the upstream velocity, only the outbreak of 2D rolls is present when LSB (Laminar Separation Bubble) burst is highly 3D. This is attributed to a lack of perturbation in the flow : the instabilities development is too long to be seen in this time window. The idea is to introduce perturbations to speed up the instabilities development. To be close to experimental conditions and to be able to compute leading edge noise, realistic synthetic turbulence is desired. New compressible simulations with RFM built inflow perturbations are currently carried out.

B. Incompressible Simulations with Inflow Turbulence

Several preliminary incompressible simulations with Code_Saturne and using the SEM to generate inflow turbulence are now presented. The simulation parameters are based on the experiments from Zhou et al. [START_REF] Zhou | Fluid forces on a very low Reynolds number airfoil and their prediction[END_REF] and Wang et al. [START_REF] Wang | Turbulent intensity and Reynolds number effects on an airfoil at low Reynolds numbers[END_REF] that were carried out in a closed-loop water tunnel to study the behavior of a LSB on a NACA 0012 with different inflow turbulence intensities. The chosen cases here are for Re c = 5.3 × 10 3 and M x = 3.3 × 10 -5 . Due to the very low Mach number, incompressible simulations are more suitable for these configurations.

Table 4 gathers information about the experimental cases. It indicates whether a turbulence grid is present in the experiments or not, the turbulent intensity, the integral length scale and the turbulence parameters used in simulations corresponding to Cases 1 and 3. Simulations were run at α = [5 • , 10 • , 15 • ]. The synthetic turbulence was built with the Gaussian shape function and 200 structures. Unfortunately, for Case 1, there is no integral length scale available in the experiments since no turbulence generator grid was present, which can also lead to anisotropy. Thus, it was chosen such that Λ f = c = 0.1 m in order to have a turbulent scales separation without creating too large structures. As seen in Table 4, the Taylor based Reynolds numbers are Re λ g = 18 or 28. These values are too low to be able to carry out a LES, it means that a quasi-DNS must be done to carry the inflow perturbations up to the airfoil. The airfoil was surrounded by a O-type mesh with desired maximal non dimensional mesh size ξ + = 10 in the streamwise direction, η + = 1 in the normal direction and ζ + = 10 in the spanwise direction. These values are low enough to efficiently capture the LSB. After the simulations were finished, the value of η + is checked and is between 0.2 and 5 on the suction side. The values of 5 are reached only on 10 cells near the trailing edge, probably due to distorted cells. Elsewhere, η + is around 1. To carry the inflow turbulence up to the airfoil with a quasi-DNS, the inlet mesh size was set to ∆ I nlet = η K /4 = 3.1 × 10 -3 m where η K is the Kolmogorov scale of the injected turbulence for Case 1. For Case 3, this gives ∆ I nlet = 1.72η K , which is too high to obtain a quasi-DNS. The inlet is located at 6c from the leading edge and the outlet at 5c from it. The up and down sides are at 5c from the leading edge and a symmetry condition is applied on it. The spanwise length is 0.5c which is enough to have an accurate description of 3D structures as seen by Thomareis and Papadakis [START_REF] Thomareis | Effect of trailing edge shape on the separated flow characteristics around an airfoil at low Reynolds number: A numerical study[END_REF]. Periodic boundary conditions are applied in the spanwise direction. The meshes are made up of 11×10 6 elements. This set-up gives a distance between the inlet and the airfoil leading edge of 6Λ f for Case 1 and 24Λ f for Case 3, which prevents any energy decay as seen in Section IV.A.4. This means that the turbulent kinetic energy set at the inlet will be the same as the one reaching the airfoil. The subgrid dynamic model is used and C s is restrained from 0 to 0.065, which is the theoretical value for a turbulent channel flow. Simulations are run for t * c = 40 with a time step ∆ t = 2.5 × 10 -3 s which gives 3.1×10 4 iterations. The time step is chosen as ∆ t = min τ W all , τ inlet /4 where τ W all = η + ν/(0.1U x ) 2 as the estimated minimum time scale present near the airfoil in the boundary layer, and τ I nlet = ∆ I nlet /U 0 as the minimum time scale introduced with the SEM. The factor 1/4 ensures the quality of the results based on what has been said in the section IV.A.3 about aliasing. This gives CF L = 0.0765 at the inlet. The lift and drag coefficients, C L and C D , are averaged on the last t * c = 20. First, it is important to note that for Cases 1 and 2, there is experimentally no LSB as the C L curves show no stall and that no particular increase is visible for α = 10 • , unlike Case 3. For α = 15 • and Case 3, the stall shows that the LSB has burst. Then, simulations results for C L and C D are shown in Fig. 15 along with experimental results from Zhou et al. [START_REF] Zhou | Fluid forces on a very low Reynolds number airfoil and their prediction[END_REF] and Wang et al. [START_REF] Wang | Turbulent intensity and Reynolds number effects on an airfoil at low Reynolds numbers[END_REF]. First, Fig. 15(a) shows agreement of the lift coefficient for α = 5 • between simulations and the three experiments. Then, for α = 10 • , there is a discrepancy between the simulation and the results from Zhou et al. [START_REF] Zhou | Fluid forces on a very low Reynolds number airfoil and their prediction[END_REF] for Case 1 but there is a good agreement with the results from Wang et al. [START_REF] Wang | Turbulent intensity and Reynolds number effects on an airfoil at low Reynolds numbers[END_REF] for Case 2. There is also a discrepancy between the simulation and the experiment for I t = 2.6% (Case 3). For this angle of attack, there is a strong influence of the turbulent intensity on the lift coefficient, even when there is no LSB on the airfoil for Cases 1 and 2. It shows that the airfoil in simulations does not see the same turbulence as in the experiments or does not not react to it the same way, since the simulation result for I t = 0.45% is closer to the experimenal results with I t = 0.6%. Furthermore, the Case 3 is not well reproduced by the simulation. In this case, an LSB is experimentally present and it seems that the simulation is not able to correctly capture it. Finally for α = 15 • , the difference between simulation and the experiments in Case 1 is even more pronounced but the error for Case 3 is very small. Again, the simulation in this case is closer to the Case 2 experimental result. This tends to show that, whatever the angle of attack, the simulations in Case 1 are closer to the experiments of Wang et al. [START_REF] Wang | Turbulent intensity and Reynolds number effects on an airfoil at low Reynolds numbers[END_REF] in Case 2 than the one from Zhou et al. [START_REF] Zhou | Fluid forces on a very low Reynolds number airfoil and their prediction[END_REF] in Case 1, especially at high angles of attack where the flow separates.

The discrepancies between simulations and experiments for Case 1 can be attributed to the fact that the experimental turbulence is not generated by a grid and is the residual turbulence of the water tunnel. Its properties are not controlled and the isotropy might not be respected. Another reason for the discrepancies is that the integral length scale in Zhou et 4 for references.

al. [START_REF] Zhou | Fluid forces on a very low Reynolds number airfoil and their prediction[END_REF] was not specified so that a choice has been made for it in the simulation, which is maybe different from the experimental one. Regarding Case 3, It can be concluded that the simulations are not able to correctly reproduce the way the airfoil reacts to the inflow turbulence or the turbulence itself, since there is no aerodynamic stall in the simulations while there is one in experiments. The LSB might not be accurately captured at α = 10 • . A reason could be that the grid carrying the turbulence from the inlet up to the airfoil is too coarse to compute the Kolmogorov scale, which means that the airfoil encounters only the largest scales of the turbulence. This should be improved in future simulations. Another reason could be the averaging time window which is only on the last t * c = 20 = 40 s in simulations while it is 5 minutes in experiments. When the LSB is a highly instationnary phenomenon and could be subject to intermittency as seen in Eljack [START_REF] Eljack | High-fidelity numerical simulation of the flow field around a NACA-0012 aerofoil from the laminar separation bubble to a full stall[END_REF], the averaging time is critical and longer simulations should be carried out in the future. Figure 15(b) shows that the simulated drag coefficient is in good agreement with all the experiments for each angle of attack. This is not surprising as this coefficient is not dependent on the turbulence intensity as shown in Wang et al. [START_REF] Wang | Turbulent intensity and Reynolds number effects on an airfoil at low Reynolds numbers[END_REF]. The drag coefficient is quasi-stationnary compared to the lift coefficient.

Finally, some flow visualisations for Cases 1 and 3 at α = 10 • are shown in Fig. 16 where the streamwise velocity is plotted for t * c = 30.7. It clearly shows the separation of the laminar boundary layer on both visualisations. While the Case 1 velocity seems to show 2D well-organised structures, Case 3 velocity is more disorganised with high negative values closer to the leading edge. It could be the LSB reattaching earlier on the aifoil as the turbulent intensity increases. This would explain the slight increase of the lift coefficient between Case 1 and Case 3. The difference between experiment and simulation for Case 3 is still to be analysed but could be a delayed reattachment of the boundary layer in the simulation.

VI. Conclusion

The present work aimed to generate inflow synthetic turbulence in Large Eddy Simulations (LES) with realistic properties such as isotropy, homogeneity or energy density spectrum, in order to accurately compute Laminar Separation Bubble (LSB) on airfoils. For future noise computations, another objective was to reduce spurious noise level due to spanwise periodicity conditions when introducing turbulence in compressible simulations. Synthetic turbulence has first been studied through spatial decaying homogeneous isotropic turbulence cases with LES. Both the Synthetic Eddy Method (SEM), included in incompressible LES with Code_Saturne, and the Random Fourier Modes (RFM) method, implemented in compressible LES with Code_Safari, show a perfect isotropy and homogeneity at every spatial position. There is a small decay of the turbulent kinetic energy on a long distance for both methods. For the SEM in incompressible simulations, it has been attributed to the prescribed Gaussian spectrum, far from correctly representing the small turbulent scales and the energy exchange from the large to small structures. This distance is responsible for a lack of dissipation up to 50Λ f . However, 150Λ f are necessary before the energy transfer mechanisms are fully established. A way to improve the method by reducing the adaptation distance will be to sum structures, each with a different size. For the RFM method in compressible simulations, the long adaptation distance has been attributed to the relatively high Mach number used, which makes it more difficult to see a decay during the considered time window. By nature, the RFM method adaptation distance should be shorter than the SEM one, since a real spectrum is prescribed, making the correct 3rd order moments easier to retrieve but no conclusion on it can be drawn based on the present simulations. The only way to rigorously compare both methods would be to do it the same code. Additionally for the SEM, it has been shown that the subgrid model and the shape function have a minor effect on the turbulence decay, but the dynamic model is best for its adaptation to the flow. Furthermore, the CFL number has a large influence on the results. When it is too high, an overestimation of the energy is present and a high-energy spurious peak appears in frequency spectra. This is here attributed to possible aliasing effects due to non-linearities of the Navier-Stokes equation. A factor of 4 between the lowest time scale we want to solve and the time step simulation is necessary to get rid of the aliasing effects.

The RFM-P method, stemmed from the RFM method, has been developped to make the inflow turbulence compatible with aerodynamic periodic boundary conditions, thus making less spurious noise when introduced in the calculation domain. Both methods are quite similar in terms of turbulent properties. The homogeneity is respected and spectra are retrieved. The RFM-P method shows exactly the same properties as the RFM one, except that the isotropy is broken, which is a natural consequence of the method construction. It is thought that this loss in isotropy is not necessarily a drawback since numerical calculations will use only a small airfoil span to limit the computational cost. The slice will then see mostly 2D structures in the spanwise direction making the isotropy not relevant. The spurious noise is reduced by at least 10 dB by using the RFM-P instead of the RFM.

Preliminary calculations with inflow turbulence on a NACA 0012 airfoil at a Reynolds number Re c = 5.3 × 10 3

have been carried out in Code_Saturne with the SEM only. While the drag coefficient is in very well agreement with experimental results, the lift coefficient, highly dependent on the turbulence characteristics, shows some discrepancies with the experiments, especially at high angles of attacks. This is attributed to a lack of information about the integral length scale and isotropy in the experiments. Furthermore, the inflow turbulence could not be accurately carried from the inlet up to the airfoil in simulations because of a too coarse mesh. Finally, the LSB and its intermittency could be badly represented as too short simulations were carried out. New improved simulations will be run to check these hypothesis.

k t = N n=1 ũ2 n . (24)
Finally, the combination of Eqs. ( 22) and [START_REF] Bogey | Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases[END_REF] gives

ũn = E(k n )∆k n . (25)
The mode amplitude ũn is determined by the spectrum. To have both the inertial and dissipation parts, the von Kármán -Saffman [START_REF] Bailly | Turbulence[END_REF] has been chosen, such as

E(k) = α s 2k t 3 (k/k e ) 4
k e 1 + (k/k e ) 2 17/6 exp -2

k k η 2 , ( 26
)
where k η is the Kolmogorov wavenumber and α s and k e are constants. Those are computed in Bailly and Juvé [START_REF] Bailly | A stochastic Approach to compute subsonic noise using linearized Euler's equations[END_REF]

where Λ f is the longitudinal integral length scale to set. The α s value represents a normalisation condition to preserve the kinetic energy of Eq. ( 22).

•

  = time or spatial average • = vector notation E (•) = integer part C s = Smagorinsky constant, notation extended to the dynamic model CF L = Courant number C L , C D = lift and drag coefficients E = turbulent kinetic energy spectrum (m 3 .s -2 ) E vv y , E ww z = one-dimensional longitudinal autospectrum in the y, respectively z, directions (m 3 .s -2 ) L x , L y , L z = length of the computational domain at the inlet in the x,y,z directions (m) L Λ f = adaptation distance for synthetic turbulence (m) M = Comte-Bellot and Corrsin's experimental mesh size of the turbulence grid (m) M x = upstream Mach number N = number of synthetic eddies or Fourier modes N E xact = desired modes number N ite = simulation iteration number N z M odes

  L x = non-dimensional time to cross a distance L x u N I = normalised and isotropic velocity fluctuation field u = velocity field (m.s -1 ) u = final velocity fluctuation field imposed at the inlet (m.s -1 ) ũn = amplitude of the n th mode (m.s -1 ) x = spatial position (m) x min = x-coordinate of the computational domain inlet plane (m) α, β = parameters for the shape functions (m -1/2 , m) α n = angular position of σ n in the (e θ n , e φ n ) plane α s = constant in the von Kármán -Saffman spectrum ∆ = simulation mesh size (m) ∆ t = simulation time step (s) = turbulent dissipation rate (m 2 .s -3 ) i j = random value (1 or -1) Λ f = longitudinal integral length scale of the turbulence (m) λ g = transverse Taylor micro-scale of the turbulence (m) ψ n = phase of the n th mode σ n = direction of the n th mode τ η = Kolmogorov time scale (s) ν

Fig. 1 a

 1 Fig. 1 a) Isocontours of v of a turbulent velocity field generated by the SEM, with S the inlet plane, B the virtual box, U c the convection velocity and v the vertical fluctuating velocity. Illustration taken from the Code_Saturne user guide [27]. b) Geometry for the n th RFM mode, σ n is in the (e θ n , e φ n ) plane, orthogonal to k n by incompressibility. Figure inspired from Bechara et al. [28].

  or less from ln k E xact p , p ∈ {1, ..., N E xact } -Increase N k z 0 from the corresponding number of modes -Keep all the k ni -Delete k E xact m from the search base such that k E xact m ∼ k ni following the previous criterion ELSE -IterationNumber = IterationNumber + 1 ENDIF ENDWHILE IF (IterationNumber = NiteMax) STOP ENDFOR For the remaining

Fig. 2 Fig. 3

 23 Fig. 2 Pseudo-algorithm for one iteration of the RFM-P method

Fig. 4

 4 Fig. 4 Turbulent kinetic energy decay for SEM1, SEM2 and SEM3, averaged on each yz plane and compared to theoretical slopes.

Fig. 5

 5 Fig. 5 One-dimensional longitudinal autospectra in the y direction a) in the inlet plane x = 0 and b) in the plane x = 71.12M = 150Λ f , compared to Comte-Bellot and Corrsin's (CBC) experiments, for the three simulations SEM1, SEM2 and SEM3

Fig. 6 Fig. 7

 67 Fig. 6 Spatial turbulent kinetic energy decay k t for different CFL numbers

  Fig 6 and the spurious peaks on Fig 7 could be explained by aliasing, disappearing only when the time step is low enough compared to the lowest time-scale to be solved. In our case, a factor 4 between the both is necessary.

Fig. 8 a

 8 Fig. 8 a) Frequency spectra computed from time-domain x velocity fluctuations data on a probe located in x = 25Λ f , y = L y /2 and z = L z /2 for different CFL numbers and b) longitudinal autospectra in the y and z directions for CF L = 0.25 at x = 0 and x = 150Λ f

Fig. 9

 9 Fig. 9 Spatial decay of the resolved turbulent kinetic energy averaged in each y-z plane for the RFM, RFM-P and RFM-P -k z = 0 simulations, for two different y axis scales

Fig. 10

 10 Fig. 10 One-dimensional longitudinal autospectra along y for the three simulations RFM, RFM-P and RFM-P -k z = 0 at a) x = 0 and b) x = 120Λ f

Fig. 11 Fig. 12

 1112 Fig. 11 Sound pressure level calculated on the two last non-dimensional times at probes located at L y /2, L z /2 and at different positions along x for the three simulations RFM, RFM-P and RFM-P -k z = 0

Fig. 13

 13 Fig. 13 One-dimensional longitudinal autospectra in the y directions for the SEM and RFM-P methods in the plane x = 0 of each incompressible and compressible simulation respectively. The spectra are compared to the theoretical ones. CSF stands for Code_Safari, CST for Code_Saturne and CBC for the Comte-Bellot and Corrsin's experiments.

Fig. 14

 14 Fig. 14 Two-dimensional slices of the vorticity amplitude (s -1 ) for a compressible simulation at α = 15 • and t * c = 336, for Re c = 10 4 and M x = 0.2.

Fig. 15

 15 Fig. 15 Averaged a) lift and b) drag coefficients for the incompressible LES simulations with SEM, and for the experiments of Wang et al. [8] and Zhou et al. [38]. See Table4for references.

Fig. 16

 16 Fig. 16 Visualisations of the instantaneous streamwise velocity at t * c = 30.7 for α = 10 • and a) Case 1 and b) Case 3

Table 1 Code_Saturne simulations characteristics for the subgrid model and the SEM shape function effects on turbulence study, C s is the Smagorinsky constant

 1

		Subgrid model	Form function
	SEM1	Dynamic (C s ∈ [0; 0.18])	Gaussian
	SEM2 Smagorinsky (C s = 0.18)	Gaussian
	SEM3	Dynamic (C s ∈ [0; 0.18])	Hat

Table 2 Code_Saturne simulations characteristics for the CFL effect on turbulence study, N ite is the number of iterations

 2

	CFL	∆ t (s)	N ite
	1	3.8819×10 -4	22260
	0.5	1.9409×10 -4	44520
	0.25 9.7047×10 -5	89030
	0.1	3.8819×10 -5 222600

Table 4 Information about some experiments of Zhou et al. [38] and Wang et al. [8], and corresponding simulations, where I t is the turbulent intensity, Λ f is the longitudinal integral length scale and Re Λ g is the Taylor transverse micro-scale based Reynolds number

 4

	Case	1	2	3
	Experiments	Zhou et al. [38] Wang et al. [8] Wang et al. [8]
	Turbulence Grid	No	No	Yes
	I t (%)	0.4-0.5	0.6	2.6
	Λ f (m)	/	/	0.025
	I t in simulations (%)	0.45	/	2.6
	Λ f (m) in simulations	0.1	/	0.025
	Re λ g	18	/	28

Acknowledgments

The author would like to acknowledge the French supercomputing centre CINES (Centre Informatique National de l'Enseignement Supérieur) for the A0052A10311 allocation in the GENCI (Grand Equipement National de Calcul Intensif) framework. The authors would also like to thank the DGA (Direction Générale de l'Armement) for their financial support. Finally, the author would like to acknowledge Sofiane Benhamadouche for his advice on the SEM in incompressible simulations, Erwan Le Coupanec and Jean-François Wald for their help with the Code_Saturne code, and Chai Koren for his assistance with the NACA mesh for incompressible simulations.

VII. Appendix

A. SEM Theoretical and Practical Considerations

This section describes the virtual box for the SEM implementation in Code_Saturne. It details also the way to compute the shape functions parameters α, β and l f used in Eqs. ( 3) and (4). In practice, a virtual box is constructed around the computational domain inlet which has dimensions [x minl f ; x min + l f ] along x, [y minl f ; y max + l f ] along y and [z minl f ; z max + l f ] along z where x min , y min , y max , z min and z max are the size of the inlet plane located in x min . The box half depth around x min is defined by the length l f which depends on the shape function used and on the integral length scale, as shown at the end of this Appendix. The fluctuating velocity is constructed on x min by summing the influence of every structures in the virtual box, before being added to the inlet velocity of the computational domain. The α and β shape function parameters, used in the Eqs. ( 3) and ( 4), are determined using the known averaged kinetic turbulent energy k t and longitudinal integral length scale Λ f . First, the constraint u i = 2 3 k t , ∀i, must be respected. When applied to Eq. ( 2), it becomes

Then, Λ f is linked to the velocity longitudinal autocorrelation normalised function along x, R N or m u x u x , by

where e x is the unit vector along the x direction. Applied to Eq. ( 2), R N or m u x u x can be written as

Together, Eq. ( 17), Eq. ( 18) and Eq. ( 19) are used to determine the shape functions parameters α and β of the Eqs. ( 3) and ( 4). For the hat function, these parameters follow the equations

with l f = 4 3 Λ f to lead the shape function value down to f Λ f = 0 on each structure border. For the Gaussian function, α and β are determined by

with l f = 4Λ f √ π = 4β, chosen to lead the shape function value down to f Λ f = 2 × 10 -7 × max f Λ f on each structure border.

B. RFM Theoretical Considerations

This section details the way to compute the mode amplitudes ũn of Eq. ( 5) for the RFM and RFM-P methods. Equation ( 11) can be found by making a link between average turbulent kinetic energy k t , Eq. ( 5) and turbulent energy spectrum E. A first expression for the kinetic energy is

and a second one is

which becomes, using Eq. ( 5),