Characterization of pressure fluctuations within a controlled-diffusion blade boundary layer using the equilibrium wall-modelled LES

In this study, the generation of airfoil trailing edge broadband noise that arises from the interaction of turbulent boundary layer with the airfoil trailing edge is investigated. The primary objectives of this work are: (i) to apply a wall-modelled large-eddy simulation (WMLES) approach to predict the flow of air passing a controlled-diffusion blade, and (ii) to study the blade broadband noise that is generated from the interaction of a turbulent boundary layer with a lifting surface trailing edge. This study is carried out for two values of the Mach number, Ma∞=0.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm Ma}}_{\infty } = 0.3$$\end{document} and 0.5, two values of the chord Reynolds number, Re=8.30×105\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm Re}}=8.30 \times 10^5$$\end{document} and 2.29×106\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.29 \times 10^6$$\end{document}, and two angles of attack, AoA =4∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=4^\circ$$\end{document} and 5∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5^\circ$$\end{document}. To examine the influence of the grid resolution on aerodynamic and aeroacoustic quantities, we compare our results with experimental data available in the literature. We also compare our results with two in-house numerical solutions generated from two wall-resolved LES (WRLES) calculations, one of which has a DNS-like resolution. We show that WMLES accurately predicts the mean pressure coefficient distribution, velocity statistics (including the mean velocity), and the traces of Reynolds tensor components. Furthermore, we observe that the instantaneous flow structures computed by the WMLES resemble those found in the reference WMLES database, except near the leading edge region. Some of the differences observed in these structures are associated with tripping and the transition to a turbulence mechanism near the leading edge, which are significantly affected by the grid resolution. The aeroacoustic noise calculations indicate that the power spectral density profiles obtained using the WMLES compare well with the experimental data.


