Advances in the Prediction of the Statistical Properties of Wall-Pressure Fluctuations under Turbulent Boundary Layers

Abstract

Analytical or empirical models of the wall-pressure power spectral density under a turbulent boundary layer are often validated on test cases in an incompressible flow regime. In this work, an analytical model based on the compressible Poisson equation for the unsteady pressure in a turbulent boundary layer is developed. The Large Eddy Simulation of the flow over a controlled-diffusion airfoil at Mach 0.5 is used to validate the assumptions made on the statistical properties of the boundary layer turbulence and to validate the prediction of the statistics of the wall-pressure fluctuations. The predicted wall-pressure spectrum also compares favorably with experimental data.

1. Introduction

The quantification of the statistical properties of the pressure fluctuations on a rigid surface under a turbulent boundary layer is a topic related to many industrial applications and has been a matter of investigation for several decades (see, for instance, the early review of Willmarth [1]). Turbulent boundary layers develop on the surfaces of all sorts of means of transportation, generating wall-pressure fluctuations which cause structural vibrations and cabin noise. From the acoustic point of view, turbulent wall-pressure fluctuations are scattered as sound waves at the trailing edge of any lifting surface. This is the dominant noise generation mechanism in fans operating with a laminar inlet flow and without interaction with other lifting surfaces [2].

Several empirical models of the single-point wall-pressure power spectral density (PSD) are available in literature [3]. The extensive review made by Lee [4] shows that they all assume the same universal wall-pressure shape and then fit the coefficients of this generic function to a given database. The limitation of these kinds of models is that their scope of application depends on the characteristics of the database on which they have been calibrated.

Another category of models is based on the Poisson equation linking the unsteady pressure to the characteristics of the turbulent flow in the boundary layer. Thanks to this link to the flow physics, these models have a potentially broader scope of application than empirical models. Some models use an exact analytical solution to the Poisson equation, such as [5,6,7]. Other models are based on the Fourier transform of the Poisson equation in the frequency-wavenumber domain, which yields an inhomogeneous Helmholtz equation. The TNO–Blake model (see the reviews by [8,9]) and the model of Lysak [10] use an approximate solution of the Helmholtz equation. This approximation has been demonstrated to be based on a questionable hypothesis concerning the statistics of the boundary-layer turbulence. The exact analytical solution to the Helmholtz equation can be used as demonstrated in the works of Panton and Linebarger [11], Remmler et al. [12] and Fischer et al. [13], at the price of a slight increase of computational time. An extensive review of the application of empirical or analytical models of wall-pressure statistics to the prediction of airfoil trailing-edge noise can be found in the work of Lee et al. [14].

The predictive approach based on the exact solution of the Helmholtz equation has been investigated in [15] with the help of direct numerical simulations of the Navier–Stokes equations [16] and experiments [17] conducted on a controlled diffusion (CD) aifoil [18] in the incompressible flow regime. In particular, the dependence of wall-pressure statistics on those of the turbulent flow in the boundary layer has been elucidated. The prediction method presented in [15] has been applied in [19] to the direct numerical simulation of the same test case with the Lattice–Boltzmann method, in order to study its sensitivity to input data variations, and it has proved its robustness. The present work aims at extending the method to the subsonic compressible flow regime by re-deriving the analytical formulation starting from the compressible form of the Poisson equation. The database selected for the validation of the model presented in this work consists of the Large Eddy Simulation (LES) of the flow around a CD airfoil, analogous to that studied in [15,16,17], at Mach $0.5$, $R{e}_c=2.29\times {10}^6$ and $4^{\circ }$ angle of attack performed by ISAE-Supaéro (see [20,21] for an exhaustive description of the computational technique and in-depth analysis of the results). The CFD data will be used to extract the information about the boundary layer turbulence that feeds the wall-pressure power spectral density model and it will also be used to compare the predicted versus the directly computed wall-pressure spectra.

The model presented in this work may be used in the context of shape optimization of blades, for instance. In fact, due to the high number of possible geometries that are analyzed, it is most common to use steady-state simulations, which do not compute the wall-pressure fluctuations directly. This information can be recovered by applying the predictive model to the steady-state velocity and turbulent kinetic energy profiles, while making a few assumptions on the statistical properties of the boundary layer turbulence, especially concerning the shape of its correlation function, its integral length scale and its anisotropy. The analysis of the LES simulation data presented in this work will guide the reader in making these assumptions.

