Analytical models of the wall-pressure spectrum under a turbulent boundary layer with adverse pressure gradient

This paper presents a comprehensive analytical approach to the modelling of wall-pressure fluctuations under a turbulent boundary layer, unifying and expanding the analytical models that have been proposed over many decades. The Poisson equation governing pressure fluctuations is Fourier transformed in the wavenumber domain to obtain a modified Helmholtz equation, which is solved with a Green’s function technique. The source term of the differential equations is composed of turbulence–mean shear and turbulence–turbulence interaction terms, which are modelled separately within the hypothesis of a joint normal probability distribution of the turbulent field. The functional expression of the turbulence statistics is shown to be the most critical point for a correct representation of the wall-pressure spectrum. The effect of various assumptions on the shape of the longitudinal correlation function of turbulence is assessed in the first place with purely analytical considerations using an idealised flow model. Then, the effect of the hypothesis on the spectral distribution of boundary-layer turbulence on the resulting wall-pressure spectrum is compared with the results of direct numerical simulation computations and pressure measurements on a controlled-diffusion aerofoil. The boundary layer developing over the suction side of this aerofoil in test conditions is characterised by an adverse pressure gradient. The final part of the paper discusses the numerical aspect of wall-pressure spectrum computation. A Monte Carlo technique is used for a fast evaluation of the multi-dimensional integral formulation developed in the theoretical part.

(where the curvature will typically generate adverse free-stream pressure gradients) are the cause of structural vibrations and noise that propagate inside the vehicle. Aircraft fuselage excitation due to pressure fluctuations is a significant source of cabin noise. Also, wall-pressure fluctuations can interfere with the signal of underwater sonar domes. More importantly for the purpose of this work, the pressure disturbances generated on the surface of wing profiles are scattered into acoustic waves at the trailing edge. This is a source of noise in low-speed fans and wind turbines. Moreover, it is the only remaining source of broadband noise for a subsonic fan operating in a homogeneous stationary flow. Due to the random nature of these disturbances, the generated noise is broadband and its intensity is directly proportional to that of the wall-pressure fluctuations, as expressed in analytical theories, among others Amiet (1976) re-addressed by .
The experimental investigation of wall-pressure fluctuations dates back to the 50s (see the reviews of Bull (1996) and Cohen & Gloerfelt (2018)) and still proves challenging nowadays. For instance, Salze et al. (2015) discussed the separation of the hydrodynamic and acoustic components of the spectrum of pressure fluctuations from the raw experimental data. In particular, the acoustic component has its energy peak at low wavenumbers for a given frequency and thus it is more likely to excite the structural modes of the underlying surface. It must also be remembered that the measured frequency-wavenumber spectrum is a convolution of the true spectrum of pressure fluctuations with the response function of the experimental apparatus. Prigent et al. (2018) discussed the application of various numerical techniques for the deconvolution of the physical spectrum from the response function, although this remains an open problem. Another question that has been long debated is the choice of the relevant boundary-layer parameters for the scaling of measured wall-pressure spectra, which is the first step towards building empirical predictive models. This question has been treated in the works of Keith, Hurdis & Abraham (1992) and Cipolla & Keith (2000).
Direct numerical simulation (DNS) of the turbulent boundary layer can provide plenty of information concerning the statistical properties of turbulence and of the wall-pressure fluctuations that it generates. Earlier attempts at DNS date back to the eighties with the work of Spalart & Leonard (1987), using a spectral method. Reynolds numbers representative of industrial applications have long been out of the reach of DNS codes until the advent of supercomputers, as seen in Borrell, Sillero & Jiménez (2013). However, DNS is still far from being considered a viable design tool. A compromise in terms of computational cost for the generation of broadband noise sources is the fast random particle mesh method, which superimposes a synthetic field of turbulent fluctuations to the mean flow calculated by a Reynolds-averaged Navier-Stokes (RANS) simulation (see Proskurov, Darbyshire & Karabasov 2017). However, this method depends on the assumed functional expression of the turbulence energy spectrum, as shown by Wohlbrandt et al. (2016), similar to the family of methods that is presented in this paper.
The fastest way to estimate the wall-pressure power spectral density (PSD) for noise prediction is to use a model, either empirical or based on the analytical solution of the Poisson equation. Most empirical models are developed by selecting several databases of wall-pressure spectra in similar flow conditions and normalising them in order to highlight the dependence on some boundary-layer parameters. Then, a mathematical expression is developed to fit the collapsed curves. Corcos (1964) made one of the first attempts to build a universal statistical model of the wall-pressure spectrum, which has been widely used to this day due to its simple mathematical formulation but is known to greatly overestimate the contribution to a given radian frequency of the wavenumber range below the convective ridge. Corrections to the low-frequency behaviour of Corcos' model have been proposed, for instance, by Ffowcs Williams (1982) and Caiazzo, D'Amico & Desmet (2016). In the work of Schlinker & Amiet (1981), the pressure PSD is normalised by the outer boundary-layer parameters, the displacement thickness δ * and the external velocity U e , and plotted as a function of the corresponding reduced frequencyω = ω δ * /U e . This model is suitable for flat-plate, zero-gradient boundary layers, but it is not able to represent the quadratic rise with ω in the low-frequency range, predicted by Kraichnan (1956b) and Phillips (1956), or the increase of the PSD magnitude under the effect of an adverse pressure gradient. The model of Chase (1980), expressed in terms of both outer and inner boundary-layer variables (with the introduction of the wall shear stress, τ w ) represents a step forward because it accounts for the low-frequency rise. More recently, the model of Goody (2004), also expressed in terms of mixed inner and outer generalised variables, has included the effect of the Reynolds number on the extension of the overlap zone observed in the middle frequency range, where the normalised pressure PSD decreases approximately as ω −0.7 . Still, these models are calibrated on zero-pressure-gradient flow databases. The model of Rozenberg, Roger & Moreau (2008), on the contrary, preserves the main features of Goody's but also takes into account the effect of the adverse pressure gradient, which can increase the wall-pressure spectrum of up to 10 dB in the low-frequency range. Initially this model was formulated with reference to τ w , which can be difficult to estimate in some cases. However, its latest formulation published in Rozenberg, Robert & Moreau (2012) takes into account additional features of the boundary layer and substitutes the wall shear stress with the maximum shear stress in the flow. Catlett et al. (2015) generalised Goody's model in order to capture adverse-pressure-gradient effects as well. Hu's model (Hu 2018) aims at representing the wall-pressure frequency spectrum under arbitrary non-equilibrium flow conditions by using the shape factor H = δ * /θ for non-dimensionalisation, a parameter that is also representative of boundary-layer development history. The most recent and comprehensive comparison of all published empirical models is given by Lee (2018), who propose a formulation that is valid for a wide range of flow conditions. The Poisson equation governing pressure fluctuations in a turbulent boundary layer can be used to build models that relate pressure and velocity statistics through relatively simple analytical formulations. Heisenberg (1948) and Batchelor (1951) addressed the problem of pressure fluctuations in homogeneous and unbounded turbulence, developing equivalent formulations. Subsequently, Kraichnan (1956a) added an infinite plane boundary condition to the homogeneous turbulence problem and then included the effect of anisotropy of the length scales of turbulence as a first step towards the understanding of a real boundary layer. This theory was then extended to the case of non-homogeneity in the direction normal to the wall by Kraichnan (1956b) and Hodgson (1961). Whereas Kraichnan's theory is formulated in the space domain, Hodgson adopted a wavenumber-domain approach. The results of these studies were essentially analytical and based on idealised flow models. Panton & Linebarger (1974) coupled the analytical wavenumber-domain solution of the Poisson equation with a realistic flow model, from which arose the necessity to develop a numerical integration method to compute the wall-pressure spectrum. Their approach was later taken up by Remmler et al. (2010), who used a RANS flow solution to provide flow input data to the model. Peltier & Hambric (2007) and Slama, Leblond & Sagaut (2018) coupled the space-domain solution of the Poisson equation with steady-state flow simulations. These different approaches will be discussed in the first part of this paper. Unlike the previous authors, Parchen (1998) developed a prediction method based on an approximate solution of the Poisson equation, following the theory developed by Blake (1986), which gave rise to the so-called TNO-Blake family of models, (see Kamruzzaman et al. 2012;Bertagnolio, Fischer & Zhu 2014;Stalnov, Paruchuri & Joseph 2016;Fischer, Bertagnolio & Madsen 2017). The difference between the exact and the approximate solutions in terms of accuracy of the predicted spectrum will be discussed as well in this work. The purpose is to provide the most comprehensive analytical approach to the prediction of the wall-pressure PSD, unifying all other approaches in a common framework. In particular, it will be shown that the fundamental problem consists of the representation of the two-point statistics of the turbulent velocity across the boundary layer, taking into account inhomogeneity and anisotropy effects. The application case selected for this study is a controlled-diffusion aerofoil at a free-stream Mach number of M = 0.05 and a chord-based Reynolds number of Re c = 1.5 × 10 5 , which presents an equilibrium turbulent boundary layer with an adverse pressure gradient on its suction side. This also represents an improvement, since most currently used models are calibrated on zero-pressure-gradient boundary layers (see Slama et al. 2018). Turbulence and wall-pressure statistics are extracted from a state-of-the-art DNS in order to validate the different steps in the construction of the model. Flow and wall-pressure statistics extracted from the DNS are compared for validation with measurements performed in the anechoic wind tunnel of the Université de Sherbrooke.
The computation of wall-pressure spectra with formulations based on the Poisson equation may present difficulties from the numerical point of view due to the presence of multi-dimensional integrals. For this reason, the last section of this work presents a Monte Carlo integration technique which can tackle this problem efficiently and at a low computational cost.

