# Wind turbine noise synthesis from numerical simulations

Andrea P. C. Bresciani* and Julien Maillard†

*Centre Scientifique et Technique du Bâtiment (CSTB), Saint-Martin-d'Hères, 38400, France.*

Sophie Le Bras‡

Siemens Industry Software SAS, 107 Avenue de la République, 92320 Châtillon, France

Leandro D. de Santana

This paper proposes a methodology for physics-based auralization of wind turbine noise with potential application to complex urban and rural scenes. Time-dependent aerodynamic broadband noise levels are computed with cost-efficient semi-analytical models that account for the presence of ground and atmosphere. The algorithm for an accurate reconstruction of the time domain sound signal is provided, with particular attention to the calculation of the Doppler effect in the time domain. Thanks to its low computational cost, this methodology can be applied simultaneously to multiple wind turbines in a complex environment. The scenario can potentially include several ground types, buildings, and other noise sources such as road or railway traffic. The numerical workflow is first validated in free-field with Schlinker and Amiet's analytical theory. The accuracy of the predicted noise levels is assessed by comparing them with measurements available in the literature for a 2.3 MW wind turbine. The spectrograms and the audio signals show the "swishing" amplitude modulation for receivers close to the cross-wind direction. Besides wind energy application, the workflow presented in this work can be applied to auralize the broadband noise generated by any rotating machine, such as a helicopter and fans.

## I. Introduction

Nowadays acoustic annoyance limits the deployment of wind turbines near dwellings. Wind turbine noise has been recognized as the most annoying sound source [1] despite the fact of having lower noise levels compared to other environmental sound sources such as wind, road, and rail. Furthermore, evidence relates wind turbine noise with sleep disturbances and, indirectly, with stress [2]. For these reasons, reliable wind turbine noise predictions are crucial to evaluate their impact on the surrounding environment, ideally prior to the installation of the wind farm.

In this context, the auralization of wind turbine noise is a valuable tool for the evaluation of dynamic noise indices, for noise acceptance studies, as a demonstrator for non-experts, or as a tool for the optimal placement of wind turbines. Different auralization techniques have been already applied to wind turbine noise. They fall into two categories: sampled-based auralization [3] and physics-based auralization [4, 5]. Whereas sample-based auralization relies on audio recordings for the construction of the sound pressure signals, physics-based auralization uses a numerical model. The advantage of physics-based auralization is that several wind farm layouts, operating and atmospheric conditions can be studied by changing the parameters of the simulation. However, a model for the prediction of wind turbine noise sources is necessary. Two aspects have to be modeled and coupled for calculating the noise levels, i.e., the aerodynamic noise emitted by the turbine and the propagation of the sound to the listener placed in the far field. Various methodologies presented by the literature focus on the BPM model [6] or Amiet's theory [7, 8] for the noise source modeling, and on ray-tracing [9–11] or parabolic equation [12, 13] methods to simulate the atmospheric propagation effects. For instance, Cotté [12] coupled Amiet's theory for the prediction of the sources with the parabolic equation method to study long-range wind turbine noise propagation.

Mascarenhas et al. [5] recently proposed a physics-based methodology to auralize wind turbine noise in free-field. Their approach relies on cross-fading window functions to smooth the transition between two consecutive time signals (*grains*) corresponding to different wind turbine blade positions. One limitation of this approach is the overlap between two grains which has to be tuned to avoid audible artifacts switching from one grain to the next one. Pieren et al. [3]

<footnote>*Ph.D. candidate, University of Twente, a.p.c.bresciani@utwente.nl.</footnote>

<footnote>†Research engineer, Centre Scientifique et Technique du Bâtiment, julien.maillard@cstb.fr.</footnote>

<footnote>‡Research engineer, Siemens, sophie.le_bras@siemens.com.</footnote>
proposed a methodology based on recorded samples where the synthesized sound is composed of a tonal component and an amplitude-modulated broadband part. The broadband part is obtained with the spectral shaping technique, which is also applied in the present work. The results of the listening test provided by Pieren et al. [3] showed that the synthesized sound can be confused with the real audio recording, thus indicating an acceptable level of realism of the synthetic noise.

The present work focuses on the physics-based auralization of wind turbines. The methodology presented here is based on a continuous-time signal, which is amplitude modulated and delayed based on the computed time-dependent noise levels. One of the advantages of this approach is that it does not require the definition of any tuning parameter, ensuring an accurate reconstruction of the noise signal. Furthermore, this approach allows an efficient decoupling between the calculation of the sound sources and the far-field propagation. A RANS-based Amiet's theory for leading- and trailing-edge noise [7, 8] is used to compute the aerodynamic sound sources of the wind turbine. It is coupled with engineering ray-based methods for outdoor sound propagation. In particular, the Harmonoise model [14, 15] is used in this work, but the same methodology can be applied to any other engineering model, such as ISO9613-2, NMPB-2008, CONCAWE, or Nord2000. These mid-fidelity models are able to provide results with a computational cost much lower than the parabolic equation or standard ray tracing for a refracting medium. In this way, this methodology can be applied simultaneously to multiple wind turbines in a complex environment, that can potentially include several ground types, buildings, and other noise sources such as road or railway traffic. The computed time-varying noise levels are used for the auralization of wind turbine noise which will be included in an existing framework for noise annoyance evaluation of complex outdoor sound scenes [16–18].

The paper is organized as follows. The modeling of the aerodynamic noise sources and the far field propagation method are described in Section II. The auralization method is then explained in Section III. Finally, the complete numerical workflow is applied to a SWT 2.3-93 wind turbine and the results are presented in Section IV.

## II. Noise source modelling and far-field propagation method

This section presents the methodology for the calculation of the broadband aerodynamic noise of a wind turbine. The methodology to obtain the noise sources from the 3D CAD of the wind turbine blade is outlined in Section II.A. The Harmonoise model used for the far-field noise propagation is briefly recalled in Section II.B. The difference with respect to the Harmonoise model described in previous literature is explained in detail.