The single and multi-point statistics of the wall-pressure fluctuations calculated in this work can be used directly in the prediction of airfoil broadband trailing-edge noise with Amiet’s statistical model and its extensions, for instance [2,22,23,24,25,26]. Indeed, within this approach, the PSD of radiated sound is demonstrated to be directly proportional to the product of the single-point wall-pressure PSD by the transverse coherence length, with both quantities taken in the close proximity of the trailing edge. Moreover, the convective speed of wall-pressure fluctuations, $U_c$, is used in Amiet’s theory to calculate the convective wavenumber, $K_c=\omega /U_c$, which is the sole responsible for the noise emitted at frequency $\omega$ in the hypothesis of frozen turbulence convection.

2. Theoretical Background

This section summarizes the demonstration of the expression of the unsteady pressure beneath a turbulent boundary layer in compressible flow conditions, laid out in the paper of Gerolymos et al. [7]. The compressible Navier–Stokes equations form the basis for this analytical model. The rationale is the same as in the incompressible problem developed in [15], but with the additional density variable. The first step consists in taking the divergence of the compressible momentum equations, which yields

$$\frac{{\partial }^2}{\partial x_i\partial t}\left(\rho u_i\right)+\frac{{\partial }^2}{\partial x_i\partial x_j}\left(\rho u_i{u}_j\right)=-{\nabla }^{2}p+\frac{{\partial }^2{\tau }_{ij}}{\partial x_i\partial x_j}+\frac{\partial }{\partial x_i}\left(\rho f_{Vi}\right)$$

(1)

where $f_{Vi}$ is the body acceleration. Subtracting from Equation (1), its Reynolds average yields the following form of the Poisson equation for the fluctuating pressure

$$\begin{array}{cc}\hfill {\nabla }^2{p}^{\prime }=& \underset{\mathrm{Quasi}-\mathrm{incompressible}\mathrm{terms}}{\underbrace{\left[-2\overline{\rho }\left(\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial {\overline{u}}_j}{\partial x_i}\right)\right]+\left[-\overline{\rho }\left(\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial u_j^{\prime }}{\partial x_i}-\overline{\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial u_j^{\prime }}{\partial x_i}}\right)\right]}}\hfill \\ \hfill +& \left[-\left({\rho }^{\prime }\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}-\overline{{\rho }^{\prime }\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}}\right)\right]+\left[-2\left({\rho }^{\prime }\frac{\partial u_i}{\partial x_j}-\overline{{\rho }^{\prime }\frac{\partial u_i}{\partial x_j}}\right)\frac{\partial {\overline{u}}_j}{\partial x_i}\right]\hfill \\ \hfill +& \left[-{\rho }^{\prime }\frac{\partial {\overline{u}}_i}{\partial x_j}\frac{\partial {\overline{u}}_j}{\partial x_i}\right]+\frac{{\delta }^2{\tau }_{ij}^{\prime }}{\partial x_i\partial x_j}+\frac{\partial }{\partial x_i}(\rho f_{Vi}-\overline{\rho f_{Vi}})\hfill \\ \hfill +& {\left[\rho \frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{1}{\rho }\frac{\mathrm{D}\rho }{\mathrm{D}t}\right)\right]}^{\prime }+{\left[-\frac{\mathrm{D}\mathbf{V}}{\mathrm{D}t}·\nabla \rho \right]}^{\prime }\hfill \end{array}$$

(2)

with the following definitions:

$$\overline{\frac{\mathrm{D}(·)}{\mathrm{D}t}}=\overline{\frac{\partial (·)}{\partial t}+u_i\frac{\partial (·)}{\partial x_i}}=\frac{\partial \overline{(·)}}{\partial t}+\overline{u_i}\frac{\partial \overline{(·)}}{\partial x_i}+\overline{u_i^{\prime }\frac{\partial {(·)}^{\prime }}{\partial x_i}}$$

and

$${\left[\frac{\mathrm{D}(·)}{\mathrm{D}t}\right]}^{\prime }=\frac{\mathrm{D}(·)}{\mathrm{D}t}-\overline{\frac{\mathrm{D}(·)}{\mathrm{D}t}}.$$

The subscripts $i,j\in \left\{1,2,3\right\}$, where the indices correspond to the coordinate system depicted in Figure 1. The first two source terms on the right-hand side of Equation (2) are known as quasi-incompressible terms and correspond respectively to the turbulence-mean shear and turbulence–turbulence interaction source terms of the incompressible Poisson equation, reproduced below for simplicity:

$${\nabla }^2{p}^{\prime }=\underset{\mathrm{Turbulence}-\mathrm{mean}\mathrm{shear}}{\underbrace{-2{\rho }_0\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial {\overline{u}}_j}{\partial x_i}}}-\underset{\mathrm{Turbulence}-\mathrm{turbulence}}{\underbrace{{\rho }_0\left(\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial u_j^{\prime }}{\partial x_i}-\overline{\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial u_j^{\prime }}{\partial x_i}}\right)}}.$$