Poisson equation solution
The Poisson equation governing the pressure fluctuations in a turbulent boundary layer derives from the divergence of the incompressible momentum equation, introducing Reynolds decomposition into mean and fluctuating quantities, then subtracting the time-averaged equation. As a result, the Laplacian of pressure fluctuations is equal to the sum of two source terms, where i, j = 1, 2, 3, with 1 the streamwise, 2 the wall-normal and 3 the transverse direction. According to (2.1) the pressure fluctuation, p, can be considered as generated by two separate mechanisms (see the review by Gerolymos, Sénéchal & Vallet (2013)). One mechanism is the turbulence-mean shear interaction, expressing the fact that the corresponding pressure fluctuations react immediately to a change of the mean flow gradient and are linear with respect to velocity fluctuations. The other mechanism is the nonlinear interaction of different components of the velocity fluctuations, which is not immediately affected by a change in mean flow gradient. Equation (2.1) is solved with the Dirichlet boundary condition lim whereby the pressure fluctuation vanishes at an infinite distance from the wall and the overall pressure tends to its mean free-stream value, and the Neumann boundary condition ∂p ∂x 2 x 2 =0 = 0 (2.3) following from the rigid wall assumption. The space-domain solution of (2.1) is expressed within a Green's function approach as where the function satisfies the boundary conditions for an impulse source. The partial derivatives applied to the velocity variables in (2.4) can be transferred to the Green's function by repeated applications of integration by parts and of the divergence theorem, yielding the equivalent formulations Equations (2.4), (2.6) and (2.7) are mathematically equivalent, so the choice is motivated by considerations of numerical stability of the integration algorithm. Peltier & Hambric (2007) and Slama et al. (2018) selected equation (2.7). However, Hu, Reiche & Ewert (2017) argued that (2.6) provides the best numerical realisation of the turbulence-turbulence interaction contribution. In any case, the spacedomain approach prescribes a three-dimensional integral for the computation of p, whereby the cross-correlation of wall-pressure fluctuations corresponds to a six-dimensional integral, which must still be Fourier transformed numerically to obtain the wavenumber-frequency spectrum that is needed in acoustic predictions. Performing these numerical integrations with limited time and computational resources would inevitably lead to a loss of accuracy. For this reason, it is preferable to derive an analytical solution for the unsteady pressure directly in the wavenumber domain.

Modified Helmholtz equation 2.2.1. Turbulence-mean shear interaction term
In this section, the wavenumber-domain expression of the unsteady pressure at the wall is derived assuming that the only significant contribution is given by the turbulence-mean shear interaction component. This approach has its roots in the analytical considerations of Kraichnan (1956a) and Hodgson (1961). More recently, however, Slama et al. (2018) argued that the two components have the same order of magnitude, based on the numerical integration of (2.7). The derivation of a wavenumber solution for the turbulence-turbulence interaction term is more cumbersome and requires a set of more restrictive hypotheses in order to get a simple closed-form formulation. For this reason, the turbulence-turbulence interaction component will be presented separately in § 2.4.
Assuming homogeneous turbulence in planes parallel to the wall, equation (2.1) can be Fourier transformed in space in the (x 1 , x 3 ) directions, yielding the following modified Helmholtz equation: where the following convention for the Fourier transform is adopted: Furthermore, assuming that the only non-vanishing mean velocity gradient is that of the streamwise component, U 1 , in the wall-normal direction, x 2 , the source term is reduced to It is demonstrated in appendix A that the wavenumber-domain solution for the pressure at the wall iŝ (2.14) The PSD of wall-pressure fluctuations is equal to the ensemble average of the product of (2.14) by its complex conjugate, The ensemble average of vertical velocity fluctuations appearing in (2.15) can be expressed as The main point in (2.17) is the correct analytical representation of the vertical velocity fluctuation cross-spectral density, ϕ 22 . The remaining terms, namely the mean velocity gradient and the vertical velocity root-mean-square (r.m.s.), are relatively easier to calculate. They can be obtained from empirical models, as in Panton & Linebarger (1974), from hot-wire measurements or from steady-state simulations, as in Remmler et al. (2010). Although the time variable has been omitted in the derivation of the solution of the modified Helmholtz equation, it is understood that the variablep is dependent on time. The wavenumber spectrum of (2.17) can be converted to a frequency spectrum first by integrating over the k 3 wavenumber and then by adopting Taylor's hypothesis (see Taylor 1938) of frozen convection of the turbulent eddies in the flow direction. Assuming constant convection speed of wall-pressure fluctuations, U c , in a close range around the point where the frequency spectrum is sought, the sole streamwise wavenumber contributing to a given radian frequency, ω, is K c = ω/U c . Therefore, the frequency spectrum of wall-pressure fluctuations is (2.18) 2.3. Turbulence modelling The analytical expression of the cross-spectral density of vertical velocity fluctuations, which closes equation (2.17), is based on a theory of isotropic turbulence which, following Panton & Linebarger (1974), is modified in two ways. In the first place, the integral length scale, Λ, is allowed to vary according to the distance from the wall. Second, the stretching of the turbulent structures in the flow direction is taken into account with a coefficient α representing the ratio of streamwise to transverse integral length scales (see also Schlinker & Amiet 1981). The cross-spectral density of anisotropic turbulence, ϕ NI 22 , is related to the isotropic spectrum, ϕ I 22 , as The ratio of length scales, α, varies across the boundary layer with the distance from the wall, x 2 . Remmler et al. (2010) postulate the existence of a function relating each wavenumber k 1 to a given value of α. However, the wall-pressure PSD at a given wavenumber k 1 results from a double integration over the distance from the wall (equation (2.17)) so that the result at each k 1 does not depend on a single value of α, but rather on the variation of α with x 2 . A characteristic of boundary-layer turbulence that is not included in this model is the 45 • orientation of the vorticity vectors due to mean straining, the effect of which was deemed by Kraichnan (1956b) negligible with respect to that of length scale anisotropy. However, this aspect has been taken into account in the models of Peltier & Hambric (2007) and Slama et al. (2018).
The expression of ϕ I 22 follows from the choice of the analytical expression of the longitudinal correlation function of turbulence, F(r) (see (B 13) and (B 17)). Two analytical expressions are considered in this work: the Gaussian and the generalised von Kármán function described in the following. The complete derivation of the corresponding cross-spectral density expressions can be found in appendix B.

Gaussian energy spectrum
The Gaussian longitudinal correlation function is defined as (2.20) with the characteristic length complying with the definition of Λ given in (B 3). The corresponding cross-spectral density of vertical velocity fluctuations is This is possibly the simplest mathematical formulation available in the literature. However, its accuracy has been questioned already by Panton & Linebarger (1974) and, more recently, by Kamruzzaman et al. (2011), for instance.

Generalised von Kármán energy spectrum
The generalised von Kármán correlation function is defined by Wilson (1997) as with the characteristic length (2.24) The corresponding cross-spectral density of vertical velocity fluctuations is with the parameter (2.26) The limit of (2.25) for x 2 → x 2 is the generalised auto-spectrum (2.27) This formulation is particularly interesting for its flexibility: different values of the parameter ν represent some of the most widely used turbulence models. First, ν = 1/3 corresponds to the original model of von Kármán (1948). This turbulence model is also at the core of the TNO-Blake family of models, as will be discussed in the next section. It can be easily verified that, for ν = 1/2, equation (2.23) corresponds to (2.28) This exponential longitudinal correlation function is the basis of Liepmann's turbulence model (see Liepmann, Laufer & Liepmann 1951). Interestingly, also Panton & Linebarger (1974) develop their formulation of the cross-spectral density of vertical velocity fluctuations from (2.28). However, instead of simply using (2.25) with ν = 1/2, they chose to formulate ϕ 22 as a double integral (see (B 10)). As a result, equation (2.18) becomes a five-dimensional integral, which can be more challenging from a numerical point of view. A possible method for the fast and accurate integration of this kind of high-dimensional formulation is presented in § 5. Finally, for ν = 7/6 the generalised formulation corresponds to the rapid distortion theory (RDT) proposed by Hunt (1973). Although this model was originally developed to account for the modification of the energy spectrum when homogeneous turbulent structures are distorted and elongated by impingement on a bluff body, it has recently been applied by Goldstein, Leib & Afsar (2017) to the study of the noise scattered by the trailing edge of a flat plate interacting with a jet parallel to its surface. Furthermore, the hypothesis that the turbulence-turbulence interaction is negligible with respect to turbulence-mean shear interaction is at the basis of RDT as well (see Bailly & Comte-Bellot 2015). It has been suggested by Christophe (2011) that the RDT formulation may be more accurate than the classical von Kármán spectrum in the low-wavenumber range (see also De Santana et al. 2016). In his dissertation de la Riva (2001) has also shown the capabilities of RDT in predicting the development of turbulence convected through a highly staggered cascade propulsor configuration formed by non-symmetric aerofoils (see also de la Riva et al. 2004).
The longitudinal correlation functions presented in this section are depicted in figure 1. The correlation curves tend to drop faster to zero going from the Gaussian to the von Kármán function. Furthermore, according to (B 3), the integral of every curve for r going from zero to infinity is equal to Λ, so that the curves that go faster to zero show higher correlation values for small r.