### A. Noise source model

The wind turbine sound sources are computed based on Amiet's theory for leading- and trailing-edge noise. The approach presented here is based on blade segmentation. In practice, the blade geometry is divided into segments in the spanwise direction as shown in Fig. 1. The aspect ratio of each segment is defined as $L_{span}/c_{mean}$, where $L_{span}$ is the spanwise extension of the segment and $c_{mean}$ is the chord length of the airfoil in the mid-span of the segment. The mid-span airfoil of each segment is used in a 2D Reynolds-averaged Navier–Stokes (RANS) simulation to compute the boundary layer parameters required by Amiet's theory for trailing-edge noise. The angle of attack $\alpha$ of the RANS simulation is computed accounting for the geometrical twist angle of the blade $\beta$, the collective pitch angle $\theta$, and the velocity triangle formed by the wind speed $w$, the induced velocity $v_i$, and the rotational speed $\Omega r$, i.e.,

$$ \alpha = \arctan\left(\frac{w - v_i}{\Omega r}\right) - \beta - \theta \quad (1) $$

where $\Omega$ and $r$ denote the angular velocity of the turbine and the radial position of the airfoil, respectively. The induced velocity $v_i$ is computed with the momentum theory [19] as

$$ v_i = w \left( 1 - \sqrt{1 - C_t} \right) / 2 \quad (2) $$

where $C_t = 2T_h/(\rho Aw^2)$ is the thrust coefficient, $T_h$ is the thrust, $A$ denotes the swept area and $\rho$ represents the density of the air. The wind speed is assumed constant over the rotor disk resulting in axis-symmetric conditions. This hypothesis allows us to drastically reduce the number of RANS and Amiet's theory simulations necessary. However, this approach can be applied also for non-uniform inflow velocity [20].

The 2D RANS simulations are performed with CFD software Simcenter STAR-CCM+ [21]. The simulations are wall-resolved, i.e. with $y^+ < 1$. The turbulence model used is k-ω SST, while the $\gamma$-transition model [22] is used to model the laminar-to-turbulent flow transition. All the RANS simulation setup, including the meshing process and the computation of the boundary layer parameters, is fully automated in STAR-CCM+.
<figure><img src="image_1.png" /><figcaption>Fig. 1 Application of the blade segmentation approach to a reconstructed SWT 2.3-93 wind turbine blade. The airfoils used for the 2D RANS simulations are indicated in blue.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2 Reference frame used in Amiet's theory for trailing-edge noise.</figcaption></figure>

For each blade segment, Amiet's theory for leading- and trailing-edge noise [7, 8] is used to compute the far-field sound pressure on the surface of a sphere with arbitrary radius $r$ in the far field. Each sphere is centered in the middle of the corresponding blade segment.

In Amiet's theory, the airfoil is modeled as an infinitely thin flat plate with chord $c = 2b$, span $L = 2d$, and at a zero angle of attack with respect to the freestream velocity $U$. If the observer is located in the acoustic and geometrical far-field and the aspect ratio $L/c$ of the flat plate is large, the power spectral density for trailing-edge noise $S_{pp}^{TE}$ is computed as [23]:

$$S_{pp}^{TE}(\mathbf{x}, \omega) = \left( \frac{\omega x_3 b}{2\pi c_0 \sigma_0^2} \right)^2 d |\mathcal{I}(x_1, \omega/U_c, k_{20})|^2 l_y(\omega/U_c, k_{20}) \Phi_{pp}(\omega), \quad (3)$$

where $\mathbf{x} = (x_1, x_2, x_3)$ is the position of the observer in the reference frame depicted in Fig. 2; $\omega$ is the pulsation; $c_0$ denotes the speed of sound; $\sigma_0 = \sqrt{x_1^2 + \beta^2(x_2^2 + x_3^2)}$; $\mathcal{I} = \mathcal{I}_1 + \mathcal{I}_2$ represents the aeroacoustic transfer function for trailing-edge noise, given in [23, 24]; $U_c = 0.7U$; $k_{20} = kx_2/\sigma_0$ with $k = \omega/c_0$; $l_y$ is the spanwise correlation length obtained from Corcos model [25] with constant $b_c = 1.47$. The wall pressure spectrum $\Phi_{pp}$ is computed with empirical models for computational efficiency. Lee's wall pressure spectrum model [26] for adverse pressure gradient flows and Goody's [27] for zero or favorable gradients are used.

Amiet's formulation for leading-edge noise reads

$$S_{pp}^{LE}(\mathbf{x}, \omega) = \left( \frac{\rho \omega c x_3}{2 c_0 \sigma_0^2} \right)^2 \pi U d |\mathcal{L}(x_1, \omega/U, k_{20})|^2 \Phi_{ww}(\omega/U, k_{20}) \quad (4)$$

where $\mathcal{L} = \mathcal{L}_1 + \mathcal{L}_2$ is the aeroacoustic transfer function for leading-edge noise [24] and $\Phi_{ww}$ is the inflow turbulence spectrum modeled using the von Kármán spectrum in this work.
<figure><img src="image_3.png" /><figcaption>Fig. 3 Directivity calculated with Amiet’s theory for leading- and trailing-edge noise for the tip segment airfoil of the SWT 2.3 MW wind turbine. The 1/3-octave power spectrum is shown for the bands with center frequencies (a) 100 Hz, (b) 1000 Hz, and (c) 8000 Hz.</figcaption></figure>

The total aerodynamic broadband noise emitted by each blade segment is the sum of the trailing- and leading-edge contributions, $S_{pp} = S_{pp}^{TE} + S_{pp}^{LE}$. For each of the N blade segments, a directivity similar to the one shown in Fig. 3 is obtained. The spherical divergence $10 \log_{10}(4\pi(r/r_0)^2)$ dB is then added to the noise levels to obtain the sound pressure levels of an equivalent point source at $r_0 = 1$ m. These spheres are used as point sources in the Harmonoise model [14, 15] to compute the time-dependent noise levels in 1/3-octave bands from 50 Hz to 10 kHz. Hence, the wind turbine is modeled as an extended source instead of a single-point source. This is necessary to correctly capture the interference dips and to predict the amplitude modulation [12].