(3)

The quasi-incompressible terms scale with the mean density component, $\overline{\rho }$, whereas the remaining terms scale with the density fluctuations, ${\rho }^{\prime }$, as demonstrated by the dimensional analysis of ([7] Appendix A). Morkovin’s hypothesis [27] implies that for a low-Mach-number flow, the effects of compressibility on turbulence are due to the gradient of the mean density, the influence of the density fluctuations being a higher-order effect. This holds true in the limit of ${\rho }^{\prime }/\overline{\rho }\ll 1$. Within this hypothesis, it is justified to retain only the quasi-incompressible terms in the solution to the compressible Poisson equation governing pressure fluctuations in a compressible boundary layer, thus obtaining an equation of the same structure as Equation (3). Furthermore, it is possible to discard in first approximation the turbulence–turbulence interaction quasi-incompressible term. This is a widely used assumption, following the work of Kraichnan [28] and Hodgson [29], who estimated that the turbulence–turbulence term makes up 4–6% of the mean-square unsteady pressure. In particular, the analytical solution of Equation (3) with the turbulence–turbulence component alone on the right-hand side, developed by Hodgson, was re-derived and generalized in [15] and applied to the CD airfoil configuration in the incompressible flow regime. It was thus confirmed that the turbulence–turbulence component is negligible with respect to the turbulence–mean shear over the low-to-middle frequency range. Therefore, the quasi-incompressible Poisson equation is

$${\nabla }^2{p}^{\prime }=-2\overline{\rho }\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial {\overline{u}}_j}{\partial x_i}.$$

(4)

The mean density, $\overline{\rho }$, at the right-hand side of Equation (4), varies in space, which is the mean difference between this equation and the incompressible one. Performing the Fourier transform of Equation (4) in space along the wall directions $(x_1,x_3)$ in time and applying the classical boundary-layer approximation, for which the mean flow speed is aligned with the $x_1$ axis, and the only non-vanishing mean velocity gradient is that of the stream-wise velocity ${\overline{u}}_1$ in the wall-normal direction $x_2$ (as depicted in Figure 1), we obtain the following inhomogeneous Helmholtz equation:

$$\frac{\partial \widehat{p}(\omega ,\mathit{k},x_2)}{\partial x_2^2}-k^2\widehat{p}(\omega ,\mathit{k},x_2)=-2\overline{\rho }\left(x_2\right)\mathrm{i}{k}_1\frac{\partial {\overline{u}}_1\left(x_2\right)}{\partial x_2}{\widehat{u}}_2(\omega ,\mathit{k},x_2),$$

(5)

where $\mathit{k}=(k_1,k_3)$. The boundary-layer model used in the derivation of the previous equation is valid for actual boundary layers as long as their thickness is much smaller than the local radius of curvature of the surface (i.e., quasi-parallel flow). Assuming that the unsteady pressure vanishes at an infinite distance from the wall, that is

$$\underset{x_2\to \infty }{lim}{p}^{\prime }=0,$$

and adopting the rigid-wall model whereby

$${\left.\frac{\partial p^{\prime }}{\partial x_2}\right|}_{x_2=0}=0$$

we can write the well-known Green’s function solution (see for instance [7]) to the differential equation above as

$$\widehat{p}(\omega ,\mathit{k},x_2=0)=2\mathrm{i}\frac{k_1}{k}{\int }_0^{+\infty }\overline{\rho }\left(X_2\right)\widehat{u}(\omega ,\mathit{k},X_2)\frac{\partial {\overline{u}}_1\left(X_2\right)}{\partial x_2}{\mathrm{e}}^{-k{X}_2}\mathrm{d}{X}_2.$$

(6)

Multiplying Equation (6) by its complex conjugate and taking the ensemble average of the result yields the general expression of the wall-pressure power spectral density:

$$\begin{array}{cc}\hfill {\Phi }_{pp}^{TM}(\omega ,k_1,k_3)=& 4\frac{k_1^2}{k^2}{\int }_{{\omega }^{\prime }}{\int }_{{\mathit{k}}^{\mathbf{\prime }}}\int {\int }_0^{\infty }\overline{\rho }\left(X_2\right)\overline{\rho }\left(X_2^{\prime }\right){\mathrm{e}}^{-(X_2+X_2^{\prime })k}\frac{\partial U_1\left(X_2\right)}{\partial x_2}\frac{\partial U_1\left(X_2^{\prime }\right)}{\partial x_2}\hfill \\ \hfill & \times \langle \widehat{u}(\omega ,\mathit{k},X_2)\widehat{u}({\omega }^{\prime },{\mathit{k}}^{\mathbf{\prime }},X_2^{\prime })\rangle \mathrm{d}{X}_2\mathrm{d}{X}_2^{\prime }\mathrm{d}{\mathit{k}}^{\mathbf{\prime }}\mathrm{d}{\omega }^{\prime }.\hfill \end{array}$$