TNO-Blake model
The family of TNO-Blake models makes use of the von Kármán turbulence spectrum in order to compute the wall-pressure spectrum. However, the original TNO-Blake model (see Parchen 1998;Bertagnolio et al. 2014) approximates the cross-spectrum of vertical velocity fluctuations as where L 2 (x 2 ) is a vertical correlation length scale, and Φ m is the so-called moving-axis spectrum, representing the contribution of the wavenumber k 1 to a given frequency ω. In most cases, Taylor's hypothesis is applied and the moving-axis spectrum is expressed as (2.30) The fundamental assumption behind the development of (2.29) is that which can be hardly justified from a physical point of view. In fact, as will be shown later in § 4.3, the vertical velocity cross-correlation coefficient is not negligible over a vertical distance corresponding to at least two integral length scales. For this reason, Fischer et al. (2017) presented an improved version of the TNO-Blake model which corresponds to the substitution of the exact equation (2.25) (with ν = 1/3) in (2.17). However, this hypothesis has the only advantage of reducing the number of dimensions of integration of the wall-pressure spectrum formulation. Then, the substitution of (2.29) and (2.30) in (2.17) yields the wavenumber wall-pressure spectrum, (2.32) Finally, assuming equation (2.30), the corresponding frequency spectrum is calculated by coupling equation (2.32) with (2.18).
The treatment of the anisotropy of turbulence length scales in the TNO-Blake family of models is similar to that presented in (2.19). However, in this case there are three coefficients, representing the anisotropy of the length scales in the streamwise, wallnormal and transverse directions, respectively. The anisotropic three-dimensional autospectrum of vertical velocity fluctuations is ϕ a,NI 22 (k 1 , k 2 , k 3 ) = β 1 β 2 β 3 ϕ a,I 22 (β 1 k 1 , β 2 k 2 , β 3 k 3 ). (2.33) Integrating the previous equation in dk 2 yields the two-dimensional auto-spectrum used in (2.29) ϕ a,NI Various values for these coefficients have been proposed over time, as can be seen in table 1. The reason for the discrepancy in the values is that each set of parameters has been tuned empirically to a different experimental database. It can be noticed that Bertagnolio et al. (2014) expressed the vertical and transverse anisotropy coefficients as a function of the free-stream pressure gradient parameter where dP/dx 1 is the gradient of the mean wall static pressure with respect to the abscissa along the aerofoil suction side. However, this dependence has been dropped in the later work of Fischer et al. (2017). None of the sets of coefficients presented in table 1 takes into account the variation of the anisotropy characteristics with the distance from the wall, a fact that is not negligible as will be shown in § 4.3. The original TNO-Blake model is also characterised by the fact that the only flow quantity extracted from computational fluid dynamics (CFD) data is the mean shear, whereas the other quantities are defined with analytic expressions. These expressions can be found in Kamruzzaman et al. (2012) and are here summarised. The vertical component of turbulent kinetic energy is

2.4.
A model of the turbulence-turbulence interaction wall-pressure spectrum The derivation of a wavenumber-domain solution of the unsteady pressure generated by the turbulence-turbulence interaction source term, originally proposed by Hodgson (1961), is thoroughly reviewed and amended in some points in appendix C. In this section only the main hypotheses and the final expression are briefly recalled, referring to the appendix C for further considerations. The following assumptions are made: (i) joint normal probability distribution of the turbulence at any two points across the boundary layer; (ii) isotropic distribution of the turbulent kinetic energy in the three spatial coordinates (2.40) (iii) cross-spectral density of velocity fluctuations having a Gaussian distribution in the wavenumber plane, k: where F 2 represents the vertical velocity correlation function: it can be assumed as Gaussian (C 32) or an exponential function (C 33); and (iv) negligible contribution of turbulence cross-spectral densities such that ϕ i2 , i ∈ {1, 3}.
Appendix C then shows that the following expression of the turbulence-turbulence wall-pressure PSD can be derived:  (2007) and Slama et al. (2018) also made the assumption of a joint normal probability distribution of turbulence in order to treat the turbulence-mean shear and turbulence-turbulence interaction terms separately. Consequently, multiplying equation (2.7) by its complex conjugate in order to get the wall-pressure correlation coefficient yields a six-dimensional integral which must still be Fourier transformed in two dimensions in order to get the wavenumber PSD. However, the solution proposed in this section is derived at the price of some restrictive assumptions and should be seen as a means to better understand the relative importance of the two wall-pressure source components.