## B. The Harmonoise model

The Harmonoise model [14] is an engineering ray-based method for outdoor noise prediction. In addition to spherical divergence, it calculates separately the main effects affecting the sound propagating in the atmosphere. These effects are atmospheric absorption, turbulence scattering, ground reflection, and diffraction. Their combined attenuation is commonly referred to as the excess attenuation. For brevity, the Harmonoise model is not presented here as the details can be found in references [14, 15]. We present here only the difference with respect to the model presented by Salomons et al. [14]. The difference concerns the calculation of the effective structure function parameter $\gamma_T$. The effective structure function parameter $\gamma_T$ affects the calculation of the excess attenuation due to turbulence scattering and the coherence factor used for the calculation of the excess attenuation due to the ground reflection. In the reference paper, it is suggested to use $\gamma_T = 5 \cdot 10^{-6}$, as a representative value for moderate turbulence. In the present work, the value of $\gamma_T$ is calculated instead using the Monin–Obukhov similarity theory.

The effective structure function parameter is defined as [14, 28]

$$ \gamma_T = \left( \frac{C_T}{T_r} \right)^2 + \frac{22}{3} \left( \frac{C_W}{c_0} \right)^2 \qquad (5) $$

where $C_T$ and $C_W$ are the structure function parameters due to temperature and velocity fluctuations, respectively, and $T_r$ is the reference air temperature. The structure function parameters $C_T$ and $C_W$ are calculated as [28]

$$ C_T^2 = \frac{3\Gamma(5/6) \sigma_T^2}{\sqrt{\pi} L_T^{2/3}} \qquad (6) $$

$$ C_W^2 = \frac{3\Gamma(5/6)}{\sqrt{\pi}} \left( \frac{\sigma_{v,s}^2}{L_{v,s}^{2/3}} + \frac{\sigma_{v,b}^2}{L_{v,b}^{2/3}} \right) \qquad (7) $$

where $\Gamma$ denotes the gamma function, $\sigma^2$ and $L$ denote the variance and length scale of the temperature fluctuations, subscript T, velocity fluctuations due to shear, subscript v, s, and due to buoyancy forces v, b. Finally, the variances and
length scales are modeled with the Monin–Obukhov similarity theory as [28]

$$ \sigma_T^2(z) = \frac{4.0T_*^2}{(1 - 10z/L_o)^{2/3}} \qquad \sigma_{v,s}^2 = 3.0u_*^2 \qquad \sigma_{v,b}^2 = 0.35w_*^2 \qquad (8) $$

$$ L_T(z) = 2.0z \frac{1 - 7z/L_o}{1 - 10z/L_o} \qquad L_{v,s}(z) = 1.8z \qquad L_{v,b} = 0.23z_i, \qquad (9) $$

where $T_* = -Q_H/(\rho c_p u_*)$ represents the surface layer temperature scale, $L_o = -u_*^3 T_s \rho c_p / (g \kappa_v Q_H)$ denotes the Obukhov length, and $w_* = (z_i g Q_H / (\rho c_p T_s))^{1/3}$ is the velocity scale for the buoyancy-induced turbulence. In the definitions of $T_*$, $L_o$ and $w_*$, $Q_H$ represents the surface heat flux, $u_*$, the friction velocity, $z_i$, the boundary layer height, $c_p = 1004$ J/(kg K), the specific heat at constant pressure, $T_s$, the surface temperature, $g = 9.81$ N/kg, the gravitational force per unit mass and $\kappa_v = 0.4$, the von Kármán constant. A typical value of the boundary layer height, $z_i$, during daytime and with a relatively flat terrain is 1000 m [28, 29]. The temperature scale, $T_*$, the friction velocity, $u_*$, and the Obukhov length are determined using the Harmonoise stability classes and the wind speed at 10 m above ground [30]. Furthermore, since $\sigma_T^2$, $L_T$ and $L_{v,s}$ are functions of the height above ground, $z$, we calculate their values at $z = (h_s + h_r)/2$, with $h_s$ and $h_r$, source and receiver heights, respectively.

To apply the Harmonoise model to a wind turbine, the rotor disk is discretized in the azimuthal direction in 360 stations, resulting in a total of 360xN elementary sources. The convective effects included in Amiet's theory are considered by evaluating the directivity and spherical divergence for the blade segment located in the present source position [31, 32]. This is necessary because of the relatively high Mach number of the tip segment of the blade ($\approx 0.2$). The Harmonoise model calculates the sound pressure levels $\tilde{L}_{p,if}$, where $if$ is the frequency index of the 1/3-octave frequency band (from 50 Hz to 10 kHz). These noise levels are used for the auralization methodology described in Section III.

### III. Auralization method

The sound pressure level at the receivers is auralized using a time domain spectral shaping synthesis. This technique transforms the 1/3-octave band frequency spectra into time domain signals. It is explained in Section III.A. Section III.B then presents the methodology to include the Doppler effect in the time domain signal as a time-varying time delay.

#### A. Time domain spectral shaping synthesis

The key idea is that stationary broadband noise can be obtained as the sum over the fractional 1/p-octave frequency bands of a band-filtered pink (or white) noise multiplied by the 1/p-octave amplitude of the noise levels. In this way, the spectrum of the pink (or white) noise is shaped to match the spectrum of the noise signal. Since the noise emission signal of a wind turbine is, in general, non-stationary, the target spectra are time-varying and updated for each time sample. The same technique has already been applied for the auralization of wind turbine noise by Pieren et al. [3]. In the present work, pink noise rather than white noise is chosen as it better matches the typical shape of wind turbine noise spectra which decrease with frequency.