(7)

The superscript $TM$ reminds one that this solution only takes into account the turbulence–mean shear interaction source term.

Due to the statistical orthogonality of the turbulence wave-vectors [30], the ensemble average of the product of vertical velocity fluctuations can be written as

$$\langle \widehat{u}(\omega ,\mathit{k},X_2){\widehat{u}}^{*}({\omega }^{\prime },{\mathit{k}}^{\mathbf{\prime }},X_2^{\prime })\rangle ={\phi }_{22}(\omega ,\mathit{k},X_2,X_2^{\prime })\delta (\mathit{k}-{\mathit{k}}^{\mathbf{\prime }})\delta (\omega -{\omega }^{\prime })$$

and then the wall-pressure PSD becomes

$$\begin{array}{cc}\hfill {\Phi }_{pp}^{TM}(\omega ,k_1,k_3)=& 4\frac{k_1^2}{k^2}\int {\int }_0^{\infty }\overline{\rho }\left(X_2\right)\overline{\rho }\left(X_2^{\prime }\right){\mathrm{e}}^{-(X_2+X_2^{\prime })k}\frac{\partial U_1\left(X_2\right)}{\partial x_2}\frac{\partial U_1\left(X_2^{\prime }\right)}{\partial x_2}\hfill \\ \hfill & \times {\phi }_{22}(\omega ,\mathit{k},X_2,X_2^{\prime })\mathrm{d}{X}_2\mathrm{d}{X}_2^{\prime }.\hfill \end{array}$$

(8)

According to Equation (8), the frequency-wavenumber PSD of the wall-pressure flucutation is directly proportional to the normal-to-wall variation of mean density, mean shear and to the cross-spectral density of vertical velocity fluctuations, ${\phi }_{22}$. The first two quantities can be computed by means of a steady-state compressible flow simulation. The latter, however, needs to be modeled for practical applications. The turbulent cross-spectral density can be expressed as the product of a wavenumber cross-spectrum and a moving-axis spectrum, where the latter expresses the frequency dependency:

$${\phi }_{22}(\omega ,\mathit{k},X_2,X_2^{\prime })={\phi }_{22}(\mathit{k},X_2,X_2^{\prime }){\Phi }_m(\omega -U_c{k}_1).$$

(9)

The wavenumber term can be represented with various analytical formulations according to different theories. The Generalized von Kármán Spectrum, as defined in [31,32], is a useful formulation of the wavenumber term because, depending on the value of one parameter, it reduces to either the classical von Kármán, the Liepmann or the Rapid Distortion Theory turbulence spectrum. Assuming homogeneous and isotropic turbulence, the Generalized von Kármán spectrum has the following expression:

$${\phi }_{22}\left(\mathit{k},x_2,x_2^{\prime }\right)=\frac{\sqrt{\overline{u_2^2}\left(x_2\right)\overline{u_2^2}\left(x_2^{\prime }\right)}{l}^2{\tilde{k}}^2{\zeta }^{\nu +2}}{\Gamma \left(\nu \right)\pi 2^{\nu +1}{\left(1+{\tilde{k}}^2\right)}^{\nu +2}}{\mathrm{K}}_{\nu +2}\left(\zeta \right)$$

(10)

with the parameter

$$\zeta =\frac{\parallel x_2-x_2^{\prime }\parallel }{l}\sqrt{1+{\tilde{k}}^2}.$$

(11)

and the characteristic length scale

$$l=\frac{\Gamma \left(\nu \right)}{\sqrt{\pi }\Gamma (\nu +1/2)}\Lambda ,$$

(12)

where $\Lambda$ is the turbulent longitudinal integral length scale. If $\nu =1/3$, Equation (10) corresponds to the model of von Kármán [33], if $\nu =1/2$ to Liepmann’s model [34] and if $\nu =7/6$ to the rapid distortion theory (RDT) [35]. In order to close the model, in earlier works [15,19] Taylor’s hypothesis of frozen turbulence convection has been applied, whereby the moving axis spectrum is a Dirac delta function centered at the convective frequency $\omega =U_c{k}_1$:

$${\Phi }_m^D(\omega -U_c{k}_1)=\delta (\omega -U_c{k}_1),$$

(13)

