1. Introduction
Global climate change has been a significant topic on the political agenda for many years (Lindzen Reference Lindzen1994). Renewable energy sources are considered the most suitable alternative to fossil fuels; among them, wind energy has become increasingly popular and the total installed wind energy capacity is expected to substantially increase and scale up in size and power rating in the next decades (Wagner, Bareiss & Guidati Reference Wagner, Bareiss and Guidati1996; Njiri & Söffker Reference Njiri and Söffker2016; Bošnjaković et al. Reference Bošnjaković, Katinić, Santa and Marić2022). Among others, acoustic pollution generated by wind turbines constitutes a serious issue for the population and wildlife species. The coexistence of several complex flow phenomena contributes to the noise generation mechanism. Therefore, a more comprehensive understanding of noise generation mechanisms and noise propagation is of crucial importance to act towards the mitigation of environmental pollution (Stöber & Thomsen Reference Stöber and Thomsen2021; Teff-Seker et al. Reference Teff-Seker, Berger-Tal, Lehnardt and Teschner2022).
Here, we first review relevant literature regarding the fluid dynamic field generated by rotors and, then, we review literature dedicated to rotors’ flow-induced noise.
1.1. Rotor fluid dynamics
The fluid dynamic field generated by a rotor immersed in a flow is composed of two main features, namely the flow over the rotating blade surfaces and the wake. Typically, a three-dimensional turbulent boundary layer (TBL) together with flow separation develops over the blades, given the variation in Reynolds number along the blade axis. They give rise to chaotic pressure fluctuations, thereby promoting conversion of turbulent kinetic energy into acoustic energy (Lee et al. Reference Lee, Ayton, Bertagnolio, Moreau, Chong and Joseph2021). In addition, a complex wake develops, characterised by the presence of helicoidal vortices, specifically the tip and root vortices. These structures evolve and interact with each other, producing wake destabilisation related to instability mechanisms (Sørensen Reference Sørensen2011). The wake is essentially composed of three regions: the near wake, the intermediate wake and the far wake (see, among others, Porté-Agel et al. (Reference Porté-Agel, Bastankhah and Shamsoddin2020) for the case of wind turbines).
The near wake is influenced by the presence of the rotor itself and its geometrical properties; a three-dimensional complex flow is present, characterised by coherent periodic helicoidal vortex structures, namely the root and tip vortices; in the case of an isolated marine propeller, a strong rectilinear vortex is evident, namely the hub vortex. These vortex structures are connected through a thin vortex sheet released by the blades. The breakdown and instability of the hub and tip vortices, occurring in the intermediate wake, are responsible for the transition from the near to the far wake, with a considerable amount of turbulent kinetic energy detected at the outer radii. Finally, a fully composite turbulent wake develops in the far wake. Further downstream, a self-similar behaviour is detected with the mean streamwise velocity deficit profile collapsing into a single Gaussian distribution (see, for example, Abkar & Porté-Agel Reference Abkar and Porté-Agel2015; Xie & Archer Reference Xie and Archer2015; Dar et al. Reference Dar, Berg, Troldborg and Patton2019).
The predominant mechanism of rotor wake instabilities has been discussed in several studies considering the instability of rotating helical vortices. Since the seminal work of Joukowski (Reference Joukowski1912), theoretical studies focused on the stability of a single helical vortex filament (see, for example, Levy & Forsdyke Reference Levy and Forsdyke1928; Widnall Reference Widnall1972), finding that three mechanisms, namely short wave and long wave instability and mutual inductance, influence the stability of a single helical vortex. In the experiments of Felli, Camussi & Di Felice (Reference Felli, Camussi and Di Felice2011), the mechanisms of evolution of propeller tip and hub vortex in the transition region and the far field were investigated. They suggested that the main tip and hub vortex destabilisation mechanism is the mutual induction mode, promoting leap-frogging and vortex grouping, which is also supported by the experimental analysis performed by Nemes et al. (Reference Nemes, Jacono, Blackburn and Sheridan2015).
Several experimental studies have been conducted on the wakes generated by wind turbines. For example, Sarmast et al. (Reference Sarmast, Dadfar, Mikkelsen, Schlatter, Ivanell, Sørensen and Henningson2014) studied the mutual-inductance instability of the tip vortex. Sherry et al. (Reference Sherry, Nemes, Lo Jacono, Blackburn and Sheridan2013a ,Reference Sherry, Sheridan and Jacono b ) analysed the interaction between the root and tip vortices, and Iungo et al. (Reference Iungo, Viola, Camarri, Porté-Agel and Gallaire2013) performed a linear stability analysis of a wind turbine wake. The tip vortex instability mechanism and its relationship with turbulent mixing were analysed by Lignarolo et al. (Reference Lignarolo, Ragni, Scarano, Ferreira and van Bussel2015), considering a model-scale wind turbine. Recently, the experimental analysis of a full-scale 2.5 MW wind turbine (Dasari et al. Reference Dasari, Wu, Liu and Hong2019) has given a better comprehension of the near-wake characteristics.
From a modelling point of view, the simplest way to study wind turbine aerodynamics is using the actuator disc model. Over the years, this methodology has been improved with the blade element momentum theory (Glauert Reference Glauert1935) to simulate rotor turbulent wake. These simplified models represent the wind turbine as a volume force placed into the fluid flow. Therefore, they are computationally inexpensive when compared with geometry-resolving simulations. The coupling between the actuator disc or actuator line model with a Navier–Stokes finite volume solver was first developed by Sorensen & Shen (Reference Sorensen and Shen2002) and successively by Mikkelsen (Reference Mikkelsen2004). The authors evaluated the effectiveness of this numerical technique, showing a high degree of correlation with measurements for power generation. More recent studies assessed the applicability of large eddy simulation (LES) coupled with the actuator model to study the dynamics of turbulent structures as well as the instability mechanism involved in wind turbine wake (Jimenez et al. Reference Jimenez, Crespo, Migoya and García2008; Ivanell et al. Reference Ivanell, Sørensen, Mikkelsen and Henningson2009; Ivanell et al. Reference Ivanell, Mikkelsen, Sørensen and Henningson2010; Sørensen Reference Sørensen2011; Sarmast et al. Reference Sarmast, Dadfar, Mikkelsen, Schlatter, Ivanell, Sørensen and Henningson2014; Sarlak, Meneveau & Sørensen Reference Sarlak, Meneveau and Sørensen2015; Martinez-Tossas et al. Reference Martinez-Tossas, Churchfield, Yilmaz, Sarlak, Johnson, Sørensen, Meyers and Meneveau2018). The instability mechanisms producing turbulence were studied numerically by Ivanell et al. (Reference Ivanell, Mikkelsen, Sørensen and Henningson2010), Sarmast et al. (Reference Sarmast, Dadfar, Mikkelsen, Schlatter, Ivanell, Sørensen and Henningson2014), and Sørensen et al. (Reference Sørensen, Mikkelsen, Henningson, Ivanell, Sarmast and Andersen2015). They confirmed that the mutual inductance mode is the main mechanism of destabilisation and loss of coherence of the tip vortices system, first suggested by the inviscid studies of Widnall (Reference Widnall1972) and successively corroborated by the experiments of Felli et al. (Reference Felli, Camussi and Di Felice2011). However, these numerical analyses are restricted by the use of the actuator line model, which does not resolve the rotor geometry, and often the instability mechanism is triggered by the use of artificial fluctuations in the fluid flow.
Recently, the immersed boundary (IB) methodology in conjunction with LES was used to take into account the rotor geometry with the aim of analysing the rotor device wake dynamics (Posa, Broglia & Balaras Reference Posa, Broglia and Balaras2021; Posa Reference Posa2022; Posa Reference Posa2023). Among them, the work of Posa et al. (Reference Posa, Broglia and Balaras2021) is of particular interest in wind turbine aerodynamics. The authors considered a very fine spatio-temporal resolution of the fluid flow, considering the same model wind turbine as in the round robin experimental campaign (Gaurier et al. Reference Gaurier, Germain, Facq, Johnstone, Day, Nixon, Di Felice and Costanzo2015). They provided details on the development of wind turbine wake instability, such as short- and long-wave instability, as well as mutual inductance. Their results are in accordance with previous experiments of Lignarolo et al. (Reference Lignarolo, Ragni, Scarano, Ferreira and van Bussel2015), and numerical analysis of Ivanell et al. (Reference Ivanell, Mikkelsen, Sørensen and Henningson2010), Sørensen (Reference Sørensen2011) and Sarmast et al. (Reference Sarmast, Dadfar, Mikkelsen, Schlatter, Ivanell, Sørensen and Henningson2014). In addition, they highlighted the process of momentum recovery and wake contraction associated with turbulent mixing. They reported that the wake instability is triggered by the shear between tip vortices and the neighbouring blade’s wake, associated with increasing levels of turbulent shear stresses and turbulent production in the region within consequent tip vortices where the mutual inductance develops.
To summarise, several experimental and numerical investigations are present in the literature, exploiting the dynamics of rotor wake, as well as the instability mechanism responsible for the transition from near to far wake. Most of them used simplified methods such as actuator technique, which may give some information about the characteristics of the wake. Geometry-resolving simulations are more recent, and few of them have been carried out using eddy-resolving techniques, like LES.
1.2. Rotor acoustics
The noise produced by rotors, especially wind turbines, is classified into three main groups, namely low-frequency noise, inflow turbulence noise and aerofoil self-noise (Wagner et al. Reference Wagner, Bareiss and Guidati1996). The rotation of blades produces the so-called thickness noise, characterised by a tonal spectrum, with peaks at the blade passing frequency (BPF) and its super-harmonics. The rotating lifting surfaces, namely the loading noise source, contribute to BPF and broadband noise, the latter in the high-frequency range; this is because the transition from laminar to turbulent boundary layer and flow separation induce fluctuating pressures on the blade surface, as well as a thin vortex sheet radiating from the blade trailing edge. This noise-generation mechanism is known as trailing-edge noise and contributes to the acoustic signature over a wide range of frequencies (Lee et al. Reference Lee, Ayton, Bertagnolio, Moreau, Chong and Joseph2021). The noise emitted by the tip vortex influences both the low- and high-frequency regions of the spectrum, and is not yet fully understood. The incoming turbulence, associated with the atmospheric boundary layer (ABL) or upwind wind turbines (in wind farm applications), produces a broadband noise. Depending on the characteristics of the turbulence, such as the time and length scales of the disturbances, the inflow turbulence noise mechanism contributes both to high- and low-frequency noise (see, among others, Sevik Reference Sevik1974; Alexander, Devenport & Glegg Reference Alexander, Devenport and Glegg2017; Wang, Wang & Wang Reference Wang, Wang and Wang2021; Zhou, Wang & Wang Reference Zhou, Wang and Wang2024).
Experimental studies have contributed to the knowledge of noise generated by rotors. Specifically in the field of wind turbines, Debruijn, Stam & Dewolf (Reference Debruijn, Stam and Dewolf1984) performed acoustic measurements in small and medium wind energy conversion (WEC) systems, finding that the radiated noise is more affected by the rotational speed than the wind inflow velocity. Van der Borg & Vink (Reference Van der Borg and Vink1994) carried out noise measurements on wind turbines, evaluating separately the contribution of tonal noise and that obtained with different blades on the same rotor, finding that changing blade geometry has a marginal effect on noise production. Oerlemans, Sijtsma & López (Reference Oerlemans, Sijtsma and López2007) performed acoustic field measurements for a three-blade wind turbine with a rotor diameter of 58 m, to characterise noise sources and to verify the dominance of trailing edge noise over other mechanisms. Ramachandran, Raman & Dougherty (Reference Ramachandran, Raman and Dougherty2014) exploited the beamforming method to measure the aerodynamic noise sources in wind turbines, confirming that the trailing-edge noise source is dominant in the noise produced by the blade.
For engineering purposes, semi-empirical and analytical methods are of common use for the quantification of rotor noise generation and propagation (see, among others, Zhu et al. Reference Zhu, Shen, Barlas, Bertagnolio and Sørensen2018; Lee et al. Reference Lee, Ayton, Bertagnolio, Moreau, Chong and Joseph2021 and references therein).
In recent years, computational aero/hydro-acoustics (CAA/CHA) has been continuously developed and employed for the prediction of noise generated by fluid flow in a wide class of engineering applications. The acoustic characterisation of rotor devices requires high-fidelity numerical simulations, and over the years, it has been exploited by several researchers. Among others, a widely used modelling approach relies on the acoustic analogy first developed in the seminal paper of Lighthill (Reference Lighthill1952), where the author derived an inhomogeneous wave equation for the acoustic pressure from the compressible Navier–Stokes equations. Under the assumption that the acoustic signal does not modify the flow field, it is possible to decompose the CAA/CHA problem into two sequential steps, first solving the fluid dynamic field and then successively the acoustic one, using fluid dynamic data. This approach is particularly advantageous in low-speed flow processes (
${\textit{Ma}} \lt 0.3$
, where
${\textit{Ma}} = U\!/c$
is the Mach number defined as the ratio between the velocity scale of the flow and the speed of sound), where the fluid-dynamic field can be considered in the incompressible regime. This is the case with wind turbines or ship propellers.
The Lighthill (Reference Lighthill1952) theory was successively extended by Ffowcs Williams & Hawkings (Reference Ffowcs Williams and Hawkings1969), considering a body in relative motion in a fluid flow. They applied the free-space Green function to solve the wave equation in an infinite homogeneous fluid medium to obtain the acoustic pressure starting from a fluid dynamic dataset. A main advantage of the use of the Ffowcs Williams–Hawkings (FWH) equation is that the total noise emitted by a body moving in a fluid flow is decomposed as the sum of two main contributions: the noise generated by the surface of the body (linear part) and that due to the three-dimensional fluid flow (nonlinear part). In particular, the linear noise component is due to the fluid displaced by the body itself (thickness term) and to the pressure fluctuations acting on it (loading term); the nonlinear term is generated by the fluctuating velocity field around the body.
Historically, the wake contribution to the noise has been omitted from the solution in many applications to simplify calculation and since its contribution to the total noise has been considered negligible, particularly in the far-field. This is because the loading term, associated with pressure variations over the solid surfaces, is the dominant noise source in many situations. For example, in the field of wind turbines, the acoustics of the complete machine under operational conditions was studied by Nelson et al. (Reference Nelson, Cain, Raman, Chan, Saunders, Noble, Engeln, Dougherty, Brentner and Morris2012) using a delayed detached eddy simulation and the linear formulation proposed by Farassat (Reference Farassat2007).
However, Brentner & Farassat (Reference Brentner and Farassat2003) first mentioned that under certain circumstances, it may be a good practice to quantify the signal produced by the wake; further, the theoretical study of Ianniello (Reference Ianniello2016) revealed that in rotor acoustics, the sound level decay rate with the distance from the source
$r$
, associated with the linear terms, in the near field is more pronounced than the expected
$ r^{-1} - r^{-2}$
laws, due to destructive kinematic auto-interaction. This phenomenon is more evident for low-Mach-number flows and as the number of rotating sources increases. Therefore, in the near-to-intermediate field, the nonlinear noise source may play a non-trivial role. The study carried out by Ianniello (Reference Ianniello2016) was confirmed in the works done by Cianferra et al. (Reference Cianferra, Petronio and Armenio2019b
) and Cianferra & Armenio (Reference Cianferra and Armenio2021) considering a marine propeller, reporting that the relative importance between dipole (linear) and quadrupole (nonlinear) terms is quite complex and that the latter may get dominant at an intermediate distance from the source. Even in aero-acoustics, the nonlinear source term may play a very important role (see, among others, Turner & Kim Reference Turner and Kim2022 and Zamponi et al. Reference Zamponi, Avallone, Ragni, Schram and Van Der Zwaag2024).
It’s worth mentioning that the characterisation of the nonlinear noise sources needs the adequate resolution of the wake, requiring expensive eddy-resolving numerical simulation (see the review paper of Wang, Freund & Lele Reference Wang, Freund and Lele2006). Moreover, the problem of compressibility delays (for a detailed discussion, see Cianferra et al. Reference Cianferra, Ianniello and Armenio2019a ) is more pronounced in the evaluation of the volume terms rather than the linear ones. These issues pose significant complications in the computational evaluation of the nonlinear terms, as addressed later in this paper. Several authors analysed the sound radiated by rotor devices, accounting for both the linear and nonlinear noise sources. Among others, Bensow & Liefvendahl (Reference Bensow and Liefvendahl2016), Stark & Shi (Reference Stark and Shi2021) and Kimmerl, Mertes & Abdel-Maksoud (Reference Kimmerl, Mertes and Abdel-Maksoud2021) used the porous formulation (Di Francescantonio Reference Di Francescantonio1997) and successive modifications (Zhou et al. Reference Zhou, Wang, Wang and He2021; Zhou et al. Reference Zhou, Zang, Wang and Wang2022) to study the hydroacoustics of marine propellers. Recently, high-fidelity numerical simulations based on the use of LES and the FWH equation were performed to analyse the noise emitted by rotors, especially for hydroacoustic purposes (see, for example, Posa, Broglia & Felli Reference Posa, Broglia and Felli2022b and successive analysis).
In contrast, the aeroacoustics characterisation of wind turbines, despite its own importance for practical implications, has not been exploited much considering high-fidelity numerical simulations.
In the present paper, we perform the acoustic analysis of a model-scale wind turbine. The study is carried out using wall-resolving LES in conjunction with the FWH equation. Specifically, we consider the rotor of the wind turbine analysed in the laboratory experiment of Gambuzza, Huck & Ganapathisubramani (Reference Gambuzza, Huck and Ganapathisubramani2023), at a Reynolds number of
$96\, 000$
, based on the wind turbine diameter and inflow velocity. We consider a numerical set-up conceptually similar to that used by Kumar & Mahesh (Reference Kumar and Mahesh2017) to reproduce the relevant features of the wake. To our knowledge, the most accurate spatio-temporal resolution of a wind turbine wake was undertaken by Posa et al. (Reference Posa, Broglia and Balaras2021), considering the immersed boundary method paired with LES to bring the wake evolution to
$7D$
. Here, we perform geometry-resolved simulations in conjunction with wall-resolved LES to reproduce the evolution of wind turbine wake up to
$8D$
and to perform acoustic analysis, at a level of accuracy large enough to correctly reproduce the underlying physics.
First, we perform the numerical simulation and analyse the fluid-dynamic field, also in view of the relevant literature. Successively, we evaluate the noise produced by the rotor using the acoustic analogy and the advective form of the FWH equation. The LES of the model-scale wind turbine allows detection of the main features of the wind turbine wake, namely tip and root vortices, as well as the thin trailing edge vortex sheet released by the blade trailing edge, and their time–space evolution from the near to the far wake. This level of accuracy can highlight the instability mechanism that involves the wake of the wind turbine. In addition, the resolution of the boundary layer that develops near the wind turbine blades provides us with the generation mechanism and transport of the thin vortex sheet released by the trailing edge of the blade, which is not correctly rated by the actuator disc method. This methodology gives us detailed information on the pressure fluctuations on wind turbine blades and velocity fluctuations along the wake, which is fundamental for acoustic purposes.
The use of such methodology allows us to understand better the mechanism of generation and propagation of sound generated by a wind turbine.
The article is structured as follows: § 2 describes the mathematical model; § 3 contains the results of the fluid dynamic field whereas the acoustic analysis is in § 4; and § 5 contains concluding remarks.
2. Methodology
2.1. Fluid dynamic model
The incompressible regime generally applies to the flow conditions of wind turbine aerodynamics. The fluid flow is governed by the mass and momentum conservation laws, and the flow variables are pressure and velocity. The fluid dynamic simulation is here performed using LES; in this context, we consider the filtered incompressible Navier–Stokes equations, in which the large, energetic and noisy scales of motion are directly resolved, and the small isotropic and dissipative structures are modelled using a subgrid-scale closure. The multiple reference frame (MRF) is used to take into account the rotational effects of the isolated wind turbine. The equation set reads as
\begin{align} \frac {\partial \overline {u}_i}{\partial x_i} = 0,& \end{align}
\begin{align} \frac {\partial \overline {u}_i}{\partial t} + \frac {\partial (\overline {u}_i\, \overline {u}_{\!j})}{\partial x_{\!j}} + 2 \epsilon _{\textit{ijk}} \omega _{\!j} \overline {u}_k + \epsilon _{\textit{ijk}} \omega _{\!j} (\epsilon _{\textit{ijk}} \omega _i r_{\!j}) &= - \frac {\partial \overline p}{\partial x_i} + \nu \frac {\partial ^2 \overline {u}_i}{\partial x_{\!j} \partial x_{\!j}} + \frac {\partial \tau _{\textit{ij}}}{\partial x_{\!j}}, \end{align}
where
$x_i$
is the spatial coordinate in the
$i$
-direction,
$\overline {u}_i$
and
$\overline {p}$
are the filtered inertial velocity component in the inertial frame of reference (as in Kumar & Mahesh Reference Kumar and Mahesh2017) and kinematic pressure, respectively,
$\omega _{\!j}$
is the angular velocity component of the rotating frame,
$\nu$
is the kinematic viscosity, and
$\epsilon _{\textit{ijk}}$
is the Levi–Civita permutation tensor. The latter two terms on the left-hand side are the Coriolis and the centrifugal force, respectively. Finally,
$\tau _{\textit{ij}}$
is the subgrid stress tensor (SGS) herein modelled using the wall-adaptive local eddy viscosity (WALE) turbulence model developed by Ducros, Nicoud & Poinsot (Reference Ducros, Nicoud and Poinsot1998). The equations are solved using the Pressure Implicit with Splitting Operators (PISO) numerical algorithm implemented in the open-source framework OpenFOAM (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998) for an unstructured finite-volume grid. We use the backward numerical scheme for time derivatives, which is implicit and second-order accurate. For spatial discretisation, we use the second-order Gauss integration scheme with linear interpolation (central differencing) for all terms. The momentum equation is solved using the bi-conjugate gradient preconditioner (Van der Vorst Reference Van der Vorst1992), while the predictor approach is used to ensure the continuity, and the Poisson equation is solved using the geometric agglomerated algebraic multigrid (GAMG) method. We use two prediction correctors and two non-orthogonal correctors in the overall algorithm. In this way, we ensure that the solution is second-order accurate in both space and time. The computation is performed considering a constant time step
$\Delta t = 1.0 \times 10^{-6} \,{\mathrm{s}}$
, which corresponds to a maximum Courant number of approximately
$Co = 10$
; this value is reached in a few grid cells in the region of the tip of the blades, whereas on average, the Courant number is of the order of
$0.017$
.
We first run a potential solver to obtain a zero divergence velocity field; after that, we employ the Reynolds-averaged Navier–Stokes (RANS) equations to develop a mean turbulent fluid dynamic field; finally, we run the large eddy simulation starting from the RANS field. To ensure that our flow field is statistically steady, we placed probes into the wake and recorded the time series of the pressure and velocity field. We analyse the time records and the relative time spectra to ensure that our numerical grid is fine enough to resolve the most energetic scales of motion. The computation of the flow field statistics is based on a total of
$10$
periods of revolution. For the present analysis, data were collected after
$31$
periods of revolution, after which the fluid flow reached a statistically steady-state condition.
2.2. Acoustic model
The acoustic pressure is evaluated using the advective form of the FWH equation (see Najafi-Yazdi, Brès & Mongeau Reference Najafi-Yazdi, Brès and Mongeau2011; Cianferra et al. Reference Cianferra, Ianniello and Armenio2019a
for the linear and nonlinear term, respectively), suited for wind tunnel-like numerical experiments, considering microphones at a fixed distance from the source. This approach allows us to account for the fact that acoustic wavefronts may deform from their spherical shape due to advection velocity
$U_\infty$
. The advective FWH equation, with
$i=1$
the direction of advection, reads as
\begin{align} 4 \pi \hat p(\boldsymbol{x} , t) = \frac {\partial }{\partial t}\int _S \left [(1-M_0) \frac {Q_i n_i}{R^*}\right ]_\tau {\rm d}S - \int _S \left [U_{\infty } \frac {\tilde {R_1 ^*} Q_i n_i}{{R^*}^2}\right ]_\tau {\rm d}S &\nonumber\\ + \frac {1}{c_0} \frac {\partial }{\partial t}\int _S \left [\frac {L_{\textit{ij}} n_{\!j} \tilde {R_i} }{R^*} \right ]_\tau {\rm d}S + \int _S \left [\frac {L_{\textit{ij}} n_{\!j} \tilde {{R_i}^*} }{{R^*}^2} \right ]_\tau {\rm d}S &\nonumber\\ + \frac {1}{c_0^2}\frac {\partial ^2}{\partial t^2} \int _V \left [T_{\textit{ij}}\frac { \tilde {{R_i}} \tilde {{R_{\!j}}}}{R^*}\right ]_\tau {\rm d}V &\nonumber\\ + \frac {1}{c_0}\frac {\partial }{\partial t} \int _V \left [T_{\textit{ij}}\left (\frac {2\tilde {{R_i}}\tilde {R_{\!j} ^*}}{{R^*}^2} + \frac {\tilde {R_i ^*}\tilde {R_{\!j} ^*} - \hat {R_{\textit{ij}}}}{{R^*}^2}\right )\right ]_\tau {\rm d}V &\nonumber\\ + \int _V \left [T_{\textit{ij}}\frac {3\tilde {R_i ^*}\tilde {R_{\!j} ^*} - \hat {R_{\textit{ij}}}}{{R^*}^3}\right ]_\tau {\rm d}V,\end{align}
where
$M_0=U_\infty /c_0$
is the advective Mach number, in which
$U_\infty$
is the free stream velocity and
$c_0$
is the speed of sound. Additionally,
$n_i$
is the unit vector normal to the surface,
$Q_i = \rho v_i$
with
$v_i$
the surface velocity,
$L_{\textit{ij}} = (p-p_0)\delta _{\textit{ij}} + \sigma _{\textit{ij}}$
, where
$\delta _{\textit{ij}}$
and
$\sigma _{\textit{ij}}$
are respectively the second-order isotropic tensor and the viscous stress tensor, where it is considered
$p_0 = 0$
Pa and
$\rho _0 = 1.0 \,{\mathrm{kg} \, \mathrm{m}^{-3}}$
, and
$T_{\textit{ij}}= \rho u_i u_{\!j} +[(p-p_0)-c_0 ^2 (\rho -\rho _0)]\delta _{\textit{ij}}-\sigma _{\textit{ij}}$
is the Lighthill tensor, here approximated as
$T_{\textit{ij}}= \rho u_i u_{\!j}$
. This holds for low-Mach-number flows (Lighthill Reference Lighthill1952) since viscous effects are negligible and in the incompressible barotropic flow regime, the term
$(p-p_0)-c_0 ^2 (\rho -\rho _0)$
is zero. Therefore, the significant contribution to the Lighthill term is given by the velocity field only. Furthermore,
$S$
is the integration surface, here coincident with the wind turbine surface and
$V$
is the volume of integration containing the turbine wake. For an observer at position
$(x_i)$
and a source at position
$(y_i)$
, the radiation vectors are given by
\begin{align}&\quad\,\,\, {R^*} = \sqrt {(x_1 - y_1)^2 + \beta ^2 \left[(x_2-y_2)^2 + (x_3-y_3)^2\right]}, \end{align}
\begin{align}&\,\,\,\, \tilde {R_1 ^*} = \frac {x_1-y_1}{R^*},\quad \tilde {R_2 ^*} = \beta ^2 \frac {x_2-y_2}{R^*}, \quad \tilde {R_3 ^*} = \beta ^2 \frac {x_3-y_3}{R^*}, \end{align}
\begin{align}& \tilde {{R}_1} = \frac {1}{\beta ^2}\big(-M_0 + \tilde {R^*}_1\big), \quad \tilde {{R}_2} = \frac {x_2-y_2}{R^*}, \quad \tilde {{R}_3} = \frac {x_3-y_3}{R^*}, \end{align}
with
$\beta = \sqrt {1-M_0 ^2}$
and
$\hat R_{\textit{ij}}$
defined as
\begin{align} \begin{bmatrix} 1, & 0, & 0\\ 0, & \beta ^2, & 0 \\ 0, & 0, & \beta ^2 \end{bmatrix} \!.\end{align}
All integral kernels must be evaluated at the emission time
$\tau$
, and this may result in complications when dealing with a large amount of data. A discussion on this matter is in § 4. In (2.3), the first and the second integral terms (thickness noise component) represent the noise generated by the displacement of the fluid mass caused by the movement of the body, whereas the third and fourth terms depend directly on the dynamic pressure over the solid surfaces (loading noise component). The three-volume integrals of (2.3) are known as quadrupole noise terms and consider all possible nonlinear sources occurring in the flow field.
2.3. Numerical set-up
We consider the model-scale wind turbine used by Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023) in their experimental analysis. It is a three-blade, fixed-pitch, controlled-speed wind turbine rotating in the clockwise direction. The diameter is
$D=0.18 \ {\mathrm{m}}$
. Figure 1 shows the blade and hub geometry, and in table 1, we report geometric features in terms of aerofoil profile, radial coordinate (
$r/R$
) with
$R = D / 2$
, non-dimensional cord (
$c/R$
) and twist (
$\beta$
). The hub is not perfectly the same as that of the experiment because the geometry specification is not reported in the paper describing the aforementioned experiment. The uniform free stream velocity is
$U_{\infty } = 8$
m s−1 and the diameter-based Reynolds number is equal to
\begin{equation} {\textit{Re}}_D = \frac {U_{\infty } D} {\nu } = {9.6\times 10^4}, \end{equation}
which corresponds to a maximum cord-based Reynolds number:
\begin{equation} {\textit{Re}}_{c,\\max} = \left .\frac { c \sqrt {U_{\infty }^2+(\varOmega r)^2} } {\nu } \right \vert _{\textit{max}} = {3.4\times 10^4}, \end{equation}
where
$\varOmega = 417.17$
rad s−1 is the wind turbine angular velocity. As reported by Chamorro et al. (Reference Chamorro, Guala, Arndt and Sotiropoulos2012) and Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023), this value of the Reynolds number is large enough to attain Reynolds-independent results in terms of wake velocity signals. The tip speed ratio is
$\textit{TSR} = \varOmega R/U_{\infty }=4.7$
, which is the value at which the maximum thrust is generated (Gambuzza et al. Reference Gambuzza, Huck and Ganapathisubramani2023); this gives a high load condition that generates strong root and tip vortices.
Wind turbine geometry, defined as aerofoil type, radius position, cord and twist.
(a) Wind turbine blade geometry; (b) zoom-in visualisation of the hub geometry.
To minimise the blockage effect caused by the lateral boundaries of the computational domain, we use the same conceptual framework as Kumar & Mahesh (Reference Kumar and Mahesh2017), presented in figure 2. Actually, we use a cylindrical domain with a radius equal to
$7D$
and a length of
$10D$
; the wind turbine is positioned at
$2D$
from the inlet section. Hereafter, we use interchangeably
$x_1,x_2,x_3$
or
$x,y,z$
for the Cartesian frame of reference. The reference coordinate system is placed at the centre of the wind turbine, and the free stream fluid flow moves along the
$x_1$
or
$x$
direction. The domain extends up to
$8D$
downstream the wind turbine, which is sufficient to develop the wake, from the near to the far field.
(a) Sketch of the three-dimensional numerical domain together with the Cartesian frame of reference; (b) schematic of the domain extension; (c) schematic of the velocity boundary condition on the wind turbine (
$\varOmega R /U_{\infty }$
).
The overall computational domain is composed of approximately 144 million hexa-dominant control volumes. The wind turbine surface is discretised with approximately
$5 \times 10^5$
quad and tria surface elements; detailed visualisation of the surface grid is reported in figure 3. The surface elements have homogeneous dimensions in the two directions, with a maximum value of
$\varDelta \kern-0.5pt/\kern-2pt D \approx 9 \times 10^{-4}$
and a minimum value equal to
$\varDelta \kern-0.5pt/\kern-2pt D \approx 2.8 \times 10^{-4}$
in the regions of large curvature, such as leading-edge, the tip of the blade and in correspondence of the thin trailing edge.
Detail of the surface mesh: (a) the root zone; (b) the blade’s leading edge; (c) the blade’s tip.
The wind turbine surface mesh is extruded with 28 layers composed of hexahedral elements, as depicted in figure 4. After defining a local frame of reference, with
$x_b$
aligned with the cord of the blade,
$y_b$
normal to the blade surface and
$z_b$
spanwise along the blade axis, the non-dimensional height of the first cell off-the-wall is
$\Delta y_b \kern-0.5pt/\kern-2pt D \approx 4.4 \times 10^{-5}$
, which, in wall units, gives the first centroid off the wall at a distance
$y_b^+\lt 1$
. The wall unit is
$\nu / u_{\tau }$
, with
$u_{\tau } = \sqrt {\tau _w / \rho }$
, where the wall shear stress is estimated as
$\tau _w =0.5 C_{\!f} U_{\infty }^2$
and the friction coefficient
$C_{\!f} = [2\log_{10} {\textit{Re}}_{c,max} -0.65]^{-2.3}$
is calculated using the formula for a turbulent flow over a semi-infinite smooth flat plate. To verify fulfilment of the condition
$y_b^+\lt 1$
at the wall, the
$y_b^+$
value was the calculated run-time during the simulation and it has been observed that it is
$O(1)$
in the trailing edge region of the blades at the outer radii only, whereas, in the rest of the surfaces, the values of
$y_b^+$
are below unity (not shown). The maximum values of
$\Delta x_b^+$
and
$\Delta z_b^+$
over the blades were verified to be smaller than
$15$
in the regions of low curvature, and approximately
$5$
in the tip and trailing edge regions. For the first layer of cells around the blades (red area of figure 4), the cell-to-cell expansion ratio is set equal to
$1.1$
. In this way, the latest cell dimension of the first extruded layer is
$\Delta x_b \kern-0.5pt/\kern-2pt D \approx \Delta y_b \kern-0.5pt/\kern-2pt D \approx \Delta z_b \kern-0.5pt/\kern-2pt D \approx 4.4 \times 10^{-4}$
(at the border between the red and grey regions of figure 4).
Representation of the near-wall layers.
The rest of the grid is illustrated in figure 5. To ensure a smooth transition between the wall region and the wake region, a refinement box of radius
$R_3 = 1.1 R$
was set, extending from
$-0.07 D$
to
$0.07 D$
with a cell dimension of
$\Delta x / D \approx \Delta y / D \approx \Delta z / D \approx 0.0027$
. Finally, the external region enclosing the turbine and containing the wake comprises a conical box with a grid resolution of
$\Delta x / D \approx \Delta y / D \approx \Delta z / D \approx 0.0055$
. The conical box has an upstream radius
$R_1 = 1.33 R$
(located in front of the wind turbine) and ends with a downstream radius of
$R_2 = 1.78 R$
at
$3D$
; afterward, the box becomes cylindrical with diameter
$R_2$
throughout the domain extension.
Refinement boxes used to resolve the turbulent structures of the wind turbine wake.
The boundary conditions are the following: at the inlet, we consider a uniform velocity in the streamwise direction
$U_\infty$
; at the outlet, we set the zero gradient condition; the slip condition is used at the lateral boundary, while the rotational velocity
$\boldsymbol{u} = {\varOmega }\boldsymbol{r}$
is set over the surfaces of the turbine (figure 2), such as blades and nacelle. A reference pressure value
$p_0 = 0 \, {\mathrm{Pa}}$
is set at the outlet while the zero gradient condition is imposed at the remaining boundaries.
For the acoustic analysis, we set a value of the speed of sound equal to
$c_0 = 330 \, {\mathrm{m}}\,{\mathrm{s}}^{-1}$
, which gives an inlet Mach number
$U_\infty /c_0 \sim 0.02$
and a tip Mach number defined as
$M_{\textit{tip}} = \sqrt {(\varOmega R)^2+U_{\infty }^2} /c_0 \sim 0.11$
; this ensures Mach similarity (see Cianferra & Armenio Reference Cianferra and Armenio2021 for details) considering, for instance, a medium-size turbine with diameter of
$ 100\, {\mathrm{m}}$
rotating with a period of
$ 8.42\, {\mathrm{s}}$
, in the presence of a wind speed equal to
$U_\infty$
.
For the analysis of sound propagation in the axial direction, we consider nine microphones placed at radial distance
$r / D = 1.0$
and at streamwise positions equal to
$x / D = 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0$
. The radial evolution of the acoustic field is evaluated considering microphones placed in the plane of the rotor (
$x / D = 0.0$
) at four radial positions,
$r / D = 1.0, 2.0, 4.0, 8.0$
. The integration volume for the nonlinear terms is a cylinder with radius
$r_v/D = 0.9$
and extending from
$x / D = -0.5$
to
$x / D = 7.0$
. The integration volume, together with the position of the microphones, is illustrated in figure 6. A sensitivity analysis on the acoustic results varying the radius of the integration volume has been conducted considering
$r_v/D = 0.7,\,0.8, \,0.9$
;
$r_v/D = 0.7$
is the minimum radius enclosing all the fluid dynamics sources. The analysis (not reported here) has shown that the results are independent of the radius of the integration volume.
Microphone positions, integration volume for the nonlinear terms (green cylinder) and wind turbine surface (red surface).
3. Results and discussion: fluid dynamic analysis
We validate our results with the experimental data of Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023) for the same TSR; to the best of our knowledge, no additional literature data are available for the wind turbine considered herein. In their laboratory experiments, the wind turbine was connected to a permanent-magnet brushed DC generator. An aluminium mast with a diameter equal to
$15.75\, {\mathrm{m}}{\mathrm{m}}$
and a height of approximately
$1.5D$
sustained the wind turbine nacelle. The tests were carried out in a wind tunnel with a blockage coefficient
$A_d/C$
equal to
$ 4.7 \,\%$
(
$A_d$
is the disk area of the wind turbine model;
$C$
is the wind tunnel cross-section). In their analysis, the authors considered the laminar inflow case as well as the turbulent ones. Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023) reported the power coefficient for the laminar inflow measurements, while the thrust coefficient is provided for the turbulent inflow cases only. The latter coefficient was found to be less affected by the free stream turbulence, compared with the former. Here, we calculate the thrust and power coefficients and compare them with the data of Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023). The coefficients are defined as
\begin{align}& C_T = \frac {T_r}{\dfrac {1}{2}\rho U_{\infty }^2 \pi R^2}, \\[-16pt] \nonumber \end{align}
\begin{align}& C_{\!P} = \frac {Q \omega }{\dfrac {1}{2}\rho U_{\infty }^3 \pi R^2}, \end{align}
in which
$T_r$
is the thrust and
$Q$
is the torque. The quantities
$T_r$
and
$Q$
are calculated considering the mean value over
$10 T$
, where
$T$
is the wind turbine revolution period.
The relative error, reported in table 2, is
$11.25\, \%$
and
$8.50 \,\%$
for the thrust and power coefficients, respectively. The percentage errors are computed as follows:
\begin{align} \epsilon _T\,\% &= \frac {\left| C_T^{\textit{LES}}-C_T^{\textit{EXP}}\right| } {C_T^{\textit{EXP}}}\times 100, \end{align}
\begin{align} \epsilon _{\!P}\,\% &= \frac {\left| C_{\!P}^{\textit{LES}}-C_{\!P}^{\textit{EXP}}\right| }{C_{\!P}^{\textit{EXP}}}\times 100, \\[6pt] \nonumber \end{align}
where the indexes
$\textit{LES}$
and
$\textit{EXP}$
refer to LES data and those of the laboratory experiment of Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023). The values of the axial load and the power generated in the physical experiment are slightly higher than those obtained from our simulation. There are possible sources of disagreement since our numerical experiment differs in some aspects from that performed by Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023). In fact, the blockage coefficient of our simulation (
$A_d/C=0.005\, \%$
) is much smaller than that of the laboratory experiment. In the latter case, some small, still not negligible, effects on the forces developed by the wind turbine may be present. However, the most significant difference consists of the presence of the mast in conjunction with the measurement apparatus in the experiments, which generates significant interactions that are not easily quantifiable. In addition, the presence of the wall at the bottom of the wind tunnel promotes the development of a boundary layer, which affects the inflow condition and therefore the wind turbine performance. Overall, the comparison between our LES and the experiment is quite satisfactory. Future research may be devoted to the analysis of a complete wind turbine, comprehensive of the mast, to analyse the effects of the latter on the overall performance of the wind turbine. The time records of the thrust and torque coefficients and the relative spectra are reported in figure 7. It is evident that the loads acting on the wind turbine exhibit a broadband character. The resolution of the LES is fine enough to resolve the load fluctuations over a wide range of frequencies, up to
$100 f/f_{T}$
, where
$f_T$
is the revolution frequency.
Thrust and power coefficients and relative error between present results and data from experiments performed by Gambuzza et al. (Reference Gambuzza, Huck and Ganapathisubramani2023).
Time history of the (a) thrust and (b) power coefficients versus the non-dimensional time
$t/T$
, where
$T$
is the revolution period; spectra of the (c) thrust and (d) power coefficients versus the non-dimensional frequency
$f/f_{T}$
, where
$f_{T}$
is the revolution frequency.
The resolution and mesh quality allow us to resolve the details of the turbulent field. As an example, we report a snapshot of the instantaneous flow field on a longitudinal plane that passes through the turbine axis (figure 8); the instantaneous velocity and pressure coefficient fields highlight the presence of a fully developed composite turbulent wake. The near and far wake are clearly distinguishable. The near wake exhibits the classical system of tip vortices together with the thin trailing edge vortex sheet. Further downstream, the tip vortices begin to lose coherence, and a population of turbulent structures emerges and spreads over the disk, characterising the intermediate-to-far wake. In figure 8(c), we report the SGS eddy viscosity made non-dimensional with the kinematic one (
$\nu _T/\nu$
). The large SGS eddy viscosity provides information on the region where most of the turbulent dissipation takes place. The values are mostly in the range of
$0{-}10$
, indicating the high resolution used in our simulation. They peak in the core of the helicoidal tip vortex and in the thin vortex structures evolving in the wake.
Instantaneous field quantities in a meridian plane: (a) streamwise velocity component made non-dimensional with
$U_\infty$
; (b) pressure coefficient
$C_p = (p-p_0)/0.5\rho U_\infty ^2$
; (c) SGS eddy viscosity made non-dimensional with kinematic viscosity
$\nu _T / \nu$
.
3.1. Wind turbine wake
Here, we discuss the characteristics of the wake that develops from the turbine downstream to
$8D$
. To highlight the coherent structures composing the turbulent wake before their instability and loss of synchronisation with the wind turbine rotation, we compute phase-averaged quantities. The phase averaging operation over a cylindrical frame of reference
$(r, \theta , x)$
for a generic quantity
$\varPsi$
is defined as follows:
\begin{equation} \hat {\varPsi }(r,\theta ,x) = \frac {1}{N} \sum _{i=1}^{N} \varPsi (r, \theta +\varOmega t_i, x, t_i). \end{equation}
In figure 9, we show the phase-averaged axial velocity, vorticity magnitude and turbulent kinetic energy along the entire wind turbine wake in the meridian plane. In addition, in figure 10, we show the snapshot of the three-dimensional wake structures individuated through the isosurface of the non-dimensional second invariant of the velocity gradient tensor
$Q$
:
\begin{equation} \frac {{\textit{QD}}^2}{U_\infty ^2} = \frac {D^2}{2U_\infty ^2}\big(\|W\|^2-\| S\|^2\big), \end{equation}
where
$S$
and
$W$
are the symmetric and deviatoric parts of the velocity gradient tensor, respectively.
Phase-averaged quantities in the meridian plane passing through the wind turbine axis: (a) axial velocity field normalised with
$U_\infty$
; (b) vorticity magnitude normalised with
$U_\infty /R$
; (c) turbulent kinetic energy (TKE) normalised with
$U_\infty ^2$
.
Snapshot of instantaneous isosurface
${\textit{QD}}^2/U_ \infty ^2=12.66$
.
The wind turbine placed in a uniform laminar flow develops a complex wake composed of three main regions, as shown in figures 9 and 10. A momentum deficit, more evident in the region between the tip and the hub, together with wake expansion (figure 9
a) are associated with the extraction of kinetic energy by the wind turbine. Behind the hub, a small recirculation zone is present beyond which the velocity in the axial direction increases rapidly to approximately
$U_\infty$
, in a region extending up to
$5D{-}6D$
. Up to roughly
$1.0D$
, the wake exhibits robust coherence, with structures synchronised with the rotation of the wind turbine, characterising the near wake. Beyond that, a decrease in the magnitude of vorticity is observed, which corresponds to an increase of the TKE at the outer radii. This phenomenon can be interpreted as the loss of coherence and synchronisation of the tip vortex system due to instability mechanisms, leading to the formation of turbulent eddies defining the intermediate wake up to approximately
$7D$
. Finally, a significant amount of TKE is spread throughout the remaining part of the wake, reducing the momentum deficit which identifies the transition to the far wake. Figure 10 also shows that in the very far wake, few large-scale structures populate the turbulent field, especially at the centre of the wake, according to the theory of decaying isotropic turbulence. Hereafter, we separately discuss the three different regions of the wake.
3.1.1. The near wake
The near wake extends roughly up to
$1.0D$
, with the presence of organised structures, namely the sheet vortical structures released by the trailing edge of the blades, the root vortices and the helicoidal tip vortex. The Q-criterion (figure 10) shows these structures, with the maximum value of vorticity observed in the tip vortex core (figure 9
b). As expected, in the near wake, the tip vortex is stable and coherent exhibiting the typical helical path. The turbulent boundary layer that develops over the blades generates a thin trailing edge vortex sheet, detectable downstream of the wind turbine in figures 9(b) and 10 as a thin sheet of intense vorticity. The root vortices of a wind turbine are composed of several coherent structures released at the root of the blades, interacting with the nacelle boundary layer, as discussed by Sherry et al. (Reference Sherry, Nemes, Lo Jacono, Blackburn and Sheridan2013a
,Reference Sherry, Sheridan and Jacono
b
). This is evident in our simulation in figures 9(b) and 10, where it is possible to detect the interaction between the root vortex structures at approximately
$0.2R$
and the vorticity downstream the nacelle; this involved vorticity system is more susceptible to destabilisation and viscous dissipation compared with that of a marine propeller (see Kumar & Mahesh Reference Kumar and Mahesh2017 and Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2019 for a discussion).
There is a clear difference between the isolated marine propeller and the wind turbine. In the former case, a strong hub vortex is present (see, among others, Felli et al. Reference Felli, Camussi and Di Felice2011; Kumar & Mahesh Reference Kumar and Mahesh2017; Posa et al. Reference Posa, Broglia and Felli2022b ), which, moving downstream, eventually may start to oscillate and to interact with the tip vortex structures, generating a complex vortex interaction and breakdown (see, for example, figures 22 and 23 of Felli et al. Reference Felli, Camussi and Di Felice2011). In contrast, in the case of the wind turbine rotor, the hub vortex is insignificant (see figure 9), and the mechanism of interaction between the hub vortex and the tip vortex systems is clearly absent. This is probably due to significant geometry differences, among them: the ratio between the hub radius and the wind turbine radius is much smaller than in the case of a marine propeller; the wind turbine blade is very elongated with a very large aspect ratio, which is not the case for the marine propeller; the ratio between the expanded area and the disc area is typically larger than 1 for the marine propeller and much smaller than 1 for the three-blade wind turbine; finally, the wind turbine is less bulky than a marine propeller, especially in the root/hub region, promoting a central region in the wake where the mean velocity is almost undisturbed (see figure 9 a).
In general, in a wind turbine wake, the thin trailing edge vortex sheet and the root vortices dissipate rapidly and are weaker than the system of tip vortices (Lignarolo et al. Reference Lignarolo, Ragni, Scarano, Ferreira and van Bussel2015; Posa et al. Reference Posa, Broglia and Balaras2021). Figure 11 shows the contour of the phase-averaged axial velocity field, vorticity magnitude and turbulent kinetic energy over cross-sectional planes at increasing distance from the wind turbine, in the very near wake (
$x / D = 0.05,0.2,0.4,0.8$
). The tip vortex system, the coherent filaments composing the root vortices and the blade wakes are distinguishable. Near the wind turbine, the axial velocity peaks in the region of the blade tip; a high level of vorticity is present along the blades with peaks in the region of the tip due to the presence of the tip vortex. The turbulent kinetic energy is large downstream of the trailing edge of the blades due to the presence of the blade trailing edge vortex sheet and in correspondence with the tip vortex. As the distance increases downstream, the streamwise velocity maintains large values in correspondence with the tips of the blades, it keeps increasing around the hub, and the defect velocity appears more homogeneous within an annular region between the hub and the tip of the blades. The blade trailing edge vortex dissipates rapidly with a corresponding decrease of turbulent kinetic energy, whereas the tip vortex dissipates TKE slower. It should be noted that a decreased level of vorticity magnitude does not correspond to a strong distortion of the vorticity field and to an increased level of TKE. In fact, at
$x=0.6D$
, the blade wakes and the coherent helicoidal structures composing the root vortices are weak but still distinguishable, suggesting their rapid dissipation. The tip vortex exhibits gradual diffusion and loss of vorticity intensity with decreasing values of turbulent kinetic energy and remains stable up to
$1.0D$
as depicted in figures 9 and 10. This provides further confirmation that the near wake of the wind turbine is dominated by the tip vortex system, as previously discussed in the literature by, among others, Lignarolo et al. (Reference Lignarolo, Ragni, Scarano, Ferreira and van Bussel2015) and Posa et al. (Reference Posa, Broglia and Balaras2021).
Phase-averaged axial velocity (left column), vorticity magnitude (central column) and turbulent kinetic energy (right column) normalised using
$U_{\infty }$
and
$R$
. From the top to the bottom, each row of panels corresponds at a certain distance downstream of the wind turbine:
$x / D = 0.05$
,
$0.2$
,
$0.4$
,
$0.8$
.
3.1.2. The intermediate wake
Beyond
$1.0D$
, the tip vortex starts losing its own coherence and synchronisation with the wind turbine rotation (figures 9 and 10), due to instabilities as already discussed by Lignarolo et al. (Reference Lignarolo, Ragni, Scarano, Ferreira and van Bussel2015) and Posa et al. (Reference Posa, Broglia and Balaras2021). Specifically, Posa et al. (Reference Posa, Broglia and Balaras2021) showed that the deviation from the helicoidal and stable configuration starts with the short-wave instability with wavelength deviations of the order of the vortex core diameter; after that, the long-wave instability increases the wavelength deviations up to few vortex core diameters and promotes mutual inductance between two successive tip vortices. Our results are consistent with the literature (see figure 12).
Snapshots of isosurfaces
$QD^2/U_\infty ^2=25.00$
computed at four time instants:
$T/4$
,
$T/2$
,
$3T/4$
,
$T$
, and phase-averaged turbulent kinetic energy normalised using
$U_{\infty } ^2$
on the meridian plane.
Instantaneous vorticity magnitude normalised using
$U_{\infty }$
and
$R$
;
$x / D = (a)\, 1.0$
; (b)
$2.0$
; (c)
$4.0$
; (d)
$6.0$
.
As the helicoidal path of the tip vortex starts to exhibit short- and long-wave instabilities, an increase of turbulent kinetic energy is detected; the trajectory of the vortical structures starts to deviate from the helicoidal path, promoting the mutual inductance mode and leap-frogging mechanism. Figure 12 clearly shows the exchange of the relative position of successive tip vortices and vortex grouping mechanisms, as described by Felli et al. (Reference Felli, Camussi and Di Felice2011) for a ship propeller, and finally the three-dimensional breakdown of the helical vortices. Leap-frogging produces the breakdown of large vortical structures into smaller eddies, thus increasing turbulent mixing, resulting in the spread of turbulent kinetic energy over a large flow region and loss of synchronism with the wind turbine revolution. This is also highlighted in figure 13, which shows the instantaneous vorticity field over four significant cross-sectional planes. The exchange of the position between two successive tip vortices is evidenced in figure 13(b), leading to a reduction of vorticity magnitude, a significant distortion of the vorticity field and generation of a wide range of turbulent eddies populating the intermediate-to-far wake (figure 13 c,d). Tip vortex destabilisation is also observable looking at phase-averaged quantities in cross-flow planes at increasing distances from the turbine (figure 14); in particular, figure 14 shows a reduction in vorticity and an increase in turbulent kinetic energy moving downstream in the wake.
In addition, the shear layer present in the tip vortex region promotes turbulent transport of axial momentum (left column of figure 14) and TKE (right column of figure 14) from the tip vortex region along the radial direction.
Phase-averaged axial velocity (left column), vorticity magnitude (central column) and turbulent kinetic energy (right column) normalised using
$U_{\infty }$
and
$R$
. From the top to the bottom, a row corresponds to a certain distance downstream of the wind turbine,
$x / D = 1.0, 2.0, 4.0, 6.0$
.
3.1.3. The far wake
In the far wake, the phase-averaged axial velocity and TKE tend to be smooth and homogeneous over the disk. The velocity fluctuations generated by the tip vortex destabilisation and shear layer in the intermediate part of the wake are transported across the entire cross-section in the far wake. In other words, turbulence progressively propagates from the tip region to the interior (figure 9
c) with significant homogenisation over the disc of diameter roughly equal to
$D$
. As for the other regions of the wake, here, we show the phase-averaged streamwise velocity, the vorticity magnitude and the turbulent kinetic energy (figure 15) at increasing distances from the turbine (
$x / D = 7.0, 8.0$
).
The phase-averaged velocity field in the axial direction becomes increasingly more homogeneous, due to the turbulent mixing and momentum wake recovery (left column of figure 15). The low value of vorticity magnitude, compared with the other regions of the wake, shows loss of coherence in the tip vortices and the absence of large, energetic vortical structures synchronised with the wind turbine’s rotation (central column of figure 15).
Finally, TKE appears distributed across the section with a substantially smaller magnitude than in the other parts of the wake (right column of figure 15). This indicates the presence of a population of smaller-scale turbulent eddies in the far wake, covering the entire disk as depicted in figure 10.
Phase-averaged axial velocity (left column), vorticity magnitude (central column) and turbulent kinetic energy (right column) normalised using
$U_{\infty }$
and
$R$
. From the top to the bottom, a row is relative to a certain distance downstream of the wind turbine:
$x / D = 7.0, 7.5$
.
4. Results and discussion: acoustic analysis
4.1. Noise sources
As discussed in § 2.2, the FWH equation (2.3) contains seven terms that represent the acoustic sources associated with different fluid-dynamic processes. The thickness noise source is related to the rotational velocity (
$\omega \boldsymbol{r}$
), which is constant in time. However, a time dependency is introduced because of the periodic variation of the relative position between the microphone and the solid surfaces, leading to a tonal noise at the blade passing frequency (
$f_{bp}$
). As for the thickness term, the loading term peaks at the blade passing frequency, although its own magnitude also depends on the pressure fluctuations and their time derivative over the surfaces. These fluctuations contribute to the high-frequency noise, as discussed later in the paper. To highlight the noise source originated by the pressure field over the blades, we analyse the pressure fluctuation coefficient (
$C_{\!p_{\textit{rms}}} = p_{\textit{rms}} / 0.5\rho U_{\infty }^2$
) acting on the blade’s surface (figure 16). We observe that the downwind side of the blade’s surface is affected by pressure fluctuations more than the upwind blade’s surface; these fluctuations are associated with the transition from laminar to turbulent boundary layer and flow separation. Significant pressure fluctuations are present in the trailing edge and moving towards the tip of the blade; this may be ascribed to the increased value of the cord-based Reynolds number. At the tip of the blade (figure 16
c), the release of the tip vortex promotes pressure fluctuations which are weaker compared with the corresponding trailing edge fluctuations (compare figures 16
a and 16
c). As discussed by Wagner et al. (Reference Wagner, Bareiss and Guidati1996), the trailing edge noise source is among the major noise sources in wind turbine aeroacoustics. In particular, the trailing-edge noise source dominates the high-frequency spectrum, perceived as a swishing sound.
Phase-averaged pressure fluctuation coefficient on the blade’s surface: downwind side (a); upwind side (b); tip of the blade (c).
Phase-averaged r.m.s. of the Lighthill source term made non-dimensional with
$\varOmega$
: (a) meridian plane; transversal planes located at
$x / D=$
(b) 0.5; (c) 1.0; (d) 2.0; (e) 3.0; (f) 5.0; (g) 7.5.
The analysis of the kinematic Lighthill source term in the fluid region,
\begin{equation} S = \frac {\partial ^2 u_i u_{\!j}} {\partial x_i \partial x_{\!j}}, \end{equation}
allows us to quantify the noise source associated with the wake. In figure 17, we show the phase-averaged root mean square (r.m.s.) of the Lighthill source term (
$S_{\textit{rms}}$
) along the whole wind turbine wake, in the meridian plane and at cross-sectional planes at increasing distance from the turbine. This quantity gives information on the noise sources mainly at the blade passing frequency. Higher values of
$S_{\textit{rms}}$
are detected in the near wake corresponding respectively to: the thin trailing edge vortex sheet released by each blade; the coherent structures composing the root vortices; the system of tip vortices syncronised with the blade passing frequency (figure 17
a,b). Similar results were obtained by Posa et al. (Reference Posa, Broglia and Felli2022b
) for the case of marine propellers. The
$S_{\textit{rms}}$
rapidly decreases in the blade wake (figures 17
b and 17
c), consistently with a strong decay of vorticity and TKE in the same region (figure 11). In the hub and tip regions, we observe large values of
$S_{\textit{rms}}$
which decrease rapidly moving downward (figures 17
a and 17
c–
e). To be noted that, beyond
$1D$
, noise is generated mainly by the tip vortex whose dynamics is described in the previous section. Finally,
$S_{\textit{rms}}$
decreases in intensity further downstream over the whole cross-sectional area (figures 17
a and 17
e–g), indicating that the near-to-intermediate wake contributes more to the production of noise. In agreement with the results of Posa et al. (Reference Posa, Broglia and Felli2022b
) for ship propellers, large values of
$S_{\textit{rms}}$
are found in regions of large coherence of the fluid flow, such as the tip and root vortices. However, in contrast to the case of marine propellers, the region of large fluctuations persists at the outer radii more than in the inner region, suggesting that the acoustic signature of the wind turbine is dominated to a large extent by tip vortices and their instability rather than by the hub vortex (which is insignificant as already discuss in the previous section) and the root vortices as the distance downstream increases. This is due to the geometrical differences between marine propellers and wind turbines, as we already pointed out in § 3.1.1.
4.2. Acoustic spectra
This section presents the analysis of the wind turbine’s acoustic signature. We show the streamwise and radial evolution of the sound pressure level (SPL) expressed in dB. The SPL is defined as
\begin{equation} \textit{SPL} = 20 \log _{10} \left (\frac {\textit{fft}(\hat p)}{ \hat p_0}\right )\!, \end{equation}
where
$\textit{fft}$
is the Fast-Fourier transform of
$\hat p$
and
$\hat p_0 = 20 \,\unicode{x03BC} {\mathrm{Pa}}$
is the reference sound pressure level in air. As discussed in § 2.2, the integrals that compose the FWH equation must be evaluated at the emission time
$\tau$
, which may introduce a significant increase in computational efforts. Cianferra et al. (Reference Cianferra, Ianniello and Armenio2019a
) introduced the maximum frequency parameter (MFP),
\begin{equation} \textit{MFP} = \frac {1}{\varDelta _{\textit{del}} f_{\textit{max}} } = 1, \end{equation}
to evaluate the maximum frequency
$f_{\textit{max}}$
beyond which neglecting the time delays may have an impact on
$\hat p$
. The maximum difference in time delay
$\varDelta _{\textit{del}}$
is defined as
\begin{equation} \varDelta _{\textit{del}} =\frac { \max \limits _{\boldsymbol{y}\in S,V}(\boldsymbol{y}(\tau ) - \boldsymbol{x}(t)) - \min \limits _{\boldsymbol{y}\in S,V}(\boldsymbol{y}(\tau ) - \boldsymbol{x}(t))}{c_0}, \end{equation}
where
$\boldsymbol{x}(t)$
and
$\boldsymbol{y}(\tau )$
are the microphone position at the reception time
$t$
and the position of the acoustic source at the emission time
$\tau$
, respectively. The criterion
$\textit{MFP}\gt 1$
gives the range of frequencies
$f\lt f_{\textit{max}}$
, where the computation of the time delays is not required. The MFP criterion provides a quantification of the general definition of a compact source, identifying an exact frequency threshold within which computation of time delays can be neglected. Note that, in the case of sources fixed in space,
$\boldsymbol{y}(\tau )$
is constant and coincides with grid nodes, making the computation of
$\varDelta _{\textit{del}}$
straightforward. However, for bodies in motion,
$\boldsymbol{y}(\tau )$
is unknown and can be computed by solving an implicit equation
\begin{equation} t - \tau - \frac {\|\boldsymbol{x}(t) - \boldsymbol{y}(\tau )\|}{c_0} = 0. \end{equation}
This equation is also adopted to compute the retarded surface defined by Ianniello, Muscari & Di Mascio (Reference Ianniello, Muscari and Di Mascio2013) and based on the retarded time approach of Brentner & Farassat (Reference Brentner and Farassat2003). We solve it numerically using the bisection method to obtain the emission times
$\tau$
, the related positions
$\boldsymbol{y}(\tau )$
and finally
$\varDelta _{\textit{del}}$
. The values of
$f_{\textit{max}}$
are then computed for each microphone for the surface elements (table 3) and for the volume elements (table 4). Note that a rough estimation of the parameter
$\varDelta _{\textit{del}}$
is given by
$D/c_0 = 5.4 \times 10^{-4} \ {\mathrm{s}}$
for the linear part and
$7D/c_0 = 0.0038 \ {\mathrm{s}}$
for the nonlinear part, for the less favourable condition (corresponding to the microphone closest to the turbine); this thumb rule estimation is consistent with the results obtained with the use of (4.4) and (4.5).
$\varDelta _{\textit{del}}$
and non-dimensional maximum frequency
$f_{\textit{max}}/\!f_T$
, where
$f_T$
is the revolution frequency, that can be resolved correctly neglecting the time delays for the linear terms of the FWH equation at different microphones.
$\varDelta _{\textit{del}}$
and non-dimensional maximum frequency
$f_{\textit{max}}/\!f_T$
, where
$f_T$
is the revolution frequency, that can be resolved correctly neglecting the time delays for the nonlinear terms of the FWH equation at different microphones.
The minimum value of
$f_{\textit{max}}$
that can be solved without introducing significant errors is obtained for the microphone placed in the very near field at
$(0,0,D)$
; its value is
$30 f_T$
for the linear terms and
$3.85 f_T$
for the nonlinear ones, where
$f_T$
is the revolution frequency. Then, the value of
$f_{\textit{max}}$
increases with distance from the source. Note that the compressibility delay is more significant for the quadrupole terms than for the linear contribution, as expected. Here, we do not consider the compressibility delay in the computation of the acoustic pressure, since it is computationally expensive, in particular, for the quadrupole terms of the FWH equation. However, in our specific case, we show that the nonlinear terms are significant in the lower frequency range only (
$f \lt 4.0 f/f_T$
), which makes our assumption acceptable.
We first calculate separately the linear and nonlinear contributions to the SPL of acoustic pressure
$\hat {p}$
(figure 18) at microphones placed at the radial distance
$r / D = 1.0$
at four streamwise positions, that is,
$x / D = 0.0, 0.5, 1.0, 2.0$
(figure 6).
Linear and nonlinear noise contributions and total noise in terms of SPL of acoustic pressure at microphones placed at
$r / D = 1.0$
along different axial positions
$x / D=$
(a)
$0.0$
; (b)
$0.5$
; (c)
$1.0$
; (d)
$2.0$
. The shaded grey area is the range of frequencies in which the compressibility delay may play a role in the evaluation of the nonlinear terms.
As expected, the linear part dominates the acoustic signature of the wind turbine, especially in the intermediate- to high-frequency range
$f \geqslant 3 f_{T}$
. However, the lowest frequency part of the spectra is dominated by the nonlinear source terms. A similar outcome was found by several authors, among the others, Cianferra & Armenio (Reference Cianferra and Armenio2021), Cianferra et al. (Reference Cianferra, Petronio and Armenio2019b
), Posa et al. (Reference Posa, Broglia and Felli2022b
), Turner & Kim (Reference Turner and Kim2022), Zamponi et al. (Reference Zamponi, Avallone, Ragni, Schram and Van Der Zwaag2024), showing the predominance of the nonlinear terms in the near-to-intermediate field. The linear contribution (figure 18
a) exhibits a fundamental peak at the blade passing frequency and is broadband in the high frequency range; the former aspect is associated with rigid rotation, while the latter is associated with pressure fluctuations over the blades (see figure 7). Moving downstream (figure 18
b–d), the peak decreases in amplitude practically to disappear at the farthest microphone, resulting in a broadband profile. Regarding the nonlinear contribution to the noise, figure 18(a,b) show several peaks at the blade passing frequency and its own multiples, in the low range of frequency; they are related to the large and organised inertial vortical structures, such as the tip vortices, trailing edge vortex sheet and root vortices. In particular, the near-wake coherent structures of the flow, associated with high vorticity and high levels of the Lighthill noise source term (figures 9
b and 17) are synchronised with the turbine rotation and thus produce tonal noise. Moving downstream, the nonlinear tonal peaks tend to be smoothed out and spread, leading to a low-frequency broadband profile. The latter is associated with the vortex instability mechanism and loss of synchronisation, occurring in the intermediate wake (see figure 9
b–c). In other words, the loss of synchronisation, at the expense of the growth of disorganised turbulence structures, results in a transition from a tonal spectrum to a purely broadband one, with an increase in amplitude at low frequencies. Similar behaviour was highlighted by Posa et al. (Reference Posa, Broglia and Felli2022b
) for the case of a marine propeller. The authors suggested that the main mechanisms that come into play are tip-vortex instability and successive grouping in the intermediate wake (see Felli et al. Reference Felli, Camussi and Di Felice2011).
The problem of the time delays in the range of frequencies
$3.85 \lt f/f_T\lt 30$
is related to the evaluation of the nonlinear terms; however, the results of figure 18 clearly show that neglecting the time delays does not introduce significant errors in the evaluation of the overall noise. In fact, for
$ f/f_T\gt 3.85$
, the nonlinear contribution to the noise is in general negligible compared with the linear one, and the noise spectrum practically coincides with that of the linear part. An exception is present in figure 18(c,d) in which an increased level of the nonlinear noise is observable in the range of high frequencies (
$ f/f_T\gt 20$
). Although in this range of frequencies the time delays may play a role, a similar behaviour was found by Wolf, Azevedo & Lele (Reference Wolf, Azevedo and Lele2012) and Turner & Kim (Reference Turner and Kim2022), when analysing the quadrupole part of the noise in the wake of an aerofoil and it was attributed to the Kelvin–Helmholtz instability developing in the separated shear layer.
The use of the FWH equation allows us to discern among the single source contributions to the overall noise spectrum. Figure 19 shows the contribution of each term of the FWH equation at the four microphones, the same as in the previous figure. The thickness term exhibits the main tonal peaks at the blade passing frequency and a few super-harmonics (figure 19
a). The main peak at
$f/f_T=3$
remains well visible moving downstream in the wake, although the signal decays rapidly over the entire frequency range. The loading terms (figure 19
b,c) exhibit a tonal behaviour at the blade passing frequency; super-harmonics are also present, as well as a broadband contribution in the high-frequency range. The amplitude of the acoustic pressure fluctuations exhibits non-monotonic decay; specifically, it increases moving from
$x / D = 0.0$
to
$x / D = 0.5$
and decreases afterward; the directivity of the loading term explains this characteristic and will be discussed in § 4.3. The three nonlinear terms (figure 19
d,e,f) exhibit a similar behaviour, with the third term (namely the seventh one of the FWH equation, panel f) being predominant over the others. This is due to the position of the microphones. Indeed, the FWH terms have different directivity, magnitude and decay rates. In the very near field (for small values of
$r$
, see (2.3)), the third nonlinear term may play a major role. As the distance
$r$
increases, the contribution of the second and first terms becomes more significant. This will be further discussed in § 4.3. All nonlinear terms (figure 19
d,e,f) exhibit a tonal behaviour for the microphones placed in the region of the near wake, while a broadband character is detected for the microphones placed in the region of the intermediate wake, as already discussed previously.
SPL of the six terms composing the FWH equation at microphones placed at the radial coordinate
$r / D = 1.0$
along different axial positions:
$x / D = 0.0, 0.5, 1.0, 2.0$
. (a) First and second linear term; (b) third linear term; (c) fourth linear term; (d) first nonlinear term; (e) second nonlinear term; (f) third nonlinear term. The shaded area of grey is the range of frequencies in which the compressibility delay may play a role in the evaluation of the nonlinear terms.
The direct formulation of the FWH equation allows us to assess the contributions of different regions of the domain to the overall acoustic spectrum. With this purpose, we perform two sets of computations, considering different shapes of subdomains to detect which fluid-dynamic structures contribute the most to the noise: the first set considers annular subdomains comprehensive of the surface of the rotor; the second set considers different cylindrical regions of the wake, roughly corresponding to the near wake and two regions of the intermediate wake.
For the first set of computations, the linear and nonlinear parts are evaluated in the inner region (
$0.0 \lt r / D \lt 0.2$
), in the middle region (
$0.2 \lt r / D \lt 0.45$
) and in the outer region (
$0.45 \lt r / D \lt 0.9$
). The three integration volumes are concentrical cylinders that extend from
$x / D = -0.5$
to
$x / D = 7.0$
. In figure 20, we show the linear contribution to the noise, considering the inner, middle and outer regions of the turbine surface. At the nearest microphone (figure 20
a), peaks at blade passage frequency
$f/f_{T}=3$
and its first super-harmonic
$f/f_{T}=6$
are visible, mainly generated by the thickness term (see figure 19
a; here, the outer and middle sections contribute equally). In the high-frequency range, pressure fluctuations in the middle region of the blades (see figure 16), feeding the loading term, are found to contribute the most.
Linear contribution of the inner (
$0\lt r / D\lt 0.2$
), middle (
$0.2\lt r / D\lt 0.45$
) and outer (
$0.45\lt r / D\lt 0.9$
) surfaces in terms of SPL at microphones placed at
$r / D = 1.0$
along different axial positions
$x / D=$
(a)
$0.0$
; (b)
$0.5$
; (c)
$1.0$
; (d)
$2.0$
.
Nonlinear contribution of the inner (
$0\lt r / D\lt 0.2$
), middle (
$0.2\lt r / D\lt 0.45$
) and outer (
$0.45\lt r / D\lt 0.9$
) volume in terms SPL at microphones placed at
$r / D = 1.0$
across different axial positions
$x / D=$
: (a)
$0.0$
; (b)
$0.5$
; (c)
$1.0$
; (d)
$2.0$
.
At downstream locations (figure 20 b–d), we observe that the middle section of the blades still dominates the linear noise, being mainly generated by the loading term (figure 19 b,c). This is consistent with the analysis of pressure fluctuations over the blades (figure 16). The analysis suggests that trailing edge pressure fluctuations are the main source of noise associated with the wind turbine surface.
In figure 21, we report the nonlinear term, considering the inner, middle and outer regions of the turbine wake. In the region of the near wake (figure 21
a,b), the outer and the middle regions contribute equally to the nonlinear noise at the blade passing frequency peak and at its super-harmonics. This originates from intense sources, namely the tip vortex system and the trailing edge vortex sheets (see figure 17
a–c). Moving downstream, the outer volume governs the noise generated in the low-frequency range (figure 21), and we observe the progressive disappearance of the peaks at blade passing frequency and its super-harmonics. The analysis of the fluid-dynamic field (figure 14) has shown that in the intermediate wake, TKE is primarily produced in the outer section due to tip vortex instability and grouping. This is the primary reason for the broadband character of the nonlinear noise. The inner volume contribution to the overall noise is negligible at all locations, highlighting that the root vortices dynamics is, by far, less important than the tip and blade vortices. For the sake of completeness, we computed the acoustic signal at additional microphones positioned at
$x / D = 3, 4, 5, 6, 7$
at the same radial distance
$r / D = 1$
(not shown). The spectral component of the linear term remains similar to that depicted in figure 20(d), with amplitude decaying with the distance from the turbine. The behaviour of the nonlinear term with increasing distance from the turbine, still in the region of the intermediate-to-far wake, remains consistent with that shown in figure 21(d) and the sound level remains relatively high. Overall, the turbulent fluctuations occurring in the outer region are found to be the primary source for all the microphones analysed in this section. A detailed analysis of the decay of linear and nonlinear terms is in § 4.3.
For the second set of computations, we evaluate the nonlinear terms of the noise associated with different regions of the wake, namely the near wake (
$-0.5 \lt x / D \lt 1.0$
), the first part of the intermediate wake (FPIW) (
$1.0 \lt x / D \lt 5.0$
) and the downstream part of the intermediate wake (DPIW) (
$5.0 \lt x / D \lt 7$
). The microphones are placed at
$r / D = 1.0$
along different axial positions
$x / D=0.5$
,
$3.0$
,
$5.0$
,
$7.0$
, thus spanning over the three regions herein considered.
As the microphone moves along the streamwise direction, the dominant acoustic contribution clearly comes from the region directly beneath it; this trend is evident across all panels. Specifically, close to the rotor (figure 22 a), the near wake is dominant, exhibiting the peak at the blade passing frequency and super-harmonics.
As the distance from the source increases, the intermediate wake contributes the most (figure 22 b–d). The relative importance of the two sub-regions of the intermediate wake is inverted moving further downstream, with the DPIW dominant in the farther microphone. This highlights the necessity to consider the whole wake in the computation of the noise in the near field.
In figure 23, we report the SPL computed at microphones placed in the turbine plane
$x / D = 0.0$
, at three radial distances,
$r / D = 2.0, 4.0, 8.0$
. Note that the spectra at
$r / D =1.0$
are already shown in figure 18(a).
For the linear part, the behaviour with increasing distance is similar to what was already observed for the microphones positioned along the x-axis. Specifically, as the radial distance increases, the tonal peaks gradually vanish, along with the broadband part of the spectrum in the high-frequency range.
For the nonlinear part, in contrast to what was observed for the microphones placed along the wake, the broadband noise at low frequencies associated with intermediate-to-far wake turbulence is not detected. Instead, the signal decays uniformly across the whole frequency range.
Nonlinear contributions to the noise associated with different regions of the wake, at four microphones placed along the wake: near wake (
$-0.5\lt x / D\lt 1.0$
), blue lines; first part of the intermediate wake (FPIW) (
$1.0\lt x / D\lt 5.0$
), yellow lines; downstream part of the intermediate wake (DPIW) (
$5.0\lt x / D\lt 7.0$
), violet lines. The total nonlinear term is shown with a red line. The microphones are placed at
$r / D = 1.0$
along different axial positions
$x / D=$
(a)
$0.5$
; (b)
$3.0$
; (c)
$5.0$
; (d)
$7.0$
.
Linear and nonlinear noise contributions in terms of SPL of acoustic pressure at microphones placed at
$x / D = 0.0$
along different radial positions
$r / D=$
(a) 2.0; (b) 4.0; (c) 8.0.
4.3. Directivity and decay rate
In the previous section, we discussed spectra at specific points along the streamwise direction out of the wake and along the radial direction. The analysis along two directions is not exhaustive since the acoustic signal typically exhibits directivity depending on the source characteristics (monopole, dipole, quadrupole). In addition, the decay rate of each FWH term (2.3) plays a role in the propagation of the signal in the near-to-far field.
Here, we analyse the directivity, which gives the preferential direction of noise propagation in the meridian plane
$z=0$
. We consider
$20$
microphones placed in concentric circles around the source, at different radii. To choose the radii, we perform a rough estimation of the acoustic near-to-far field transition, in accordance with the literature. The choice depends on a reference wavelength, for example,
$\lambda = c_0 /3f_T$
, which is the wavelength associated with the blade passage frequency. The condition
$K_{\lambda } = d/\lambda \lt 1$
, where
$d$
is the distance from the source, may be used to delimit the acoustic near field from the far field. Based on this discussion, we consider the radii
$2D$
,
$10D$
, and
$100D$
, which respectively give
$K_{\lambda } \sim 0.22, 1.1, 10.8$
. We consider separately the linear and nonlinear contributions, and the overall noise, in figure 24.
Directivity of the root mean square of the acoustic pressure fluctuations
$C_{\hat p_{\textit{rms}}}$
at microphones placed in concentric circles on the
$x{-}y$
plane, at: (a–c)
$2D$
; (d
–f)
$10D$
; (g–i)
$100D$
. (a,d,g) Linear contribution; (b,e,h) nonlinear contribution; (c,f,i) overall contribution. The arrow indicates the streamwise direction. Note that at the radial distance of 2D the microphones between
$30^{\circ}$
–
$330^{\circ}$
are not shown since they lie inside the FWH integration volume.
The linear part exhibits the typical dipole shape with lobes aligned with the flow direction and with decreasing amplitude with increasing distance from the source (figure 24 a,d,g). The dipole-type pattern aligns with expectations and is consistent with findings from other rotor studies (see, for example, Wang et al. (Reference Wang, Wang and Wang2021), in spite of differences in test case and operating conditions). The directivity of the linear term, which exhibits a dipole-like shape, explains the non-monotonic decay of the acoustic intensity moving downstream, as already observed in figure 19(b,c) about the spectral analysis.
The nonlinear part exhibits a quadrupole shape (figure 24 b,e,h). In the near field, the observed asymmetry is due to the wake, the lobe in the streamwise direction being much more pronounced than the others. With increasing distance from the source, the magnitude of the signal decays rapidly and a more symmetric quadrupole shape tends to recover.
The overall noise level in the near-field (
$K_{\lambda } = 0.21$
), figure 24(c), exhibits a non-symmetric dipole shape, given by the summation of the two contributions. This highlights the importance of the nonlinear terms in the direction of the wake, near the source. To be noted is that the directivity is expected to be symmetric with respect to the
$0^\circ {-}180^\circ$
axis. The slight asymmetry has to be ascribed to imperfect statistical convergence.
The directivity in the rotor plane is not reported because the axial symmetry of the wind turbine motion results in a homogeneous noise propagation over the circumferential directions. In this case, the decay rate along the radial direction is more suited to evaluate the relative importance of the linear and nonlinear noise source contributions. For this purpose, we calculate the root mean square of the acoustic pressure coefficient (
$C_{\hat p_{\textit{rms}}} = \hat p_{\textit{rms}} /0.5 \rho U_{\infty } ^2$
) considering three revolution periods, up to a distance of
$1000D (K_{\lambda } \sim 108)$
. In figure 25, for the sake of clarity, we show the directivity up to
$100D\, (K_{\lambda } \sim 10.8)$
.
Decay of the root mean square of the acoustic pressure coefficient: (a) radial decay; (b) decay along the streamwise direction out of the wake.
In figure 25(a), we report the decay rate of the linear and nonlinear terms together with that of the overall noise along the radial direction in the turbine plane.
The linear and nonlinear terms decay as
$r^{-1}$
in the far field, with the linear term dominating over the other; this finding is consistent with the literature, where the linear term is recognised as the principal contributor to the overall noise in the far field. In the near field, linear and nonlinear contributions decay much faster than
$r^{-1}$
; specifically, the nonlinear terms decay roughly as
$r^{-4}$
; the term-by-term analysis (not reported) shows that the third nonlinear term contributes the most, in accordance with results shown in figure 19 (compare panels d,e,f).
The distance at which the decay rate changes to
$r^{-1}$
, in this case, around
$z=4D$
, may depend on various features, such as the rotational Mach number and the number of blades. In the near-field, the linear term decays as
$r^{-4}$
.
This pronounced decay would seem anomalous, given (2.3), which indicates
$r^{-2}$
in the near-field and
$r^{-1}$
in the far-field. It can be explained by the destructive interference associated with the rotational motion of the noise sources, as discussed in detail by Ianniello (Reference Ianniello2016). In that paper, the author found that the kinematic auto-interaction of the noise sources produces destructive interference, leading to a decay rate higher than the expected
$r^{-1}$
law. In particular, the author found that considering
$N$
rotating sources, the decay rate of the linear term approaches the
$r^{-(N+1)}$
law in the rotor plane in the near field; our results are in accordance with that of the rotpole model. This applies to the linear terms, due to geometrical reasons. Interestingly, our results are consistent with those of Posa et al. (Reference Posa, Broglia and Felli2022a
) for a system composed of a seven-blade marine propeller placed upstream of a hydrofoil and Posa et al. (Reference Posa, Broglia, Shi and Felli2025) for conventional and tubercled five-blade marine propellers. Indeed, along the radial direction, the authors observed decay trends and relative contribution of the different FWH terms with increasing distance similar to those herein reported. However, such comparisons should be approached with caution due to the influence of several key parameters, such as the rotational Mach number, the rotation frequency and the number of blades.
Figure 25(b) reports the decay rates along a line parallel to the wind turbine axis, at a distance of
$r / D=1$
, out of the wake. The linear term decays as
$r^{-1}$
as the distance increases. Up to
$9D$
, the nonlinear term dominates; this is expected since in the region up to
$8D$
, the selected microphones are placed close to the turbulent wake, leading to a non-monotonic behaviour peaking at
$4.5 D$
. Beyond that, the nonlinear part drops very rapidly, decaying asymptotically as
$r^{-1}$
.
To summarise, beyond
$10D$
, the decay rate of the overall noise fits the
$r^{-1}$
law, due to the linear part, indicating that the latter becomes dominant in the far field.
When comparing the decay rates (panel a versus panel b of figure 25), mutual interaction between rotating noise sources is evident in the turbine disk since it gives rise to a destructive effect (higher decay rate), while the mutual interaction is not clearly identified out of the rotor plane, along the streamwise direction. The destructive effect is described using the rotpole model by Ianniello (Reference Ianniello2016), where the author showed decay rates steeper than
$r^{-1}$
along axial and radial directions, by analysing a multi-rotpole configuration. In contrast, in the present study, we found that the decay along the wake axis adheres consistently to the
$r^{-1}$
law. One possible explanation for this discrepancy, as previously asserted by Ianniello (Reference Ianniello2016), is that the rotpole model relies on a constant pressure point-like source, while in the present case, realistic unsteady pressure loads on the blade surface are present. The broadband pressure distribution over the wind turbine surface tends to break up the synchronisation of the signals, being more similar to a set of time-dependent distributed sources along the surface. This is particularly true when investigating turbulent flow-induced pressure disturbances, such as the turbulent boundary layer at the trailing edge of a blade. Conversely, the concept of multiple rotpoles is fully representative of the phenomenon when the observer is placed on the turbine disk and produces accurately the mechanisms of kinematic auto-interaction. This can be attributed to the nearly constant pressure distribution in the tip region resulting from the reduced instability induced by the tip vortex.
To show the importance of unsteady pressure fluctuations on the decay rate of the linear terms, we extend the rotpole model considering a stochastic pressure signal (stochastic rotpole). In particular, we consider a tri-rotpole model with the rotating point-sources characterised by pressure
$\tilde {p}$
being a stochastic quantity; we use the same geometrical and rotational properties of the model-scale wind turbine herein considered (radius
$R$
, rotational velocity
$\varOmega$
, speed of sound
$c_0$
) to have the same rotational Mach number. We consider the classical formulation of the FWH equation for the loading component in terms of multiple rotating noise sources:
\begin{equation} 4\pi\! \hat p(\boldsymbol{x},t) = \frac {1}{c_0} \frac {\partial }{\partial {t}} \sum _{l=1} ^{N} \frac {\tilde {p}_l \boldsymbol{n}_{l} \boldsymbol{\cdot }\hat {\boldsymbol{r}}_{l}}{r_{l} |1-M_{r_l} |} + \sum _{l=1} ^N \frac {\tilde p_l {\boldsymbol{n}}_{l} \boldsymbol{\cdot }\hat {\boldsymbol{r}}_{l}}{{r_{l}}^2 |1 - M_{r_l} |} ,\end{equation}
where
$\tilde {p}_l = p_s + p_{rand}$
is the stochastic (normal distribution) pressure fluctuation, in which
$p_s = 1 \,{\mathrm{Pa}}$
and
$p_{rand}$
is a random distribution with zero mean value and assigned standard deviation,
$\boldsymbol{n}_l$
is the normal vector associated with the
$l$
th rotating source,
$\hat {\boldsymbol{r}_{l}}$
is the versor connecting the source,
$r_l$
is the source-observer distance,
$M_{r_l} = \boldsymbol{v}_l\boldsymbol{\cdot }\boldsymbol{r}_l/c_0$
is the rotational Mach number, where
$\boldsymbol{v}_l$
is the velocity vector of the rotating source (for more details, see Ianniello Reference Ianniello2016) and
$N$
is the number of rotating sources, which is equal to 3 in this case. For clarity, we consider the normal versors
$ \boldsymbol{n}_1$
,
$\boldsymbol{n}_2$
and
$\boldsymbol{n}_3$
, the first pointing in the radial direction, the second is tangent to the circumferential trajectory of the rotating sources and the last is normal to the rotor plane (for more details, see figure 12 of Ianniello Reference Ianniello2016). Three case are here considered, namely
$p_{\textit{rms}}/\!p_s=0$
,
$p_{\textit{rms}}/\!p_s=0.01$
,
$p_{\textit{rms}}/\!p_s=0.1$
. The first one coincides with the original rotopole model.
Acoustic pressure coefficient decay rate in the radial direction, considering the tri-rotpole model:
$\tilde {p}_{\textit{rms}} /p_s = 0.0$
yellow solid line;
$\tilde {p}_{\textit{rms}} /p_s = 0.01$
red-dot line;
$\tilde {p}_{\textit{rms}} /p_s= 0.1$
blue-diamond marker line. Contribution associated to: (a) the normal
$\boldsymbol{n}_1$
l; (b) the normal
$\boldsymbol{n}_2$
. Note that the contribution associated to the normal
$\boldsymbol{n}_3$
is zero.
Acoustic pressure coefficient decay rate in the streamwise direction, considering the tri-rotpole model:
$\tilde {p}_{\textit{rms}} /p_s = 0.0$
, yellow solid line;
$\tilde {p}_{\textit{rms}} /p_s = 0.01$
, red-dot line;
$\tilde {p}_{\textit{rms}} /p_s= 0.1$
, blue-diamond marker line. Contribution associated to: (a) the normal
$\boldsymbol{n}_1$
; (b) the normal
$\boldsymbol{n}_2$
; (c) the normal
$\boldsymbol{n}_3$
.
In figures 26 and 27, we show the acoustic pressure coefficient decay rate along the radial and streamwise direction, respectively. The microphones are placed along the same lines as shown in the previous figure 6. The original tri-rotpole model, in the near field, exhibits the
$r^{-4}$
and
$r^{-7}$
decay rate in the radial and streamwise direction, respectively; moving towards the far field, an asymptotic decay is evident (compare our figures 26 and 27 with figures
$27\;\mathrm{and}\;28$
of Ianniello Reference Ianniello2016). However, we observe that pressure fluctuations induce a faster transition from the near-field high-order decay rate to the asymptotic ones. Specifically, the increase in the level of fluctuations anticipates this transition. As also argued by Ianniello (Reference Ianniello2016), this indicates that the presence of pressure fluctuations, always present in a turbulent flow field, reduces or even eliminates the destructive kinematic auto-interaction present in the case of rotating constant-pressure signals. This stochastic rotpole model highlights the role of fluid dynamic pressure fluctuations in noise propagation, providing additional conceptual support to our results.
To conclude, it is evident that in the rotor plane, where the tip of the blade contributes the most to loading noise generation, pressure fluctuations given by the released tip vortex are weak (as depicted in figure 16 c), thereby promoting a higher decay rate due to the phenomena of kinematic auto-interaction. However, along the wake of the rotor, the downwind face of the blade where the fluid flow exhibits separation and turbulent transition induces high turbulent pressure fluctuations, so that the phenomenon of kinematic auto-interaction is practically absent and the decay is linear.
5. Conclusions
In the present paper, we performed the acoustic characterisation of a model-scale wind turbine, using LES and the acoustic analogy. Considering perfect Mach similarity (Cianferra & Armenio Reference Cianferra and Armenio2021), the present analysis is representative of medium-size wind turbines with diameter of the order of
$100 \, {\mathrm{m}}$
rotating with a period of approximately
$8 \, {\mathrm{s}}$
and, consequently, characterised by low rotational Mach number (
$M_{\textit{tip}} \sim 0.10$
).
The model-scale wind turbine herein considered is that of the experimental campaign of Gambuzza & Ganapathisubramani (Reference Gambuzza and Ganapathisubramani2021), except for the absence of the mast. The bulk flow quantities were in good agreement with the experimental data, also in light of differences in geometrical and flow configurations.
The fluid dynamic analysis showed the following features: a turbulent boundary layer develops over the blades together with a composite system of tip vortex, blade wake and root vortices in the near wake; the intermediate wake is characterised by destabilisation of the tip vortex and the mutual inductance mode triggering the leap-frogging mechanism; the fully turbulent wake populated by a number of turbulent structures is present in the far field. Our results agree with those of the relevant literature. In particular, these fluid dynamic features are also shared with those observed in the case of an isolated ship propeller, in spite of the different working conditions (turbine versus thruster) and substantially different geometry. However, the wake of a marine propeller is characterised by the presence of a strong rectilinear hub vortex, in contrast with the case of a wind turbine in which the hub vortex seems to be almost insignificant.
The noise generation mechanisms are discussed in light of the fluid dynamic field. The analysis showed the presence of two main noise source mechanisms, one associated with fluctuating pressure loads on the wind turbine surface and the other associated with large, coherent and inertial structures of motion populating the wake. The trailing edge noise source is associated with the fluctuating pressure field in the turbulent boundary layer. It exhibits broadband characteristics and rules the high-frequency range. The pressure loads, together with wind turbine rotation, are also responsible for the tonal peaks related to the blade passing frequency and its super-harmonics. These peaks are more evident near the source, while, as the distance from the source increases, they tend to vanish. The contribution of the wake to the noise is significant in the low-frequency range
$f \lt 4f_T$
. Among the vortical structures present in the wake, the role of vortex grouping has been highlighted: it causes the loss of tonal peaks, at the blade passing frequency, and generates low-frequency broadband noise.
The contribution to the noise originating from different regions has been investigated, considering two different configurations, the first composed of three different annular concentric regions, and the other composed of different sub-regions of the wake. The analysis of the first configuration confirmed that the trailing edge vortex system is the most important noise source, as regards the linear part of the noise; while the tip vortex system, together with its instability and evolution in the intermediate wake, is the dominant mechanism of the nonlinear noise source. This analysis puts in evidence similarities and differences between wind turbines and marine propellers; both devices share the importance of the tip vortex system and of its own evolution as a nonlinear noise generator; however, in the case of a ship propeller, the acoustic response is influenced by a close interaction between the systems of tip and hub vortex, which is not the case for the wind turbine. The analysis of the second configuration highlights the importance of the entire wake for the evaluation of the nonlinear contribution to the noise in acoustic near field.
To evaluate the preferential direction of noise propagation, we have carried out the analysis of the directivity on the plane orthogonal to the rotor passing through its centre. The linear contribution exhibits a nearly symmetrical dipole shape in spite of the presence of advection; the nonlinear contribution exhibits a non-symmetrical and elongated quadrupole-like shape in the near field, recovering symmetry moving away from the rotor. The directivity of the overall noise in the near field exhibits a non-symmetric dipole shape, due to the combination of the linear and nonlinear contributions. The larger lobe is elongated downstream, highlighting that the nonlinear part contributes the most in the wake direction. Moving in the acoustic far-field, a nearly symmetrical dipole shape is recovered, which is due to the linear part of the noise, the nonlinear one becoming marginal. From a practical point of view, of particular interest for industrial applications, we may argue that the directivity of an isolated wind turbine field is substantially represented by a three-dimensional dipole oriented along the wind turbine axis, with symmetry between the lobes achieved moving towards the far field.
Finally, we analysed the decay rate of the noise level along the radial and axial directions. Along the radial direction, both the linear and nonlinear terms decay as
$r^{-4}$
in the near field and as
$r^{-1}$
in the far field. However, in the direction perpendicular to the rotor plane, the linear contribution decays as
$r^{-1}$
, whereas the nonlinear part exhibits a non-monotonic behaviour caused by the presence of the elongated wake. These results are in close alignment with those reported by Posa et al. (Reference Posa, Broglia and Felli2022a
) and Posa et al. (Reference Posa, Broglia, Shi and Felli2025) for the case of a marine propeller.
The near-field decay of the linear term is well explained by the rotpole model developed by Ianniello (Reference Ianniello2016) and its stochastic extension reported in the present paper. In the absence of energetic fluctuations, the anomalous decay
$r^{-4}$
along the radial direction is associated with the destructive kinematic auto-interaction of the three rotating sources. The phenomenon is evident in the rotor plane, caused by the low level of pressure fluctuations at the tip of the blades, associated with a stationary vortex structure. However, in the direction of the axis of the turbine, the high level of pressure fluctuations present over the blades surface inhibits the kinematic auto-interaction and the
$r^{-1}$
decay is recovered.
Overall, the analysis of the sound level decay rate and the directivity suggest that for the rotational Mach number herein considered, the volume terms may play a role near the source (
$K_{\lambda} \lt 1$
).
In conclusion, these outcomes reveal important characteristics of the noise radiated by wind turbines. These characteristics can be used for the acoustic characterisation and mitigation of wind turbines, and may be useful for propagation modelling in industrial applications, although a more comprehensive study should consider the presence of the mast and the presence of the terrain. This is the topic of successive studies.
Acknowledgements
We acknowledge Valorem SAS (213 cours Victor Hugo 33 323 BEGLES CEDEX – FRANCE) for the support and the exchange of information and knowledge in the wind turbine field (https://www.valorem-energie.com/en/).
Funding
The authors acknowledge the support by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA). This research is part of the project BluEcho – From science to policy: assessing impacts and developing solutions for ship traffic and offshore wind farms through detailed soundmaps – CUP J93C23002010007 – funded by the Sustainable Blue Economy Partnership. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support (project name: IsB29-WTBAR, OriginID: HP10BR95C2).
Declaration of interests
The authors report no conflict of interest.
Open access