For each segment of each blade, the pressure time signal $\tilde{p}$ (no Doppler effect is considered at this stage) is computed as

$$ \tilde{p}(t) = \sum_{i_f=1}^{N_b} p_{\text{rms},i_f}(t) \frac{n_{i_f}(t)}{n_{\text{rms}}} \qquad (10) $$

where $N_b$ is the number of 1/3-octave frequency bands, $n_{i_f}(t)/n_{\text{rms}}$ is the band-filtered normalized pink noise, and $p_{\text{rms},i_f}(t)$ is the root-mean-square (RMS) amplitude of the signal for each 1/3-octave band. The band-filtered pink noise is obtained from pink noise with zero-phase band filtering using 12th order Butterworth filters. The RMS amplitude of the signal is computed from the sound pressure levels $\tilde{L}_{p,if}$ as

$$ p_{\text{rms},i_f}(t) = p_0 \sqrt{10^{\tilde{L}_{p,if}(t)/10}} \qquad (11) $$

with $p_0 = 2 \cdot 10^{-5}$ Pa. Since $p_{\text{rms},i_f}(t)$ is sampled at a frequency related to the azimuthal discretization and $n_{i_f}(t)$ is sampled at the sampling frequency of the audio signal (e.g., $f_s = 44.1$kHz), $p_{\text{rms},i_f}(t)$ is re-sampled at $f_s$ by linear interpolation.
As an example, Fig. 4 shows the spectral representation of the spectral shaping synthesis applied to a given analytical power spectral density. In this figure, only the 1/3-octave bands from 400 to 5000 Hz are considered for better readability. The power spectral density of the final pressure time signal follows closely the analytical spectrum over the frequency range of interest.

<figure><img src="image_4.png" /><figcaption>Fig. 4 Spectral representation of the spectral shaping synthesis. The input power spectral density (PSD) is depicted in red; the PSDs of the band-filtered pink noise for the 1/3-octave bands of interest are depicted with black dashed lines; the black continuous line represents the PSD of the final pressure time signal.</figcaption></figure>

The procedure described above is applied to each blade segment for an arbitrary number of revolutions. The contributions from each segment are not summed yet because the Doppler effect is considered for each segment separately.

## B. Doppler effect and fractional delay interpolation

In the time domain, the Doppler effect is the result of the time-varying propagation time $\tau(t)$ due to the source-receiver relative motion. So, the pressure time signal $\tilde{p}(t)$ has to be delayed by an amount $\tau$ to include the Doppler effect. The Doppler-shifted pressure time signal $p(t)$ is computed as

$$p(t) = \tilde{p}(t - \tau(t)) \qquad (12)$$

where $\tau(t) = d(t - \tau(t))/c_0$; $d(t - \tau(t)) = \|\mathbf{x}_r - \mathbf{x}_c(t)\|$ is the distance between convected source position $\mathbf{x}_c(t)$ and the receiver $\mathbf{x}_r$. The convected source position is obtained by convecting the emission position $\mathbf{x}_e(t)$ at the Mach number of the fluid $\mathbf{M}_f$ (the wind speed) for the propagation time $\tau$, i.e. $\mathbf{x}_c(t) = \mathbf{x}_e(t) + \mathbf{M}_f c_0 \tau(t)$. It follows a second-order equation in $\tau$:

$$(c_0 \tau)^2 = \|\mathbf{x}_r - \mathbf{x}_e - \mathbf{M}_f c_0 \tau\|^2 \qquad (13)$$

which has the solution

$$c_0 \tau = \frac{-M_f R_f + \sqrt{M_f^2 R_f^2 + (1 - M_f^2) R^2}}{1 - M_f^2} \qquad (14)$$

with $R(t) = \|\mathbf{R}(t)\| = \|\mathbf{x}_r - \mathbf{x}_e(t)\|$ and $R_f(t) = \mathbf{R}(t) \cdot \frac{\mathbf{M}_f}{\|\mathbf{M}_f\|}$. A similar derivation for the propagation time can be found in the work of Sinayoko et al. [32], where the authors supposed $\mathbf{x}_c(t) \approx \mathbf{M}_f c_0 \tau(t)$ for an observer in the far-field and a reference frame centered in the rotor center. It is worth noting that the convected source position can be approximated with the emission position for a low Mach number of the wind. With this approximation, the propagation time is computed as

$$\tau(t) = \frac{\|\mathbf{x}_r - \mathbf{x}_e(t)\|}{c_0} \qquad (15)$$
It should be mentioned that $\mathbf{x}_e$ in Eqs. (14) and (15) is known only at emission time $t_e$ and not at the reception time $t$. Indeed, $\mathbf{x}_e = \mathbf{x}_e(t_e)$ is simply the position of the blade segment at $t_e = \Omega\psi$, with $\psi$ azimuthal position of the blade segment. Hence, we can re-write Eq. (12) as

$$p(t_e) = \tilde{p}(t_e - \tau(t_e))$$

where $\tau(t_e)$ is a known value that can be computed directly with Eq. (14).

As done previously for $p_{rms,if}(t)$, the propagation time $\tau$ is re-sampled at the audio frequency, but, in this case, with 1/3-order polynomials. A third-order polynomial is used to ensure the second derivative of the delay $\tau(t)$ to be continuous as suggested by Maillard [33]. This yields a continuous Doppler shift and avoids audible artifacts.

In practice, Eq. (12) cannot be applied in that form because $\tilde{p}$ is known only at the discrete time samples $t_n$ and, in general, $\tau(t_n)$ is not a multiple of the time step $\Delta t = 1/f_s$, where $f_s = 44.1$ kHz is the sampling frequency. So, the signal must be interpolated. Following the work of Laakso et al. [34], a third-order Lagrange interpolator is used to obtain sufficient accuracy at high frequency. The first step is to separate the delay $\tau$ in his integer and fractional part, $D$ and $d$, respectively:

$$\tau(t_n) = D + d \qquad (16)$$