Methods
Governing equations. In this paper, the following notation is adopted: x 1 or x is the streamwise coordinate; x 2 or y is the wall-normal coordinate; and x 3 or z is the spanwise coordinate. The governing equations used in this study are the full Favre-filtered compressible Navier-Stokes which, for a calorically perfect gas, read where ρ is the density, U i is the velocity component in the x i -direction, P is the pressure, E = P /[ρ(γ − 1)] + U i U i /2 is the total energy, R is the gas constant, T is the temperature, and γ = c P /c V is the ratio of specific heats, which is kept constant at 1.4. The overline (resp. tilde) denotes the filtered (resp. Favre filtered) value. The variables are either spatially filtered or ensemble-averaged quantities depending on the use of LES or RANS equations. The latter model is used for the wall-model part. It is assumed that both the filtered stress tensor τ ij and the filtered heat flux Q j can be expressed in a way similar to their instantaneous counterparts, but applied to filtered quantities where µ is the molecular viscosity that has a power-law dependence on temperature The parameter in (2) is the molecular thermal conductivity, Pr = µc P / is the molecular Prandtl number kept constant at 0.7, and S ij d = S ij − δ ij S ij /3 is the deviatoric part of the rate-of-strain tensor, which is defined as S ij = (∂ U i /∂x j + ∂ U j /∂x i )/2 , where δ ij is the Kronecker delta symbol. The other parameters µ t and t in (2) are the turbulent eddy viscosity and conductivity, respectively. The LES solution is formally defined everywhere in the computational domain, and the wall-model equations are solved on a separate embedded grid near the wall. The only difference between the LES and the wall-model equations is in the computation of µ t and t . In the LES framework, µ t = µ SGS and t = SGS are the SGS eddy viscosity and conductivity, respectively. Note that the isotropic part of the modeled turbulent stress tensor is neglected in (2). This is often the case in the LES approach for low Mach number flows. The same assumption is made in the zero-equation RANS modeling context, but is less justified. equilibrium wall-model. The equilibrium wall-model is derived from the compressible Reynolds-averaged Navier-Stokes equations with the boundary layer scaling approximations and neglecting unsteady and convective terms such as pressure gradient. The equilibrium model essentially relies on a grid resolution based on the outer layer requirements using the largest scales of the boundary layer thickness. The inner layer lies within the first cells normal to the wall and its behavior is modeled through the momentum wall normal flux. A typical resolution is about 20 cells per turbulent boundary layer thickness, δ , in each spatial direction 10 . However, in this study, a fundamental question arises when considering the pressure fluctuations that are responsible for noise generation. In fact, to the best of our knowledge, there are no studies in the literature that guarantee that using a WMLES approach with a resolution of 20 cell per δ results in an accurate prediction of the turbulent structures responsible for the pressure fluctuations and hence, that can achieve an overall discretization that can predict the noise sources. However, because the location of the maximum pressure fluctuations is within the outer layer, this resolution appears to generate reasonable and satisfying results. The verification of this assumption is one of the primary goals of this paper.
In the following, we briefly review the equilibrium model used herein combined with the WMLES approach. Hereafter, we use the subscript " wm " to denote a quantity at the wall. Thus, given the instantaneous magnitude of the wall-parallel velocity U || and the instantaneous temperature T in the LES at height x 2 = l wm perpendicular to the wall, we estimate the instantaneous wall shear stress vector τ w,i by solving the following system using two ordinary differential equations (ODEs) 10 where the wall-model eddy viscosity is taken as  (5) and (6) are κ = 0.41 , A + = 17 and Pr t,wm = 0.9 . For compressible flows, the scaled distance from the wall in the Van Driest damping factor ( x + 2 in (6)) is computed using a semi-local scaling given by x + 2 = x 2 √ ρ|τ w | 17 . These two coupled ODEs are solved using the tridiagonal matrix or Thomas algorithm (TDMA) applied in a segregated manner, i.e., by alternating TDMA sweeps of the momentum and energy equations with updated eddy-viscosity µ t,wm in between. The pressure is assumed to be constant in the x 2 direction and is obtained from the LES solution at the exchange location, while the density is computed using the temperature via the equation of state.
One major difficulty in the parallel implementation of the model for unstructured grids is when the boundary condition for the wall-stress model (exchange location) and the wall location are located on different processors in terms of the domain decomposition associated with parallel computations. This issue is fixed by linking each wall face with its exchange location during the pre-processing step using parallel communicators, which are then used during the computation. This last point is important in terms of numerical errors that can arise within the first cell of the LES mesh. Indeed, these errors may significantly decrease the accuracy of the solution even though the model is correct. The recommended choice of the exchange location above the wall boundary is typically ∼ 0.1δ 14 . The grid spacing in the wall normal direction must be smaller than the resolution needed for fully turbulent simulations, which is a limitation of the WMLES approach for transitional flow. In fact, the typical wall normal grid spacings for WMLES based on turbulent boundary layer thickness, δ , might be too coarse to resolve the pre-transitional laminar boundary layer and its instability. Park and Moin 18 reported that, for WMLES in transitional cases, at least the size of the first grid cell near the wall �x + 2 < 20 is required to marginally resolve the integral length scales of the pressure-producing eddies near the wall. However, when using such a fine wall normal spacing, it may not always be possible to have the LES input taken at a distance even close to ∼ 0.1δ 19 . In the present study, because the flows are often transitional, the exchange location is taken at least three cells away from the wall, which lies at approximately 0.08-0.1δ. computational setup. The compressible Navier-Stokes equations are solved numerically using the massively parallel CharLES X solver. This solver implements a cell-centered finite volume scheme that is minimally dissipative. The Euler flux is computed by a blend of a non-dissipative central scheme and a dissipative upwind scheme where the blending parameter 0 ≤ f αf ≤ is precomputed based on the local grid quality. To avoid numerical instabilities, the dissipative upwind-flux contribution is significant only in the region of relatively poor grid quality 20 . Indeed, the proportion of the upwind flux f αf is not a static parameter, but it scales with the local departure of the global advection matrix (constructed with the central scheme) from a skew-symmetric matrix. For all the grids used in this study, the upwind proportion is less than 1.5% ( f αf < 0.015 ) in the regions with non-zero Reynolds stresses. Thus, the small-scale near-wall eddies and the large eddies passing through the separated shear layer, which are important in the dynamics of the separating and reattaching mechanisms, are largely unaffected by the numerical dissipation.
The numerical method is formally second-order accurate in space, although it achieves fourth-order accuracy on a uniform mesh containing only hexahedral cells. Time integration is performed using a third-order lowstorage Rung-Kutta-Wray scheme 21 . Unresolved turbulent scales are modeled using the constant-coefficient Vreman subgrid scale model 22 .
The computational grid around the controlled-diffusion airfoil is topologically an O-type mesh due to the round leading and trailing edges, with boundary layer clustering at the airfoil walls. The boundary layer clustering blocks are structured blocks, whereas the rest are unstructured to allow for a quick outward coarsening. The computational grid is also equipped with wake blocks with a large angle of opening. Note that independently of the angle of attack to be simulated, the blocks are not rotated; hence, the mesh is not changed, as the inlet velocity angle is instead adapted. The computational domain extends 20c in both the streamwise, x 1 , and wall-normal, x 2 , directions, to allow for a velocity inlet boundary condition to act as free-stream (see Fig. 1). Unphysical numerical reflections at the computational boundaries are avoided by the choice of appropriate boundary conditions. Characteristic boundary conditions are used at all inflow and outflow boundary conditions 23,24 . At the airfoil surface, an adiabatic, no-slip condition is applied. Periodic boundary conditions are imposed in the spanwise direction.
Although there are no experimental data dealing with the spanwise correlation length, Grebert et al. 25 found in their WRLES that a spanwise extent L x 3 of at least 2δ max showed reasonable agreement with experimental data for an accurate estimate of the correlation length ℓ z .
Three flow conditions from the matrix of computation proposed by the SCONE project 16 are studied here, and their details are summarized in Table 1. For C * ,c 1 cases, the resolution level is very similar to the standard DNS resolution criteria 26,27 ; see Table 2. For these simulations, the mesh resolution quality was also verified by Pope's criterion 27 , which monitors the quantity IQ k = K/(K + K SGS ) , where K and K SGS are the resolved and the SGS turbulent kinetic energies, respectively. On the one hand, the resolved turbulent kinetic energy is evaluated as K = (U where ν t is the kinematic turbulent viscosity, C M a constant set to 0.069 28 , and is estimated as the cubic root of   Table 2, we found that IQ k ≥ 0.94 for the three WRLES cases, which means that the resolution level is close to that of a DNS mesh. With the wall-model configuration, no near-wall streaks can be captured. In fact, the friction momentum flux is imposed by the model. Therefore, there is no need to follow any resolution requirement imposed by these structures. Note that when considering the x 2 direction, for example, the number of points remain essentially unchanged compared to WRLES. Instead, the first cell close to the wall increases in size. This ensures a proper resolution in the outer layer and drastically increases the time step because of the Courant-Friedrichs-Lewy (CFL) limit, which is directly dependent on this length scale. The overall cost ratio between the WRLES and the WMLES simulations is close to 50. The grid resolution used for the WMLES is chosen to resolve the flow scales in the outer layer and thus, the grid spacing is scaled by the local boundary layer thickness. As a consequence, the mesh does not strongly depend on the Reynolds number. As indicated in Table 2, the same grid is used for all WMLES computations. The cases labeled with C ⋆,c 3 all have the same grid resolution but the wall-model is not activated to investigate its effect at an iso-mesh resolution. The grid distributions in wall units for these computations along the chord are shown in Fig. 2. The grid sizes for all the WMLES cases have x 1 ≈ x 3 . In order to capture the boundary layer transition on the upper surface of the blade and achieve as smooth a flow as possible at the airfoil trailing edge, the grid spacing in the wall normal direction is fine, i.e., y + is less than 20 for the three cases. To allow errors due to the subgrid modeling and numerics to be made arbitrarily small, as shown in the study by Kawai and Larsson 14 , five grid points ( l wm = y 5 ) off the wall in the LES mesh are matched to the wall-model top boundary in this work. The variations of the height of the exchange location in viscous units are included in Fig. 2d. In the present study, the local wall-model layer thickness is set to a maximum of 4 times the local wall-normal grid-spacing and a user-defined minimum thickness l ud . The effects of the exchange location on the flow field, especially on the velocity profiles, will be investigated in future simulations.
Finally, it is important to emphasize that for the WMLES calculations, the non-dimensional time-step size �t = �t ⋆ U ∞ /c (where the superscripts ⋆ denote dimensional quantities) is larger due to the coarse grid resolution near the wall, and approximately two orders of magnitude bigger than that used for the WRLES computations. Furthermore, the maximum CFL number is ∼ 0.4.

