1. Introduction
Leading-edge noise, also known as aerofoil–turbulence interaction noise, is produced by the scattering of turbulent velocity fluctuations of a given incoming flow by the aerofoil's leading edge. It is a dominant noise mechanism for many applications, such as wind turbines, cooling systems and turbofan engines. Regarding wind turbines, it is known that the interaction of the blades with atmospheric turbulence causes unwanted noise, particularly at low frequencies (Buck, Oerlemans & Palo Reference Buck, Oerlemans and Palo2018). For engines with multi-row rotor systems, the wake from turbulence interacting with rotor blades impinges on the downstream blade and is a dominant noise source (Lyu, Ayton & Chaitanya Reference Lyu, Ayton and Chaitanya2019). Alongside these industrial reasons, we may also consider social reasons, such as counteracting harmful noise pollution.
In the literature, two methods of particular interest have been developed experimentally and analytically to reduce leading-edge noise. The first is using serrations on the leading edge (Roger, Schram & De Santana Reference Roger, Schram and De Santana2013; Ayton & Kim Reference Ayton and Kim2018; Ayton & Paruchuri Reference Ayton and Paruchuri2018; Lyu & Ayton Reference Lyu and Ayton2020; Teruna et al. Reference Teruna, Avallone, Casalino and Ragni2021). Serrations have been shown to reduce noise across various frequencies – however, the mathematical modelling can be more involved and introduce aerodynamic difficulties such as drag and loss of lift.
The other popular approach is to use porous/impedance plates, in which either full or partial porosity is introduced and shown to have acoustic benefits (Geyer, Sarradj & Giesler Reference Geyer, Sarradj and Giesler2012; Roger et al. Reference Roger, Schram and De Santana2013; Ayton et al. Reference Ayton, Colbrook, Geyer, Chaitanya and Sarradj2021a; Teruna et al. Reference Teruna, Avallone, Casalino and Ragni2021). Currently, there are limitations on how we model porosity mathematically and the need to balance the aeroacoustic benefits of porosity with aerodynamic costs. The focus of this paper is somewhat different; rather than modifying the aerofoil to reduce unwanted noise (thereby impacting the aerodynamics), we introduce a novel methodology to understand how altering the flow can affect the noise from an aerofoil.
We achieve our aim through a theoretical model of aerofoil–turbulence interaction noise. Early work by Amiet describes how to derive a cross-power spectral density (cross-PSD) of the surface pressure for an aerofoil in turbulent flow by considering first a single ‘gust’ component of the turbulent flow and then integrating it over a wavenumber–frequency spectrum (Amiet Reference Amiet1975, Reference Amiet1976). The solution relates the turbulent velocity experienced by the plate to the pressure jump along the plate via a transfer function 
$g$. It comprises two parts: a gust-scattering solution and a turbulence spectrum. The first part requires details of the aerofoil and its rigid boundary condition, whilst the latter requires details of the incoming turbulence and any anisotropy.
Our modelling follows the same path as Amiet's, and we first focus on the gust scattering. Rather than employing Amiet's method exactly, we use a more direct approach to explicitly calculating the far field noise from a single gust without requiring Curle's integral. This is achieved using the Wiener–Hopf technique instead of Amiet's approach via Schwarzschild's solution (Amiet Reference Amiet1975). To further assist theoretical progress, we assume, as done so by Amiet, that the aerofoil is an infinitesimally thin leading edge extending infinitely downstream, thus occupying 
$x>0$, 
$y=0$. We suppose the plate extends infinitely in the spanwise 
$z$-direction. This geometric assumption is commonly known as thin aerofoil theory. It has been shown to have good agreement with experimental data in the case of thin aerofoils in (isotropic) turbulent streams.
Turning to the second part required for modelling leading-edge noise, we need a turbulence spectrum detailing anisotropic fluctuations. Most prior theoretical work on leading-edge noise only considers an isotropic spectrum, which may be constructed via Gaussian decomposition (Wohlbrandt et al. Reference Wohlbrandt, Hu, Guérin and Ewert2016). Anisotropic turbulence has been addressed via the axisymmetric spectrum (Kerschen & Gliebe Reference Kerschen and Gliebe1981), but this has never been incorporated into a theoretical prediction model. This paper combines these two turbulence modelling approaches to develop an anisotropic Gaussian decomposition model for upstream turbulence while addressing other considerations to mathematically and physically model the turbulence. The literature is comparatively scarce on investigating the effects of anisotropic turbulence in the context of leading-edge noise, with the key research on this being Gea-Aguilera et al. (Reference Gea-Aguilera, Zhang, Chen, Gill and Node-Langlois2015, Reference Gea-Aguilera, Gill, Zhang, Chen and Node-Langlois2016); Gea-Aguilera, Gill & Zhang (Reference Gea-Aguilera, Gill and Zhang2017); Zhong & Zhang (Reference Zhong and Zhang2019) and Gea-Aguilera et al. (Reference Gea-Aguilera, Karve, Gill, Zhang and Angland2021). In these papers, it is found that modifying anisotropy has a noticeable effect on turbulence spectra due to the redistribution of turbulence kinetic energy. Therefore, by modelling this, we may better understand the physical phenomena causing changes in aerofoil noise. In practice, turbulence is rarely isotropic, so accounting for anisotropy reflects the physical system of interest. Although the previous authors have used anisotropic Gaussian kernels to generate synthetic anisotropic turbulence, this approach did not consider an analytical weighting function. This was developed in a prior study by the authors (Hales et al. Reference Hales, Ayton, Kisler, Mahgoub, Jiang, Dixon, de Silva, Moreau and Doolan2022). Therefore, in this paper, we further develop the model to represent the underlying physics of the problem better. Experimentally, we observed that the turbulent behaviour varied significantly between the wall-normal direction and the streamwise and spanwise directions; however, the differences in the features of the anisotropic turbulence were far less significant between the streamwise and spanwise directions.
We will take a different approach in our turbulence modelling and assume the ‘axial’ direction (using terminology from the stated papers) is now the wall-normal direction, whereas the ‘transverse’ direction consists of the streamwise and spanwise plane containing the aerofoil. We will then apply the method of Gaussian decomposition introduced by Wohlbrandt et al. (Reference Wohlbrandt, Hu, Guérin and Ewert2016). However, we will adapt this to axisymmetric turbulence by considering anisotropic kernels and verifying the analytical representations. With this description of the turbulence, we can analyse the impact of our parameters representing the ratio of turbulence parameters. We discuss our methodology for obtaining data for the interaction of axisymmetric turbulence generated by a cylinder wake with a flat-plate model. With this data, we present an initial test of the Gaussian decomposition method by comparing various models for incident turbulence spectra with spectra calculated from experimental predictions. Finally, we verify our model against this data and outline more utility for our Gaussian decomposition, which can be included simply in our overall PSD model and will demonstrate an even better agreement with the experimental noise data.
2. Review of leading-edge noise
We first review the interaction of a general turbulent disturbance with a thin leading edge, as depicted in figure 1 and solved by Amiet (Reference Amiet1975, Reference Amiet1976), under the assumption that nonlinearities such as turbulence–turbulence interactions can be neglected and the mean flow Mach number is small. We invoke rapid distortion theory (RDT) to express the total fluid velocity as
\begin{equation} \boldsymbol{{u}}=\boldsymbol{{U}}_{c}+\boldsymbol{\nabla}\phi_s + \boldsymbol{{u}}^{(I)}, \end{equation}
where 
$\boldsymbol { {U}}_{c}$ is the steady mean flow velocity as it impinges upon the leading edge (the convection velocity), 
$\boldsymbol { {u}}^{(I)}$ denotes the turbulent distribution and 
$\phi _s$ is a velocity potential containing only scattered contributions from the leading edge. Without an edge 
$\phi _s=0$ (Goldstein Reference Goldstein1978).
Representation of the semi-infinite plate problem.
We can relate this velocity potential to pressure fluctuations via
\begin{equation} p_s={-}\rho_0\frac{{\rm D}\phi_s}{{\rm D}t}, \end{equation}where
\begin{equation} \frac{{\rm D}}{{\rm D}t}=\frac{\partial}{\partial t}+\boldsymbol{U_c}\boldsymbol{\cdot}\boldsymbol{\nabla}, \end{equation}
$\rho _0$ is the fluid density and 
$t$ denotes time.
Knowing 
$p_s$, we may calculate the PSD
\begin{equation} \varPsi(\omega)=\lim_{T\to\infty}\frac{\rm \pi}{T}p_s(\omega)p_s^*(\omega), \end{equation}
where 
$p_s^*$ denotes the complex conjugate, 
$\omega$ the angular frequency and 
$T$ is the time over which we measure our signal (Blake Reference Blake1970). This quantity describes how the power of the pressure jump is distributed across the frequency range when taking into account two separate signals and averaging over time. Our model will be of the form
\begin{equation} \varPsi(\omega,\theta)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} {g}(\boldsymbol{{k}},\theta)\varPhi_{22}(\boldsymbol{{k}})\delta\left(k_1-\frac{\omega}{U_{\infty}}\right) \text{d}k_1\,\text{d}k_2\,\text{d}k_3. \end{equation}
The details of deriving this are left to Appendix A. Our model (2.5) consists of two key components that need careful mathematical analysis. The turbulence spectrum 
$\varPhi _{22}(\boldsymbol { {k}})$ and the transfer function 
${g}$ that accounts for scattering effects. This function 
${g}$ can be considered a transfer function (Amiet Reference Amiet1976) and will be the focus of this section.
We will find 
${g}$ by analytically calculating the scattered field in § 2 before integrating this over the whole of wavenumber space using our spectrum 
$\varPhi _{22}$. While Amiet used Schwarzschild's solution, which will only give the solution on the surface, we will apply the Wiener–Hopf technique to obtain the solution in 
$\mathbb {R}^2$. This avoids propagating the surface pressure to the far field via an additional integral and can provide the far field nose directly.
2.1. Gust-scattering solution
For our simplified model, we will assume that the blade is a semi-infinite plate with streamwise variable 
$x$, the spanwise variable 
$z$ and the plate at 
$y=0$. Figure 1 demonstrates the physical problem. We will assume the plate is infinitely thin with an infinite span and will seek the scattered solution of a single frequency of turbulence. We aim to find the total far field scattered pressure occurring as our incoming gust impinges on the aerofoil, 
${g}(\boldsymbol { {k}},\theta ),$ for use in (2.5).
Given our disturbance convects within a uniform steady mean flow with Mach number 
$M={U_{c}}/{c_0}$, the three-dimensional convective wave equation for our scattered velocity potential 
$\phi _s$ is
\begin{equation} \left( \frac{{\rm D}^2}{{\rm D}t^2} - c_0^2\nabla^2 \right )\phi_s=0. \end{equation}Our incoming disturbance for a single gust is given by
\begin{equation} \phi_i=w_2\exp({{\rm i}\boldsymbol{{k}}\boldsymbol{\cdot}\boldsymbol{{x}}-{\rm i}\omega t}). \end{equation}
Here, we have intrinsically non-dimensionalised our length and velocity by a length scale 
$L$ relating to our experimental set-up and by the convection velocity 
$U_{c}$, the velocity at which turbulent eddies convect when impinging upon the leading edge. We assume our scattered wave is time harmonic with a specific exponential dependence in the spanwise direction; that is, we have a 
$\exp ({{\rm i}(k_3z-\omega t)})$ dependence on 
$\phi _s$. Furthermore, we assume our turbulence is frozen so that the acoustic wavenumber 
$k={\omega }/{c_0}$ can be replaced with the streamwise wavenumber 
$k_1$ via 
$k=Mk_1$ as per Taylor's hypothesis.
Using these assumptions, we can rearrange (2.6) and obtain a convected Helmholtz equation
\begin{equation} (1-M^2)\frac{\partial^2 \phi_s}{\partial x^2} +\frac{\partial^2 \phi_s}{\partial y^2} +2{\rm i}k_1M^2\frac{\partial \phi_s}{\partial x} +(M^2k_1^2-k_3^2)\phi_s=0. \end{equation}
We also require the pressure jump 
$\varDelta [p]=p(x,0^+)-p(x,0^-)$ for the total field upstream of the plate to be zero by continuity,
\begin{equation} \varDelta [\phi_s]= 0 \quad y=0,\ x<0, \end{equation}and zero through flow along the plate corresponding to
\begin{equation} \frac{\partial\phi_s}{\partial y}={-}{\rm i}k_2w_2\text{e} ^{{\rm i}k_1x} \quad y=0,\ x>0. \end{equation}
To eliminate the 
${\partial \phi _s}/{\partial x}$ term from our governing equation, we define 
$\beta =\sqrt {1-M^2}$ and use the convective transform
\begin{equation} \widetilde{\phi_s}(x,y)= \phi_s(x,y)\exp\left({\frac{{\rm i}k_1M^2x}{\beta^2}}\right), \end{equation}followed by the Prandtl–Glauert transformation
\begin{equation} \tilde{\tilde{\phi}}_s(x,y)=\widetilde{\phi_s}\left(\frac{x}{\beta},y\right). \end{equation}Applying these transformations, we obtain a standard Helmholtz equation and retain the Neumann boundary condition, but with slightly shifted wavenumbers
$$\begin{gather} \frac{{\partial^2 \tilde{\tilde{\phi}}_s}}{{\partial x^2}} + \frac{\partial^2 \tilde{\tilde{\phi}}_s}{\partial y^2} +{k^{*}}^2\tilde{\tilde{\phi}}_s=0, \end{gather}$$
$$\begin{gather}\varDelta{\left[\tilde{\tilde{\phi}}_s\right] }= 0,\quad y=0,\ x<0, \end{gather}$$
$$\begin{gather}\frac{\partial\tilde{\tilde{\phi}}_s}{\partial y}={-}{\rm i}k_2w_2\text{e} ^{{\rm i}\tilde{k}x}, \quad y=0,\ x>0, \end{gather}$$where we have denoted the ‘Prandtl–Glauert wavenumbers’ as
\begin{equation} \tilde{k}=\frac{k_1}{\beta^3}, \quad k^{*}=\frac{\sqrt{M^2k_1^2-\beta^2k_3^2}}{\beta}. \end{equation}
These wavenumbers account for our convective flow and our temporal and spanwise assumptions. We choose branch cuts for 
$k^*$ to ensure 
${\rm Im}(k^*)\geq 0$. For the rest of the paper, we drop the tildes and consider our problem in the new coordinates.
2.2. The Wiener–Hopf solution
Following the standard Wiener–Hopf procedure and the application of Liouville's theorem via the edge conditions (see Appendix B for more details), we obtain the solution to (2.13) as
\begin{equation} \phi_s(x,y)=\frac{-w_2k_2\,{\rm sgn}(y)}{2{\rm \pi}\gamma^-\left(-\tilde{k}\right)}\int_{-\infty}^{\infty}\frac{1}{\gamma^+(\alpha)\left(\alpha+\tilde{k}\right)} \exp({-{\rm i}\alpha x-\gamma|y|})\,\text{d}\alpha, \end{equation}
where 
$\gamma =\sqrt {\alpha ^2-{k^*}^2}$ with branch cuts taken so that 
$\gamma (0)=-{\rm i}k^*$, and 
$\gamma =\gamma ^+\gamma ^-$ is a multiplicative Wiener–Hopf factorisation of this function, that is, 
$\gamma ^+(\alpha )=\sqrt {\alpha +k^*}$ is analytic in the upper half-complex plane and 
$\gamma ^-(\alpha )=\sqrt {\alpha -k^*}$ is analytic in the lower half-complex plane.
To relate this to the far field pressure, we first perform the steepest descent approximation on our integral solution to find the velocity potential in the far field
\begin{equation} D(\boldsymbol{{k}},\theta)=\lim_{r\to\infty}{\sqrt{r}}{\phi_s(x,y,\boldsymbol{{k}})}, \end{equation}
where 
$r=\sqrt {x^2+y^2}$ is the radial distance from the source. Omitting the details for brevity, we find
\begin{equation} D(k_1,k_3,\theta)=\frac{\cos{\left(\dfrac{\theta}{2}\right)}}{\sqrt{k^*+\tilde{k}}(\tilde{k}-k^*\cos\theta)}. \end{equation}We then undo our convective transform to return our solution back to being in terms of pressure, that is we use the relation
\begin{align} p_s&={-}\rho_0\left(U_c\frac{\partial}{\partial x}-{\rm i}\omega\right)\phi_s\nonumber\\ &={-}\rho_0U_c\left(\frac{\partial}{\partial x}-{\rm i}k_1\right)\phi_s, \end{align}however, we disregard the dimensional constants that scale the result since they will be absorbed when we dimensionalise the PSD later in our results section.
After undoing our convective transform and applying (2.18), we obtain the transfer function that we denote by 
${g}$
\begin{align} {g}&=\frac{M^2k_2^2}{\rm \pi} \left|k_1\frac{1-2M^2}{1-M^2}+k^*\cos\theta\right|^2 \, \big|\text{e}^{{\rm i}k^*r}\big|^2 \, \big|D(k_1,k_3,\theta)\big|^2\nonumber\\ &=\frac{M^2k_2^2}{\rm \pi}\tilde{g}(k_1,k_3,\theta), \end{align}from which our theoretical model for the spectral sound pressure level (SPL) would be
\begin{equation} \varPsi(\omega,\theta)=\frac{M^2}{\rm \pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}k_2^2\, \tilde{g}\left(\frac{\omega}{U_c},k_3,\theta\right) \varPhi_{22}\left(\frac{\omega}{U_{c}}, k_2,k_3\right) \text{d}k_2 \,\text{d}k_3, \end{equation}
here, we are implicitly assuming our shifted wavenumbers 
$k^*$, 
$\tilde {k}$ are evaluated at 
$k_1={\omega }/{U_{c}}$ due to frozen convection, where 
$U_c$ is the convection velocity of eddies at the leading edge of the plate.
The 
$k_2$ dependence in (2.20) occurs only through 
$\varPhi _{22}(\boldsymbol { {k}})$, so we can reduce the double integral to a single integral by pre-integrating the spectrum.
We retain the term 
$|\text {e}^{{\rm i}k^*r}|$ because some 
$k^*$ may be complex. These subcritical gusts are known to have little (but potentially some) contribution to the overall sound; therefore, when evaluating (2.19), we pick 
$r=100$ and still include them (Roger & Moreau Reference Roger and Moreau2005). Since these subcritical gusts have 
${\rm Im}(k^*)>0$, the gusts act as evanescent waves in the far field.
3. Turbulence spectrum models
Having derived a solution to the transfer function 
$\tilde {g}$, we now need a fast and accurate turbulence spectrum 
$\varPhi _{22}(\boldsymbol { {k}})$ to account for the properties of the incoming turbulence.
3.1. The isotropic vertical velocity turbulence spectrum
To give an understanding of the turbulence itself, we will briefly review the fundamental statistical derivations behind the spectral functions. Many popular turbulence textbooks will cover this in detail (Batchelor Reference Batchelor1953; Tennekes & Lumley Reference Tennekes and Lumley1972; Hinze Reference Hinze1975). Appendix C of Grasso et al. (Reference Grasso, Jaiswal, Wu, Moreau and Roger2019) contains a succinct review within the context of noise reduction.
The correlation tensor 
$\boldsymbol{{R}}_{ij}$ is defined by
\begin{equation} \boldsymbol{{R}}_{ij}(\boldsymbol{{r}})=\frac{\overline{u_i(\boldsymbol{x},t)u_j(\boldsymbol{x}+\boldsymbol{r},t)}}{\overline{u_i (\boldsymbol{x},t)u_j(\boldsymbol{x},t)}}, \end{equation}which describes average spatial velocity fluctuations over the spatial domain.
We define the full three-dimensional energy spectral density, 
$\varPhi _{ij}(\boldsymbol {k})$ to be the full spatial Fourier transform of 
$\boldsymbol{{R}}_{ij}(\boldsymbol {x})$. This function splits the correlation function's behaviour across the wavenumbers we are concerned with, 
$\boldsymbol { {k}}$. Thus, we consider three-dimensional turbulence with our gust solutions when integrating our model over all these wavenumbers.
The two most popular turbulence models are the von Kármán model (von Kármán Reference von Kármán1948) and the Liepmann model (Liepmann, Laufer & Liepmann Reference Liepmann, Laufer and Liepmann1951), which differ by the assumed energy scaling of the inertial subrange and are derived from empirical results and dimensional arguments. An alternative model represents Hunt's RDT (Hunt Reference Hunt1973), in which the rapid distortion of turbulence is considered and studied, leading to a new 
$\varPhi _{22}$ with a much steeper decrease in energy within the inertial subrange.
All three significant models are encompassed by the generalised von Kármán model (Wilson Reference Wilson1997) which begins with the correlation function given as
\begin{equation} \boldsymbol{{R}}_{22}(r;p)=\frac{1}{2^{p-{7}/{2}}\varGamma(p-\frac{5}{2})} \left(\frac{r}{L}\right)^{p-{5}/{2}}{{\rm K}}_{p-{5}/{2}}\left(\frac{r}{L}\right), \end{equation}
where 
${{\rm K}}_{p-{5}/{2}}$ is the modified Bessel function of the second kind of order 
$p-\frac {5}{2}$ and 
$L$ is a characteristic length scale to be defined in terms of the integral length scale 
$\varLambda$. The von Kármán model corresponds to 
$p=17/6$ in (3.2), while the Liepmann model uses 
$p=3$ and the RDT model 
$p=11/3$. At this point, unless specified, the von Kármán model directly refers to the specific (
$p=\frac {17}{6}$) example of the umbrella term ‘generalised von Kármán model’.
We define the energy density spectrum to be
\begin{equation} {E}(\kappa)=\tfrac{1}{2}\int\varPhi_{ii}(\boldsymbol{{k}})\,\text{d}A(\kappa), \end{equation}
where we integrate over the surface of spheres of which 
$\text {d}A$ is an element and tells us the density of turbulence kinetic energy for each 
$\kappa =|\boldsymbol { {k}}|$ (Batchelor Reference Batchelor1953).
For isotropic turbulence, we obtain the energy density spectrum for the generalised von Kármán model to be
\begin{equation} E(\kappa;p)=2k_T L(p) C(p) \frac{(\kappa L(p))^4}{(1+(\kappa L(p))^2)^p}, \end{equation}
where 
$C(p)={\varGamma (p)}/{\varGamma (\frac {5}{2}) \varGamma (p-\frac {5}{2})}$ is a normalising constant that ensures the total integral over 
$|\boldsymbol { {k}}|$ is equal to the turbulent kinetic energy 
${k_T}$.
The turbulent kinetic energy is defined as
\begin{equation} k_T=\tfrac{1}{2}(u_{rms}^2+v^2_{rms}+w^2_{rms}) \end{equation}
where we use the notation 
$(u_{rms},v_{rms},w_{rms})$ to be the root mean square (r.m.s.) velocity in each corresponding axial direction. For isotropic turbulence, 
$u=v=w$ and so we can set 
$k_T={3u_{rms}^2}/{2}$.
Using Batchelor's formula relating the energy density spectrum to 
$\varPhi _{ij}$, which holds specifically for the case of isotropic turbulence, we find that the spectrum required for the PSD is
\begin{equation} \varPhi_{22}(\boldsymbol{{k}};p)=\frac{E(\kappa;p)}{4{\rm \pi}\kappa^4} (k_1^2+k_3^2). \end{equation}
The integral length scale 
$\varLambda$ is defined via
\begin{align} \varLambda&=\int_0^{\infty}\boldsymbol{{R}}_{11}(x)\,\text{d}\kern0.7pt x\nonumber\\ &=\frac{3}{2k_T}\int_0^{\infty}\left( \int_{\mathbb{R}^3}\varPhi_{11}(\boldsymbol{{k}}) \text{e}^{-{\rm i}k_1x}\,\text{d}\boldsymbol{{k}} \right), \end{align}and the characteristic length scale used in (3.4) is given by (Durbin & Petterson Reference Durbin and Petterson2001)
\begin{equation} L(p)=\frac{\varLambda \varGamma(p-\frac{5}{2})}{\sqrt{\rm \pi}\varGamma(p-2)}=\varLambda L^*(p). \end{equation}
Putting this together gives the generalised 
$\varPhi _{22}$ for isotropic turbulence
\begin{equation} \varPhi_{22}(\boldsymbol{{k}};p)=\frac{3u_{rms}^2CL^5}{4{\rm \pi}}\frac{k_1^2+k_3^2}{\left(1+L^2|\boldsymbol{{k}}|^2\right)^p}. \end{equation}We compare the three models alongside the standard Gaussian model
\begin{equation} \varPhi_{22}^G(\boldsymbol{{k}})=\frac{u_{rms}^2L^5}{{\rm \pi}^4}(k_1^2+k_3^2)\exp{\left(-\frac{L^2\left|\boldsymbol{{k}}\right|^2}{\rm \pi}\right)}, \end{equation}
in figure 2, where we plot the one dimensional turbulence spectrum 
$\varTheta _{22}(k_1)$, given by
\begin{equation} \varTheta_{22}(k_1;p)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\varPhi_{22}(\boldsymbol{{k}};p)\,\text{d}k_2\, \text{d}k_3. \end{equation}Comparisons of the 
$\varTheta _{22}(k_1)$ model (3.11) while using the generalised von Kármán model (3.9) and varying 
$p$ alongside the standard Gaussian model (3.10). Here, we use 
$\varLambda =0.0133\ \text {m}$, 
$u_{rms}=3.92\ \text {m}\ \text {s}^{-1}$ as the length scale and r.m.s. velocity used for each model.
3.2. An axisymmetric anisotropic model
We now turn to extend these models to anisotropic turbulence. A useful model for anisotropic turbulence, more specifically axisymmetric turbulence, is presented in Kerschen & Gliebe (Reference Kerschen and Gliebe1981) for modelling turbulent jet noise
\begin{align} \varPhi_{ij}(\boldsymbol{{k}})&=[k^2\delta_{ij}-k_ik_j]F+[(k^2-(k_m\lambda_m)^2)\delta_{ij}\nonumber\\ &\quad -k_ik_j-k^2\lambda_i\lambda_j +k_m\lambda_m(\lambda_ik_j+k_i\lambda_j)]G, \end{align}
where the unit vector 
$\boldsymbol {\lambda }$ gives the direction of the mean flow. The functions 
$F$ and 
$G$ are given by
\begin{equation} F=\frac{2l_a l_t^{4}u_a^2} {{\rm \pi}^{2}[1+l_a^2k_1^2+l_t^2(k_2^2+k_3^2)]^{{3}/{2}}} \quad G=\left(2\frac{u_t^2}{u_a^2}-\frac{l_t^2}{l_a^2}-1 \right)F. \end{equation}
The subscripts 
$a$, 
$t$, distinguish parameters based on axial and transverse directions, chosen to represent the anisotropic behaviour in different spatial directions. Here, the 
$l_{a,t}$ values refer to the length scale equivalents to 
$L$ and 
$u_{a,t}$ values refer to velocities equivalent to 
$u$ similar to the isotropic case.
Kerschen and Gliebe's model was explicitly derived for jet turbulence (Kerschen & Gliebe Reference Kerschen and Gliebe1981), thereby the axial direction is taken to be the direction of the mean flow, while the transverse direction is that of the plane normal to this, where it is assumed the turbulence statistics are spatially uniform.
This model was adapted to leading-edge noise by Gea-Aguilera et al., in which the four parameters 
$l_{a,t}$, 
$u_{a,t}$ are discussed in depth and formally defined (Gea-Aguilera et al. Reference Gea-Aguilera, Karve, Gill, Zhang and Angland2021). The length scales are equal to integral length scales in the contextual direction. In contrast, the axial velocity is set to be the r.m.s. velocity in the direction of the mean flow. However, this paper notes that the ‘transverse velocity’ is defined arbitrarily in the original literature. It is defined as some velocity 
$u_t$ required to satisfy the condition
\begin{equation} 2\frac{u_t^2}{u_a^2}-\frac{l_t^2}{l_a^2}\geq0. \end{equation}
This inequality is forced in Kerschen & Gliebe (Reference Kerschen and Gliebe1981) to ensure that the velocity correlation tensor 
$\boldsymbol{{R}}_{ij}$ is non-negative.
Following these definitions, the integral length scale conditions would require equality between the two velocities within the transverse plane.
Our experimental study can be represented by design as axisymmetric around the plane perpendicular to the cylinder. Thus, we take the axial direction 
$\boldsymbol {\lambda }$ as the 
$y$-direction (wall normal), while we assume the plane containing our plate (
$x,z$)-plane is the transverse direction. For completeness, the three wavenumber-frequency spectra 
$\varPhi _{ii}, \ i=1,2,3$ can be written as
\begin{equation} \left.\begin{gathered} \varPhi_{11}(\boldsymbol{{k}})=\left(k_2^2+\left(2u_r^2-l_r^2\right)k_3^2\right)F\\ \varPhi_{22}(\boldsymbol{{k}})=\left(k_1^2+k_3^2\right)F\\ \varPhi_{33}(\boldsymbol{{k}})=\left(k_2^2+\left(2u_r^2-l_r^2\right)k_1^2\right)F \end{gathered}\right\},\end{equation}
where 
$F$ is given in (3.13a,b) and 
$u_r={u_t}/{u_a}$, 
$l_r={l_t}/{l_a}$.
We now formally define
\begin{equation} u_t^2=\frac{u_{rms}^2+w_{rms}^2}{2}, \end{equation}which satisfies the inequality (3.14) and ensures our model makes physical sense.
In individual measurements during a separate set of experiments, we used stereo-particle image velocimetry (PIV) to obtain three velocity components in a cylinder wake. Data are obtained at Reynolds numbers (based on cylinder diameter) 
$Re_D=20\times 10^3$ and 
$30\times 10^3$ to and downstream locations 
$x/D$ of 5 and 14. The experimental set-up and measurement techniques are described in Dixon et al. (Reference Dixon, Jiang, de Silva, Moreau and Doolan2022). The PIV data provide the turbulence characteristics of the incoming anisotropic turbulent flow in the 
$x$ and 
$z$ directions, which we need as input for our mathematical model. The results at different Reynolds numbers and downstream locations both indicated that 
$u_{rms}\approx 1.2 w_{rms}$, which we would expect; therefore, we will make the simplification 
$u_{rms}=w_{rms}$ so that 
$u_t=u_{rms}$ and 
$u_a=v_{rms}$.
This approach is inspired by the contextual design of our experiment, which was chosen to ensure significant ratios can be tested. This varies from the numerical study (Gea-Aguilera et al. Reference Gea-Aguilera, Karve, Gill, Zhang and Angland2021) in which the ratios are inverted due to the axial direction being the 
$x$-axis.
Implementing these initial changes into the axisymmetric framework gives the model
\begin{equation} \varPhi_{22}^a(\boldsymbol{{k}};p)= \frac{u_a^2l_t^4l_a\varGamma(p)}{{\rm \pi}^{3/2}\varGamma\left(p-\frac{5}{2}\right)} \frac{k_1^2+ k_3^2}{[1+l_t^2k_1^2+l_a^2k_2^2+l_t^2k_3^2]^p}. \end{equation}
Here, we have normalised 
$\varPhi _{22}$ for the turbulence kinetic energy requirement for anisotropic spectra. Note that, in doing so, we will be using 
$\sigma ={2}/{u_r^2}-{l_t^2}/{l_a^2}$, to be compared with the form 
$\sigma =2u_r^2-\varLambda _r^2$ that would be an equivalent definition for the model in Kerschen & Gliebe (Reference Kerschen and Gliebe1981) and Hales et al. (Reference Hales, Ayton, Kisler, Mahgoub, Jiang, Dixon, de Silva, Moreau and Doolan2022). The equation that our normalised model must satisfy is
\begin{equation} \int_{\mathbb{R}^3} \varPhi_{ii}(\boldsymbol{{k}};p)\,\text{d}\boldsymbol{{k}}=2k_T=u_a^2(2+u_r^2), \end{equation}
where the sum 
$\varPhi _ii$ can be found using (3.15). It remains to relate our length scale parameters 
$l_{a,t}$ to turbulent integral length scales that may be measured in experiments. We use the two equations
\begin{equation} \varLambda_{11}^{(1)}=\frac{\rm \pi}{u_{{ rms}}^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\varPhi_{11}(k_1=0,k_2,k_3;p)\,\text{d}k_2\,\text{d}k_3 \end{equation}and
\begin{equation} \varLambda_{22}^{(2)}=\frac{\rm \pi}{v_{{ rms}}^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\varPhi_{22}(k_1,k_2=0,k_3;p)\,\text{d}k_1\,\text{d}k_3. \end{equation}
Note that for simplicity we define 
$\varLambda _1\equiv \varLambda _{11}^{(1)}$ and 
$\varLambda _r\equiv {\varLambda _{22}^{(2)}}/{\varLambda _{11}^{(1)}}$, so that the isotropic limit is 
$\varLambda _r\to 1$. Solving these to find our parameters 
$l_{a,t}$ in terms of 
$\varLambda _{1,r}$ gives
\begin{equation} l_a=\varLambda_1\varLambda_rL^*(p), \quad l_t=\varLambda_1L^*(p), \end{equation}
with 
$L^*$ as in (3.8), while the condition (3.14) now becomes
\begin{equation} \frac{2}{u_r^2}-\frac{1}{\varLambda_r^2}\geqslant 0, \end{equation}which will be reasonable for all cases of anisotropy studied experimentally in the paper.
Putting all this together, we obtain the final anisotropic spectrum
\begin{equation} \varPhi_{22}^a(\boldsymbol{{k}};p)= \frac{u_a^2(L^*(p)\varLambda_1)^5\varLambda_r\varGamma(p)}{{\rm \pi}^{3/2}\varGamma\left(p-\frac{5}{2}\right)} \frac{k_1^2+ k_3^2}{[1+{L^*}^2(p)\varLambda_1^2\left(k_1^2+\varLambda_r^2k_2^2+k_3^2\right)]^p}. \end{equation}As the introduction mentions, wall-normal velocity fluctuations generate pressure fluctuations that can be attributed to the sound we hear. Our new approach, in which anisotropic spectra are combined with an analytic transfer function as in (2.20), indicates how leading-edge noise is generated physically. We aim to create a clearer model, adapted to the context of leading-edge noise, in which the parameters of interest now relate to the direction of flow for consistency with previous literature; however, all ratios compare the wall-normal with the streamwise direction, and it is these ratios we aim to understand.
3.3. Approximation by Gaussian decomposition
Evaluating (3.23) is often impractical. Any subsequent integral it appears in may not have closed forms or may be challenging to compute efficiently and accurately. An example of the former issue is 
${E}(\kappa )$, which is discussed in Appendix D. We can thus use a Gaussian transform method to approximate turbulence spectra with Gaussian kernel filters (Alecu, Voloshynovskiy & Pun Reference Alecu, Voloshynovskiy and Pun2005; Wohlbrandt et al. Reference Wohlbrandt, Hu, Guérin and Ewert2016). This paper will use the terminology Gaussian decomposition when referring to this model.
To the authors’ knowledge, the method of Gaussian decomposition has been investigated solely for isotropic turbulence, such as the models in § 3.1. Here, we extend the model to the anisotropic case with two length scales 
$l_a$ and 
$l_t$.
The idea of Gaussian decomposition in the anisotropic case is to find an analytical weighting function 
$f$ that satisfies the integral equation
\begin{equation} \varPhi_{22}^a(\boldsymbol{k};p) =\int_0^\infty\int_0^\infty f(l_a,l_t)\tilde{\varPhi}_{22}^{a,G}(\boldsymbol{{k}})\,\text{d}l_t\,\text{d}l_a, \end{equation}
where 
$\tilde {\varPhi }_{22}^{a,G}$ is the anisotropic Gaussian kernel function
\begin{equation} \tilde{\varPhi}_{22}^{a,G}(\boldsymbol{k})=\frac{l_a^3l_t^4u_t^2(2+u_r^2)(k_1^2+ k_3^2)}{{\rm \pi}^4\left((1+\sigma)l_a^2+l_t^2\right)}\exp\left(-\frac{l_a^2k_2^2}{\rm \pi}\right) \exp\left({-\frac{l_t^2(k_1^2+k_3^2)}{\rm \pi}}\right), \end{equation}
which is a turbulence spectrum in its own right but scaled depending on the arbitrary length scales 
$l_a$, 
$l_t$, which we integrate over.
After finding 
$f$, we approximate the resulting integral equation by summation
\begin{equation} \varPhi_{22}^a(\boldsymbol{{k}};p)\approx\sum_{m=0}^Nf(p;l_m)\Delta l_m\varPhi_{22}^{a,G}(\boldsymbol{{k}};p;l_m), \end{equation}
where the 
$l_m$ values are chosen to be within a sensible range of length scales representing the flow itself, while the 
$\Delta l_m$ values are weightings distributed with the trapezoid rule between some scalar multiple of the largest scale used in our summation and the smallest (Wohlbrandt et al. Reference Wohlbrandt, Hu, Guérin and Ewert2016). The former length scale will be the integral length scale of turbulence, while the latter is the Kolmogorov length scale
\begin{equation} l_k=\left(\frac{\nu^3}{\epsilon} \right)^{{1}/{4}}, \end{equation}
with 
$\nu$ the dynamic viscosity and 
$\epsilon$ the dissipation rate. This ensures we capture the lengths of the largest and smallest scales of turbulence that contribute to our far field noise.
To find 
$f$, we perform similar changes of variables in each direction to Wohlbrandt et al. (Reference Wohlbrandt, Hu, Guérin and Ewert2016). Details of this are found in Appendix B and reduce the right-hand side to two Laplace transforms. Taking the inverse transforms, in turn, gives us
\begin{equation} \left.\begin{gathered} f(p;l_a,l_t)=\frac{4{\rm \pi} u_r^2(L\varLambda_1)^5\varLambda_r}{(2+u_r^2)\varGamma \left(p-\frac{5}{2}\right)}\tilde{f}(l_a)\tilde{g}(l_a,l_t),\\ \tilde{f}(p;l_a)=\left(\frac{l_a^2}{{\rm \pi}(L\varLambda_r\varLambda_1)^2}\right)^p\exp \left(-\frac{l_a^2}{{\rm \pi}(L\varLambda_r\varLambda_1)^2}\right),\\ \tilde{g}(l_a,l_t;p)=\frac{(1+\sigma)l_a^2+l_t^2}{l_a^2l_t^3}\delta\left(l_t^2-\frac{l_a^2}{\varLambda_r^2}\right). \end{gathered}\right\} \end{equation}
However, we note that the 
$l_t$ integral of 
$g$ in (3.24) may be evaluated immediately to obtain our spectrum as an approximation dependent on only one length scale, 
$l$
\begin{equation} {\varPhi}_{22}^a(\boldsymbol{{k}};p)=\int_0^\infty f(p;l)\varPhi_{22}^{a,G}(\boldsymbol{{k}})\,\text{d}l, \end{equation}with
\begin{equation} \varPhi_{22}^{a,G}(\boldsymbol{{k}})=\frac{u_a^2l^5\varLambda_r(k_1^2+k_3^2)}{{\rm \pi}^4} \exp\left(-\frac{l^2(k_1^2+\varLambda_r^2k_2^2+k_3^2)}{\rm \pi}\right). \end{equation}After this integration, we obtain the weighting function
\begin{equation} \hat{f}(p;l)=\frac{2}{l\varGamma\left(p-\frac{5}{2}\right)}\left(\frac{l}{\sqrt{\rm \pi} \varLambda_1L}\right)^{2p-5}\exp\left({-\frac{l^2}{{\rm \pi}{L}^2\varLambda_1^2}}\right). \end{equation}Therefore, our normalised approximated anisotropic spectrum can be written in the form
\begin{equation} \overline{\varPhi_{22}}(\boldsymbol{{k}};p)=\sum_{m=0}^{M} \frac{\hat{f}(p;l_m)\Delta l_m}{\sum_{m=0}^{M}\hat{f}(p;l_m)\Delta l_m}\varPhi_{22}^{a,G}(\boldsymbol{{k}};l_m). \end{equation}
An interesting observation is that there are no anisotropy parameters in our weighting function 
$f$; thus, the approximation method generalises from isotropic to anisotropic turbulence solely by using anisotropic Gaussian kernels with the same turbulent properties instead. We shall discuss this feature in more detail later in § 5.
4. Experimental methods
This section introduces the experimental methodology for collecting turbulence and noise data to validate the mathematical model. Two experiments were conducted. The first characterised a cylinder wake as a source of anisotropic turbulence in an open jet anechoic wind tunnel. The second measured the acoustic radiation generated by a flat plate placed in the turbulent wake of the cylinder. The first set of experiments used PIV to determine the characteristics of the anisotropic turbulent inflow. In the second set of experiments, the radiated sound from the leading edge of the flat plate was captured with a phased microphone array. All experiments were performed in the UNSW anechoic wind tunnel (UAT). The two sets of experiments are used to validate the mathematical model. The PIV data can provide turbulence characteristics, which are needed as inputs to the model (3.32). Moreover, the far field noise measurement taken at the same flow conditions will be used to validate the far field noise predictions from the mathematical model. The UAT consists of an open jet placed inside an anechoic room, as illustrated in figure 3. The test section is in the potential core of the open jet. The cross-section of the outlet is 
$0.455 \times 0.455\ \textrm {m}^2$. The overall volume of the anechoic chamber is 
$3 \times 4.17 \times 2.15\ \textrm {m}^3$. The UAT has a free-stream turbulence intensity of 0.7 % at a mean flow speed of 
$20\ \textrm {m}\ \textrm {s}^{-1}$. Further details about the facility are available in Doolan et al. (Reference Doolan, Moreau, Awasthi and Jiang2019). For the acoustic measurements, a flat-plate aerofoil with a chord length of 352 mm was used. The geometry of the test model is illustrated in a schematic in figure 4. The airfoil's trailing edge is deliberately serrated to mitigate trailing-edge noise, thus allowing for better separation of noise sources when conducting acoustic beamforming, especially at lower frequencies. Two Reynolds numbers (based on chord length) were investigated: 
$Re_c=4.65\times 10^5$ and 
$Re_c=6.52\times 10^5$. The flat plate model was positioned at an angle of attack of 
$0^{\circ }$ in all cases. For the generation of anisotropic turbulence, a cylinder was placed upstream of the flat plate parallel to its leading edge as depicted in figure 5. Six different distances between the centre of the cylinder and the plate's leading edge (
$\Delta x = 8{D}, 9.5{D}, 11{D}, 12.5{D}, 14{D}$ and 
$16D$, where 
$D$ is the cylinder diameter) were investigated separately. In addition, baseline measurements were performed without the cylinder in place. Acoustic measurements were performed with a 64 channel microphone array described in § 4.2.
Schematic of UNSW anechoic wind tunnel.
Schematic of the flat-plate aerofoil test model. The spatial dimensions shown are in millimetres. Figure adapted and used, with the permission of the authors, from Ayton et al. (Reference Ayton, Karapiperis, Awasthi, Moreau and Doolan2021b). LE: leading edge; TE: trailing edge.
Schematic for (a) the planar PIV measurement for the cylinder wake and (b) the measurements of radiated leading-edge noise in an anisotropic turbulent flow.
4.1. Particle image velocimetry
The PIV technique characterised the anisotropic turbulent flow field. The PIV system used is a planar PIV system that can obtain two-dimensional velocity components. In this set of experiments, the velocities measured are 
$u$ and 
$v$, representing the 
$x$ and 
$y$ velocity components. A cylinder of diameter 
$D=22$ mm was placed 
$6D$ upstream of the measurement field of view (FOV) to generate anisotropic turbulence in its wake region. The FOV for the PIV measurement is on an (
$x,y$) plane and has dimensions of 
$286\ \textrm {mm}\times 143\ \textrm {mm}$ in the streamwise (
$x$) and vertical (
$y$) directions, respectively. The origin of the coordinate system is located at the centre of the cylinder. The measurement FOV covers a streamwise distance between 
$\Delta x = 6D$ and 
$19D$ and a vertical distance between 
$\Delta y = - 3.25D$ and 
$3.25D$. Figure 5(a) shows a schematic of the planar PIV experiment. For the PIV measurement, a set of 3000 double-frame images was obtained for each case. A double-pulse Litron Nd:YAG PIV laser generated a laser sheet to illuminate the tracer particles. The lasers have a wavelength of 532 nm and an output energy of 
$250\ \textrm {mJ}\ \textrm {pulse}^{-1}$. Tracer particles were generated using a diethylhexyl sebacate fluid and a LaVision aerosol generator with a mean diameter of 
$1\ \mathrm {\mu } \textrm {m}$. A PCO.panda 26 DS CMOS camera with a maximum resolution of 
$5120\times 5120$ pixels was used to acquire images at a frequency of 5 Hz. The laser and camera synchronisation was controlled using a LaVision Programmable Timing Unit (PTU X). The software used for PIV data processing was Davis 8.4. The raw images have a 
$5120\times 2560$ pixel resolution and were processed using final interrogation windows of 
$32\times 32$ pixels with 50 % overlap. Therefore, each instantaneous velocity field had 
$322\times 161$ vectors. Spurious vectors were removed from the resulting vector fields using a median test during post-processing (Westerweel & Scarano Reference Westerweel and Scarano2005). The uncertainty of the PIV displacements is approximately 0.1 pixels (Adrian & Westerweel Reference Adrian and Westerweel2011), which corresponds to 
$\pm 5\ \mathrm {\mu }\textrm {m}$ in the present work.
For the present study, the anisotropic turbulent flow field is characterised by its flow statistics and the integral length scale. The Reynolds’ stresses can provide turbulence kinetic energy and dissipation rates. The turbulence dissipation rate obtained from two-dimensional planar PIV measurements was estimated using the following equation (Wang et al. Reference Wang, Yang, Wu, Ma, Peng, Liu and Wang2021):
\begin{equation} \varepsilon=\nu\left(\overline{2\left(\frac{\partial u}{\partial x}\right)^2} + \overline{2\left(\frac{\partial v}{\partial y}\right)^2} + \overline{3\left(\frac{\partial v}{\partial x}\right)^2} + \overline{3\left(\frac{\partial u}{\partial y}\right)^2} + \overline{2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}}\right). \end{equation}
Figure 6 shows the variation of the velocity r.m.s. values and the dissipation rate in the streamwise (
$x$) direction. The anisotropy in the cylinder wake can be quantified by calculating the ratio between the fluctuating streamwise (
$x$), 
$u$, and vertical (
$y$), 
$v$, velocity components. Note that 
$u$ and 
$v$ in experiments correspond to 
$u_a$ and 
$u_t$ in the analytical model, respectively. The longitudinal integral length scales in the streamwise (
$\varLambda _{11}^{(1)}$) and vertical (
$\varLambda _{22}^{(2)}$) directions can be obtained by integrating the spatial auto-correlation function 
$\boldsymbol{{R}}_{ii}$ for the corresponding fluctuating velocity component (Pope Reference Pope2000)
\begin{equation} \varLambda_{ii}^{(i)} = \int_{0}^{\infty} \boldsymbol{{R}}_{ii} \,\text{d} r_i, \end{equation}with
\begin{equation} \boldsymbol{{R}}_{ii}(r_i) = \frac{\overline{u_i(\boldsymbol{x},t)u_i(\boldsymbol{x}+r_i,t)}}{\overline{u_i(\boldsymbol{x},t)^2}}, \end{equation}
where 
$\boldsymbol {x}$ denotes the location of the velocity vector, 
$u_i$ is the fluctuating velocity in the 
$i$-direction and 
$r_i$ is the distance between two locations in the 
$i$-direction. The auto-correlation function for each 
$\Delta x$ location was calculated using the data in the range 
$\pm 2D$. The integral length scales were then obtained by integrating the Laplacian fitting of the auto-correlation function. Tables 1 and 2 summarise the integral length scale results.
Variation of the (a) r.m.s. of the fluctuating velocity and (b) dissipation rate along the centreline (
$y=0$ mm) of the cylinder wake.
Model parameters for 
$U_{\infty }=20\ \textrm {m}\ \textrm {s}^{-1}$.
Model parameters for 
$U_\infty =28\ \textrm {m}\ \textrm {s}^{-1}$.
4.2. Acoustic measurement
A 64 channel array of GRAS type 40 PH 1/4 inch microphones was used to measure noise in the UAT (figure 5b). The uncertainty of the sound pressure measurement is 
$\pm 1$ dB over 0.01 to 20 kHz. The microphone array samples acoustic data simultaneously on each channel (microphone). Subsequent signal processing uses these data to produce a ‘beamforming map’ over a ‘scanning plane’. Here, the scanning plane is the plane coincident with the flat-plate aerofoil. The beamforming map is a spatial representation of acoustic source strength on the scanning plane. It is produced during data processing by altering the phase on each channel, assuming point sources are at discrete points on the scanning plane. Sources are identified when these phase shifts reinforce, creating local maxima in the beamforming map.
The microphones are distributed in a spiral shape on the array. The chosen spatial distribution of the microphones provides an optimised localisation and quantification of sound sources (Prime, Doolan & Zajamsek Reference Prime, Doolan and Zajamsek2014). The array plane is located 1.15 m from the leading edge to ensure it is in the acoustic and geometric far field relative to the sound sources on the flat plate's leading edge. The acoustic signals were analysed in a frequency range of 0.2–10 kHz. The array plane is parallel to the flat plate, and the array centre microphone is aligned with the centre of the flat plate's leading edge. The recorded acoustic signals are sampled at 65 536 Hz. The time signals are Fourier transformed into cross-spectral matrices using Welch's overlapped segment averaging estimator. The overlap is 50 %, the averaging block size is 8192, and the Hann window is used, resulting in a frequency resolution of 8 Hz. Complete details of the beamforming algorithm can be found in Sarradj (Reference Sarradj2012).
To obtain an acoustic spectrum for leading-edge noise that is not contaminated with extraneous noise, such as noise emanating from the junction of the flat plate and the top and bottom plate, beamforming results are spatially integrated to contain sources from the leading-edge region only. To account for the influence of the point spread function (PSF) of the microphone array, the source power integration method is applied (Brooks & Humphreys Reference Brooks and Humphreys1999), which significantly reduces the effect of the PSF in its contaminating contribution to the overall integrated frequency spectrum of the sound pressure. The integration process first sums the source power estimates of the scanning grids in the integration region. Then it scales the result to the sound pressure spectrum obtained from the array centre microphone for an equivalent monopole source in the centre of all scanning grids.
The spatial integration region is rectangular and centred on the flat plate's leading edge, spanning over 0.2 m. The streamwise length of the integration region varies with frequency to ensure that acoustic sources within a 6 dB dynamic range are considered. Figure 7 shows a sample beamforming map for the case at 
${U_\infty = 28\ \textrm {m}\ \textrm {s}^{-1}}$ and 
$\Delta x = 12.5D$, where the dominant noise source is located at the airfoil leading edge and all noise sources within the 6 dB dynamic range are included in the integration region. The final frequency spectrum obtained from this post-processing sequence is presented in 
$1/12\textrm {th}$ octave bands and represents the overall integrated spectral density of the sound pressure in 
$\textrm {Pa}^2\ \textrm {Hz}^{-1}$.
Beamforming map at 2 kHz for the test case at 
${U_\infty =28\ \textrm {m}\ \textrm {s}^{-1}}$ and 
$\Delta x = 12.5D$. The locations of the array microphones, upstream cylinder and airfoil are indicated by grey circles, a grey box and a green box, respectively. A blue dashed box marks the integration region.
5. Results
We have derived both components of our analytical model (2.5) and discussed the methodology behind an experimental campaign to validate this model; therefore, we present our findings comparing the two. Our focus will be on the benefits of Gaussian decomposition as a method to better represent essential features of our experimental set-up, whilst highlighting the versatility and efficiency of the method. Tables 1 and 2 give the empirical constants used in our model obtained experimentally as described in the previous section. The value of 
$\Delta x$ represents the streamwise location of the leading edge at which the parameters in the other columns are determined (see figure 4). Assuming Taylor's frozen turbulence hypothesis, the convection velocity 
$U_c$ is set as the local mean flow speed in the vicinity of the leading edge. An experimental investigation into the convection velocity by Lin & Hsieh (Reference Lin and Hsieh2003) determined 
$U_c$ to converge around 
$0.9 U_\infty$ at 
$\Delta x = 10D$. However, their results showed a moderate Reynolds number dependency and were obtained for lower Reynolds numbers than those investigated in the current study. The value of 
$u_{rms}$ depicts the root mean square of the streamwise turbulent velocity fluctuations, while 
$u_r$ corresponds to the ratio between 
$v_{rms}$ and 
$u_{rms}$, obtained via PIV. Similarly, 
$\varLambda _1 = \varLambda ^{(1)}_{11}$ is the streamwise integral length scale, and 
$\varLambda _r$ represents the ratio between the integral length scales 
$\varLambda ^{(2)}_{22}$ and 
$\varLambda ^{(1)}_{11}$. They are determined by applying (4.2) and (4.3) to the fluctuating velocity vectors from the PIV measurements. Those are also used to calculate the Kolmogorov length scale 
$l_k$ using (4.1).
5.1. Modelling anisotropic turbulence
The Gaussian decomposition model is beneficial for a variety of reasons. First, we find that a small number of individual kernels are required to give an approximation for the total PSD, 
$\varPsi (\omega,\theta )$, in which the axisymmetric von Kármán type model is used (2.20), as seen in figure 8. Second, there are both numerical and analytical benefits to using exponential kernels instead of the standard von Kármán type models. An example of this is shown explicitly in Appendix D for the energy density spectrum 
$E(\kappa )$.
Convergence of the Gaussian decomposition method when used to approximate the leading-edge noise within our model 
$\mathcal {P}(f)$ (5.8).
Often, having exponentials in our spectral functions will allow for easily calculable closed forms after integration, such as when spectral functions such as 
$\varPhi _{ij}$ are integrated to obtain the spherically symmetric energy density function (D5). Although the details are omitted for brevity, one can derive that the weightings required for a Gaussian decomposition model of a spectral function like 
$E$ will also be the same, and will converge just as quickly. Therefore, we can propose that it may be possible to fit a Gaussian model to a spectral representation of the data and then use this within the scattering model.
More specifically, in the Gaussian decomposition model, we are summing 
$N$ kernels with length scales distributed within the range in which we expect our eddies to carry almost all the turbulent kinetic energy. However, instead of decomposing a typical von Kármán type model into weighted Gaussians, we explore summing up Gaussians that are weighted to better match a given physical example of incoming turbulence. The generalised von Kármán model has a parameter 
$p$, which is used to tune the decay rate in the inertial subrange. Thus, we take a new approach to reflect better the differing power laws in different regions of 
$k_1$ to add kernels that would be used to approximate the corresponding von Kármán type model with 
$p$ scaling.
As an initial test, we will look at model comparisons for the spectral quantity
\begin{equation} \varTheta_{22}(k_1)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\varPhi_{22}(k_1,k_2,k_3)\,\text{d}k_2\,\text{d}k_3. \end{equation}
We will compare the incident PIV turbulent measurements, converted to spectral quantities using the fast Fourier transform, with the isotropic von Kármán model, the anisotropic von Kármán model and the anisotropic Gaussian decomposition model with 
$N=35$. Using the notation
\begin{equation} k_e=\frac{\sqrt{\rm \pi}\varGamma\left(p-2\right)}{\varLambda_1\varGamma\left(p-\frac{5}{2}\right)} \quad \hat{k}_1=\frac{k_1}{k_e}. \end{equation}The following formulae are used:
$$\begin{gather} \varTheta^{VK}_{22}(k_1)=\frac{u_a^2\varLambda_1}{6{\rm \pi}}\frac{3+8\hat{k}_1^2}{\left(1+\hat{k}_1^2\right)^{11/6}}, \end{gather}$$
$$\begin{gather}\varTheta^{VK,a}_{22}(k_1)=\frac{u_a^2u_r^2\varLambda_1\varLambda_r}{6{\rm \pi}}\frac{3+8\hat{k}_1^2}{\left(1+\varLambda_r^2\hat{k}_1^2\right)^{11/6}}, \end{gather}$$
$$\begin{gather}\varTheta^{Exp}_{22}(k_1;l)=\frac{u_a^2u_r^2l}{2{\rm \pi}^2}\left({\rm \pi}+2l^2{k_1}^2\right)\exp\left({-\frac{l^2k_1^2}{\rm \pi}}\right). \end{gather}$$
Equation (5.5) refers to the spectrum for each Gaussian kernel, where 
$p$ may be varied. A Gaussian decomposition model which demonstrates an improved fit to the data is
\begin{equation} \varTheta^{GD}_{22}(k_1;l;p)=\sum_{m=1}^{35}\varTheta^{Exp}_{22}(k_1;l_m;p^*_m), \end{equation}where
\begin{equation} \boldsymbol{p^*}=\left[\boldsymbol{p}_{{24}/{6}}, \ \boldsymbol{p}_{{16}/{6}}, \ \boldsymbol{p}_{{22}/{6}}, \boldsymbol{p}_{{16}/{6}}, \ \boldsymbol{p}_{{16}/{6}}, \ \boldsymbol{p}_{{18}/{6}}, \ \boldsymbol{p}_{{16}/{6}}\right], \end{equation}
is a vector with individual components 
$\boldsymbol {p}_n$ defined as a block of five consecutive values of 
$n$.
Although we tailored this value for a good agreement, it was solely for one data set at 
$U_{\infty }=20$ with 
$\Delta x=8D$ (figure 9a). We plot the isotropic and anisotropic von Kármán models alongside our anisotropic Gaussian decomposition model for all data sets in figures 9 and 10.
Comparisons of 
$\varTheta _{22}$ models: anisotropic Gaussian decomposition (GD) as per (5.6), isotropic von Kármán (VK), (5.3), and anisotropic von Kármán (VK), (5.4), with experimental measurements at 
$U_{\infty }=20\ \textrm {m}\ \textrm {s}^{-1}$. Panels show: (a) 
$\Delta x=8D$; (b) 
$\Delta x=9.5D$; (c) 
$\Delta x=11D$; (d) 
$\Delta x=12.5D$; (e) 
$\Delta x=14D$; ( f) 
$\Delta x=16D$.
Comparisons of 
$\varTheta _{22}$ models: anisotropic Gaussian decomposition (GD) as per (5.6), isotropic von Kármán (VK), (5.3), and anisotropic von Kármán (VK), (5.4), with experimental measurements at 
$U_{\infty }=28\ \textrm {m}\ \textrm {s}^{-1}$. Panels show: (a) 
$\Delta x=8D$; (b) 
$\Delta x=9.5D$; (c) 
$\Delta x=11D$; (d) 
$\Delta x=12.5D$; (e) 
$\Delta x=14D$; ( f) 
$\Delta x=16D$.
In every example, we see that the isotropic von Kármán model underestimates the fitted curve, particularly when 
$\Delta x$ is smaller since there is less distance between the cylinder and the leading edge; thus the turbulence at the PIV window can be considered ‘more anisotropic’. Although our model was constructed to give the best fit in figure 9(a), the chosen parameter values give an improved fit compared with both von Kármán models for low and mid-frequencies across the entire parameter range. The agreement for the different flow speed 
$U_{\infty }=28\ \textrm {m}\ \textrm {s}^{-1}$ is particularly good for the 
$\Delta x=12.5D, \ 14D, \ 16D$ simulations. In contrast, all models under-predict in figure 9( f). This may be pinpointed to a surprisingly small 
$u_{rms}$ value in table 1. We reiterate that the only matching we have done for the Gaussian decomposition model is for the 
$U_{\infty }=20\ \textrm {m}\ \textrm {s}^{-1}$, 
$\varDelta _x=8D$ case.
From this initial investigation, we can conclude that including anisotropy is crucial to proper data modelling. Different behaviours in different wavenumber regions of the turbulence can be modelled well with Gaussian decomposition, generalising to different parameter values and data sets. The key takeaway is that only one experiment appears necessary to calibrate the Gaussian decomposition model, allowing the user to reproduce improved representations of velocity spectra with different characteristics.
5.2. Leading-edge noise model validation
To estimate the leading-edge noise, we use a non-dimensionalised model by considering a turbulence spectrum and gust solution independent of the length scales of the model. The result will then be dimensionalised by multiplying by 
${U^2}/{L}$ having set the dimensional constants 
$L$ as the chord of the plate and 
$U$ as the convection velocity 
$U_c$, the mean flow velocity in the vicinity of the plate.
When comparing the convection velocity 
$U_c$ with the inlet mean velocity 
$U_{\infty }$ we observe that the effect of the cylinder decreases the mean flow considerably, as seen in tables 1 and 2. Incorporating this detail into our model is crucial while accurately representing our scattering set-up by focusing on what is happening near the edge.
Finally, we divide the result by a reference pressure 
$p_a=2\times 10^{-5} \ \textrm {Pa}$, then calibrate with a vertical shift by constant 
$C=179$ (which is independent of the test case being used)
\begin{equation} \mathcal{P}(\,{f})=20 \log_{10} \left(\frac{L \varPsi\left(f,\dfrac{\rm \pi}{2}\right)}{U_c^2p_a}\right) + 171. \end{equation} We begin by plotting the experimental acoustic results (figure 11) and observe an unusual feature. Whilst the PIV incident turbulence exhibited a von Kármán type profile, the scattered noise data exhibit an RDT type decay profile of 
$f^{-{11}/{3}}$. We attribute this to the nature of the experiment, whereby a cylinder generates anisotropy. This induces a strong shear profile in the mean flow (Wu & Zhang Reference Wu and Zhang2019). Thus, scattered surface fluctuations undergo distortion en route to the far field measurement location. Again, this is a feature of the experimental set-up and would not be present in other contexts.
Experimental SPL data plotted to measure the power scaling 
$f^p$.
This key feature is omitted in our scattering model, but could be accounted for, such as in Ayton & Peake (Reference Ayton and Peake2015); however, this is unnecessarily intricate and difficult to compute. This goes against the primary motivation of creating an efficient, simple and versatile model. In summary, the use of RDT as the initial spectrum may be thought of as pre-distorting the turbulence on the approach of the edge rather than post-distorting the sound on its way to the far field. Since the scattering model provides a transfer function between the incident and scattered turbulence, we believe this is a robust and effective way to manage the background shear profile rather than overcomplicating the scattering model itself.
We hypothesise that the velocity spectrum best suited to account for this would be the RDT model (Hunt Reference Hunt1973) with length scale parameters deduced from the PIV data.
In figures 12(a) and 12(b), we plot comparisons at different 
$U_{\infty }$ of the experimental data and RDT model predictions, demonstrating the overall agreement of the trend in which noise decreases as the anisotropic ratios decrease. There is a consistent under-prediction at low frequencies, possibly due to vortex shedding in the cylinder wake. We see better agreement across the mid and high frequencies for all anisotropic ratios. The change in velocity and length scales shifts the overall spectrum to the right, creating a broader peak and matching the observed behaviour, albeit not capturing the previously mentioned spike due to the cylinder vortex shedding.
Comparisons of experimental far field PSD data with the analytic model (5.8) at different anisotropic ratios and mean flow velocities. In each figure the left panel indicates the experimental data while the right is our prediction. Panels show: (a) 
$U_{\infty }=20\ \textrm {m}\ \textrm {s}^{-1}$, empirical constants taken from table 1; (b) 
$U_{\infty }=28\ \textrm {m}\ \textrm {s}^{-1}$, empirical constants taken from table 2.
We further improve the model by taking 
$k_1={2{\rm \pi} L\varrho f}/{U_c}$ with 
$\varrho =3/2$ a constant used to effectively shift the whole spectrum to the left, allowing the peak to demonstrate a better agreement with the data. That is, we consider
\begin{equation} \mathcal{P}(f)=20\log_{10}\left(\frac{L\varPsi\left(\varrho f,\dfrac{\rm \pi}{2}\right)}{U_c^2p_a}\right)+179. \end{equation}
The explanation for why this is necessary may lie with the convection velocity being lower than measured. Right at the plate itself, our eddies convent slower than first anticipated. We plot the 
$\Delta x=8D$ comparisons for both mean flow velocities in figure 13 with this new shift included to demonstrate the two agreements more carefully. Although the scaling concerning mean flow speed 
$U_{c}$ is accurate at high frequencies, the impact of the aforementioned horizontal shift 
$\varrho$ is a crucial difference at low frequencies since it better fits the peak, suggesting future work may need to consider this difference with more care. An improved agreement is also observed for the remaining anisotropic ratios at the two studied velocities. It is worth reiterating that these two alterations are implemented for better agreement with only one simulation (
$U_{\infty }=20, \Delta x=8D$). This agreement carries over from those choices to the second data set.
Comparison of improved analytical model (5.9) with experimental data for both the 
$U_{\infty }=20$, 
$28\ \textrm {m}\ \textrm {s}^{-1}$ cases when 
$\Delta x=8D$.
5.3. Improvements with Gaussian decomposition
The primary focus of this study is verifying that anisotropy is an essential feature of turbulence modelling that can be included in a theoretical mathematical model for leading-edge noise with ease. Having verified the model, we turn to the secondary goal, that this theoretical model can be adapted via Gaussian decomposition to produce a versatile and efficient alternative that can be used in experimental campaigns to make predictions and better agreements with data.
We first formally define our Gaussian decomposition PSD model as
\begin{equation} \mathcal{P}^{GD}(f)=\sum_{m=1}^M \mathcal{P}^{GD}_m(f;l_m;p_m) ,\end{equation}
where each 
$\mathcal {P}^{GD}_m(f;l_m;p_m)$ can be thought of as a PSD model of the form (5.8) in which we have selected a Gaussian kernel such as (3.10) with length scale 
$l_m$ and scaling parameter 
$p_m$ used within 
$\varPsi$.
From the previous subsection, we deduce that we have to tailor our Gaussian decomposition model around kernels with larger 
$p_m$ values than those for the velocity field to recalculate the RDT type effects due to the pressure distortions in the sheared flow.
In other words, the best choice representing the incoming turbulence is not necessarily the best choice for our scattering model (if there were no induced shear flow, it would be). The only way to ensure we account for these unknown distortions is to match them to a measurement that includes them, namely the scattered SPL. We, therefore, take one SPL measurement to generate our turbulence decomposition, which can then be reused for all other predictive cases.
We alter our chosen 
$\boldsymbol {p}^*$ from (5.7) to
\begin{equation} \boldsymbol{p^*}=\left[\boldsymbol{p}_{{43}/{12}}, \ \boldsymbol{p}_{{22}/{6}}, \ \boldsymbol{p}_{{45}/{12}}, \boldsymbol{p}_{{22}/{6}}, \ \boldsymbol{p}_{{23}/{6}}, \ \boldsymbol{p}_{{24}/{6}}, \ \boldsymbol{p}_{{{60}/{6}}}\right], \end{equation}
which was chosen before to fit the 
$U_{\infty }=20$, 
$\Delta x = 8D$ data set. We retain all other PIV data values from tables 1 and 2. We define
\begin{equation} \mathcal{P}(f)=\sum_{m=1}^{35}\mathcal{P}^{GD}(f;l_m;p_m), \end{equation}
to be the corresponding PSD model and plot the resulting PSD in figure 14. We find that, by using these 
$p$-values, we can better reflect the sharp peak in the data, and then we smoothly return to the gradient achieved using the RDT model.
Finally, we test our hypothesis of how versatile this approach can be when applying the 
$\boldsymbol {p^*}$ decomposition model to each data set from tables 1 and 2, despite it being determined from one data set alone.
We consider each velocity separately and average over every 
$\Delta x$ value, then investigate the change in noise against our model while repeating the above procedure. That is, we plot
\begin{equation} \bar{\mathcal{P}}(\,{f})=\frac{1} {6}\sum_{i=1}^6{\mathcal{P}}_i(\,{f}), \end{equation}
with each 
${\mathcal {P}}_i(\,{f})$ referring to the model output when using empirical data from the 
$i$th data set from the respective table in this section.
As seen in figure 15, we find an even better agreement when using our Gaussian decomposition model instead of the RDT model at almost all frequencies. There are notable improvements at low and high frequencies for 
$U_{\infty }=20$ and at mid-range frequencies for 
$U_{\infty }=28$. Since we have used averages, this new approach will improve the data fitting for various parameter values. Despite fitting our model to only one data set, the improvement carries on for all 11 other data sets, which is unsurprising since the same level of distortion will occur as it is induced by the same cylinder each time. This correction is solely due to our experimental set-up, which was required to generate suitable anisotropic ratios.
Comparisons of the change in noise 
$\Delta \bar {\mathcal {P}}(f)$ (5.13) when compared with experimental data for both the 
$p$-varied Gaussian decomposition model and the previously used RDT model (5.8) having averaged over all six experimental data sets for 
$U_{\infty }=20\ \textrm {m}\ \textrm {s}^{-1}$ (a) and 
$U_{\infty }=28\ \textrm {m}\ \textrm {s}^{-1}$ (b).
From this, we begin to see how this particular model is reproducible and that the anisotropy modelling within the turbulence spectrum ensures the mathematical model is suitable for various styles of incidental turbulence. By utilising the Gaussian decomposition model, we can generate an arbitrary turbulence spectrum 
$\varPhi _{22}$ for our model or otherwise. In particular, we can take one measurement to establish the behaviour of the turbulence in one case and then generalise this to a variety of other types and contexts.
With this, we find that, by looking at one set of results (
$\Delta x=8D, U_{\infty }=20$) and using the model to improve the agreement with the data, we can apply this to all our other data sets and see similar improvements. Thus, we can take a minimal data approach to our modelling and, having reproduced experimental data in this case, we can potentially apply the approach to other cases.
6. Conclusion
In this paper, we have derived a model for leading-edge noise that implements changes in turbulent properties between different directions in a simple and numerically efficient manner. Our model shows that introducing anisotropy to the turbulent mean flow can affect the PSD shape and height. The model was matched for one simulation and successfully demonstrated improvements for all other cases in which the velocity or anisotropic ratios varied.
We have further outlined an experimental campaign which measures leading-edge noise in a turbulent stream. These experimental results inform both the turbulence spectrum modelling and the gust scattering.
The duality of the work between the gust scattering and the turbulence spectrum is beneficial, and the approach should be amenable to other useful areas in aeroacoustics.
Comparisons of the experimental data and theoretical predictions have shown good agreement. When constructing our model, we accounted for shear flow distorting the turbulence that scatters off the leading edge. We anticipate that anisotropic turbulence created by a turbulence grid that avoids shear generation would be modelled similarly with a Gaussian decomposition model that shares more features with the von Kármán model, as seen when modelling the incident velocity.
The method of Gaussian decomposition was applied to obtain an even better agreement with the experimental data. We first tested the ability to use a tailored decomposition of kernels to fit the one-dimensional spectrum 
$\varTheta _{22}(k_1)$ to experimental measurements. After demonstrating that this was possible and providing better agreement with all data sets, we discussed how this approach would work for the SPL after accounting for the near-edge effects.
This led to an altered Gaussian decomposition model that demonstrated an even better fit than the anisotropic RDT model. The significant application of this theory in future work will be to use experimental data that are easy to obtain, such as from hot-wire measurements, and utilise these within our analytical model via Gaussian decomposition to find the correct weightings. With these weightings, we may expect a good agreement with our far field noise data.
Funding
A.D.G.H. acknowledges support from EPSRC studentship EP/T517847/1 and would like to thank M.J. Priddin for his helpful discussions regarding the turbulence spectrum and help with proofreading. L.J.A. acknowledges support from EPSRC Early Career Fellowship EP/P015980/1. This research was supported partially by the Australian Government through the Australian Research Council's Discovery Projects funding scheme (project DP210102288).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Model derivations
Amiet presumes the incoming turbulence 
$\boldsymbol { {u}}^{(I)}=\boldsymbol {\nabla }\phi ^{(I)}$ takes the form
\begin{equation} \phi^{(I)}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}w_2(\boldsymbol{{k}})\exp({{\rm i}\boldsymbol{{k}} \boldsymbol{\cdot}\boldsymbol{{x}}-{\rm i}\omega t})\,\text{d}k_2\,\text{d}k_3, \end{equation}
where 
$\boldsymbol { {k}}=(k_1,k_2,k_3)$ and convects with the mean flow, such that 
$k_1={\omega }/{U_{c}}$.
The function 
$w_2(\boldsymbol { {k}})$ is defined as the Fourier transform of our incident vertical velocity field, 
$v$. We calculate the pressure scattered from a single gust at fixed 
$\boldsymbol { {k}}$ as some
\begin{equation} p_g(\boldsymbol{{k}},\theta)=w_2(\boldsymbol{{k}})P(\boldsymbol{{k}},\theta), \end{equation}
where 
$\theta$ is a far field observer angle.
Then, by linearity, our total far field pressure due to the scattering of the turbulent flow will be the integral of all these Fourier components of the incoming turbulence
\begin{equation} p_t(\omega,\theta)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}w_2(\boldsymbol{{k}}) P(\boldsymbol{{k}},\theta)\delta\left(k_1-\frac{\omega}{U_{c}}\right) \text{d}k_1\,\text{d}k_2\,\text{d}k_3, \end{equation}
we have also applied the frozen turbulence assumption to relate 
$k_1$ to 
$\omega$.
Then, the sound pressure spectral density is given by the time-averaged statistical variable
\begin{equation} \varPsi(\omega,\theta)=\lim_{T\to\infty}\frac{\rm \pi}{T}p_t(\omega,\theta)p_t^*(\omega,\theta). \end{equation}Noting that
\begin{equation} \varPhi_{22}(\boldsymbol{{k}})=\lim_{T\to\infty}\frac{\rm \pi}{T}w_2(\boldsymbol{{k}})w_2^*(\boldsymbol{{k}}), \end{equation}the scattered PSD can be written as
\begin{equation} \varPsi(\omega,\theta)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} |P(\boldsymbol{{k}},\theta)|^2\varPhi_{22}(\boldsymbol{{k}})\delta\left(k_1-\frac{\omega}{U_{\infty}}\right) \text{d}k_1\,\text{d}k_2\,\text{d}k_3. \end{equation}Appendix B. The Wiener–Hopf solution for the gust-scattering solution
We Fourier transform our continuity conditions on the wake, having assumed our solution is of the form
\begin{equation} \phi_s(x,y)=\frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty} A(\alpha) \exp({-{\rm i}\alpha x-\gamma |y|})\,\text{d}\alpha, \end{equation}
where 
$\gamma =\sqrt {\alpha ^2-{k^*}^2}$, with the branch cut of this square root taken to ensure 
$\gamma$ has a positive real part and 
$A(\alpha )$ is the unknown function of 
$\alpha$ we solve for
\begin{equation} \left.\begin{gathered} A(\alpha)-B(\alpha)=u_1^+,\\ -\gamma \left(A(\alpha)+B(\alpha) \right) =u_2^+. \end{gathered}\right\} \end{equation}
Each 
$u_i^+$ represents an unknown analytic function in the upper-half complex plane. Due to the symmetry of our other boundary condition along the plate, we can assume 
$A(\alpha )=-B(\alpha )={u_1^+}/{2}={U^+}/{2}$. We Fourier transform the Neumann boundary condition along the plate to obtain
\begin{equation} -\gamma A(\alpha)=\frac{k_2}{\alpha+\tilde{k}}+L^-, \end{equation}
where 
$L^-$ is analytic in the lower-half complex plane.
Putting this together, our Wiener–Hopf equation is
\begin{equation} \frac{-\gamma U^+}{2}=\frac{k_2}{\alpha+\tilde{k}}+L^-. \end{equation}
To solve this, we rearrange (B4) so the left-hand side is upper analytic and the right-hand side is lower analytic; therefore, each side is equal to some entire function 
$E(\alpha )$
\begin{equation} \underbrace{ \gamma^+(\alpha)U^+(\alpha)+\frac{2k_2}{\left(\alpha+\tilde{k}\right)}\frac{1}{\gamma^-\left(-\tilde{k}\right)}}_{\text{Upper analytic}} =\underbrace{ \frac{-2k_2}{\alpha+\tilde{k}}\left(\frac{1}{\gamma^-(\alpha)}-\frac{1}{\gamma^-\left(-\tilde{k}\right)}\right)-\frac{L^-(\alpha)}{\gamma^-(\alpha)}}_{\text{Lower analytic}}=E(\alpha). \end{equation}
Here, we have used the fact that 
${1}/({\alpha +\tilde {k}})$ is upper analytic and then applied the method of pole decomposition to additively factorise 
${1}/{(\alpha +\tilde {k})\gamma ^-(\alpha )}$.
To find 
$U^+$, we argue that 
$E(\alpha )=0$ by showing that each side of the equation tends to zero as 
$\alpha \to \infty$ and then applying Liouville's theorem. This follows naturally as a consequence of the edge conditions
\begin{equation} \left.\begin{gathered} \phi_s={O}(1) \quad\text{as}\ x\to 0,\\ \frac{\partial \phi_s}{\partial y}={O}\left(x^{-{1}/{2}}\right) \quad \text{as} \ x\to 0, \end{gathered}\right\} \end{equation}
which must be Fourier transformed to provide behaviour as 
$\alpha \to \infty$.
Having found 
$U^+$, we can then deduce 
$A$ and thus 
$\phi _s$
\begin{equation} \phi_s(x,y)=\frac{-k_2\,\textrm{sgn}(y)}{2{\rm \pi}}\int_{-\infty}^{\infty} \frac{1}{\gamma^-\left(-\tilde{k}\right)\gamma^+(\alpha)\left(\alpha+\tilde{k}\right)}\exp({-{\rm i}\alpha x-\gamma|y|})\,\text{d}\alpha. \end{equation}Our above result noticeably does not undo the Prandtl–Glauert transform, which would be unnecessarily complex when doing the steepest descent procedure. Instead, we leave our answer in ‘Prandtl–Glauert space’. This discrepancy will have very little impact on our approximation since, in practical applications, our Mach number will be small. As stated in Myers & Kerschen (Reference Myers and Kerschen1995), we can relate the Prandtl–Glauert polar coordinates back to physical space via
\begin{equation} \left.\begin{gathered} \tan\theta_{PG}=\beta\tan\theta\\ r_{PG}=r\sqrt{1-M^2\sin^2\theta}. \end{gathered}\right\} \end{equation}Appendix C. Deriving the anisotropic spectrum weighting function
Recall that we are aiming to find a function 
$f(l_a,l_t)$ solving
\begin{equation} \frac{u_a^2u_r^2(L\varLambda_1)^5\varLambda_r\varGamma(p)(k_1^2+ k_3^2)}{{\rm \pi}^{3/2}\varGamma\left(p-\frac{5}{2}\right)[1+L^2\varLambda_1^2\left(k_1^2+\varLambda_r^2k_2^2+k_3^2\right)]^p}= \int_0^\infty\int_0^\infty f(l_a,l_t)\tilde{\varPhi}_{22}^{a,G}(\boldsymbol{{k}})\,\text{d}l_t\,\text{d}l_a. \end{equation}We define
\begin{equation} S_a=k_2^2, \quad S_t=k_1^2+k_3^2, \quad T_a=\frac{l_a^2}{\rm \pi}, \quad T_t=\frac{l_t^2}{\rm \pi}. \end{equation}After using these coordinate transformations and rearranging, we obtain
\begin{align} &\frac{4u_r^2(L\varLambda_1)^5\varLambda_r\varGamma(p))}{{\rm \pi}(2+u_r^2)\varGamma\left(p-\frac{5}{2}\right)}[1+L^2\varLambda_1^2\left(\varLambda_r^2S_a+S_t\right)]^{{-}p} \nonumber\\ &\quad =\int_0^\infty\int_0^\infty f(\sqrt{{\rm \pi} T_a},\sqrt{{\rm \pi} T_t})\frac{T_a T_t^{3/2}}{(1+\sigma) T_a+T_t}\text{e}^{- T_aS_a}\text{e}^{- T_tS_t}\,\text{d}T_t\,\text{d}T_a\nonumber\\ &\quad =\mathcal{L}\left[\mathcal{L}\left[ f\left(\sqrt{{\rm \pi} T_a},\sqrt{{\rm \pi} T_t}\right)\frac{T_aT_t^{3/2}}{(1+\sigma) T_a+T_t}\text{e}^{{-}T_aS_a}\text{e}^{{-}T_tS_t} \right]\left(\frac{l_t^2}{\rm \pi}\right)\right]\left(\frac{l_a^2}{\rm \pi}\right). \end{align}Inverting the Laplace transforms and simplifying, we obtain
\begin{align} f(l_a,l_t)&=\frac{4{\rm \pi} u_r^2(L\varLambda_1)^5\varLambda_r}{(2+u_r^2)\varGamma(p-\frac{5}{2})}\frac{(1+\sigma)l_a^2+l_t^2}{l_a^2l_t^3} \left(\frac{l_a^2}{{\rm \pi}(L\varLambda_r\varLambda_1)^2}\right)^p\delta\left(l_t^2-\frac{l_a^2}{\varLambda_r^2}\right)\nonumber\\ &\quad \times\exp\left(-\frac{l_a^2}{{\rm \pi}(L\varLambda_r\varLambda_1)^2}\right), \end{align}which we can re-write into the form
\begin{equation} \left.\begin{gathered} f(l_a,l_t)=\frac{4{\rm \pi} u_r^2(L\varLambda_1)^5\varLambda_r}{(2+u_r^2)\varGamma\left(p-\frac{5}{2}\right)}\tilde{f}(l_a)\tilde{g} (l_a,l_t),\\ \tilde{f}(l_a)=\left(\frac{l_a^2}{{\rm \pi}(L\varLambda_r\varLambda_1)^2}\right)^p\exp\left(-\frac{l_a^2}{{\rm \pi}(L\varLambda_r\varLambda_1)^2}\right),\\ \tilde{g}(l_a,l_t)=\frac{(1+\sigma)l_a^2+l_t^2}{l_a^2l_t^3}\delta\left(l_t^2-\frac{l_a^2}{\varLambda_r^2}\right), \end{gathered}\right\} \end{equation}as per (3.28).
Appendix D. Approximating 
$E(\kappa )$ for anisotropic turbulence
Recall that we introduced the energy density spectrum
\begin{equation} {E}(\kappa;p)=\tfrac{1}{2}\int\varPhi_{ii}(\boldsymbol{{k}};p)\,\text{d}A(\kappa). \end{equation}This appendix aims to study an essential application of the Gaussian decomposition method. We derive an analytical model for the energy density function for anisotropic turbulence.
From (3.3), we may see how the procedure for Gaussian decomposition would also apply to 
$\varPhi _{11}$ and 
$\varPhi _{33}$. The value of 
$E(\kappa )$ can be found experimentally with less difficulty than other turbulence statistics. It is a scalar measure of the turbulent kinetic energy held by eddies with wavenumber 
$|k|$; however, there are no representations of this for anisotropic turbulence. Analytically calculating the anisotropic energy density function 
$E^a(\kappa ;p)$ using (3.3) is impractical without the Gaussian approximation. Integrating with exponential functions is analytically tidier since our methods often use Fourier or Laplace transforms to solve for unknowns or switch from physical to wavenumber space.
Therefore, since our weighting function only depends on the length scale parameter and not on our integration variables in (3.3), we can see that the corresponding anisotropic energy density Gaussian, which we denote 
$E^a(\kappa )$, can be approximated using our already derived Gaussian decomposition model
\begin{equation} E^{a}(\kappa;p)=\sum_{m=0}^{M}\frac{\hat{f}(p;l_m)\Delta l_m}{\sum_{m=0}^{M}\hat{f}(p;l_m)\Delta l_m}E^{a,G}(\kappa;p;l_m). \end{equation}Thus, we only need to calculate the energy density spectrum for an arbitrary anisotropic Gaussian kernel
\begin{equation} E^{a,G}(\kappa)=\frac{\kappa^2}{2}\int_0^{2{\rm \pi}}\int_{0}^{\rm \pi} \varPhi_{ii}\left(\kappa\sin\theta\cos\phi,\kappa\cos\theta,\kappa\sin\theta\sin\phi\right)\sin\theta\,\text{d}\theta\,\text{d}\phi. \end{equation}Inserting
\begin{equation} \varPhi_{ii}=\frac{u_a^2u_r^2\kappa^2l^5\varLambda_r(2\cos^2\theta+(1+\gamma)\sin^2\theta)}{{\rm \pi}^4}\exp\left(-\frac{(l\kappa)^2}{\rm \pi}(\cos^2\theta+\varLambda_r\sin^2\theta)\right), \end{equation}and integrating gives us
\begin{equation} E^{a,G}(\kappa)=\frac{u_a^2u_r^2\varLambda_r l^5\kappa^4}{{\rm \pi}^3}\exp\left({-\frac{(\kappa l)^2}{\rm \pi}}\right)\underbrace{\left[\frac{\gamma-1}{z^2}{\rm e}^{{-}z^2}-\frac{\sqrt{\rm \pi}}{2}\left(\frac{\gamma-1}{ z^3}-2\frac{\gamma+1}{z}\right)\text{erf}(z) \right]}_{\varphi(z)}, \end{equation}
where 
$z={\kappa l \sqrt {\varLambda _r^2 -1}}/{\sqrt {{\rm \pi} }}$.
When plotting the function numerically, we circumvent issues caused by the presence of the error function in 
$\varphi (z)$ by using expansions for small 
$z$ (accounting for the limit 
$\kappa \to 0$ as well as 
$\varLambda _r\to 1$) while also taking the limit 
$\kappa \to \infty$, or equivalently 
$z\to \infty$, 
$\varLambda _r\neq 1$
\begin{align}
E^{a,G}(\kappa)&= \frac{u_a^2u_r^2\varLambda_r
l^5\kappa^4}{{\rm \pi}^3}\times \exp\left({-\frac{(\kappa
l)^2}{\rm \pi}}\right)\nonumber\\
&\quad \times \left(\frac{4(\gamma+2)}{3}+\frac{4(\gamma+4)}{15}z^2+\frac{2(\gamma+6)}{35}z^4
+{O}(z^6) \right) \quad \text{as} \ z\to 0,
\end{align}
\begin{align} E^{a,G}(\kappa)&= \frac{u_a^2u_r^2\varLambda_r l^5\kappa^4}{{\rm \pi}^3}\nonumber\\ &\quad \times \exp\left({-\frac{(\kappa l \varLambda_r)^2}{\rm \pi}}\right) \left(\frac{2}{z^2}+\frac{\gamma}{z^4}+{O}\left(\frac{1}{z^6}\right) \right) \quad \text{as} \ z\to \infty, \ \varLambda_r<1. \end{align} A demonstration of the convergence in 
$N$ of the model is given in figure 16.
Convergence of the Gaussian decomposition method when used to approximate the von Kármán 
$E(\kappa )$ model (D2).
Open access