The stencil of points used for the interpolation is determined using the integer part:

$$\tilde{p}_1 = \tilde{p}(t_n - D - 3) \quad \tilde{p}_2 = \tilde{p}(t_n - D - 2) \quad \tilde{p}_3 = \tilde{p}(t_n - D - 1) \quad \tilde{p}_4 = \tilde{p}(t_n - D)$$

The interpolated sample $p(t_n) = \tilde{p}(t_n - \tau(t_n))$ is computed as

$$p(t_n) = \sum_{i=1}^{4} \tilde{p}_i h_i \qquad (17)$$

where the weights $h_i$ are computed from the fractional part as:

$$h_1 = d(d-1)(d-2)/6 \quad h_2 = -d(d-1)(d-3)/2 \quad h_3 = d(d-2)(d-3)/2 \quad h_4 = -(d-1)(d-2)(d-3)/6$$

Finally, the contributions from each segment are summed to obtain the pressure time signal of one blade. The other two blades are considered by delaying the pressure time signal of one blade by 1/3 and 2/3 of the rotation period. The procedure explained for a time-dependent delay is applied also for these two constant time delays. The noise of the complete wind turbine is obtained by summing the contribution from the three blades.

It might be argued that frequency-dependent propagation effects must be applied to the Doppler-shifted signal, which means that the Doppler effect must be considered first and applied to the source signal. However, thanks to the fact that the methodology is based on 1/3-octave bands, the frequency-dependent propagation effects, and the Doppler shift are commutative operations: in the present formulation, the propagation effects $h_{rms,if}$ are simple multiplication factors. Indeed, if the Doppler effect is considered first, the source signal in each frequency band, $\tilde{y}_{if}(t)$, has to be delayed,

$$y_{if}(t) = \tilde{y}_{if}(t - \tau) = A_{rms,if}(t - \tau) \frac{n(t - \tau)}{n_{rms}}$$

with $A_{rms,if}$, the gain due to the sound power and the directivity of the source. Then, the propagation effects are considered at emission time

$$p(t) = \sum_{i_f=1}^{N_b} h_{rms,if}(t - \tau) y_{if}(t)$$

This approach is equivalent to the approach adopted in the present work, that is computing the signal at reception time

$$\tilde{p}(t) = \sum_{i_f=1}^{N_b} h_{rms,if}(t) A_{rms,if}(t) \frac{n_{if}(t)}{n_{rms}}$$

and delaying it afterward, i.e. Eq. (12). An alternative approach for which the Doppler effect must be considered first is the case in which the impulse response at emission position, $h(t)$, is convolved with the Doppler shifted source signal,
<figure><img src="image_5.png" /><figcaption>Fig. 5 Mesh used for the tip segment airfoil of the SWT 2.3-93 wind turbine. (a) View of the entire airfoil. (b) Close-up view of the trailing-edge. The red lines close to the trailing-edge correspond to the extraction locations for the boundary layer parameters, i.e. $x/c = 0.99$.</figcaption></figure>

y(t):

$$p(t) = \int_{-\infty}^{\infty} h(t)y(t-t')dt'$$
$$y(t) = \tilde{y}(t-\tau) = A(t-\tau)\frac{n(t-\tau)}{n_{\text{rms}}}$$

The approach suggested in this work simplifies the algorithm since the Doppler shift is applied to the signal generated by each blade segment, instead of to each frequency band of each blade segment. Furthermore, the source power and directivity, $A_{rms,if}$, and the propagation effect, $h_{rms,if}$, can be multiplied and stored as a single variable until the application of the Doppler effect.

## IV. Application to a 2.3 MW wind turbine

The blade of the SWT 2.3-93 wind turbine has been reconstructed using the geometry and the airfoils from references [35, 36]. As shown in Fig. 1, the blade is divided into six segments. The five outer segments have an aspect ratio of 3 to satisfy the large aspect ratio assumption of Amiet's theory. The blade segment closer to the root has an aspect ratio of 4.2 to include the blade up to the radial position $R = 9$ m, corresponding to the first cross-section with an airfoil shape. The operating conditions are chosen to match the acoustic measurements of [37]: the rotational speed (RPM) is 14, the wind speed is 8 m/s and the collective pitch angle is $-2^\circ$. An integral length scale of 312 m and a turbulence intensity of 10.7% of the wind speed have been chosen for the leading-edge noise model.

Figure 5 shows the 2D mesh generated with Simcenter STAR-CCM+ for the airfoil of the tip segment. Each mesh counts approximately 67 000 elements. The structured layer counts 40 cells in the wall-normal direction resulting in a $y^+ < 0.5$. The mesh is refined close to the trailing-edge to calculate accurately the boundary layer parameters. A grid sensitivity study has been done in previous works [38], showing a convergence of the boundary layer parameters for even coarser meshes and a low dependency of the boundary layer parameters on the mesh element size.

Figure 6 shows the boundary layer parameters extracted from the RANS simulation and used for the calculation of the wall pressure spectra, $\Phi_{pp}$. The boundary layer thickness $\delta$, the displacement thickness $\delta^*$, the momentum thickness $\Theta$, the friction coefficient $C_f$, and the dimensionless pressure gradient $dC_p/dx$ are displayed as a function of the radial position $r$, along the blade span, where $r = 0$ corresponds to the root of the blade and $r = 46.5$ m to the tip. $C_p$ denotes the pressure coefficient. On the suction side, the pressure gradient is positive, i.e., adverse, and Lee's wall pressure spectrum model [26] is applied, whereas, on the pressure side it is negative, i.e., favorable, and Goody's model [27] is used.

For far-field noise propagation, the wind turbine is placed on flat ground with constant acoustic properties. The flat ground is representative of the wind turbine test site [37]. The results are presented for two ground types. Following the norm NF S31-133, we will refer to ground type D to indicate a ground with flow resistivity $\sigma = 200$ kPa s/m$^2$, which is representative of an uncompacted, loose ground, such as grass. We will use ground type H for a flow resistivity
<figure><img src="image_6.png" /><figcaption>Fig. 6 Boundary layer parameters extracted at $x/c = 0.99$ from the 2D RANS simulations. (a) Boundary layer thickness $\delta$. (b) Displacement thickness $\delta^*$ and momentum thickness $\Theta$. (c) Friction coefficient $C_f$. (d) Dimensionless pressure gradient $dC_p/dx$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 7 (a) OASPL and (b) 1/3-octave noise levels in free-field conditions. The levels calculated with the present auralization methodology are depicted with a continuous curve, whereas the analytical results computed with Schlinker and Amiet [31]’s theory are shown with a dashed curve.</figcaption></figure>

$\sigma = 20000 \text{ kPa s/m}^2$, representative for the hard measurement board used in the experimental campaign [37]. Results will be shown in terms of 1/3-octave spectra averaged over one revolution or as a function of time. The average levels are computed as follows. After the A-filtering (when performed), the signal $p(t)$ is filtered with 8th order bandpass Butterworth filters to obtain the pressure time signal $p_{if}$ for each 1/3-octave band. The noise levels are then computed for each 1/3-octave band as

$$L_{p,if} = 10 \log_{10} \left( \frac{\frac{1}{T} \int_{0}^{T} p_{if}^{2}(t) dt}{p_{0}^{2}} \right) \quad (18)$$

where $T$ is the length of the time signal. The time-dependent levels in 1/3-octave (i.e. the 1/3-octave spectrogram) are computed as follows: 1) the pressure signal for each 1/3-octave band is computed band-pass filtering the time domain signal with 8th order Butterworth filters; 2) the time-dependent noise levels are computed as in Eq. (18) with $T = T_p/100$ and $T_p$, the time for one rotation; 3) the noise levels are phase-averaged over 40 rotations. For the calculation of the overall sound pressure level (OASPL), step 1) is skipped and Eq. (18) is applied to the complete time signal. The signal processing is performed in Python with the `acoustics` package (version 0.2.4).