where $U_c$ is the convective speed of velocity fluctuations. In this case, the wall-pressure spectral density at a given frequency $\omega$ is uniquely determined by turbulent fluctuations at the streamwise wavenumber $K_1=\omega /U_c$. However, the Dirac formulation of the moving-axis spectrum does not ensure the decay of the wall-pressure coherence function in the streamwise direction, $x_1$. Furthermore, it can be interesting to keep the functional dependency of the wall-pressure PSD expressed in Equation (8) on both $\omega$ and $k_1$. To solve both problems, the formulation of the moving-axis spectrum given in the TNO-Blake model (see [8]) as a Gaussian function centered around the convective frequency, can be used:

$${\Phi }_m^G(\omega -U_c{k}_1)=\frac{1}{{\alpha }_2\sqrt{2\pi }}exp\left[-\frac{{(\omega -U_c{k}_1)}^2}{2{\alpha }_2^2}\right]$$

(14)

with

$${\alpha }_2=\frac{c_6{U}_c}{{\Lambda }_2}.$$

In the original TNO model, the coefficient $c_6$ is assumed constant and equal to 0.05.

Finally, this model needs to account for the anisotropy of turbulent structures, defined as the ratio of longitudinal-to-transverse integral length scales, $\alpha$. Following Schlinker and Amiet [36], the anisotropic turbulence vertical cross-spectral density, ${\phi }_{22}^{NI}$, is defined with respect to the isotropic formulation of Equation (10), ${\phi }_{22}^I$, as

$${\phi }_{22}^{NI}(k_1,x_2,x_2^{\prime },k_3)=\alpha {\phi }_{22}^I(\alpha k_1,x_2,x_2^{\prime },k_3).$$

(15)

In principle, the anisotropic coefficient $\alpha$ varies in the normal-to-wall direction; however, it has been shown in [15,19] that a constant representative value of $\alpha$ may suffice for the prediction of the wall-pressure spectral density in the frequency range of interest.

3. Study of the Controlled-Diffusion Airfoil Test Case

3.1. Large Eddy Simulation of the Flow around a Controlled-Diffusion Airfoil

This section presents the post-processing of the LES of the flow around the CD airfoil, which was in turn performed at ISAE-Supaéro [20,21] in the framework of the European project SCONE. The wall-resolved LES (WRLES) has been performed by means of the massively parallel compressible solver CharLES^X, which offers a low dissipation cell-centered finite-volume scheme. The numerical solution method is second-order accurate in space and third-order accurate in time, using a low-storage Runge–Kutta–Wray scheme for explicit time advancement. The sub-grid turbulent scales (SGS) are represented by Vreman’s model. The unstructured computational grid has an O-grid topology and is made up of 225.5 M cells. The computational domain extends $20c$, where c is the airfoil chord, in both the streamwise and normal-to-wall directions and $0.1c$ in the transverse direction. The dimensionless wall distance, $y^{+}$, on the CD airfoil surface is below $1.12$. The normalized time step is $2.68\times {10}^{-5}$. The mesh resolution quality has been verified by means of Pope’s criterion [37], which measures the ratio of resolved to total (i.e., sum of resolved and SGS) turbulent kinetic energy. This ratio is found to be at least $0.94$ for all the WRLES simulations performed in the SCONE project, meaning that the resolution level is close to that of a DNS mesh.

3.2. Boundary Layer Profiles Close to the Trailing Edge of the Airfoil

As observed in the previous section, the boundary layer profiles and the statistical properties of the turbulent flow are the necessary input data to the wall-pressure PSD analytical model. These quantities are extracted from the LES simulation close to the trailing edge of the CD airfoil, at $2\%$ of chord upstream of the trailing edge, corresponding to the position of pressure sensor number 26 in the experimental setup (see Figure 2). It is interesting to predict the wall-pressure statistics as close as possible to the edge of the profile, because they determine the amplitude and frequency distribution of the corresponding radiated noise (see [26] for an investigation of this dependence). Figure 3 and Figure 4 show the mean stream-wise and the RMS normal-to-wall velocities, respectively. The mean density, $\overline{\rho }$, has a moderate $4\%$ variation from the wall to the edge of the boundary layer, as can be observed in Figure 5.

3.3. Turbulence Statistics

3.3.1. Streamwise and Transverse Velocity Correlation Coefficient

Unsteady velocity time series have been extracted from the LES only on three mesh layers parallel to the airfoil surface, depicted in Figure 2. Furthermore, as can be observed in Figure 2, the upper mesh layer lays outside the boundary layer; therefore, the study of turbulence statistics will be focused on the two inner mesh layers.

The first statistical property of the turbulent field that is calculated is the correlation coefficient. Considering a homogeneous turbulence field, this quantity is defined as