Results
In this section, we simultaneously: (i) investigate the generation of airfoil trailing edge broadband noise that arises from the interaction of a turbulent boundary layer with an airfoil trailing edge and (ii) study the performance of the equilibrium boundary layer wall-model presented in "Equilibrium wall-model". Describing a flow as it passes a controlled-diffusion blade is challenging because it involves complex physical processes: laminar separation, turbulent transition, turbulent reattachment, and turbulent separation near the trailing edge on the suction side. We compare our numerical results with the experimental database generated within the framework of the CRORTET Clean Sky project 29 .  A fully laminar boundary layer is present on the lower (pressure) side up to the trailing edge of the airfoil, as well as a transitional and turbulent boundary layer on the upper (suction) side. A transition from a laminar to a turbulent state occurs in the shear layer resulting from the flow separation in the leading edge region, which results in the massive generation of vorticity downstream of the leading edge, with large vortices shed from the suction side of the airfoil. However, some smaller vortices still remain attached to the wall and roll over the airfoil suction side, grazing the trailing edge. As the curvature of the controlled-diffusion airfoil changes, the adverse pressure gradient leads to an increase of the boundary layer thickness. At the trailing edge, the laminar boundary layer coming from the pressure side destabilizes into a small vortex-shedding. This Von Kármán street then interacts with the fully turbulent vortical structures issuing from the upper side. Figure 4 depicts the level of vorticity and the size of turbulent structures in the flow at the same instant for all three C 1 simulations. We observe that the vorticity and the vortex shedding is more intense for the C 1,c 1 case (i.e., the WRLES) than for C 1,c 2 and C 1,c 3 cases. Furthermore, in the C 1,c 1 case, many more structures are resolved. The nature of the developed turbulent structures is broadly influenced by the size and the structure of the recirculation bubble near the leading edge (see Fig. 4). Furthermore, the WMLES (i.e., the C 1,c 2 case) does not show a laminar separation, and the flow becomes turbulent without the clear 2D vortex breakdown, as is shown in the C 1,c 1 and C 1,c 3 cases. Indeed, the vortices are first generated very close to the wall due to the flow instability induced by the sudden increase in the wall shear stress τ w (which is fed into the LES as flux boundary conditions) by activating µ t,wm in the wall-model at the leading edge. As shown in Fig. 5, by increasing both the Reynolds and the Mach numbers, the flow topology observed in the former case C 1 is also reproduced in the C 2 case, and we can see that, even when the flow seems to be similar at x 1 /c > 0.2 , the processes that reach the turbulent state are completely different. The comparison of the WMLES results with the reference WRLES configuration (the C 2,c 1 case) indicates that the equilibrium dynamic wall-model does not adequately capture the formation of two-dimensional intermittent laminar separation vortices, which are induced by the gradual adverse pressure gradient near the leading edge. Therefore, the wall-model failure directly translated into an incorrect estimate of the momentum flux and thus a different boundary layer thickening. This interesting feature is presented in Fig. 6, which shows that the wall-model simulation cannot capture the leading edge bubble properly. In fact, the laminar part of the boundary layer is treated with a wall-model, which enhances the friction and hence leads to an artificial growth of the boundary layer. Therefore, because the turbulent boundary layer has a weaker resistance against separation, the separation occurs earlier. For a very small recirculation region as in the current cases, the wall-model pushes the flow to recover a turbulent boundary layer. This numerical behavior is problematic and is the cause of the very small size of the recirculation region. However, the vortical structures at the trailing edge and the wake are similar to the reference WRLES case; see Figs. 4 and 5.
Mean flow profiles. In Fig. 7, we present the time-and spanwise-averaged pressure coefficient, C P , computed as  www.nature.com/scientificreports/ where P is the time-and spanwise-averaged wall pressure and P ∞ is the reference pressure taken at the outlet boundary. The whole pressure side of the blade is subjected to a favorable pressure gradient and the boundary layer is completely attached. On the suction side, the separation bubble creates a strong pressure drop close to the leading edge and large fluctuations of the wall shear stress are observed in the reattachment region. Downstream of the separation bubble, the mean pressure increases slightly up to mid-chord, then an adverse pressure gradient is observed up to the trailing edge. The mean flow features described in the previous section are recovered quantitatively with the numerical results, with a slight mismatch for C P at the beginning of the separation zone. The results obtained with the WMLES approach are very close to the results of the reference WRLES configuration. However, note that the discrepancies observed at the leading edge are not surprising since our numerical simulations do not account for the jet shear layers and the possible introduction of the free-stream turbulence that are present in the experiments 30 . We observe a small difference in the leading edge region between the WRLES and WMLES cases. However, we found very good agreement between the WMLES calculation and the reference WRLES for the surface pressure coefficient near the trailing edge region, and we found significant improvements between the iso-resolution configurations C 3,c 2 (with the wall-model) and C 3,c 3 (without the wall-model). Therefore, we can conclude that the WMLES approach is particularly suitable for the estimation of the leading edge turbulent boundary layer characteristics. At this stage, we must acknowledge that the accurate prediction of C P is a common finding in most external aerodynamics calculations. This result is most probably attributed to the outer-layer nature of the mean pressure, which appears to be less sensitive to the flow details of the near-wall turbulence 31 .
The profiles of the mean velocity normal to wall U ⊥ are shown in Fig. 8 for the cases C 1 and C 2 . The profiles of the case C 3 are roughly similar to the case C 2 ; unless otherwise stated, only the results of C 2 are shown. The velocity profiles show almost no deviation compared to the reference WRLES. The relative errors of the three mean velocities are less than 2% for the five locations considered. The turbulent boundary layer over the airfoil is about 10-20 mm thick, which, with our grid resolution, results in roughly 20-40 points per boundary later thickness, δ . Assuming that the flow at a given station can be approximated by a local canonical zero-pressuregradient turbulent boundary layer, the expected error in the mean velocities can be estimated as ε m = 0.4�/δ (with an isotropic grid size), which yields values of 2% for the current grid resolutions. The previous estimation www.nature.com/scientificreports/ is strictly valid for the zero-pressure-gradient turbulent boundary layer and, as such, it should be understood only as representative of the errors in complex geometries. However, this estimation provides a useful reference for the expected performance of WMLES in the absence of 3D and non-equilibrium effects for the current grid resolutions. In Fig. 9, the Reynolds stress component at five different locations are also plotted. The trend followed by all the stress components is correctly captured, although their magnitudes are systematically underpredicted by ∼ 18% . In fact, they present a slight deviation consistent with the modeling approach, which supports a limited fraction of the turbulent fluctuations. This behavior is typical for WMLES, which usually tends to overestimate the streamwise fluctuations while underestimating the transverse components. This is a direct effect of the grid resolution; the turbulent flow cannot carry identical turbulent structures in the near-wall region, which affects the Reynolds stress structures. However, the difference in the Reynolds stresses is mainly localized in a narrow  www.nature.com/scientificreports/ region very close to the wall. Thus, it should not significantly affect the pressure fluctuations whose maximum value is located further away from the solid wall. The prediction capability of the WMLES reduces close to the trailing edge for the streamwise component of Reynolds stress. Note that at iso-resolution, the resolved Reynolds stress computed by WMLES are in good agreement with the reference WRLES simulations compared to the case where the wall-model is not activated. In Fig. 10, the mean wall-parallel velocity profiles on the suction side are presented in log scale. The most basic wall-models are based on analytical laws that provide a direct link between the velocity normal to the wall U ⊥ at a certain wall distance x 2 and the wall-shear stress τ w = ρU 2 τ . The most widely used law is the Reichardt law-of-the-wall 32 : with U + ⊥ = U ⊥ /U τ and y + = yU τ /ν . Although this law is dedicated to equilibrium flows, it should be applicable to out-of-equilibrium flows as long as the input position is chosen below the end of the log-layer region. The constant κ and C are the same than those of the typical von Kármán classical law of the wall 33 Both analytical approaches are compared in Fig. 10. Notice that the wall-parallel velocity component profiles show a fast increasing boundary layer thickness after 50% of the chord, due to a strong adverse pressure gradient on the blade. At x 1 /c = 0.9 , the flow is nearly at the edge of the separation, but is still attached ( ∂P /∂x 2 is positive). All WMLES profiles show quite a large log-region, which collapse on the Reichardt law-of-the-wall for y + ≤ 300 . If the wall-model input is chosen below that threshold, the equilibrium wall-model in this study appears to be sufficient. Nagib and Chauhan 34 have emphasized the non-universality of the Kármán constant κ through various experiments on canonical flows and showed that FPG leads to higher values of κ while adverse  www.nature.com/scientificreports/ pressure gradient gives lower κ . They also showed that the constant C can be negative. Therefore, the comparability between our results and the two proposed analytical approaches could be improved by finding an optimal choice of κ and C values. In Fig. 11, we show the variation of the shape factor ( H 12 = δ ⋆ /θ ), the momentum thickness θ , and the displacement thickness δ ⋆ against the chordwise distance x 1 /c along the aerofoil upper surface. In the case of aerofoil flows, the traditional definition of boundary layer parameters as used for flat plates must be adapted to be consistent with the spatially varying potential flow. In the present study, the definitions used for the displacement thickness δ ⋆ and the momentum thickness θ follow the recommendation of Wagner et al. 35 . Near the leading edge, the flow is laminar and the shape factor is approximately 2.6 for C * ,c 1 cases, which is a typical value of the Blasius profile, while the C * ,c 2 cases feature a higher value ∼ 2.9 (zoom not shown). The displacement thickness reaches a local maximum in the region of the separation bubble while the momentum thickness is nearly zero, which results in a very large shape factor. Downstream of the reattachment region ( x 1 /c ≃ 0.03 for all cases), the shape factor decreases abruptly, and then close to the mid-chord where there is close to a zero pressure gradient, a plateau is observed at ∼ 1.5 for C 1 , 1.35 for C 2 and ∼ 1.28 for C 3 , which is smaller than the typical value of the Klebanoff profile H 12 ≃ 1.4 observed for turbulent boundary layers). Along the second half of the chord, both the two thicknesses and the shape factor increase because of the adverse pressure gradient. A larger difference is observed between C 2 and C 3 . Those differences depend both on the Reynolds number and the pressure gradient 36 . Acoustic field. Analyzing the divergence of the velocity field is particularly useful to delineate the sources and the strength of the acoustic waves in the flow field. An instantaneous picture of the field of ∇ · U on the midspan plane is shown in Figs. 12 and 13 along with contours of iso-entropy s = log(P /ρ γ ) . The C 3 simulation is found to be qualitatively similar to the C 2 simulation. Two major acoustic source regions can be visually identified: the trailing edge and the transition/reattachment zones on the upper surface. Contrary to the study of Wu et al. 37 , the near wake does not appear to contribute significantly to the self-noise. The transition/reattachment noise source appears to be even stronger than the trailing edge noise for this flow configuration. This is likely caused by the high local Mach number due to the strong flow acceleration near the leading edge where the blade curvature is higher. An additional sound source is visible near the leading edge on the suction side of the lifting surface. This sound source is mainly caused by the massively accelerated vortex structures near the leading edge 38 . The transition/reattachment noise source radiates mostly downstream as it is shielded by the upper blade surface (with only a slight diffraction at the leading edge), whereas the trailing edge source shows an anti-symmetric pattern on each side of the blade, as expected from classical trailing edge noise theories 39 . Note Figure 11. Variations of shape factor H 12 , momentum thickness θ , and displacement thickness δ ⋆ with chordwise distance x 1 /c between WRLES computations (solid lines, cases C * ,c 1 ) and WMLES computations (dashed lines, cases C * ,c 2 ).
Scientific RepoRtS | (2020) 10:12735 | https://doi.org/10.1038/s41598-020-69671-y www.nature.com/scientificreports/ the existence of a noise source near the trailing edge. The vortical structures in the turbulent boundary layer pass the trailing edge and generate sound. Visually, the WMLES approach seems to reproduce a qualitatively similar acoustic source field to the reference WRLES simulations. For the higher Mach number cases C 2 -C 3 , we notice a stronger radiation in the presence of the separating shear layer, which generate additional noise sources and affect the propagation of acoustic waves upstream by the refraction through the separated shear layer. The secondary source, at higher frequency (closer wave fringes), appears less on the suction side, near the leading edge, close to the reattachment point of the recirculation bubble when using WMLES, which results from the fact that the current equilibrium WMLES does not show properly the laminar separation.
Frequency spectra of surface pressure fluctuations. Here we assess the capability of wall-modeled LES to predict surface pressure fluctuations on the suction side of the airfoil surface at x 1 /c = 0.976 , as shown in Fig. 14. The frequency spectra of pressure fluctuations (PSD) as well as the cross-spectral and auto-spectra have been estimated numerically based on the Welch method of periodogram 40 , which minimizes the variance of the PSDs estimator 41 . The overall pressure signals from our computations, extracted from probes that are placed in the first wall-normal cell along the airfoil, are subdivided into N s /16 sub-blocks, where N s represents the number of pressure recordings of each case. We adopt the Welch method combined with a Hanning window throughout the spectral analysis while using a FFT. The overlapping of the sub-blocks is set to 50%. The frequency spectrum is then obtained by averaging the periodograms of all the sub-blocks. The normalized sampling frequency is set to 3,000, i.e., we collect 3,000 pressure signals per each probe and convective time scale. The PSD is plotted as a function of the streamwise location and the normalized frequency (or Strouhal number St = fc/U ∞ ). The trailing edge noise displays a low frequency broadband spectrum with most of the energy of the signal below St = 9 . The PSD of each case features noticeable differences, as indicated by the slopes that are superimposed onto each dimensionless PSD. In general, the WMLES is found to underpredict the spectral level of the experimental data by 2-4 dBs in the low to intermediate range of frequency. The high range of frequency of the experimental pressure spectra is more accurately captured by the WMLES compared the reference WRLES's. The most striking finding in this figure is that at iso-resolution, WMLES cases are in better agreement than the cases without the wall-model, mainly at the low frequency range. Overall, all of the spectra in Fig. 14 demonstrate that the broadband content is predicted reasonably well by the WMLES, which provides a good baseline study to assess and quantify noise prediction. The tonal peaks that arise from the experimental spectra on Fig. 14 are attributable to acoustically untreated duct models of the wind-tunnel facility. Moreover, the observed disparities at low Strouhal number may be partly caused by the installation effects of the experiment that are not accounted for in the numerical simulations, but also to a lack of statistical convergence of the WRLES spectra as a result of the time series being short. The wall-pressure statistics computed from the time-domain analysis are deduced from the temporal pressure cross-correlation calculated as (11) R a,b (ξ 1 , 0, ξ 3 , τ ) = �P ′ (x 1 , 0, x 3 , t)P ′ (x 1 + ξ 1 , 0, x 3 + ξ 3 , t + τ )�,  www.nature.com/scientificreports/ where P ′ is the surface pressure fluctuations, a and b represent two arbitrary points on the airfoil separated by distance ξ 1 , ξ 3 in the streamwise and spanwise directions, respectively. The brackets �·� indicate an ensemble average. The time delay between two signals τ is normalized on the chord c, and the inflow velocity U ∞ . The effect of streamwise separation is depicted in Fig. 15. The cross-correlation R a,b (ξ 1 , 0, 0, τ ) is computed using the probe located at x 1 /c = 0.975 as a reference for the computation of correlation with upstream probes that are separated by ξ 1 . We note the diminution of the correlation between the two signals with increasing separation between the signals. The decreasing peak of cross-correlation is due to the change in the pressure signature in the turbulent boundary layer, since the greater the distance between the probes, the greater the distance over which the structures can change and evolve, which yields a lower correlation between the pressure signals. For the WMLES cases, we notice that the pressure signals are correlated over a larger separation distance, and thus for a higher time delay relative to the reference WRLES configurations. This can be explained by the fact that the structures evolve less rapidly when using the wall model, i.e., the structures are constrained by using the wall-model, leading to much lower correlation value of R a,b (ξ 1 , 0, 0, τ ).
Using the function R a,b (ξ 1 , 0, 0, τ ) , it is possible to estimate the value of convection velocity of the vortical structures U c (ξ 1 ) as a function of the separation distance ξ 1 in the streamwise direction between probes as where τ U edge /δ ⋆ max represents the time lag, which corresponds to the maximum of the cross-correlation R a,b (ξ 1 , 0, 0, τ ) between two pressure signals separated by a normalized distance ξ 1 /δ ⋆ .
In Fig. 16, the normalized convection velocity is extracted as a function of the longitudinal separation ξ 1 , taken as the reference probe located at x 1 /c = 0.978 , for both the reference WRLES and WMLES computation. The convective velocity ranges from 0.55 · U ∞ to 0.80 · U ∞ with an apparent dependency on the Reynolds number or angle of attack. WMLES seems to cause a higher streamwise convective velocity U c with respect to the baseline WRLES configuration. www.nature.com/scientificreports/ cross-spectral density and coherence. One of the main quantities generally used to describe wall-pressure statistics is the cross-spectral density in the spanwise direction, which is defined by where P ′ is the time FFT transform of the fluctuated wall pressure. The complex conjugate is denoted by · * . According to the coherence function, which is used to indicate the strength of the relation between the pressure fluctuations at two separate locations as the frequency varies, as the disturbed flow is convected in the streamwise (13) Ŵ P a P b (ξ 1 , 0, ξ 3 , ω) = � P ′ * (x 1 , 0, x 3 , ω), P ′ (x 1 + ξ 1 , 0, x 3 + ξ 3 , ω)�, www.nature.com/scientificreports/ or in the spanwise direction, it is possible to express the cross-spectrum in non-dimensional form. Hence, the coherence function can be expressed as where Ŵ P a P a , Ŵ P b P b are the pressure auto-spectra at the two positions a and b, respectively. The different analytical expression coherence functions proposed in the literature are a function of the angular frequency ω (or equivalently the Strouhal number St = ωδ/u τ ), and normalized with respect to ξ 3 and U c , which are computed from the cross-correlation functions. Despite its simplicity, the Corcos model 42 is still extensively used in many practical applications, since it provides a simple analytical form and clear physical significance. This model is based on the assumption that the correlation lengths in both the streamwise and spanwise directions are statistically independent, and thus the coherence takes the form where k c = ω/U c is the convective wavenumber. In the formulation, α 1 and α 3 are constants that measure the loss of coherence in the streamwise and spanwise directions, for which different values can be found in the literature.
A second semi-empirical model, which is an extension of the Corcos model, has been proposed by Efimtsov 43 and improved by Salze et al. 44 . This model explicitly considers the compressibility effects by introducing a number of additional tunable constants. The coherence function is expressed for this model, according to the multiplicative hypothesis (cf. (15)), as where 1 and 3 are the correlation lengths in the streamwise ( � 1 = U c /ωα 1 ) and spanwise direction ( � 3 = U c /ωα 3 ), respectively. For a free stream Mach number below 0.75, the spanwise correlation length is modeled as follows 44 (14) γ a,b (ξ 1 , 0, ξ 3 , ω) = |Ŵ P a P b (ξ 1 , 0, ξ 3 , ω)| Figure 16.  45 , the Efimtsov-Salze's constant a 4 constrains the amplitude of mid-to high-frequency coherence lengths, similar to the constants in the Corcos model. The frequency at which the Efimtsov-Salze model breaks away from the Corcos model is governed by the parameter a 5 , while the low-frequency roll off is controlled by the parameter a 6 . The major limitations of the Salze-Efimtsov model are the tunable constants a 4,5,6 , which generally depend on the Mach and Reynolds numbers, as well as on the separation distance. These parameters are estimated in this study at x 1 /c = 0.976 by fitting the coherence data points obtained from the reference WRLES configurations, by means of a nonlinear least-squares optimization procedure. The coefficients thus obtained in the range of Mach numbers considered are a 4 = 0.47 , a 5 = 17.80 and a 6 = 1.0 . We depict in Fig. 17 the spanwise coherence length at different locations close to the trailing edge for the three cases. We see that the collapse of the curves for both the Corcos or Salze-Efimtsov models is achieved at high frequencies for all the tested computations at x 1 /c = 0.976 . We obtain less satisfactory agreement between the theoretical predictions and the numerical data far upstream x 1 /c = 0.976 (for which the parameters of Salze-Efimtsov model are computed). Nevertheless, by comparing spectra obtained for the two locations x 1 /c = 0.909 and x 1 /c = 0.976 , it can be concluded that very close to the trailing edge, the wall-pressure statistics are well established and almost stationary justifying the use of radiation models based on a wall-pressure statistics at a single point close to the trailing-edge 46,47 . The Salze-Efimtsov model predicts the peak in coherence lengths and compares reasonably well to the Corcos model. The Salze-Efimtsov also achieves a much better agreement with the present simulated dataset at lower frequencies. The lack of agreement for the Corcos model at low frequencies is predicted, since the Corcos model is not suited for this range of frequencies.
Moreover, good agreement between the predictions of the Salze-Efimtsov model and the numerical data is obtained in all cases. Unlike the Corcos approximation, the Salze-Efimtsov model is able to provide a good agreement in the low-frequency range of the coherence functions, especially at x 1 /c = 0.976 for all the matrix computation. We note that the comparison of coherence length distribution at Re = 2.29 × 10 6 between both models and the simulated results is better than the case Re = 8.30 × 10 5 , suggesting a Reynolds number dependence of the shape of these models. We can also see that the WMLES curves fit the reference WRLES computations reasonably well, regardless of the studied case.  The prohibitive computational cost of resolving the inner region of turbulent boundary layers prevents using the WRLES approach to study high-Reynolds number turbulent flows in complex geometries. However, using the wall-modeled LES approach can reduce this computational cost. The major objective of this work was to evaluate the feasibility and accuracy of the WMLES approach for the prediction of broadband noise generated by a controlled-diffusion blade. To this end, we carried out detailed comparisons between the results obtained with the equilibrium WMLES approach, experimental measurements, and the wall-resolved LES reference solutions. LES with wall-modeling accurately predicted the pressure coefficient distribution, velocity statistics (including the mean velocity), and the trace of Reynolds tensor components. We also found that, for the wall-modeled LES cases, the instantaneous flow structures resembled those observed in the reference wall-resolved LES, except near the leading edge region. In addition, the boundary layer thickness, displacement thickness, and momentum thickness along the airfoil chord showed a convergent behavior towards those obtained using the wall-resolved LES. Surface pressure fluctuations, which act as the sources of the broadband noise at the far-field, were also extracted in the form of surface pressure spectral densities along the airfoil chord for both the wall-modeled LES and wall-resolved LES calculations. Our results indicate that the pressure spectral density profiles from the wall-modeled LES compare well with the experimental profiles. Furthermore, we found that wall-modeled LES is more comparable with the experimental data in the high-frequency region than the WRLES reference case.
The convection velocity showed an increase towards an asymptotic value as the separation distance, ξ 1 , increased, and we found that the WMLES results overestimated the WRLES computations. We compared the spanwise correlation length to the Corcos 42 model, which assumes a frequency-independent correlation length, and the Salze-Efimtsov model 44 , which, in contrast, introduces a frequency dependence. We found both WMLES and WRLES computations compare well with both models over the higher range of frequency, which accounted reasonably well for the effects of the pressure gradients. Nevertheless, the Salze-Efimtsov model was in better agreement with the computational results for the lower range of frequency, but two parameters needed to be tuned. Overall, the equilibrium WMLES presents an interesting approach to (i) investigate in more detail the turbulent pressure field induced by turbulent boundary layer over the controlled-diffusion airfoil, and (ii) advance the development of new versions of trailing edge noise models on coarser grids, which is needed for fast simulations in industrial design.