At first, the auralization methodology is validated in free-field with the analytical results computed with the frequency domain approach described by Sinayoko et al. [32] and originally proposed by Schlinker and Amiet [31]. For this purpose, the excess attenuation computed by Harmonoise is not considered. The atmospheric absorption is also neglected as it is not included in Schlinker and Amiet [31]’s theory. Figure 7 shows the OASPL and the 1/3-octave sound pressure levels for three bands phase-averaged over 40 rotations. The receiver is at ground level and close to the crosswind direction (descending side of the turbine). The horizontal distance from the turbine is 100 m and the angle with the crosswind direction is 10° towards downwind. The OASPL and the 1/3-octave levels follow closely the analytical solution. An overshoot by 0.5 dB of the lower peaks is noticed at high frequency. At low frequency (63 Hz) the 1/3-octave signal is noisier; a smoother curve could be obtained by phase-averaging over more revolutions. The difference (not shown) between the averaged 1/3-octave noise levels computed with Schlinker and Amiet [31]’s theory and with the present approach is less than 0.5 dB. This difference is not audible and is considered negligible for the scope of the present work. The amplitude modulation is noticeable in the mid and high-frequency range. Instead, it is less important at low frequencies.

Then, the results are compared against the measurements of Leloudas [37]. Figure 8 shows the 1/3-octave noise levels, 100 meters downwind of the wind turbines, on the ground. The two black dashed lines are obtained by adding and subtracting 3 dBA from the measured spectrum. The 1/3-octave noise levels referred to as *Harmonoise* are obtained from the pressure time signals of 40 revolutions of the blades. The Harmonoise meteorological model as described by Defrance et al. [30], with 5 m/s wind speed (10 m above ground) and stability class 1, representative for a sunny day, has
been used to obtain a sound speed profile $c(z) = az + b \log(z/z_0) + c_0$, with $a = -0.0109$ 1/s, $b = 0.1004$, $z_0 = 0.1$ m and $c_0 = 343$ m/s. The sound speed profile is linearized, as explained by Salomons et al. [14] for high sources. Finally, the spectrum referred to as *Free-field Amiet* solution in Fig. 8 is computed with the free-field frequency domain approach of references [31, 32]. The atmospheric absorption is not included in *Free-field Amiet* solution, whereas it is considered for the *Harmonoise* one. Furthermore, 6 dBA are added to the *Free-field Amiet* solution to account for the hard measurement board. As expected, *Harmonoise* is closer to the measurements at high frequencies where the atmospheric absorption is important whereas Amiet's and the Harmonoise solutions overlap at low frequencies. The *Harmonoise* spectrum is within the ±3 dBA error band for most of the 1/3-octave bands.

<figure><img src="image_8.png" /><figcaption>**Fig. 8** Third-octave noise spectra for an observer on the ground at 100 m downwind of the wind turbine.</figcaption></figure>

Figure 9 shows the 1/3-octave spectrogram, the OASPL, and the time signal of the synthesized sound for a receiver at 1.5 m above the ground and close to the crosswind direction. The exact position in the horizontal plane is the same as the receiver of Fig. 7. In this case, ground type D has been used as closer to the real environment. The amplitude modulation ("swishing") is audible and visible in the spectrogram due to the position of the observer close to the rotor plane and the relatively small distance from the tower. As explained by Oerlemans and Schepers [39], the swishing is caused by the trailing-edge noise directivity function which is mainly emitting sound towards the leading-edge of the blade. Furthermore, the convective amplification factor enhances this phenomenon. In Fig. 10, the angle between the crosswind direction receiver is 50°. In this case, the "swishing" is less audible and more difficult to identify in the spectrogram.

## V. Conclusions