$$R_{ij}({\eta }_1,{\eta }_2,{\eta }_3)=\frac{\langle u_i(x_1,x_2,x_3)·u_j(x_1+{\eta }_1,x_2+{\eta }_2,x_3+{\eta }_3)\rangle }{{\sigma }^2}$$

(16)

where ${\sigma }^2$ is the variance of turbulent velocity fluctuations. In the framework of the Generalized von Kármán theory, this correlation coefficient is expressed analytically as

$$\begin{array}{cc}\hfill R_{ij}({\eta }_1,{\eta }_2,{\eta }_3)=& \frac{1}{2^{\nu -1}\Gamma \left(\nu \right)}{\left(\frac{\eta }{l}\right)}^{\nu }\left\{\frac{{\eta }_i{\eta }_j}{{\eta }^2}{K}_{\nu }\left(\frac{\eta }{l}\right)\right.\hfill \\ \hfill & \left.+{\delta }_{ij}\left(1-\frac{{\eta }_i{\eta }_j}{{\eta }^2}\right)\left[(\nu +1)K_{\nu }\left(\frac{\eta }{l}\right)-\frac{1}{2}\frac{\eta }{l}{K}_{\nu +1}\left(\frac{\eta }{l}\right)\right]\right\}\hfill \end{array}$$

(17)

where $\eta =\sqrt{{\eta }_1^2+{\eta }_2^2+{\eta }_3^2}$ and $K_{\nu }$ is the modified Bessel function of the second kind of order $\nu$. The correlation coefficient $R_{22}$ is of particular interest because it is related to the vertical velocity cross-spectral density, ${\phi }_{22}$, by the double spatial Fourier transform:

$${\phi }_{22}(k_1,{\eta }_2,k_3)=\frac{{\sigma }^2}{4{\pi }^2}\int {\int }_{-\infty }^{+\infty }{R}_{22}({\eta }_1,{\eta }_2,{\eta }_3)cos\left({\eta }_1{k}_1\right)cos\left({\eta }_3{k}_3\right)\mathrm{d}{\eta }_1\mathrm{d}{\eta }_3.$$

(18)

Therefore, the choice of the analytical formulation of ${\phi }_{22}$ could be validated directly by calculating the vertical velocity correlation coefficient, $R_{22}\left({\eta }_2\right)$, from the CFD data and comparing it with Equation (17) (see [15], Figure 12 for an example of this comparison). However, as the unsteady velocity time series has been recorded only on two mesh layers inside the turbulent boundary layer, there are not enough points for comparing the analytical and the directly computed correlation. For this reason, the choice of the value of $\nu$ in Equations (10) and (17) will be based on the calculation of the streamwise ($R_{11}\left({\eta }_1\right)$) and transverse ($R_{33}\left({\eta }_3\right)$) velocity correlation coefficients, under the assumption that the vertical correlation behaves in the same way.

Figure 6 and Figure 7 depict the streamwise velocity correlation coefficient calculated on mesh layers 1 and 2, respectively. The following Figure 8 and Figure 9 depict the transverse velocity correlation coefficient on the same mesh layers parallel to the airfoil surface. It can be observed in these plots that the Rapid Distortion Theory (RDT) [35] convincingly describes both the longitudinal and the transverse correlation coefficients inside the boundary layer. This is consistent with the finding of Magnaudet (see [38]) that the RDT is the leading-order approximation describing the short-term and long-term evolution of turbulent boundary layers in the limit of the large Reynolds number. Furthermore, the fundamental assumption of the RDT is that the turbulence–turbulence interaction is negligible with respect to the turbulence–mean shear interaction (see also [30]), which is exactly the same hypothesis that has been made in the derivation of Equation (8). Therefore, the use of the RDT to model the vertical velocity cross-spectral density, ${\phi }_{22}$, is arguably the most consistent choice with the fundamental assumption of the model that has been laid out in Section 2.

3.3.2. Integral Length Scale and Anisotropy Coefficient

The longitudinal and transverse integral length scales are available only at three locations above the wall, as shown in Figure 10. However, a continuous function of $\Lambda$ versus the distance from the wall, $x_2$, is necessary in order to calculate the integral in Equation (8). For this reason, a working hypothesis is made to represent the transverse length scale with a physically-consistent model based on the hyperbolic tangent function, following the empirical considerations of Parchen [39] (see also [8] and [15], Equation (2.38)). The model is here re-scaled in order to match the transverse length scale as closely as possible. Considering the available values of the ratio between the longitudinal and transverse scales, the value of $\alpha =3.5$ will be used in the prediction of wall-pressure statistics.

4. Direct Computation and Empirical Modeling of Wall-Pressure Statistics