Heuristic investigation of the influence of the turbulent kinetic energy distribution
The previous section has highlighted the fact that the distribution of turbulent kinetic energy over the wavenumber spectrum can be represented by various analytical expressions. For this reason, the shape of a predicted wall-pressure PSD not only depends on the mean velocity and the r.m.s. vertical fluctuation profiles of the corresponding boundary layer, but also on the prescribed functional representation of the cross-spectrum of vertical fluctuations. The aim of this section is to make some heuristic considerations on the latter effect by using a simplified model of the boundary-layer profiles over a flat plate without a pressure gradient. This model, proposed by Kraichnan (1956b), assumes that the turbulent kinetic energy is constant in the vertical direction and that it is isotropic along the three axes, so that the r.m.s. velocity is Furthermore, the mean streamwise velocity gradient is assumed to decay exponentially with the wall-normal direction Kraichnan found that (3.2) is a fair approximation of the derivative of the logarithmic law of the wall if the parameter γ is of the order of magnitude of the inverse longitudinal integral length scale, Λ. This is admittedly a highly artificial flow model and serves only as a first approximation towards the understanding of a real boundary layer, as the one treated in § 4, and its main advantage is merely its analytical simplicity. However, it will allow us to make some important considerations on the relationship between the model of turbulence statistics and the shape of the corresponding wall-pressure spectrum. For this purpose, all turbulence models discussed in § 2.3 will be used. In particular, the mathematical properties of the Gaussian and Liepmann models allow us to carry out analytically the integrations of (2.17). To begin with the Gaussian spectrum, substituting equations (2.22), (3.1), (3.2) in (2.17) and integrating with respect to X 2 and X 2 yields the following idealised wall-pressure spectral distribution: Similarly, the substitution of (2.25) (with ν = 1/2) and of (3.1) and (3.2) in (2.17) yields 2). Equations (3.3) and (3.4) provide direct, if approximate, expressions of the wall-pressure energy distribution among the wavenumber spectrum which are a direct consequence of the respective assumptions on the energy distribution of the boundary-layer turbulence.
3.1. Comparison with Goody's model The empirical model developed by Goody (2004) represents the basic physical features of a wall-pressure spectrum under a zero-pressure-gradient turbulent boundary layer. At the lower frequencies, which are dominated by larger-scale turbulent structures present in the outer part of the boundary layer, the spectrum rises as ω 2 . The ω −0.7 decay corresponds to an overlap region that collapses with both inner and outer  boundary-layer scaling parameters (see Goody 2004). It has been argued that this region of the spectrum is mostly influenced by the physics of the log layer (see Bradshaw 1967). Finally, the wall-pressure spectrum decays as ω −5 for ω → ∞ (see also Moreau & Roger 2005, p. 49). These trends are represented in the equation represents the ratio of the outer layer to inner layer time scales. It can be noticed that the use of inner and outer boundary-layer parameters for normalisation provides a better overlap of different experimental curves over the whole frequency range (see Keith et al. 1992;Cipolla & Keith 2000;Goody 2004). Equation (3.5) is compared with the distributions produced by the turbulence models of § 2.3, following the assumption that the main contribution to the wall-pressure spectrum is that of the turbulence-mean shear interaction. The results are depicted in figure 2. The Gaussian and Liepmann curves are obtained by means of (3.3) and (3.4), respectively. The wavenumber spectra given by these equations are transformed into frequency spectra by integrating numerically equation (2.18). The von Kármán and RDT curves are obtained by numerical integration of (2.17) and (2.18) in the X 2 , X 2 and k 3 variables. The following parameters, representative of a subsonic turbulent boundary layer at M = 0.2 and Re c = 5.85 × 10 5 , have been used: The results show that all turbulence models considered in this study are able to represent the quadratic rise of the wall-pressure PSD at low frequencies. In particular, the Liepmann and von Kármán curves are quite close to Goody's in that range. Concerning the high-frequency range, only the Liepmann and von Kármán curves show a decay close to ω −5 , whereas the Gaussian curve drops much faster. Also the RDT curve drops slightly faster than the Liepmann and von Kármán ones. The ω −5 slope has been predicted theoretically by Blake (1986), who attributed it to the dominating influence of the viscous sub-layer. This is consistent with the fact that the Gaussian turbulence spectrum represents only large-scale turbulence structures, whereas the three turbulence spectral models classified under the generalised von Kármán spectrum take into account also the energy contained in the finer-scale structures such as those present in the viscous sub-layer. It is thus shown that the mere knowledge of the wavenumber distribution of turbulent kinetic energy, with only a simplified model of the boundary-layer mean and r.m.s. velocity profiles, is sufficient to justify the low-frequency and high-frequency slopes of the wall-pressure spectrum observed experimentally. On the contrary, the ω −0.7 slope of the overlap region cannot be retrieved by analytical considerations alone, at least not without a realistic model of the boundary-layer velocity profiles.
3.2. Effect on the mean-square wall pressure due to turbulence-mean shear interaction The mean-square wall pressure due to the turbulence-mean shear interaction is defined as the integral of the wall-pressure PSD over the wavenumber spectrum This integration can be performed analytically in order to understand how a given formulation of ϕ 22 influences the predicted mean-square wall pressure, which is in turn directly proportional to the overall noise emitted by the trailing edge of an aerofoil, for instance. However, analytical solutions can be obtained only at the price of a further simplification, that is imposing γ = 0. In this case, the mean shear would be constant across the boundary layer and the mean velocity profile would be linear. This is then representative only of the physics of the laminar sublayer of a turbulent boundary layer. The application of (3.6) to the Gaussian-based wall-pressure spectrum of (3.3) yields Similarly, the integration of the Liepmann-based wall-pressure spectrum of (3.4) yields It is thus verified that the flattening of the longitudinal correlation function of turbulence from a Gaussian to an exponential curve (see figure 1) makes the turbulence-mean shear interaction component of the mean-square wall pressure increase by a factor 2, within the assumptions of the simplified boundary-layer model and assuming that the length scales of the two turbulence models are equal. The same result was derived by Kraichnan (1956a) by integration of the wall-pressure correlation coefficient in the space domain. The study of a real boundary layer presented in § 4 will allow us to understand if this relationship between longitudinal correlation function and mean-square wall pressure holds also in the case of a real boundary layer.
3.3. Application to the turbulence-turbulence interaction component Equation (2.42) can be integrated analytically in order to estimate the order of magnitude of the mean-square wall pressure due to the turbulence-turbulence source term of the Poisson equation. For illustrative purposes, the mean-square velocity u 2 is considered constant in the x 2 direction. The result then depends on the choice of the F 2 function. If the Gaussian function of (C 32) is substituted in (2.42) and the resulting formulation is integrated with respect to X 2 and X 2 , one obtains This result is obtained with a procedure similar to that described in § D.1. Then, the integration of (3.9) with respect to dk yields If, on the contrary, the exponential function of (C 33) is substituted in (2.42), one obtains In this case, the integral in dX 2 dX 2 is solved by using the variable transformation for Lk = −2 (which is always verified because L and k are positive by definition). Then, substituting equation (3.12) in (2.42) leads to which corresponds to Hodgson's equation (B.16) except that the present coefficient 28 appears as 22 in the manuscript, which is most likely a typewriting error. Finally, the integration of (3.13) with respect to dk yields (3.14) These considerations based on a simplified boundary-layer model suggest that the shape of the vertical correlation function does not alter significantly the mean-square wall pressure for constant u 2 . Again, the application of the turbulence-turbulence interaction theory to the DNS data presented in § 4.6 will allow us to assess the validity of this result in case of a real boundary layer. It can also be noticed that the turbulence-turbulence mean-square wall pressure does not depend on the integral length scale, unlike the analogous expressions derived for the turbulence-mean shear component in § 3.2.
4. Comparison of the models with direct numerical simulation and measurements of the flow over a controlled-diffusion aerofoil 4.1. Wall-pressure measurements Measurements of the surface-pressure fluctuations on a controlled-diffusion (CD) aerofoil were performed in the anechoic open-jet wind tunnel at Université de Sherbrooke (see Padois et al. (2015) for an assessment of the capabilities of this facility for aeroacoustic investigations). The CD aerofoil is instrumented with remote microphone probes (RMP) on the suction side (see . The RMPs are positioned along the streamwise and spanwise directions, as can be seen in figure 3. Each probe comprises a Knowles FG 23329 P07 microphone remotely connected via a T junction to a pinhole of diameter of 0.5 mm on the suction side of the aerofoil and to a 2 m long anechoic termination on the other end to damp end reflections and avoid resonances (see for more details Perennes & Roger 1998). This microphone has a flat response in the frequency range of 100 to 10 kHz. For the calibration of the microphones a TMS microphone (130P10 ICP) with a diaphragm diameter equal to 6.35 mm was used. It is known to have a flat response in the frequency range of 100 to 20 kHz. The TMS microphone was calibrated using a Larson Davis 200 calibrator which has an uncertainty of approximately ±0.2 dB. The data were acquired using a NI 9234 module with a 24 bit resolution. The acquisition frequency was 25 600 Hz and the length of the signals was 30 s. Since the actual probes were not flush mounted but instead were placed at a remote location, this leads to attenuation of the surface pressure registered at the remote location. This attenuation can be attributed both to a change in section of the connecting tubes and the finite length of the tubing. Hence the actual signal measured by the remote microphone has an attenuated amplitude and is phase shifted. To account for these effects an in situ calibration procedure was realised.
Lastly, it should be kept in mind that, due to the finite size of the pinhole, some spatial averaging is anticipated especially at higher frequencies. Gravante et al. (1998) performed a set of wall-pressure measurements with pinhole microphones of several diameters and showed that spectral attenuation can be avoided until f + = 1 if 12.0 < d + < 18.0 , where f + = f ν T /u 2 τ and d + = d u τ /ν T are the frequency of interest and diameter of the pinhole respectively, both normalised by wall units. It was also shown that the diameter of the pinhole has no influence on the mean-square wall pressure up to d + = 27. Using the friction velocity calculated from the DNS, the d + for sensor 21 was estimated to be equal to 26.0, whereas d + for sensor 24 was estimated to be equal to 22.2. Likewise, f + = 1 corresponds to f = 39.1 kHz for sensor 21 and to f = 28.6 kHz for sensor 24. However, since the measurement range for the Knowles FG 23329 P07 is restricted to 10 kHz the attenuation effects are expected to be small. However, Roger (2017) has shown with analytical considerations that the attenuation due to integration effects at high frequencies depends on the statistical properties of the wall-pressure field to be measured. Then, the DNS-based wall-pressure statistics could be used in the future to understand better this effect.

DNS of the flow over the controlled-diffusion aerofoil
The three-dimensional DNS of the flow over the CD aerofoil is conducted using a multi-block structured high-order accurate compressible Navier-Stokes solver (Sandberg 2015) taking the mean installation effects ( jet width) into account, as suggested by Moreau et al. (2003). The volume data around the aerofoil were recorded during the simulation at a sampling frequency of 78 kHz for 7 flow-through times based on the reference velocity and aerofoil chord length, which was found enough to yield properly converged flow statistics in previous large-eddy simulations (LES) and DNS on this aerofoil at this flow condition (see Wang et al. 2009;Moreau et al. 2011;Wu, Moreau & Sandberg 2019). A fourth-order central standard-difference scheme with Carpenter boundary stencils (Carpenter, Nordström & Gottlieb 1999) is applied for the spatial discretisation in the streamwise and cross-wise directions. A spectral method using the FFTW3 library is used in the spanwise direction. Time marching is achieved by an ultra-low-storage five-step fourth-order Runge-Kutta scheme (Kennedy, Carpenter & Lewis 1999) with a constant time step of t = 6.4 × 10 −8 s. Characteristic-based boundary conditions (described in Sandberg & Sandham 2006;Jones 2008;Jones, Sandberg & Sandham 2008) (see also Kim & Lee 2003) are used for this simulation on domain boundaries to avoid unphysical reflections (see for more details Wu et al. 2018). An adiabatic, no-slip condition is set on the aerofoil surface. The grid refinement in the boundary layers is guided by the data obtained from experiments (Moreau et al. 2003;Wu et al. 2016). In the spanwise direction, 96 Fourier modes are employed with 100 % de-aliasing, which corresponds to 194 collocation points in physical space. These points are distributed with equal distance over a spanwise width of 12 % of the aerofoil chord length and are shown to resolve properly the Kolmogorov scale (see Wu et al. 2018Wu et al. , 2019. The O-grid around the boundary layer comprises 3341 × 279 × 194 grid points, which give finally a grid resolution of x + 1,ave ≈ 6.0, x + 2,ave ≈ 0.8 and x + 3,ave ≈ 5.5 based on the friction velocity. These average values are calculated on the rear part of the aerofoil surface corresponding to 40 % of the aerofoil chord. Further details on the validation of the DNS can be found in Wu et al. (2019) including turbulent kinetic energy budgets.
Mean velocity profiles extracted from the DNS data on the aerofoil surface and in the wake compare favourably with hot-wire measurements performed in the anechoic wind tunnel of the Université de Sherbrooke (UdeS) (see Wu et al. 2018Wu et al. , 2019. The calculated pressure coefficient, where p 0 is the free-stream static pressure, is compared in figure 4(a) with wall-pressure measurements taken in the same wind tunnel, which has a free-stream turbulence intensity below 0.4 %. A laminar separation bubble can be observed at the leading edge, whose extension is well captured by the DNS. The onset of an adverse pressure gradient is visible on the rear half of the suction side. The development of the pressure gradient on the suction side of the aerofoil is presented in figure 4(b) using the Clauser parameter, defined as, The pressure gradient with respect to the abscissa along the aerofoil, dP/dx 1 , increases continuously from the mid-chord to the trailing edge.
The wall-pressure PSD calculated from the DNS results is compared in figure 5 with the measured one at three pressure sensors on the suction side of the aerofoil close to the trailing edge, namely sensor 21 (with β C = 5.01), sensor 24 (with β C = 6.28) and sensor 26 (with β C = 9.72). The boundary-layer parameters used for the normalisation of the curves are summarised in table 2. The boundary-layer thickness corresponds to the x 2 coordinate where the total pressure reaches 99 % of its maximum constant value along the wall-normal direction. The total pressure is used instead of the velocity because the aerofoil is loaded and therefore the external flow is deflected so that the external velocity is not constant. It can be noticed that at the lower frequencies the DNS spectrum is slightly below the experimental one. This is due to the fact that, whereas the simulation does not include the jet unsteadiness and acoustics, some significant effect from the jet is inevitable in the experiment, which prevents us obtaining the theoretical positive slope in the wall-pressure spectra.
In the case of sensor 26, a high-frequency hump can be noticed, which is attributed to the contamination from the acoustic sources in the wake of the aerofoil. For this reason, the following analysis will focus on sensors 21 and 24 only. Finally, figure 6 presents the boundary-layer profiles of mean streamwise velocity and vertical r.m.s. fluctuation extracted from the DNS corresponding to sensors 21 and 24. Both plots are normalised in wall units. The effect of the increase of β C from sensors 21 to 24 can be observed in the mean velocity plot of figure 6(a), particularly in the wake zone of the boundary layer corresponding to x + 2 > 100. Also the vertical r.m.s. fluctuation of figure 6(b) is seen to increase with the Clauser parameter as expected.

Calculation of turbulence statistics
The DNS of the flow around the CD aerofoil is post-processed to gain information on the statistics of turbulence that are part of the above analytical model for the prediction of the wall-pressure spectrum. The wall-normal velocity correlation coefficient, R 22 , is analysed. In fact, the cross-spectral density of wall-normal velocity fluctuations, ϕ 22 , which appears in (2.17), is the double spatial Fourier transform of R 22 (see (B 9)). The correlation coefficients calculated from the DNS results are compared with the analytical formulations that are derived from a given longitudinal correlation function, F(r), according to (B 6). This comparison assesses the accuracy of the analytical formulation of ϕ I 22 . The R 22 profiles discussed below are extracted from planes parallel to the suction side of the CD aerofoil using seven flow-through times for averaging. The results correspond to sensor 24. The same calculations performed for sensor 21 are not reported here for brevity as they yield equivalent results. The profiles are plotted versus the distance from the reference point, r, normalised by the longitudinal integral length scale, Λ, calculated at each wall-normal position. Furthermore, R 22 is calculated for several streamwise, r 1 , and transverse r 3 separations. The curves of R 22 (r 1 ) and R 22 (r 3 ) nearly collapse if the coordinate r 1 is divided by the coefficient α that represents the ratio of streamwise to transverse longitudinal integral length scales (see (2.19)). Figure 7(a) shows the R 22 (r 1 ) and R 22 (r 3 ) profiles extracted from a plane at x + 2 = 23. In this case, it is sufficient to take α = 1.88 for the curves to coincide, at least for small separations, r. The transverse correlation R 22 (r 3 ) is higher than R 22 (r 1 ) for longer separations, whereas R 22 (r 1 ) exhibits negative lobes. These facts point to the presence of coherent turbulent structures elongated in the transverse direction (see Sillero, Jiménez & Moser 2014, figure 1). Wu et al. (2018) have shown that the onset of these structures is related to the increase of the adverse pressure gradient and that also the wall-pressure correlation coefficient increases in the transverse direction. In this case, all analytical formulations seem to overestimate the numerical correlation coefficients. Only the Liepmann formulation is in close agreement with the numerical curves for small values of r. The von Kármán formulation is not shown in this analysis as its results are very close to that of Liepmann in all calculations. Figure 8 presents the same analysis of the velocity fluctuations extracted from a plane at x + 2 = 65. The elongation of the turbulence structures farther away from the wall is represented by the coefficient α = 1.55. In this case, the Gaussian formulation overestimates the negative undershoots of the numerical curves. On the contrary both the RDT and the Liepmann formulations give good approximations to the results. In particular, the Liepmann curve is closer to the computed R 22 for r < 1, whereas for r > 1 the RDT curve is the closest. Figure 9 refers to a plane close to the outer edge of the boundary layer (x + 2 = 165), where the elongation of the turbulent structures is α = 1.35. It can be noticed that the Gaussian theory (figure 9b) reproduces well the negative parts of the curve of R 22 . Figures 9(c) and 9(d) compare the DNS curves with the RDT and Liepmann formulations, respectively. In this case, the former analytical formulation provides a slightly better representation of the numerical results than the latter.
Vertical velocity correlation coefficient calculated on a plane at x + 2 = 23 in the streamwise (r 1 ) and transverse (r 3 ) directions. Ratio of streamwise to transverse length scales: α = 1.88.
In the results presented so far it has been assumed that x 2 = 0. In order to complete this analysis, the calculation of R 22 between two planes parallel to the wall is now presented. Figure 10 shows the correlation coefficient computed between the plane at x + 2 = 65 and a second plane parallel to the wall which has a normalised distance x 2 /δ = 0.024 from the first. Due to the wall-normal separation, the maximum value of R 22 , found at r 1 = r 3 = 0, is logically lower than unity. Visibly, the three analytical formulations predict differently the decay of R 22 with x 2 . In this case, the RDT and Liepmann formulations are in better agreement with the numerical data than the Gaussian. Finally, figure 11 presents the R 22 coefficient calculated between the plane at x + 2 = 65 and another plane at x 2 /δ = 0.05 from the first. It can be noticed that the peak of R 22 (r 1 ) corresponds to a r 1 slightly below zero, whereas the plot of R 22 (r 3 ) is centred on r 3 = 0. This is due to the fact that the velocity time series has been extracted directly at the mesh cell centres, in order to avoid the numerical error introduced by interpolation. The mesh follows the curvature of the aerofoil surface, so the cell centres are not perfectly aligned in the wall-normal direction, whereas there is perfect alignment in the transverse direction. It must be emphasised that if the fundamental assumption of the TNO-Blake model, expressed in (2.31), was true, then the DNS-based statistics presented in figures 10 and 11 would vanish for any r.
The decay of the maximum R 22 with the increase of the vertical separation is plotted in figure 12. The reference point for this calculation has been taken at  FIGURE 8. Vertical velocity correlation coefficient calculated on a plane at x + 2 = 65 in the streamwise (r 1 ) and transverse (r 3 ) directions. Ratio of streamwise to transverse length scales: α = 1.55. roughly half the boundary layer thickness for both sensors 21 and 24. The second point for the calculation of the cross-correlation coefficient has been chosen moving upwards in the vertical direction. As will be shown below, the longitudinal integral length scale of turbulence is nearly constant in this region of the boundary layer. For this reason, it has been possible to normalise the abscissa of figure 12 using the averaged Λ in the outer part of the boundary layer. It can be observed that, especially for small separations, the RDT is in best agreement with the DNS data corresponding to both sensors. For higher vertical separations, the DNS data are between the RDT and Liepmann theoretical curves. The cross-correlation coefficient of vertical velocity fluctuations is significantly higher than zero on a vertical extension corresponding to two integral length scales, thus further contradicting equation (2.31).
The general conclusion that can be drawn from this study is that especially the RDT and, second, the Liepmann models provide an acceptable description of the turbulence statistics in this adverse-pressure-gradient boundary layer. The present result is also very consistent with what Magnaudet (2003) found: the rapid distortion theory is the leading-order approximation capable of describing short-and long-time evolutions of turbulent boundary layers in the limit of large Reynolds number. In any case, all analytical formulations of the cross-spectral density of wall-normal turbulent fluctuations, ϕ 22 , that have been presented in § 2.3 will be applied to the prediction of the wall-pressure spectra on the surface of the CD aerofoil in § 4.5. In this way, the

R 22 (r 1 /å) R 22 (r 3 ) Gaussian
. Vertical velocity correlation coefficient calculated on a plane at x + 2 = 165 in the streamwise (r 1 ) and transverse (r 3 ) directions. Ratio of streamwise to transverse length scales: α = 1.35. relationship between the shape of turbulence statistics and that of the corresponding wall-pressure PSD will be assessed.
In order to compute ϕ 22 with any analytical formulation, it is necessary to know the variation of the length scale Λ with the distance from the wall. This is depicted in figure 13(a) corresponding to sensors 21 and 24, for which the wall-pressure spectrum will be calculated. Since Λ is a quantity that is defined in the theory of homogeneous turbulence (see, for instance Wilson (1997Wilson ( , 1998 and appendix B), it should only be calculated with respect to spatial directions in which the turbulence can be, to a good approximation, considered homogeneous. In this case, the correlation of the streamwise velocity component, R 11 , has been calculated as a function of the streamwise separation, r 1 , on planes at given distances from the wall. According to (B 3), Λ is the integral of this correlation function. Since the tail of the function tends to be noisy, similarly to the R 22 functions previously examined, R 11 has been calculated over several parallel streamwise lines and the results have been averaged in order to obtain a smooth curve allowing a robust calculation of the integral. Interestingly, the plateau of Λ/δ corresponds to a higher value than the classical 0.14 of Prandtl theory (see Panton & Linebarger 1974). Figure 13(b) shows the trend of α along the wall-normal direction: the length scales are equal (α = 1) for x 2 /δ > 0.8 corresponding to both pressure sensors. Then, α slowly rises to 2.0 corresponding to x 2 /δ = 0.2 and reaches the value 3.0 at the inner edge of the boundary layer, but only for sensor 24. Instead, the maximum α equals 2.5 corresponding to sensor 21. Yet,

Estimation of the convective speed of wall-pressure fluctuations
The wavenumber wall-pressure spectra predicted with (2.17) will be converted into the frequency domain by using (2.18) for the sake of comparison with the measured and directly computed spectra. For this purpose, it is necessary to define the value of the convective speed of the wall-pressure fluctuations, U c . This quantity is, in principle, a function of the frequency. The empirical model of U c developed by Smol'yakov (2006) for zero-pressure-gradient boundary layers, however, predicts that this quantity will tend to a constant value with the increase of ω. Salze et al. (2014) have shown that this trend is also valid for adverse and favourable pressure gradients. The reference measured and computed spectra correspond, in the present application case, to a range of frequencies where, according to the results presented by these authors, U c can be considered to a good approximation a constant fraction of the external flow speed, U 0 . The estimation of U c is made by means of the cross-correlation of two pressure signals extracted from the DNS over seven flow-through times at a sampling frequency of 78 663 Hz. This has been done for both sensors 21 and 24. The two points on the surface of the aerofoil are aligned in the transverse direction and are separated by a given distance in the streamwise direction. The points around sensor FIGURE 11. Vertical velocity correlation coefficient calculated between a reference plane at x + 2 = 65 and a second plane parallel to the wall at x 2 /δ = 0.05 from the first. Ratio of streamwise to transverse length scales: α = 1.55. 21 are separated by 0.85δ, whereas those around sensor 24 are separated by δ. Since the largest size of the vortical structures corresponds roughly to the boundary-layer thickness, pressure signals recorded at locations more than δ apart are expected to be statistically uncorrelated. Figure 14 depicts the wall-pressure cross-correlations versus the time lag, τ . The maximum cross-correlation is indicated with a circle in each plot. The distance between the two points where the pressure signals were recorded divided by the time lag corresponding to the maximum cross-correlation gives an estimate of the convection speed. It is found that for sensor 21, U e /U c = 1.51, whereas for sensor 24, U e /U c = 1.54. These values are consistent with those found in the previous experimental investigations of Roger & Moreau (2004) and  on the CD aerofoil.

Prediction of the wall-pressure PSD
The wall-pressure PSD is predicted by using the mean shear, ∂U 1 /∂x 2 , and r.m.s. of the wall-normal velocity fluctuations, u 2 , as input to (2.17) and (2.18). The mean and r.m.s. velocity profiles have been presented in figure 6. The effect of the anisotropy of the length scales on the predicted wall-pressure PSD can be verified, as a preliminary test, by choosing one wall cross-spectral density formulation and performing the calculation of ϕ pp twice. The first calculation is made taking into account in (2.19) the variation of α with x 2 (depicted in figure 13). The second calculation is performed with α = 1. The results of this test are shown in figure 15 for sensors 21 and 24 using the RDT turbulence model. It is then shown that the wall-pressure PSD in the frequency range of interest is determined by the isotropic part of the boundary layer. As a matter of fact, the integrand of (2.17) is directly proportional to ∂U 1 /∂x 2 and u 2 , multiplied by an exponential function that decays proportionally to the product of the distance from the wall by the magnitude of the wavenumber vector. Consequently, it is possible that the wall-pressure PSD is mostly determined by the near-wall physics of the boundary layer, as suggested by the study of Panton & Linebarger (1974). The outer part of the boundary layer is expected to play an important role only at very low frequencies, which correspond to small k 1 values. However, since a comparison with measured and DNS wall-pressure PSD can be made only in the medium-low-frequency range, the following computations will be performed assuming α = 1 without any loss of accuracy in the range of frequencies of interest. Since all relevant physical quantities are known from the DNS, the prediction of wall-pressure PSD can be used to assess the effect of the different analytical formulations of ϕ 22 presented in § 2.3. The results for sensors 21 and 24 are presented in figure 16. Comparing the results in figures 16 and 1 clearly highlights the relationship between the longitudinal correlation function of turbulence and the corresponding predicted wall-pressure PSD. The flattening of F(r) from the Gaussian to the von Kármán correlation function, indicated by the arrows in figure 1, corresponds to a marked increase of the wall-pressure PSD in the middle-frequency range and a slight decrease at high frequencies, also indicated with the arrows in figure 16. Overall, the flattening of the turbulence correlation function corresponds to an increase of the mean-square pressure at the wall, a result that has been foreseen by means of purely analytical considerations in § 3.2. It can be noticed that, for ω δ/U e > 5, the RDT and Gaussian ϕ 22 expressions provide the best agreement with the directly computed curve, whereas in the lower-frequency range the Liepmann and RDT-based curves are much closer to the DNS. The good agreement between the DNS and RDT curves is particularly visible in the sensor 24 plot, whereas for the sensor 21 the Liepmann model performs slightly better at low frequencies. As a general conclusion, considering that the difference between the various formulations at high frequencies is relatively small, the RDT and possibly Liepmann turbulence spectrum provide a proper description over the whole frequency range that has been considered (0.15 < ω δ/U e < 25). This is consistent with the considerations made on the turbulence statistics earlier in this section. Furthermore, the use of the turbulence-mean shear interaction source term alone is sufficient to retrieve the correct level of the wall-pressure PSD on the explored range of frequencies. However it is interesting to apply also the turbulence-turbulence interaction formulation in order to estimate its contribution to the overall spectrum.
4.6. Application of the turbulence-turbulence interaction component formulation The formulation presented in § 2.4 can be applied to the prediction of the turbulenceturbulence interaction wall-pressure PSD in the same way as the turbulence-mean shear formulation. The results obtained using both F 2 formulations are presented in figure 17. It can be seen that the curve obtained with the Gaussian vertical correlation function is only slightly higher than the curve obtained with the exponential F 2 formulation, as suggested by the previous analytical considerations. These numerical results suggest that the turbulence-turbulence source is significant only at high frequencies but thoroughly underestimates the main low-frequency contribution to the overall mean-square pressure. The turbulence-turbulence (TT) curves presented in figure 17 exhibit a plateau for ωδ/U e < 5, a trend that can also be observed in figure  17 of Slama et al. (2018). However, this conclusion must be taken with caution since it is obtained by making such a critical hypothesis as that of a normal probability distribution and isotropic turbulence. This should be further evaluated in the DNS results in the future.

Application of the TNO-Blake model
The TNO-Blake model can be applied to the DNS data in order to assess the effect of the substitution of the turbulence cross-spectrum with the single-point auto-spectrum and of the hypothesis on the boundary-layer profiles of u 2 and Λ on the predicted wall-pressure PSD.
First, the empirical equations used in the original TNO-Blake model can be compared with the boundary-layer profiles extracted from the DNS on the CD aerofoil corresponding to the sensors 21 and 24. Figure 18 compares the directly computed u 2 with the results of the application of (2.36) to the mean-shear profiles. It appears from these results that the empirical formulation fails to reproduce the DNS curve. The small local peaks in the modelled u 2 curve are due to the numerical calculation of the velocity gradient of the mean flow interpolated on a line perpendicular to the  blade surface. On the contrary, it is shown in figure 19 that (2.39) provides a slightly underestimated approximation of the computed Λ.
The accuracy of the TNO-Blake model in the prediction of the wall-pressure PSD can be evaluated in two steps. First, equation (2.32) is computed using the profiles of u 2 2 and Λ given by the empirical equations (2.36) and (2.39), respectively. Second, and perhaps more interestingly, equation (2.32) is computed using the profiles of u 2 and Λ extracted from the DNS of the controlled-diffusion aerofoil, as done in § 4.5. Both sets of calculations are performed assuming isotropy of the length scales (α = 1), following the results presented in § 4.5.
In the first test the modelled wall-pressure PSD is underestimated with respect to the DNS result, as can be seen in figure 20, and this is a direct consequence of the inaccurate representation of boundary-layer profiles. The results of the second test are depicted in figure 21. In this case, the wall-pressure PSD obtained with the exact expression of the von Kármán cross-spectral density, already presented in § 4.5, has been reproduced for comparison. It can be seen that the approximate solution of the TNO-Blake model and the one based on the exact formulation of the turbulence crossspectral density yield similar results in terms of wall-pressure PSD at high frequency but diverge at low frequency. In the range ω δ/U e < 5 the wall-pressure PSD predicted using the vertical velocity auto-spectrum underestimates the prediction obtained with the exact formulation of the cross-spectrum but, thanks to this very fact, it also gets closer to the DNS result. Indeed, in § 4.5, the von Kármán theory was shown to overestimate the wall-pressure PSD at low frequencies. In the TNO-Blake model, this inaccuracy is compensated by substituting the cross-spectral density with the singlepoint auto-spectrum, although the latter assumption is quite wrong as shown in § 4.3. Section 2.3.3 has illustrated the treatment of the anisotropy of length scales by means of the coefficients presented in different versions of the TNO-Blake model (see table 1). The final test of this section concerns the effect that these different sets of coefficients have on the predicted wall-pressure PSD. The results, obtained for sensors 21 and 24, are depicted in figure 22, where the DNS spectrum and the TNO-Blake prediction obtained with the hypothesis of isotropy of the length scales are reported for comparison. The earliest set of coefficients proposed by Bertagnolio et al. (2014) significantly overestimates the DNS result at low frequencies. However, the later sets of coefficients proposed by Stalnov et al. (2016) and Fischer et al. (2017) are more conservative in the sense that they do not vary significantly the prediction with respect to the isotropic case.

Application of a Monte Carlo integration method to high-dimensional wallpressure spectrum formulations
The analytical formulation of the turbulence-mean shear wall-pressure spectrum component developed in § 2.2 is a three-dimensional integral, which does not pose problems of numerical computation. However, in some cases it might be desirable to . In this case, as mentioned earlier, the space-domain wall-pressure correlation coefficient is obtained by computing a six-dimensional integral, followed by a two-dimensional integration to obtain the wavenumber-domain PSD. It is not advised to perform a high-dimensional numerical integration by quadrature. In fact, taking N equally spaced samples in each of the d integration dimensions, the integrand function must be calculated at N d points in a d-dimensional hyperspace. Therefore, even for a moderate number of samples N in each dimension, the number of function evaluations rises exponentially with the number of dimensions, a problem which is known to mathematicians as the curse of dimensionality. Another issue concerns the rate of convergence of the integral in multiple dimensions. For instance, if a trapezoidal rule is used the uncertainty of the integral estimation is proportional to N −2/d for large N and similar trends can be demonstrated for other quadrature rules (see James 1980). Additional difficulties arise, in this particular integration problem, from the singularity points in the integration domain where x = x . The solution proposed by Peltier & Hambric (2007) was Gauss-Legendre integration, while Slama et al. (2018) used a Kriging-based algorithm. In order to provide a further alternative, this paper investigates the application of one of the Monte Carlo integration methods, for which the number of samples of the integrand function, N, and the rate of convergence of the estimation of the integral are independent of the number of dimensions. Different Monte Carlo algorithms applied to the problem of wall-pressure spectrum prediction have been compared and evaluated in Grasso, Jaiswal & Moreau (2018), namely quasi-random sampling, importance sampling and recursive stratified sampling. It has been shown that the importance sampling method provides the fastest convergence of the integral with respect to the number of samples (see Grasso et al. 2018, figure 3). This method has been used by Panton & Linebarger (1974) and Remmler et al. (2010) for the computation of wall-pressure PSD. Its only disadvantage is that it is tailored to the specific integrand function, so that an algorithm must be devised for every integration problem. As an alternative, this section focuses on the recursive stratified sampling which does not impose conditions on the integrand function, not even continuity, while still requiring an acceptable computational cost, as will be shown below.
5.1. Recursive stratified sampling This method has been developed by Press & Farrar (1990) and is implemented in the widely used MISER algorithm. In order to understand its advantages, it is necessary to take a step back and introduce the simplest form of Monte Carlo integration: the quasi-random sampling. Given a generic function f : R n → R, sampled with N random vectors, x i , uniformly distributed over the integration domain, the function mean estimator is defined as Then, the estimator of the integral is the product of f by the size of the integration volume. It can be shown (see James 1980) that the variance of the estimator is related to the variance of the function over the integration domain as and that the following asymptotic relation holds for the error of the estimation Thus, the convergence of the quasi-random sampling is independent of the number of dimensions of integration. The recursive stratified sampling differs from quasi-random sampling because it divides the integration region V into two equal and disjoint subvolumes, V a and V b . The sub-volumes are then sampled at N a and N b = N − N a random points, respectively, drawn from a uniform probability distribution function. In this case, the function estimator is defined as where the f a,b terms are computed by applying equation (5.1) to each sub-volume. It can also be shown that the variance of f is where σ = Var( f ). If N a satisfies equation (5.6), the variance of the function mean estimator is The advantage of this method is that Var( f ) as defined in (5.7) is never greater than the variance of the estimator computed with the quasi-random sampling, defined in (5.2). In practice, a fixed number of function evaluations N is allocated at the beginning of the execution of the algorithm. Then each dimension of integration is divided in two disjoint sub-volumes. A given fraction of the N samples is used to estimate the variance of the function on each sub-volume, so that the remaining points can be allocated according to (5.6). In order to converge to a more accurate estimation, the algorithm can be applied recursively on the sub-volumes where the integrand function has greater variance.

5.2.
Application to the turbulence-mean shear component of the wall-pressure spectrum As mentioned at the beginning of this section, the expression of the turbulence-mean shear component of the wall-pressure PSD resulting from the substitution of (2.17) in (2.18) is a three-dimensional (3-D) integral which can be calculated numerically with quadratures in a reasonable computational time. However, the same turbulence-mean shear component can be expressed by a four-dimensional (4-D) integral if the one-dimensional integral expression of ϕ 22 given in (B 13) is substituted in (2.17). This expression can also be made five-dimensional (5-D) if the two-dimensional (2-D) integral expression of ϕ 22 given in (B 10) is used instead. Thus, this problem in which the same physical quantity can be expressed by a 3-D, 4-D or 5-D integral formulation is particularly suited for the evaluation of the performance of a Monte Carlo method. The recursive stratified sampling is applied to the calculation of the turbulence-mean shear wall-pressure PSD, with the 4-D and 5-D formulations, corresponding to sensor 24 of the CD aerofoil using the boundary-layer profiles presented in § 4. The von Kármán model of turbulence statistics has been selected for this study. Analogous tests have been presented in Grasso et al. (2018) using the Gaussian turbulence model, which has a simpler mathematical formulation. The MISER algorithm used to obtain the following results is available in the Scikit-Monaco package for Python (see Buignon 2013). Figure 23 presents the convergence of the 4-D and 5-D formulations of the wall-pressure PSD calculated with the MISER algorithm for an increasing number of sample points. The spectra have been computed with N ranging from 2000 to 2 × 10 5 points. It can be seen that the integration is challenging from the numerical point of view only at the highest and lowest frequencies. However, starting from 2 × 10 4 sampling points, all curves coincide. It is thus shown that the convergence of the MISER algorithm depends on the number of sampling points but not on the number of dimensions of integration, as postulated by the theory. The calculation of the wall-pressure PSD at 25 values of ω, using 2 × 10 4 sampling points for each frequency, has been performed in parallel on 4 cores of an Intel Xeon Processor EG-1630 v3 with base frequency of 3.70 GHz. The computational time is 8.6 s for the 3-D formulation and 10.2 s for both the 4-D and the 5-D integral formulations. Increasing the number of sampling points for each frequency of one order of magnitude, that is 2 × 10 5 , the computational time corresponding to 25 values of ω integrated in parallel on 4 cores rises to 1 min for the 3-D formulation and to 1 min 13 s for both the 4-D and the 5-D integral formulations. Finally, figure 24 presents the somewhat obvious but necessary verification that the integration algorithm converges to the same prediction of the wall-pressure PSD using the 3-D, 4-D and 5-D formulations, which are equivalent from the analytical point of view. This application shows that the MISER algorithm can be recommended for the prediction of wall-pressure PSD due to its fast and robust convergence, the ease of implementation and the absence of conditions on the shape of the integrand function.

Conclusions and perspectives
This work has presented a comprehensive approach to several analytical formulations of the spectrum of wall-pressure fluctuations under a turbulent boundary layer available in the literature. The discriminant character of published models usually consists in the choice of an analytical expression of the turbulence cross-spectral density. For this reason, the most widely used turbulence spectral models have been reviewed and presented in a unified framework. The two sources of wall-pressure fluctuations, turbulence-mean shear and turbulence-turbulence interaction, have been modelled in the wavenumber domain. In particular, the influence of the different analytical formulations of turbulence statistics on the shape and level of the resulting wall-pressure PSD has been studied for the first time. The relationship between the shape of the longitudinal correlation function of turbulence and that of the corresponding predicted wall-pressure spectrum has been elucidated by applying different turbulence models to the same boundary-layer mean profiles extracted from the DNS of the flow around a controlled-diffusion aerofoil under significant adverse pressure gradients (Clauser parameter of 5.01 and 6.23 in the two test cases, respectively). The post-processing of DNS data allowed us to compute directly two quantities that are usually modelled in the case of RANS-based predictions, namely the longitudinal integral length scale of turbulence, Λ, and the anisotropy factor, α. It has been found that only the inmost part of the boundary layer, in this application, is affected by a strong anisotropy of the turbulence length scales. For this reason, a good agreement between the predicted and directly computed wall-pressure PSDs in the range of frequency of interest could be obtained by assuming α = 1. The inclusion of the variation of α with the wall-normal coordinate in the integral formulation of the wall-pressure PSD did not yield appreciably different results. Also, the turbulence-turbulence interaction component resulted to be negligible with respect to the turbulence-mean shear, except at the highest frequencies.
The turbulence spectral model that best fits the turbulence statistics extracted from the DNS is the rapid distortion theory, followed by Liepmann's model. As a consequence, the predictions based on these turbulence models are the closest to the directly computed wall-pressure PSD. The capabilities of the RDT to represent the evolution of turbulence in various configurations have also been highlighted by de la Riva et al. (2004), Santana et al. (2016) and Goldstein et al. (2017). The predictions based on the von Kármán and on the Gaussian model are found to overestimate and underestimate the PSD for ωδ/U e < 10, respectively. It has also been verified that the flattening of the longitudinal correlation function corresponds to a higher low-frequency wall-pressure PSD and, consequently, to a higher mean-squared wall-pressure fluctuation. Interestingly, this result is also suggested by purely analytical considerations on a simplified model of boundary layer originally proposed by Kraichnan (1956b).
The TNO-Blake model, which is an approximate form of the wall-pressure PSD, has also been reviewed. Its fundamental assumption of the statistical de-correlation of vertical velocity fluctuations for non-vanishing vertical separations has been proven wrong by the analysis of the turbulence statistics extracted from the DNS. Yet, provided that realistic flow parameters extracted from the DNS are used, the TNO-Blake model provides reasonable predictions at both selected adverse-pressure-gradient conditions. Among the various versions of the model, the most recent ones proposed by Stalnov et al. (2016) and Fischer et al. (2017) also yield the most accurate predictions.
Finally, a contribution to the problem of the numerical computation of wall-pressure PSD with high-dimensional integral formulations has been presented. The Monte Carlo integration with the recursive stratified sampling technique is robust and its implementation is straightforward thanks to open-source libraries.
The results presented in this paper can provide a useful reference for the improvement of existing RANS-based wall-pressure PSD prediction methods that share the same analytical foundation. The next development that can be envisaged is the extension of the theory hereby presented to the case of compressible subsonic flows. The main problem will be that of the quantification of the effect of the variation of density across the boundary layer on the wall-pressure PSD.
where the physical source terms are represented byQ(k, x 2 ). Equation (A 1) can be solved with the method of variation of parameters (see Ince 1956), whereby the complete solution to (A 1), Y, is split into two terms The term Y c is the solution of the complementary homogeneous ODE and Y p is the specific solution to the inhomogeneous equation. Substituting the ansatẑ p(k, x 2 ) = e λ x 2 in (A 3) we obtain In order to find the specific solution, Y p , let y 1 = e −kx 2 and y 2 = e +kx 2 . Then, the functions y 1 and y 2 are used to calculate the Wronskian, W(y 1 , y 2 ) = e −kx 2 e +kx 2 d d x 2 e −kx 2 d d x 2 e +kx 2 = 2 k.

(A 5)
The specific solution is Substituting equation (A 5) in (A 6) we obtain The complete solution to (A 1) iŝ where the coefficients A and B can be determined from the boundary conditions. In particular, for the solution to vanish at infinity according to (2.2), the coefficient B must be equal to 0. The coefficient A must satisfy equation (2.3), for which the x 2 derivative of the solution vanishes at the wall. Taking the derivative of (A 8), we obtain for x 2 = 0 Finally, the wavenumber-domain solution for the wall pressure iŝ The same solution is obtained by Gerolymos et al. (2013, appendix B) following a Green's function approach.

Appendix B. Derivation of turbulence cross-spectral densities
This appendix presents the concepts from the theory of homogeneous and isotropic turbulence that are necessary to understand the derivation of the turbulence crossspectral densities used in § 2.3. For a complete discussion of this topic, the reader can refer to handbooks by Batchelor (1953), Pope (2000) or Bailly & Comte-Bellot (2015). The analytical expressions derived in this appendix are summarised in the works of Wilson (1997Wilson ( , 1998.

B.1. General definitions
The normalised correlation coefficient of the i, jth components of a homogeneous vector field is defined as With the further assumption of isotropy of the vector field, the longitudinal correlation function is defined as a function of the normalised correlation coefficients as F(r) = R 11 (r, 0, 0) = R 22 (0, r, 0) = R 33 (0, 0, r).

(B 2)
From this follows the definition of the longitudinal integral length scale of turbulence as Furthermore, the lateral correlation function is defined as G(r) = R 11 (0, r, 0) = R 22 (r, 0, 0) = · · · . (B 4) Batchelor (1953) provides the following relationship between the longitudinal and the lateral correlation functions as well as the analytical formulation of the normalised correlation coefficient as In order to take into account the inhomogeneity of turbulence in the vertical direction for a boundary-layer flow, the coordinates of the two points used to compute the crosscorrelation must be written explicitly, and the variance must be written as coherently with the approach of Panton & Linebarger (1974).

B.1.1. Direct Fourier transform
One way to define the cross-spectrum of vertical velocity fluctuations is the double spatial Fourier transform of the normalised correlation coefficient in the plane defined by the wall, R 22 (r 1 , r 2 , r 3 ) cos(k 1 r 1 ) cos(k 3 r 3 ) dr 1 dr 3 .

B.1.2. Inverse Fourier transform
An alternate way to derive ϕ 22 is by performing a one-dimensional inverse Fourier transform, which may be easier to solve depending on the specific form of F(r). In first place, the cosine Fourier transform of F(r) is defined as the one-dimensional spectral function,F Then, it is shown by Batchelor (1953) that the energy spectrum of turbulence is related to the spectral function by the following equation: In turn, the three-dimensional energy spectral density is derived from the energy spectrum E(κ) according to the relation Finally, the inverse spatial Fourier transform in the vertical direction of (B 16) yields B.2. Gaussian spectrum The Gaussian energy spectrum has a rather simple mathematical formulation and therefore it can be used to illustrate the derivation of the cross-spectral density, ϕ 22 , from the longitudinal correlation function, F(r), using either the direct or the inverse Fourier transform method.

B.2.1. Direct Fourier transform
Substituting equation (2.20) in (B 6), the vertical velocity correlation coefficient is obtained as Then, substituting equation (B 18) in (B 13) yields Since the integral of (B 19) is formulated as a subtraction, it can be split into two terms, which are both re-worked using the technique of integration by parts. Using the tables of integrals involving Bessel functions of Rosenheinrich (2016), both terms can be simplified to yield

B.2.2. Inverse Fourier transform
The application of (B 14) to (2.20) yieldŝ Then, the energy spectrum follows from (B 15) The three-dimensional spectral density follows from (B 16) as and, finally, the inverse Fourier transform of (B 23) yields which corresponds to (B 20).

Appendix C. Derivation of the turbulence-turbulence interaction wall-pressure spectrum
This appendix reviews the derivation of the model proposed by Hodgson (1961) for the wavenumber wall-pressure spectrum generated by turbulence-turbulence interaction, clarifying some steps and correcting some mistakes found in the original document. Hodgson's main contribution with respect to the earlier studies of Batchelor (1951) and Kraichnan (1956a) was the introduction of the inhomogeneity of the flow in the wall-normal direction, which makes the task more cumbersome from a mathematical point of view.

C.1. Differential equation formulation
The turbulence-mean shear and turbulence-turbulence source terms of (2.1) are assumed to be independent. Therefore, the governing equation of the pressure fluctuations is reduced to The flow obeys the continuity equation, from which we derive that Using (C 3), equation (C 1) is re-formulated as In (C 1) the wall-normal direction, indicated by the subscript 2, has been singled out. The subscripts i, j ∈ {1, 3}, where the values 1 and 3 represent the streamwise and transverse directions, respectively. Since the boundary-layer flow is assumed to be homogeneous in planes parallel to the wall, equation (C 1) can be Fourier transformed in the two corresponding directions. The dependence on the wall-normal direction, on the contrary, must be kept explicit in the following derivations.
Assuming two functions f (x) and g (x), where x ∈ R m , the Fourier transform of their product is the convolution of their transforms, that is (C 5) Then, the source terms of equation (C 4) are transformed in the following way and finally The same relations hold for the mean terms of (C 4) as well. Then, the following modified Helmholtz equation is obtained: − mean terms dk , (C 9) which corresponds to (B.3) of Hodgson (1961). The mean terms will not appear explicitly in the following derivation, but it is understood that they are subject to the same transformations as the fluctuating quantities. They will be written again at the point of the derivation where they cancel out. Equation (C 9) can be expressed in a This corresponds to (B.5) of Hodgson (1961), but in the original reference one parenthesis is misplaced.
C.2. Green's function solution Assuming the boundary conditions p = 0 for x 2 = 0, dp/dX 2 → 0 for x 2 → ∞, u 2 = 0 for x 2 = 0, the general solution of the (C 15) is obtained with the following Green's function: So, the pressure at the wall can be written in the wavenumber domain aŝ The derivative by X 2 appearing in the second term of the right-hand side of (C 17) can be eliminated by performing the integration in dX 2 , The term in square brackets of (C 18) vanishes for X 2 → ∞ due to the exponential function and vanishes also for X 2 = 0 due to the boundary conditions. For this reason, equation (C 17) can be re-written aŝ which corresponds to (B.6) of Hodgson (1961).

(C 24)
Consistently with the previous definitions, some second-order statistics vanish due to the statistical orthogonality of wavenumber vectors, e.g.
Substituting equations (C 22), (C 23) and (C 24) in (C 21) and integrating with respect to dk and dk , ϕ TT pp (k) = ρ 2 0 ∞ 0 ∞ 0 e −k(X 2 +X 2 ) dX 2 dX 2 × R 2 1 k 2 (2k i k j − k i k j − k i k j )(2k l k m − k l k m − k l k m )ϕ il (X 2 , X 2 , k ) ϕ jm (X 2 , X 2 , k − k ) As a consequence of the use of (C 22), the mean terms that appeared in (C 21) have been cancelled out. Furthermore, it is possible to simplify the coefficient of the first term of the integral in the following way: (2k i k j − k i k j − k i k j )(2k l k m − k l k m − k l k m )ϕ il ϕ jm = 4k i k j k l k m − 2k i k j k l k m − 2k i k j k l k m −2k i k j k l k m + k i k j k l k m + k i k j k l k m −2k i k j k l k m + k i k j k l k m + k i k j k l k m ϕ il ϕ jm = (4k i k j k l k m − 8k i k j k l k m + 4k i k j k l k m )ϕ il ϕ jm .

(C 27)
In fact, k i k j k l k m , k i k j k l k m , k i k j k l k m , k i k j k l k m all represent the same kind of combinations of the elements of k and k , namely those that take one coefficient from k and the remaining three coefficients from k. Likewise, terms k i k j k l k m , k i k j k l k m , k i k j k l k m , k i k j k l k m represent the same combinations where two factors are taken from k and two from k . So, the first term of the integral can be re-written as 1 k 2 (2k i k j − k i k j − k i k j )(2k l k m − k l k m − k l k m ) ϕ il (X 2 , X 2 , k ) ϕ jm (X 2 , X 2 , k − k ) = 4k l k m k 2 (k i k j − 2k i k j + k i k j ) ϕ il (X 2 , X 2 , k ) ϕ jm (X 2 , X 2 , k − k ).

(C 28)
The formulation that appears in Hodgson (1961) is slightly different. There we can find in (B.11) 2k l k m k 2 (2k i k j − k i k j − k i k j ) ϕ il (X 2 , X 2 , k ) ϕ jm (X 2 , X 2 , k − k ).

(C 29)
However, this formula lacks all the fourth-order monomials where two factors are taken from k and two from k . This is most likely a typewriting error.
C.4. Turbulence modelling A possible expression of the cross-spectrum of two components of the turbulent velocity fluctuations is given by ϕ ij (X 2 , X 2 , k) = u 2 (X 2 ) 16 π L 4 (k 2 δ ij − k i k j ) e −k 2 L 2 /4 F 2 , (C 30) where i, j ∈ {1, 3} and the use of a Gaussian function indicates a kind of 'isotropic' distribution of energy in the plane of the wavenumbers (k 1 , k 3 ). Similarly, the crossspectrum of vertical velocity fluctuations can be defined as ϕ 22 (X 2 , X 2 , k) = u 2 (X 2 ) 16 π L 4 k 2 e −k 2 L 2 /4 F 2 .
(C 31) In both (C 30) and (C 31), the term that represents the distribution of turbulent kinetic energy over the planar wavenumbers (k 1 , k 3 ) is multiplied by a factor F 2 that represents the correlation of turbulence between planes parallel to the wall, having vertical coordinates X 2 and X 2 , respectively. Consistently with the theory presented in § 2.3, the correlation of velocity between planes parallel to the wall should be expressed as F G 2 = e − ( X 2 −X 2 ) 2 /L 2 (C 32) that is also a Gaussian function. In particular, using (C 32) in (C 31), equation (2.22) is immediately obtained. However, Hodgson (1961) expresses this coefficient as Since the exponential function is flatter than the Gaussian, the corresponding turbulence field is significantly correlated over a greater part of the boundary layer. The cross-spectral densities obtained by substituting equation (C 33) in (C 31) or (C 30) cannot be derived from a given longitudinal correlation function, as shown in appendix B. However, they are grounded in experimental evidence, being the Fourier transform of the correlation coefficients given by Grant (1958). In any case, the assumption behind this formula is that the turbulent kinetic energy, K, is equally distributed among the three spatial directions, so that u 2 (X 2 ) = 2 3 k T (X 2 ). (C 34) Furthermore, Hodgson (1961) argues on the basis of a mirror flow model that the contribution of the turbulence cross-spectra of the kind ϕ i2 , with i ∈ {1, 3}, is small enough to be negligible. Accepting this assumption and substituting equations (C 28), (C 30), (C 31) in (C 26), the following expression is obtained: ϕ TT pp (k) = ρ 2 0 L 8 256 π 2 ∞ 0 u 2 (X 2 ) u 2 (X 2 ) e −k(X 2 +X 2 ) F 2 2 dX 2 dX 2 × R 2 e −k 2 L 2 /4 e −|k−k | 2 L 2 /4 × 4k l k m k 2 (k i k j − 2k i k j + k i k j ) (k 2 δ il − k i k l ) × (|k − k | 2 δ jm − (k j − k j )(k m − k m )) + 4(k i − k i )(k l − k l ) (k 2 δ il − k i k l ) |k − k | 2 dk .

Nomenclature
C f friction coefficient c airfoil chord c 0 speed of sound erf error function erfc complementary error function p F q (a 1 , . . . , a p ; b 1 , . . . , b q ; x) generalised hypergeometric function F Fourier transform operator F turbulence vertical correlation function J 0 Bessel function of the first kind of order 0 k (k 1 , k 3 ), planar wavenumber vector