This paper presents a cost-efficient methodology for wind turbine noise simulations and physics-based auralization. The time-dependent aerodynamic noise levels are computed by coupling the semi-analytical RANS-informed Amiet's theory with the Harmonoise model. The approach is based on 2D RANS around the blade airfoils are used to accurately predict the boundary layer parameters necessary to calculate the free-field trailing-edge noise with Amiet's theory. Leading-edge noise is also considered. The Harmonoise model is exploited for the propagation of noise in the atmosphere by modeling the wind turbine as an extended source. Finally, an algorithm is proposed to accurately reconstruct the time domain signal, including the Doppler effect. The advantage of the described methodology is that the calculation of the sound sources with Amiet's theory is fully decoupled from the far-field propagation part. Furthermore, thanks to the low computational cost, this methodology can be applied to complex urban and rural scenes, which can potentially include multiple turbines, buildings, and other noise sources such as road or railway traffic. The implementation of the wind turbine model and of the reconstruction of the time-domain signal is validated with the analytical strip theory [31] for free-field propagation. The accuracy of the average noise levels is assessed through comparison with measurements available in the literature for a 2.3 MW wind turbine. The 1/3-octave noise levels are within ±3 dB for most of the frequency bands. This agreement is considered satisfactory considering the uncertainties related to the atmospheric conditions, and to the reconstruction of the wind turbine blade geometry. Finally, we present the spectrograms and the pressure signals for a receiver located close to the crosswind direction. As already noticed in previous works, the
<figure><img src="image_9.png" /><figcaption>Fig. 9 Third-octave spectrogram, OASPL and time signal of the synthesized noise for a receiver close to the crosswind direction. The "swishing" is clearly audible and visible in the spectrogram.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10 Third-octave spectrogram, OASPL and time signal of the synthesized noise for a receiver far from the crosswind direction. The "swishing" is barely audible.</figcaption></figure>
amplitude modulation (*swishing*) associated with the trailing-edge noise directivity is observed in the spectrograms and clearly audible in the audio signals.

## Acknowledgments

This research was supported by the European Commission through the H2020-MSCA-ITN-209 project zEPHYR (grant agreement No 860101). We thank J. Christophe and M. Roger for sharing reference [40].

## References

[1] Janssen, S. A., Vos, H., Eisses, A. R., and Pedersen, E., "A comparison between exposure-response relationships for wind turbine annoyance and annoyance due to other noise sources," *The Journal of the Acoustical Society of America*, Vol. 130, No. 6, 2011, pp. 3746-3753.

[2] The Expert Panel on Wind Turbine Noise and Human Health, *Understanding the Evidence: Wind Turbine Noise*, Council of Canadian Academies, Ottawa, Canada, 2015.

[3] Pieren, R., Heutschi, K., Müller, M., Manyoky, M., and Eggenschwiler, K., "Auralization of wind turbine noise: Emission synthesis," *Acta Acustica united with Acustica*, Vol. 100, No. 1, 2014, pp. 25-33. https://doi.org/10.3813/AAA.918683.

[4] Mascarenhas, D., Cotté, B., and Doaré, O., "Physics-based auralization of wind turbine noise," *9th International Conference on Wind Turbine Noise*, 2021, pp. 1-6.

[5] Mascarenhas, D., Cotté, B., and Doaré, O., "Synthesis of wind turbine trailing edge noise in free field," *JASA Express Letters*, Vol. 2, 2022, p. 033601. https://doi.org/10.1121/10.0009658.

[6] Brooks, T. F., Pope, S. D., and Marcolini, M. A., "Airfoil self-noise and prediction," Tech. rep., NASA, Langley Research Center, 1989.

[7] Amiet, R. K., "Acoustic radiation from an airfoil in a turbulent stream," *Journal of Sound and Vibration*, Vol. 41, No. 4, 1975, pp. 407-420. https://doi.org/10.1016/S0022-460X(75)80105-2.

[8] Amiet, R. K., "Noise due to turbulent flow past a trailing edge," *Journal of Sound and Vibration*, Vol. 47, No. 3, 1976, pp. 387-393. https://doi.org/10.1016/0022-460X(76)90948-2.

[9] McBride, S., Burdisso, R., and Parra, J. D., "An efficient noise modeling tool for wind turbines including sound propagation in arbitrary weather conditions," *ICSV 2016 - 23rd International Congress on Sound and Vibration: From Ancient to Modern Acoustics*, 2016.

[10] Ecotière, D., "Can we really predict wind turbine noise with only one point source?", 6th International Conference on Wind Turbine Noise, Glasgow, 2015.

[11] Elsen, K. M., and Schady, A., "Different sound source setups in the simulation of wind turbine sound propagation," *9th International Conference on Wind Turbine Noise*, Remote from Europe, 2021, pp. 1-6.

[12] Cotté, B., "Coupling of an aeroacoustic model and a parabolic equation code for long range wind turbine noise propagation," *Journal of Sound and Vibration*, Vol. 422, 2018, pp. 343-357. https://doi.org/10.1016/j.jsv.2018.02.026; URL http://dx.doi.org/10.1016/j.jsv.2018.02.026.

[13] Sessarego, M., and Shen, W. Z., "Noise Propagation Calculations of a Wind Turbine in Complex Terrain," *Journal of Physics: Conference Series*, Vol. 1452, No. 1, 2020. https://doi.org/10.1088/1742-6596/1452/1/012063.

[14] Salomons, E., Van Maercke, D., Defrance, J., and De Roo, F., "The Harmonoise sound propagation model," *Acta Acustica united with Acustica*, Vol. 97, No. 1, 2011, pp. 62-74. https://doi.org/10.3813/AAA.918387.

[15] van Maercke, D., and Defrance, J., "Development of an Analytical Model for Outdoor Sound Propagation Within the Harmonoise Project," *Acta Acustica united with Acustica*, Vol. 93, 2007, pp. 201-212.

[16] Maillard, J., and Jagla, J., "Auralization of Urban Traffic Noise - Quantitative and Perceptual Validation," *Proc. of Congrès Français d'Acoustique (CFA)*, 22-25 April 2014, Poitiers, France, 2014.