The time series of wall-pressure fluctuations, $p^{\prime }$, recorded on the airfoil-surface grid, can be used to compute directly various statistical properties. Figure 11 shows the coherence function of $p^{\prime }$ for various values of the Strouhal number, $St$, taking the reference point at $12\%$ of the chord upstream of the trailing-edge and shifting the second point in the stream-wise direction. This reference point has been taken more upstream in order to be able to assess the coherence function both upstream and downstream. The plot shows a faster decay of the correlation function for a higher $St$. Figure 12 represents the phase of the cross spectrum of wall-pressure fluctuations versus the stream-wise coordinate. Assuming that the phase is linear, its slope is directly proportional to the convective speed of wall-pressure fluctuations, $U_c$, for a given frequency, $\omega$. Following the method described in [40], the linear function $\omega {\eta }_1/U_c$ is fitted to the phase in order to estimate the value of $U_c$. This value is needed for the prediction of the wall-pressure PSD because it is part of the definition of the moving-axis spectrum. Furthermore, it is needed in the analytical modeling of the trailing-edge noise, according to Amiet’s model, as mentioned in the introduction. As can be observed in the plots, the linear fit is not a trivial operation, neither at low nor at high frequency, due to the scarcity of the sample points.

The variation of $U_c$ with the frequency is shown is Figure 13. Interestingly, Smol’yakov’s empirical model [41] predicts a peak at the same value of the normalized abscissa $\omega {\delta }_1/U_{\infty }$, where ${\delta }_1$ is the boundary layer displacement thickness. However, the amplitude of the peak is subject to the uncertainty of the linear interpolation. On the contrary, for normalized frequencies above $0.6$, the ratio stabilizes around $U_c/U_{\infty }=0.5$. It is therefore safer to use this constant value of $U_c$ in the predictions presented in the next section.

Finally, Figure 14 presents the transverse coherence function ${\gamma }^2(\omega ,{\eta }_3)$, as a function of the chord-based Strouhal number, for various separations ${\eta }_3$ in the span-wise direction. As expected, the coherence decreases with a separation at any frequency. It exhibits a maximum for some relatively low frequencies and a drop at very low frequencies. Yet, this drop is hard to estimate accurately because of the suspected lack of convergence in the numerical simulations. A log-normal distribution of the coherence has often been found in previous low-speed experiments, as reported for instance in [42]. Attempting such a fit on the present data is less convincing at high frequencies, as illustrated in Figure 15. The theoretical expression plotted as the blue mesh reads

$${\gamma }^2(f,{\eta }_3)=\frac{A}{\sigma \sqrt{2\pi }St}{\mathrm{e}}^{-{({\log }_{10}\left(St\right)-{\log }_{10}\left(S{t}_0\right))}^2/\left(2{\sigma }^2\right)},$$

with $S{t}_0=1.5-0.06{\eta }_3$, $A=20.5{\mathrm{e}}^{-0.19{\eta }_3}$, $\sigma =0.53-0.02{\eta }_3$, the separation being in mm. In particular, the high-frequency behavior deviates from the trend at larger separations. It is also worth noting that the decay with ${\eta }_3$ at a given frequency is not exactly exponential, according to both the fit and the extracted data. Further investigation is needed to elucidate whether the coherence includes specific effects of compressibility.

5. Prediction of Wall-Pressure Statistics with the Analytical Model Based on the Poisson Equation

5.1. Single-Point Statistics

As a first example of the application of the model based on the Poisson equation, the single-point PSD of wall-pressure fluctuations is predicted by means of Equation (8). In this case, the wavenumber component of the cross-PSD of vertical velocity fluctuations is modeled by means of the RDT. Taylor’s hypothesis completes the turbulence model and allows us to reduce the computational time by eliminating one dimension of integration. The same strategy has been followed successfully in [15,19]. The input data are the boundary layer profiles presented in Section 3.2 and the integral length scale presented in Figure 10, along with the values of $\alpha$ and $U_c$ determined empirically from the LES results. The integration of Equation (8) is performed by means of a recursive stratified sampling Monte-Carlo method (see [15,43] for further details) implemented in the Python scikit-monaco library [44]. The wall-pressure PSD prediction is shown in Figure 16, where it is compared with a measured spectrum, obtained by NLR in the framework of the CRORTET project [45], and with the spectrum calculated directly from the LES. The measured and LES spectra are taken at $2\%$ of the chord upstream of the trailing edge, consistently with the extraction of the boundary layer profiles. The convincing overall agreement confirms the relevance of the model. In particular, the high-frequency roll-off is very well captured.

5.2. Multi-Point Statistics