[17] Maillard, J., and Kacem, A., "Auralization Applied to the Evaluation of Pedestrian and Bike Paths in Urban Environments," *Proc. of Internoise 2016, Hamburg, Germany, August 21-24, 2016*.
[18] Maillard, J., Kacem, A., Martin, N., and Faure, B., "Physically-based auralization of railway rolling noise," *Proceedings of the 23rd International Congress on Acoustics, 9-13 September, Aachen, Germany*, 2019.

[19] *Aerodynamics of Wind Turbines*, John Wiley & Sons, Ltd, 2009, Chap. 3, pp. 91–155. https://doi.org/https://doi.org/10.1002/9781119994367.ch3, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119994367.ch3.

[20] Christophe, J., Le Bras, S., Chavez, R., Vidal, J., Panjwani, B., and Samseth, J., "Application of the coupled atmospheric-aeroacoustic simulation platform to a wind farm benchmark," Tech. rep., UPWARDS Deliverable D4.5, 2022.

[21] Siemens Digital Industries Software, "Simcenter STAR-CCM+, version 2021.3," Siemens 2021.

[22] Siemens Digital Industries Software, "Simcenter STAR-CCM+ User Guide v. 2021.1," Siemens 2021.

[23] Roger, M., and Moreau, S., "Back-scattering correction and further extensions of Amiet's trailing-edge noise model. Part 1: theory," *Journal of Sound and Vibration*, Vol. 286, No. 3, 2005, pp. 477–506. https://doi.org/10.1016/j.jsv.2004.10.054.

[24] Bresciani, A. P., Le Bras, S., and de Santana, L. D., "Generalization of Amiet's theory for small reduced-frequency and nearly-critical gusts," *Journal of Sound and Vibration*, Vol. 524, 2022, p. 116742. https://doi.org/https://doi.org/10.1016/j.jsv.2021.116742.

[25] Corcos, G. M., "The structure of the turbulent pressure field in boundary-layer flows," *Journal of Fluid Mechanics*, Vol. 18, No. 3, 1964, p. 353–378. https://doi.org/10.1017/S002211206400026X.

[26] Lee, S., "Empirical wall-pressure spectral modeling for zero and adverse pressure gradient flows," *AIAA Journal*, Vol. 56, No. 5, 2018, pp. 1818–1829. https://doi.org/10.2514/1.J056528.

[27] Goody, M., "Empirical spectral model of surface pressure fluctuations," *AIAA Journal*, Vol. 42, No. 9, 2004, pp. 1788–1794. https://doi.org/10.2514/1.9433.

[28] Ostashev, V., and Wilson, D., *Acoustics in Moving Inhomogeneous Media*, Taylor & Francis (2015), 2016.

[29] Kamrath, M. J., Ostashev, V. E., Wilson, D. K., White, M. J., Hart, C. R., and Finn, A., "Vertical and slanted sound propagation in the near-ground atmosphere: Amplitude and phase fluctuations," *The Journal of the Acoustical Society of America*, Vol. 149, 2021, pp. 2055–2071. https://doi.org/10.1121/10.0003820.

[30] Defrance, J., Salomons, E., Stuijt, I., Heimann, D., Plovsing, B., Watts, G., Jonasson, H., Zhang, X., Premat, E., Schmich-Yamane, I., Aballéa, F.-E., Baulac, M., and Roo, F., "Outdoor Sound Propagation Reference Model Developed in the European Harmonoise Project," *Acta Acustica united with Acustica*, Vol. 93, 2007, pp. 213–227.

[31] Schlinker, R. H., and Amiet, R. K., "Helicopter Rotor Trailing Edge Noise," Tech. rep., National Aeronautics and Space Administration, 1981.

[32] Sinayoko, S., Kingan, M. J., and Agarwal, A., "Trailing edge noise theory for rotating blades in uniform flow," *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, Vol. 469, No. 2157, 2013. https://doi.org/10.1098/rspa.2013.0065.

[33] Maillard, J., "Prediction and auralization of construction site noise," *Proc. of EURONOISE, 26-28 October 2009, Edinburgh, Scotland*, 2009.

[34] Laakso, T., Valimaki, V., Karjalainen, M., and Laine, U., "Splitting the unit delay [FIR/all pass filters design]," *IEEE Signal Processing Magazine*, Vol. 13, No. 1, 1996, pp. 30–60. https://doi.org/10.1109/79.482137.

[35] Churchfield, M., "A method for designing generic wind turbine models representative of real turbines and generic Siemens SWT-2.3-93 and Vestas V80 specifications," *National Renewable Energy Laboratory, Golden, Colorado*, 2013.

[36] Christophe, J., Buckingham, S., Schram, C., and Oerlemans, S., "zEPHYR - Large On Shore Wind Turbine Benchmark," Data set, Mar. 2022. https://doi.org/10.5281/zenodo.6380879, URL https://doi.org/10.5281/zenodo.6380879.

[37] Leloudas, G., "Optimization of Wind Turbines with respect to Noise," master's thesis, Technical University of Denmark, Technical University of Denmark, Anker Engelunds Vej 1, Building 101A, 2800 Kgs. Lyngby, Nov. 2006.

[38] Bresciani, A. P. C., Boatto, U., Bras, S. L., Bonnet, P., and de Santana, L. D., "Influence of blade deflections on wind turbine noise directivity," *Journal of Physics: Conference Series*, Vol. 2257, 2022, p. 012012. https://doi.org/10.1088/1742-6596/2257/1/012012, URL https://iopscience.iop.org/article/10.1088/1742-6596/2257/1/012012.
[39] Oerlemans, S., and Schepers, J. G., "Prediction of wind turbine noise and validation against experiment," *International Journal of Aeroacoustics*, Vol. 8, No. 6, 2009, pp. 555-584.

[40] Amiet, R. K., "Frame of Reference Consideration for the Forward Flight Noise Problem," Tech. rep., United Aircraft Research Laboratories, 1974.