Equation (8) can also be used to model the wall-pressure coherence function. Assuming homogeneity of the wall-pressure fluctuations on the wall surface, the corresponding wall-pressure coherence function can be formulated as the (normalized) inverse spatial Fourier transform of the wall-pressure PSD:

$${\gamma }^2(\omega ,{\eta }_j)=\frac{{\left|{\int \int }_0^{\infty }{\Phi }_{pp}(\omega ,k_1,k_3)cos\left(k_j{\eta }_j\right)\mathrm{d}{k}_1\mathrm{d}{k}_3\right|}^2}{{\left|{\int \int }_0^{\infty }{\Phi }_{pp}(\omega ,k_1,k_3)\mathrm{d}{k}_1\mathrm{d}{k}_3\right|}^2}$$

(19)

where $j\in \left\{1,3\right\}$. The substitution of the frequency-wavenumber wall-pressure spectrum given by Equation (8) in Equation (19) for the calculation of the coherence implies the use of the Gaussian moving-axis spectrum formulation of Equation (14) in order to separate the frequency and wavenumber variables.

The Dirac formulation of the moving-axis spectrum following the hypothesis of frozen turbulence, on the contrary, is not applicable to the calculation of the coherence function because it does not ensure the decay of the coherence with distance. For this reason, Blake [46] proposes the following empirical correction to the Dirac formulation, which introduces an exponential decay in the longitudinal direction:

$${\Phi }_m(\omega -U_c{k}_1,{\eta }_1)=\delta (\omega -U_c{k}_1){\mathrm{e}}^{-{\gamma }_B|\omega {\eta }_1/U_c|}$$

(20)

where ${\gamma }_B$ is an empirical parameter ranging between $0.2$ and $0.3$. Substituting Equation (20) in Equation (8), the longitudinal coherence function is simply

$${\gamma }^2(\omega ,{\eta }_1)={\mathrm{e}}^{-2{\gamma }_B|\omega {\eta }_1/U_c|}.$$

(21)

Figure 17 compares the longitudinal wall-pressure coherence extracted from the LES simulation and predicted by means of the numerical integration of Equation (19), using the Gaussian formulation of the moving-axis spectrum given in Equation (14). The empirical exponential decay of Equation (21) is also plotted in Figure 17, using the value ${\gamma }_B=0.3$. The reference point is the same for which the single-point statistics have been computed; the ${\eta }_1$ coordinate is directed upstream. The Gaussian formulation of the moving-axis spectrum better captures the LES-predicted trend at the lowest Strouhal number, whereas for increasing frequency the difference between the two formulations decreases. Interestingly, the Gaussian formulation predicts a second smaller local maximum of the longitudinal coherence function, which is difficult to see in the LES data, especially at high frequencies, due to the relative coarseness of the surface mesh.

Figure 18 shows the transverse correlation as a function of the Stouhal number for two values of the transverse separation, since the LES time series of pressure fluctuations have been recorded for a smaller number of points in the ${\eta }_3$ direction. The LES transverse coherence is compared with the model based on the Gaussian moving axis spectrum alone, since Equation (21) does not apply to the transverse direction. The modeled curve is close to the LES data for low Strouhal numbers, but the LES decays faster after $St=10$. This divergence may be partially explained with the limited length of the time series available for the Fourier transform.

6. Conclusions

This work has extended a model of the statistical properties of wall-pressure fluctuations, initially developed within the assumption of incompressible flow, to the subsonic compressible flow regime. A simple analytical formulation is obtained thanks to the assumption that the effect of density fluctuations in the boundary layer is negligible with respect to that of the spatial variation of mean density. The LES simulation of the flow over a CD airfoil at Mach $0.5$ has been taken as validation test case. Similarly to what was found in previous incompressible-regime test cases [15,19], the RDT model reasonably predicts the correlation function of boundary-layer turbulence; it has therefore been used for the prediction of wall-pressure PSD with good agreement with the experimental and directly computed results.

Empirical and semi-analytical models are usually developed separately for the single-point wall-pressure PSD and for the multi-point coherence function. In this work, on the contrary, a unified analytical framework for the prediction of both single and multi-point statistics has been presented. Although the comparison of the model with directly computed coherence functions suffers from the limited length of the recorded time signal, the model captures the main trends. It will be interesting to apply the same analysis presented in this work to the results of scale-resolving simulations, where the flow is extracted at more mesh layers above the wall and for longer time intervals.

The model hereby presented will be used to estimate the sources of the trailing-edge noise of airframe operating in subsonic conditions on the basis of steady-state flow simulations, thus enabling the implementation of a low-computational cost acoustic evaluation in the early stages of the design of innovative geometries.