© The Author(s), Published by EDP Sciences, 2023

1 Introduction

Wind energy is a renewable energy source and the steady rise in the number of wind farms can be seen evidently [1]. However one of the concerning downsides to the installation of wind turbines is the impact of the noise generated by the moving blades. This leads to potential causes of annoyance and disturbances in sleep patterns for people living in the vicinity of wind farms [2, 3]. Other noise sources from wind turbines such as blade-tower interaction producing tonal noise [4, 5], mechanical noise from the gearbox and generator [6] may also be disturbing in some circumstances [7]. The trade-off between the necessity for an energy source and concerns on the generated noise emphasizes the importance of the development for a synthesis tool for wind turbine noise. Characteristics of wind turbine noise due to the rotational motion of the blade such as swish, thump or more generally amplitude modulation [8] are not well described by noise level predictions. A physics-based synthesis tool would be useful to have control on the prevailing physical parameters of the wind turbine system and accordingly predict the generated noise and study its impact in the design phase of the wind farms. Furthermore, accurately synthesizing wind turbine noise signals can be used to study psychoacoustic questions such as sleep disturbance and annoyance caused by wind turbine noise.

The synthesis of wind turbine noise has lately been of increasing interest. However, the existing synthesis models are seen to be either sample-based [9] or to consider only the source [10]. The restriction arising from the sample-based synthesis is that the resultant signals are dependent on the prevailing parameters of the obtained sample. The synthesis can be difficult to extend for other settings that are not observed in the base samples [9]. Furthermore, the synthesis is directly based on the recorded noise which may be contaminated with background noise and wind noise [11, 12]. The synthesis model of Lee et al. [10] which considers only the source is based on computing the aerodynamic noise directly in the time domain based on the Ffocws Williams and Hawkings analogy. The drawback of this model is that it does not include the effect of propagation which is known to change the sound considerably. Analytical models are commonly used in the context of auralization to account for the propagation effect such as ground reflection, effect of the turbulence scattering etc. [13–15] as they are rather simple and computationally effective.

The broadband wind turbine noise is mainly due to aeroacoustic mechanisms of the rotating blade edges interacting with the turbulent flow fluctuations. Other important components that influence the received noise are the propagation of this sound through the atmosphere and its reflection off an impedance ground. The sound received from a fixed point source radiating in an outdoor environment depends on the source characteristics such as power and directivity and also on the propagation effects observed. Similarly, a wind turbine blade can be modelled as a set of moving points sources whose noise levels depend on the propagation environment. Various frequency-domain models of wind turbine noise sources have been proposed in the literature and it is known that the main sources of broadband noise generated by the wind turbine blades are turbulent inflow noise (TIN), trailing edge noise (TEN) and stall noise [16–18]. Among these noise sources TIN and TEN occur for all angles of attack while stall noise is the most intense but is produced only when the angle of attack is large. Thus, it can be said that in the simplest case scenario, TIN and TEN are observed as dominant sources of wind turbine noise.

To account for the propagation effects through the atmosphere the parabolic equation methods have been used widely [19–21] and have also been applied to wind turbine noise models [22–25]. The moving monopole approach described by Cotté [24] is reported to work well using the classical parabolic equation (PE) with the influence of the wind speed gradient accounted for by considering the effective sound speed. However, based on a new formulation proposed by Ostashev et al. [26] the effect of the wind speed gradient can be included in the PE model without using the effective sound speed approximation. This allows one to accurately account for wind speed gradients, although a recent paper has shown that the error associated with the effective sound speed approximation is small in the context of wind turbine noise [27]. Another important effect when an acoustic shadow zone is present is the scattering due to turbulence in the atmosphere. This is accounted for in the Harmonoise model [28] which we implement and study here.

In this article, we synthesize noise from a 2.3 MW wind turbine over a flat ground, while including the refraction effect due to the wind speed gradient, ground reflection and scattering due to turbulence. To study the effects of atmospheric propagation on wind turbine noise we benefit from the PE for a moving medium [26], the Harmonoise model [28] with the moving monopole model [24] coupled with a synthesis tool for wind turbine noise [29]. To our knowledge, this is the first time that wind turbine noise signals are synthesized from physical models of the aeroacoustic sources and the atmospheric propagation. We first describe briefly in Section 2.1, the source model for the wind turbine noise which is used to obtain the noise spectra in free field. The propagation modeling of the wide-angle PE in a moving medium (WAPE-MM) which is now used in the moving monopole model and its dependence on few physical parameters is then explained in Section 2.2 along with the implementation of the Harmonoise model. The synthesis method used to generate the wind turbine noise signals [29] is reviewed in Section 2.3. In Section 3, the configurations used in the study are described and the associated propagation effects are studied for three source heights. Finally, audio signals are provided for a few test cases in Section 4, in order to emphasize various source and propagation physical phenomena.

2 Physics-based synthesis of wind turbine noise

The sound pressure level (SPL) in dB due to a point source in the far-field can be estimated by considering the noise level obtained in free field that is corrected by the propagation effects. For the wind turbine system, the noise from the blades can be modelled as multiple point sources along the length of the blade, whose individual generated sound propagates through the atmosphere [16, 24]. The blades of the wind turbine are divided into a number of segments whose contributions to the total sound are summed incoherently for a receiver at (x_R, z_R) in the propagation plane that makes an angle θ with respect to the wind direction (see Fig. 1a). As a flat ground is considered, the elevation is zero. The rotation of the blades are accounted for by considering discrete angular positions γ (Fig. 1b). The SPL for a fixed receiver due to a segment at a particular angular position γ modelled as a point source is thus given as [30]:

$$\mathrm{SPL}\left(\omega ,\gamma \right)={\mathrm{SPL}}_{\mathrm{FF}}\left(\omega ,\gamma \right)+\Delta L\left(\omega ,\gamma \right)-{\alpha }_{\mathrm{abs}}\left(\omega \right)R\left(\gamma \right),$$(1)where ω = 2πf is the angular frequency, SPL_FF is the sound pressure level observed at the receiver position in free field, ΔL is the sound pressure level relative to free field, α_abs is the atmospheric absorption coefficient [30] and R is the direct distance between the point source and the receiver. SPL_FF is characteristic of the source while (ΔL − α_absR) includes the atmospheric propagation effects. The source modelling for the TEN and TIN mechanisms at each blade segment is reviewed in Section 2.1 and the propagation model is explained in Section 2.2 followed by the synthesis method in Section 2.3.

Figure 1 Schematics for (a) the modelled wind turbine blades with the receiver represented as a sphere, and (b) the segmented blade in motion represented as a set of point sources whose heights are interpolated to the nearest neighbor in the propagation calculation based on a set of Nh source heights (Nh = 5 here).

2.1 Source model

The wind turbine noise source model developed by Tian and Cotté [16] predicts the trailing edge noise (TEN) and turbulent inflow noise (TIN) generated by a segmented wind turbine blade based on Amiet’s theory. For TEN, the power spectral density (PSD) of the acoustic pressure at a far-field observer (x, y, z) produced by an airfoil with large span L to chord c aspect ratio (L > 3c) is given by [31, 32]:

$$S_{\mathrm{pp}}^{\mathrm{TEN}}(x,y,z,\omega )={\left(\frac{\mathrm{kcz}}{4\pi S_0^2}\right)}^{2}2L{\mathrm{\Phi }}_{\mathrm{pp}}\left(\omega \right)l_y\left(\omega ,\frac{\mathrm{ky}}{S_0}\right)\times {\left|{\mathcal{L}}_{\mathrm{TE}}\left(x,\frac{\omega }{U_{\mathrm{c}}},\frac{\mathrm{ky}}{S_0}\right)\right|}^2,$$(2)where k is the acoustic wavenumber, S₀ is the convection-corrected distance between the source and the observer [32], Φ_pp is the wall pressure fluctuation spectrum (WPS), l_y is the spanwise correlation length estimated by Corcos model [33], U_c is the convected velocity and ${\mathcal{L}}_{\mathrm{TE}}$ is the aeroacoustic transfer function for TEN. For the pressure side of the airfoil, the WPS is calculated using Goody’s model [34] while for the suction side of the airfoil, the WPS is calculated using Lee’s model [35], which is an improved version of Rozenberg’s model [36] for strong adverse pressure gradient. The boundary layer parameters necessary for the calculation of the WPS are obtained using XFoil calculations at 99% of the chord.

In the case of TIN, the PSD of the acoustic pressure at a far-field observer position (x, y, z) produced by a large aspect-ratio airfoil is given by [32, 37]:

$$S_{\mathrm{pp}}^{\mathrm{TIN}}(x,y,z,\omega )={\left(\frac{{\rho }_0\mathrm{kcz}}{2S_0^2}\right)}^2\pi U_0\frac{L}{2}{\mathrm{\Phi }}_{\mathrm{ww}}\left(\frac{\omega }{U_0},\frac{\mathrm{ky}}{S_0}\right)\times {\left|{\mathcal{L}}_{\mathrm{TI}}\left(x,\frac{\omega }{U_0},\frac{\mathrm{ky}}{S_0}\right)\right|}^2,$$(3)where ρ₀ is the air density, U₀ is the inflow velocity of the wind, Φ_ww is the two-dimensional energy spectrum of the velocity fluctuations and ${\mathcal{L}}_{\mathcal{TI}}$ is the aeroacoustic transfer function for TIN. For the turbulent inflow, the Kolmogorov energy spectrum is used:

$${\mathrm{\Phi }}_{\mathrm{ww}}(k_1,k_2)=\frac{4}{55\sqrt{\pi }}\frac{\mathrm{\Gamma }\left(\frac{1}{3}\right)}{\mathrm{\Gamma }\left(\frac{5}{6}\right)}{C}_K{ϵ}^{\frac{2}{3}}(k_1^2+k_2^2{)}^{-\frac{4}{3}},$$(4)where k₁ and k₂ are the turbulent wavenumbers in the chordwise and spanwise directions, C_K ≈ 2.0 is the Kolmogorov constant [38, 39], Γ is the Gamma function and ϵ is the turbulence dissipation rate (TDR). As the Kolmogorov spectrum is dependent only on one parameter it is preferred to the von Kármán spectra in the presented model. As noted by Buck et al. [39], the Kolmogorov spectrum does not differ from the von Kármán spectrum within the audible frequency range for typical atmospheric conditions.

The frequency-domain model approximates the complete rotation of the blade as a series of translations between discrete angular positions γ. The PSD obtained from equations (2) and (3) that include the noise directivity of the airfoils can be used to calculate the response of each blade segment at (x, y, z) corresponding to the angular position γ. The segmentation is done ensuring the segment span is greater than the spanwise turbulence correlation length so that the individual contribution of the segments can be summed in an incoherent manner at the receiver. However, the formulation of the PSD is done for a non-rotating airfoil. Following Sinayoko et al. [40], to account for the convective amplification and Doppler effect associated with the blade motion, the instantaneous PSD for a moving airfoil for a fixed receiver position is given by:

$$S_{\mathrm{pp}}\left({\mathbf{x}}_{\mathbf{0}},\omega \right)=\left(\frac{{\omega }_{\mathrm{e}}}{\omega }\right)S_{\mathrm{pp}}^{\mathrm{N}}\left(\mathbf{x},{\omega }_e\right),$$(5)where $S_{\mathrm{pp}}^{\mathrm{N}}(\mathbf{x},{\omega }_{\mathrm{e}})$ is the PSD of the TEN or TIN for a fixed blade, ω_e and ω are the emitted and observed frequencies, x₀ and x are the observer coordinates in the hub (stationary frame) and blade (moving frame) coordinate systems. With the reference pressure p_ref = 20 μPa, the sound pressure level in free field SPL_FF can be obtained using equation (5) and the relation:

$${\mathrm{SPL}}_{\mathrm{FF}}=10\mathrm{lo}{\mathrm{g}}_{10}\left(\frac{S_{\mathrm{pp}}\cdot 1 \mathrm{Hz}}{p_{\mathrm{ref}}^2}\right).$$(6)

Thus for a receiver with coordinates defined by x_R, z_R and θ (Fig. 1a), using the relations equations (2)–(6), the noise in free field from each segment is obtained at each discrete blade angular positions.

2.2 Propagation model

The propagation effects that have an influence over a flat impedance ground in the outdoor environment involve refraction due to the wind speed and temperature profiles, ground reflection, atmospheric absorption and scattering due to turbulence in the atmosphere [30]. Based on the implemented models, the effects of the propagation can be included in the term ΔL − α_absR in equation (1). We present in Section 2.2.1 the parabolic equation in moving medium that includes the effects of refraction and ground reflection. We consider an empirical model in Section 2.2.2 to account for the effect of scattering due to turbulence to eventually obtain the term ΔL.

2.2.1 Refraction and ground reflection using the parabolic equation in a moving medium

For a sound wave propagating through a medium that is moving with an arbitrary velocity U_x(z), Ostashev et al. [26] proposes a formulation which can be solved using the Crank-Nicholson (CN) algorithm with a Padé (1,1) approximation. In a two-dimensional vertical plane (x, z), assuming that the air density is a constant equal to ρ₀, equations (27) and (39) of reference [26] for the sound pressure $\widehat{p}(x,z)$ and the scaled velocity potential $\widehat{\varphi }(x,z)$ in the frequency domain reduce to:

$$\widehat{p}\left(x,z\right)=\left(1+\frac{i{M}_x}{k_0}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\right)\widehat{\varphi }\left(x,z\right),$$(7)

$$\left(\frac{\mathrm{\partial }}{\mathrm{\partial }x}-i{k}_0{\zeta }_x^2\sqrt{1+\eta +\widehat{\mu }}+i{k}_0\widehat{\tau }\right)\widehat{\varphi }\left(x,z\right)=0,$$(8)where M_x(z) = U_x(z)/c₀ is the Mach number of the medium moving in the x direction, k₀ = ω/c₀ is the wavenumber associated with the reference sound speed c₀, ${\zeta }_x^2=(1-M_x^2{)}^{-1}$, η = (c₀/c_z)² − 1 is the deviation of the refractive index from unity, c_z is the sound speed at height z, $\widehat{\mu }=\frac{1}{{\zeta }_x^2{k}_0^2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}$, and $\widehat{\tau }=M_x{\zeta }_x^2\sqrt{1+\eta }$. In the absence of flow, M_x = 0, ${\zeta }_x^2=1$ and $\widehat{\tau }=0$, thus the classical parabolic equation is retrieved [30]:

$$\left(\frac{\mathrm{\partial }}{\mathrm{\partial }x}-i{k}_0\sqrt{1+\eta +\frac{1}{k_0^2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}}\right)\widehat{p}=0.$$(9)

Using the Padé (1,1) approximation, and introducing the variable $\overline{\varphi }$ related to the velocity potential $\widehat{\varphi }$ by $\widehat{\varphi }(x,z)=\exp \left(i{k}_{0}x\right)\overline{\varphi }(x,z)$, equation (8) can be rewritten:

$${\mathrm{\Psi }}_1\left(x,z\right)\frac{\mathrm{\partial }\overline{\varphi }}{\mathrm{\partial }x}=i{k}_0{\mathrm{\Psi }}_2\left(x,z\right)\overline{\varphi },$$(10)where the operators Ψ₁ and Ψ₂ are given by:

$${\mathrm{\Psi }}_m=h_{m,0}+\frac{h_{m,2}}{k_0^2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2},\hspace{1em}m=\mathrm{1,2}.$$(11)

The coefficients h_m,j are written as:

$$h_{\mathrm{1,0}}=1+b_{\mathrm{1,1}}\eta ,\hspace{1em}{h}_{\mathrm{1,2}}=b_{\mathrm{1,1}}/{\zeta }_x^2,$$

$$h_{\mathrm{2,0}}=a_{\mathrm{1,1}}{\zeta }_x^2\eta -(1+b_{\mathrm{1,1}}\eta )\tilde{\tau },\hspace{1em}{h}_{\mathrm{2,2}}=a_{\mathrm{1,1}}-b_{\mathrm{1,1}}\tilde{\tau }/{\zeta }_x^2,$$with a_1,1 = 1/2, b_1,1 = 1/4, and the function $\tilde{\tau }$ is defined as:

$$\tilde{\tau }=M_x{\zeta }_x^2(\sqrt{1+\eta }-M_x)=\widehat{\tau }-M_x^2{\zeta }_x^2.$$

As in the classical wide-angle parabolic equation, the Crank-Nicholson algorithm can be used to reduce equation (10) to a matrix system that can be easily solved. The variable $\overline{\varphi }$ is discretized using a cartesian mesh of size Δx and Δz: ${\varphi }_m^n=\overline{\varphi }(m\Delta x,n\Delta z)$. The domain is bounded by a ground impedance condition at z = 0 and by an absorbing layer at the top of the domain. The details are given in the Appendix.

The acoustic pressure $\widehat{p}$ can be calculated from ${\varphi }_m^n$ at x_m = mΔx and z_n = nΔz using a second-order centered finite difference scheme (Eq. (84) of [26]):

$$\widehat{p}\left(x_m,z_n\right)=e^{i{k}_0{x}_m}\left[\left(1-M_x\right){\varphi }_m^n+\frac{i{M}_x}{2k_0\Delta x}\left[{\varphi }_{m+1}^n-{\varphi }_{m-1}^n\right]\right].$$(12)

The starting field used for the calculation of the two-dimensional sound pressure $\widehat{p}$ is defined as (see Appendix G in ref [30])

$$\overline{\varphi }\left(0,z_s\right)={\varphi }_0\left(0,z-h_S\right)+\frac{1-\beta }{1+\beta }{\varphi }_0\left(0,z+h_S\right)$$(13)where h_S is the source height, β is the normalized ground admittance and ϕ₀ is expressed as

$${\varphi }_0\left(0,z\right)=S\sqrt{i{k}_0}\left(1.3717-0.3701k_0^2{z}^2\right)e^{-\frac{k_0^2{z}^2}{3}},$$(14)with S as the pressure amplitude taken as unity. The two-dimensional sound pressure calculated is then converted into the three-dimensional pressure as $\left|{\widehat{p}}_{3\mathrm{D}}\right|=\left|\widehat{p}\right|/\sqrt{x}$. In free-field we know that $\left|{\widehat{p}}_{3\mathrm{D},\mathrm{FF}}\right|=\left|S\right|/R$ and thus the relative sound pressure level for a receiver at the position (x_R, z_R) is given by

$$\Delta L_{\mathrm{PE}}=10{\log }_{10}\left(\frac{|{\widehat{p}}_{3\mathrm{D}}{|}^2}{|{\widehat{p}}_{3\mathrm{D},\mathrm{FF}}{|}^2}\right)=10{\log }_{10}\left(\frac{\left|\widehat{p}\right(x_{\mathrm{R}},z_{\mathrm{R}}){|}^2{R}^2}{\left|S\right|x_{\mathrm{R}}}\right).$$(15)where the pressure amplitude S is taken as unity.

2.2.2 Scattering due to turbulence using the Harmonoise model

The Harmonoise engineering model is used only to consider the effect of the scattering due to turbulence [28]. In this model, the total relative sound pressure level is written as:

$$\Delta L=10{\log }_{10}\left(10^{\frac{\Delta L_{\mathrm{PE}}}{10}}+10^{\frac{\Delta L_{\mathrm{scat}}}{10}}\right),$$(16)where ΔL_PE is the relative sound pressure observed at the receiver due to the propagation effects obtained by the WAPE-MM in Section 2.2.1, and ΔL_scat is the sound level that reaches the receiver by scattering due to the turbulent fluctuations of wind speed and temperature in the atmosphere. This term is given by:

$$\Delta L_{\mathrm{scat}}=25+10{\log }_{10}\left({\gamma }_{\mathrm{T}}\right)+3{\log }_{10}\left(\frac{f}{1000}\right)+10{\log }_{10}\left(\frac{R}{100}\right),$$(17)where f is the frequency, R is the direct distance between source and the receiver and γ_T is the turbulence dependent parameter which is given by:

$${\gamma }_{\mathrm{T}}={\left(\frac{C_{\mathrm{T}}}{T_0}\right)}^2+\frac{22}{3}{\left(\frac{C_{\mathrm{v}}}{c_0}\right)}^2.$$(18)

Here, C_v and C_T are the turbulent structure parameters for wind speed and temperature fluctuations respectively and T₀ is the temperature in Kelvin. As shown by Ostashev and Keith Wilson (see Fig. 1 of Ref. [41]), the contribution from wind speed fluctuations is at least one order of magnitude higher than the contribution from temperature fluctuations, except in light wind conditions and close to the ground. As a result, the effect of turbulent temperature fluctuations is neglected in the following, and the structure parameter C_v is calculated from the Kolmogorov spectral model for consistency with the source model. Based on the expressions for structure function given in references [38, 42], we obtain:

$$C_v^2=C_K{ϵ}^{\frac{2}{3}}.$$(19)

Thus the term ΔL_scat that accounts for turbulence scattering in the Harmonoise model becomes:

$$\Delta L_{\mathrm{scat}}=25+10{\log }_{10}\left(\frac{22}{3}\frac{C_K{ϵ}^{2/3}}{c_0^2}\right)+3{\log }_{10}(f/1000)+10{\log }_{10}(R/100).$$(20)

2.3 Synthesis method based on the moving monopole model

For the calculation of ΔL, the source and the receiver positions are evidently important parameters that influence propagation effects. However, coupled with the source model of the wind turbine described in Section 2.1, the computation of the propagation effect for each segment of the blade at each height is expensive and time consuming. The moving monopole model allows the computation of ΔL for a finite number of source heights and approximates the height of each the noise source of the wind turbine model to the nearest-neighboring height as shown in Figure 1b. It has been shown in reference [24] that for a receiver distance greater than 200 m, the wind turbine noise predictions converge when seven source heights or more are considered in the moving monopole model. Following this suggestion, we compute ΔL for seven heights distributed evenly between the highest and lowest rotor heights to account for the propagation effects for each segment at each height. As the directivity of the airfoil is accounted for in the source model itself (Sect. 2.1), the synthesis of the individual point sources as monopoles is justified.

With SPL_FF computed using the source model and ΔL − α_absR using the propagation model, we obtain the frequency response for each moving segment at each discrete angular position. To synthesize a time signal from the PSD of acoustic pressure at the receiver, we use the method of Mascarenhas et al. [29]. The method consists of firstly synthesizing the contribution of each blade segment at every angular position by converting the frequency domain response between the frequencies f_min and f_max into a time signal with the help of the Inverse Discrete Fourier transform. Random phase is assigned to each frequency bin between f_min and f_max. To avoid the potential artificial perception of the sharp drop at f_max, the PSD is extrapolated beyond this frequency with the slope of f⁻². However this artifact is not perceived by a distant receiver (>200 m) for frequencies above 2 KHz when atmospheric absorption is included. The individual time signals of the segments at each angular position γ called grains are then arranged on the basis of their corresponding propagation time.

The duration of each considered grain corresponds to the time Δt_γ associated with the blade rotation over two successive angular positions γ and the propagation time difference Δt_r between two successive segment positions. The time duration of each grain is thus given by:

$$T_{\Delta \gamma }=\Delta t_{\gamma }+\Delta t_r,$$(21)where Δt_γ = Δγ/Ω depends on the rotational speed Ω of the wind turbine blade and the angular step Δγ, and Δt_r = Δr/c₀ depends on the difference Δr between the propagation distances corresponding to successive angular positions which can be positive or negative. The grain of duration T_Δγ contains N_Δγ samples depending on the sampling frequency.

For a smooth transition between two grains an adaptive window function is used, conserving the size of each grain while also maintaining the power of the respective grains [29]. The window function W[k] of N samples is defined with the overlapping functions f[χ] and g[χ] of length w_l with a unit response of length N_unit between them as can be seen in Figure 2. The overlapping functions used are given by $f\left[\chi \right]=\sin \left(\frac{\pi \chi }{2}\right)$ and $g\left[\chi \right]=\cos \left(\frac{\pi \chi }{2}\right)$, where χ ∈ [0, 1] is the normalized time index.

Figure 2 The window function W[k] of length N with the overlap functions f[χ] and g[χ]. The white area indicates the grain in consideration and the grey shaded area indicates the previous and next grains. The blue dashed lines represent the overlap functions of the adjacent grains.

As the segment is in motion, the duration of each grain is different due to Δt_r. Based on the required duration of the grain, the length of the unit response N_unit is correspondingly adjusted while maintaining the amount of the overlapping functions the same for all the grains. The amount of overlap to be used is determined on the basis of the smallest grain in the system. The overlap amount Ψ defines how much of the smallest grain is overlapped by the preceding and the following grain while the length of the unit response is set to zero. An overlap amount of Ψ = 100% which is suggested and used [29] meaning that the smallest grain would be overlapped completely by the preceding and following grains while still conserving its power and duration. The propagation time difference Δt_r is calculated with the reference sound speed c₀. In reality, the propagation time will be slightly different due to refraction effects.

Using a ray tracing model, we have compared in reference [43] (Appendix B) the propagation time of each grain Δt′_r using various wind speed profiles with the propagation time observed without wind Δt_r. It is found that the maximum difference Δt′_r − Δt_r that is observed is less than 3 ms for a receiver at x_R = 1000 m. It is thus justified to disregard the difference in the duration observed due to the propagation effects. As each of the blade segments in rotation are assumed to be uncorrelated, the grains are synthesized separately for TEN and TIN and then summed together at the corresponding time step. The synthesis of the total wind turbine noise is done with the inclusion of the propagation effects for a few test cases described in Section 4.

3 Propagation effects in various conditions

3.1 Configurations studied

In this section, we focus on the propagation effects for the same conditions which are to be tested in the wind turbine noise synthesis of Section 4. To this aim, the WAPE-MM method implemented in this study is calculated for a few source heights z_S between 35 m and 125 m that spans the entire rotor height of the wind turbine considered. The wind speed profile U_x(z) used follows a power law given by:

$$U_x\left(z\right)=U_{\mathrm{ref}}{\left(\frac{z}{z_{\mathrm{ref}}}\right)}^{\alpha },$$(22)where U_ref is the reference wind speed at the height z_ref. Since the plane of propagation is rotated by an angle θ with respect to the x-z plane (Oxz) as shown in Figure 1a, the effective wind speed profile used in the simulations is U(z) = U_x(z)cosθ. Different values of the shear exponent α are used between 0.2 (near-neutral atmosphere) and 0.5 (very stable atmosphere) [44].

The ground impedance is calculated via the Miki model [45] for a rigidly-backed layer of thickness e and flow resistivity σ. These ground impedance parameters can be estimated using the method suggested by Guillaume et al. [46]. Based on Table 2 of reference [46], a grass ground in summer is characterized by the mean parameters σ = 354 kNs/m⁴ and e = 0.0157 m, while the same grass ground in winter is characterized by the mean parameters σ = 631 kNs/m⁴ and e = 0.0060 m. To account for the turbulence in the source model as well as scattering effect, we consider three turbulence levels as shown in Table 1. The considered levels correspond to the range of the values found in the literature [38, 42].

All calculations of the WAPE-MM model are done with a grid resolution of Δx = Δz < λ/20, where λ is the wavelength of the propagated sound. The maximum height of the grid taken here as 320 m with the absorption layer beginning at 80% of the grid height. A second-order Salomons initial starter at x_S = 0 is used at the required source height z_S [30]. A constant temperature profile is set to the reference temperature of 10°C. The same reference temperature is used to calculate the absorption coefficient with the relative humidity of 80%. The number of frequencies equally distributed per third octave band used in the computation is given in Table 2.

3.2 Results

Figure 3 shows the relative sound pressure level ΔL at 500 Hz influenced by the wind flowing in the positive x direction for a receiver upwind (x < 0) and downwind (x > 0) As the WAPE-MM has an angular validity of 30° [26], the shaded area in Figure 3 shows where the calculation of ΔL is not valid. ΔL is plotted for different wind shear coefficients α, source heights and ground impedance parameter with the effect of turbulence scattering. It can be seen from Figure 3a that the wind shear α shifts the positioning of the produced interference patterns and the start of the shadow zone. The analytical solution matches well with the implemented WAPE-MM for α = 0. In Figure 3b the influence of the height of the source can be seen clearly in which the shadow zone for sources closer to the ground starts earlier. Finally in Figure 3c, the change in the ground parameters are more significant upwind rather than downwind. Here the influence of the scattering due to turbulence in the shadow zone can also be seen. High turbulence also influences the level of ΔL in the interference dips downwind.

Figure 3 Relative SPL ΔL calculated at f = 500 Hz for the receiver at height zR = 2 m and xR upwind (−x) and downwind (+x) with Uref = 8 m/s at zref = 80 m for (a) different wind shear exponents α with the source at zS = 80 m and a hard ground with no turbulence in the atmosphere, (b) α = 0.3 and with the grass ground in the summer and different source heights with no turbulence in the atmosphere, and (c) α = 0.3 with the source at zS = 80 m and different ground conditions and including turbulence in the atmosphere. The shaded area shows where the implemented WAPE-MM is not validated (solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).

Figures 4 and 5 show the spectrum of ΔL obtained between 100 Hz to 2000 Hz for the source at different heights and a receiver downwind and upwind at a distance of x_R = 500 m. The results for a grass ground in summer are compared to the results for a grass ground in winter. It can be seen from Figure 4 that the ground impedance parameters modify the level of ΔL and the positioning of the interference patterns. However, this change is more pronounced in the upwind condition as seen from Figure 5. The start of the shadow zone calculated for different heights can be seen in Figure 5 for a receiver at a distance of 500 m, in Figure 6 for a receiver at 800 m and in Figure 7 at a receiver distance of 1000 m. At 500 m, the shadow zone appears only for z_S = 35 m, while it is present for the two lowest source heights at 800 m, and for all source heights at 1000 m. Figures 6 and 7 show that the sound from sources that are higher, which corresponds to the top of the rotor plane, propagate over larger distances. The difference of ΔL(f) on the basis of the source height is what is experienced by the rotating blade of the wind turbine. It is clear from Figure 3b and Figure 4–7 that ΔL(f) cannot be modelled using a single source height and hence requires multiple heights to appropriately account for the propagation effects.

Figure 4 Narrowband spectrum of the relative SPL ΔL calculated for the receiver downwind at xR = 500 m with (a) a grass ground in the summer and (b) grass ground in the winter (α = 0.3 used with Uref = 8 m/s at zref = 80 m. Solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).

Figure 5 Narrowband spectrum of the relative SPL ΔL calculated for the receiver upwind at xR = 500 m with (a) a grass ground in the summer and (b) grass ground in the winter (α = 0.3 used with Uref = 8 m/s at zref = 80 m. Solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).

Figure 6 Narrowband spectrum of the relative SPL ΔL calculated for the receiver upwind at xR = 800 m with a grass ground in the summer (solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).

Figure 7 Narrowband spectrum of the relative SPL ΔL calculated for the receiver upwind at xR = 1000 m with a grass ground in the summer (solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).

4 Test cases

As in reference [16], the modelled wind turbine has a hub height of H = 80 m with a blade of length 45 m divided into 8 segments. The rotational speed Ω = 1.47 rad/s is set and one rotation is divided into N_γ = 36 discrete angular blade positions as this value of N_γ provides a good quality of sound synthesis [29]. The receiver is placed at a distance x_R and a constant height z_R = 2 m with the orientation θ with respect to the wind direction (Fig. 1a). The wind speed profile implemented follows the power-law profile as in equation (22) with the wind speed U_ref = 8 m/s at the reference height z_ref = 80 m and various wind shear exponents α. The SPL_FF is computed for the TEN and TIN mechanisms for the set of frequencies between 50 Hz and 5000 Hz that are given in Table 2.

As it is very time consuming to calculate ΔL_PE above 2 kHz and wind turbine noise levels are very low at frequencies beyond this value, ΔL_PE is calculated only from 50 Hz until 2300 Hz using 7 heights equally distributed between 35 m and 125 m as per the moving monopole model. The level of the relative SPL ΔL_PE for the third octave band at center frequency 2000 Hz is taken as the value of ΔL_PE for the frequencies between 2300 Hz and 5000 Hz in order to have a continuous decay of the noise level. We synthesize the signals of the wind turbine noise at the sampling frequency of f_s = 44.1 kHz between the frequencies 100 Hz and 5000 Hz for two complete rotations of the blades for various test cases which are described below. In the signals provided in the .mp3 format, 0 dBFS corresponds to 0.1 Pa. All the sound files associated to each individual case can be found in the respective table so that the reader can assess the differences as mentioned.

4.1 Free-field vs. Propagation effects

We begin by comparing the synthesis done in free field with the synthesis inclusive of the propagation effects. The signals are obtained for a receiver at distance x_R = 500 m downwind θ = 0°, close to crosswind θ = 80° and upwind θ = 180°. For the wind speed profile, the shear exponent is taken as α = 0.3. For this test case, a grass ground in summer [46] is implemented via the Miki-model. The turbulence considered in the TIN (Eq. (4)) and also in ΔL_scat for the scattering due to the atmospheric turbulence (Eq. (20)) corresponds to the medium level in Table 1. The different parameters for the cases A are given in Table 3. The spectrogram for the synthesized cases A1-1,2 and A2-1,2 are shown in Figure 8 and 9 and the difference due to the propagation effects can be seen clearly. For a quantitative comparison between the test cases, we plot the spectra of the SPL averaged over one rotation obtained from the frequency domain in Figure 10a. The amplitude modulation (AM) due to the rotation of the blades is the maximum difference in SPL(f, γ) during one rotation calculated as:

$$\mathrm{AM}\left(f\right)=\underset{\gamma }{\max }\left[\mathrm{SPL}\left(f,\gamma \right)\right]-\underset{\gamma }{\min }\left[\mathrm{SPL}\left(f,\gamma \right)\right].$$(23)

Figure 8 Spectrogram for the cases A in free field: (a) A1-1 and (b) A1-2.

Figure 9 Spectrogram for the cases A with the inclusion of the propagation effects: (a) A2-1 and (b) A2-2.

Figure 10 Narrowband spectra of the (a) time-averaged SPL and (b) AM of the test cases A which are given in Table 3.

From the time-averaged SPL in Figure 10a, it is evidently seen that noise received in free field does not contain the interference patterns due the ground reflection. The interference dips vary from 150 Hz to 250 Hz depending on the orientation of the receiver. The decrease in the SPL for cases A2 at higher frequencies starting from 1500 Hz is due to the atmospheric absorption. It can also be seen that the lower SPL in the cross-wind direction in comparison to upwind and downwind is maintained as expected [17] even with the addition of the propagation effects. The difference between the synthesized signals with and without propagation effects is clearly audible.

During the rotation of each blade, the SPL is modulated for each frequency. The maximum difference between the amplitude of the SPL in one rotation varies for each case as is plotted in Figure 10b. The AM which is absent in the free-field condition upwind is audibly present with the inclusion of the propagation effects. The level of modulation for the upwind and downwind condition is nearly the same and almost zero in free field conditions (Cases A1-1,3), which was expected as no propagation effect is considered and the propagation distance is almost the same for all grains. The AM for the case of A1-2 originates mainly from the orientation of the receiver which is close to the rotational plane of the blade. The AM is significantly changed for most frequencies with the addition of the propagation effects with the maximum change in the AM observed upwind (A2-3) for all frequencies. For the downwind (A2-1) and crosswind (A2-2) conditions, the increase in the AM is related to the ground effect that varies with the different source heights experienced by each segment during the movement of the blade as seen in Figure 4. The influence of the propagation effect is more significant in the upwind condition (A2-3) where there is an increase of around 4 dB in AM between 200 Hz and 300 Hz with a maximum increase of close to 8 dB around 1000–2000 Hz. This high increase in the AM is understood using Figure 5a in which it can be seen that due to the upward refraction caused by the wind speed profile, the lowest source height z_S = 35 m is already close to the shadow zone while the rest are within the limits of the shadow zone. Thus the blade segments close to the tip fluctuate from inside to outside the shadow zone during the course of one rotation. This inclusion of the propagation effects increases the level of the AM which is comparable to the crosswind condition. Close to the crosswind condition (θ = 80°, A2-2), the AM appears to slightly decrease at higher frequencies compared to that observed in free field.

4.2 Effect of the ground impedance variation between summer and winter

We synthesize three cases similar to A2 while replacing the grass ground in summer with the grass ground in winter and keeping all other parameters the same. The test cases are described in Table 4.

The comparative influence of the two types of grounds can be observed in Figure 11. From Figure 11a, it can be seen that the SPL difference between the two grounds for the downwind direction (θ = 0°, B1-1) and close to crosswind direction (θ = 80°, B1-2) is small. This difference is more visible between the two grounds for the upwind condition at the lower frequencies but it is difficult to audibly distinguish the difference. Similarly, the difference is small for the AM as seen in Figure 11b for cases B1-1 and B1-2. The AM is increased by around 3 dB for the downwind case B1-3 between the frequencies 250 Hz and 450 Hz. The difference observed in the propagation effects due to the two grounds can be clearly seen in Figure 4 and 5. However, the relatively small change in the SPL and AM seen in Figures 11a and 11b in the two grounds can be explained by the fact that the synthesized noise inclusive of propagation effect is a summation of all the segments passing through different source heights used for the calculation of the propagation effect. The outcome of this synthesis is an averaged out effect of the influence of the ground reflection. Another parameter which is known to influence the limits of the shadow zone is the considered wind profile. We again synthesize the noise with the parameters described in the test cases A2 while varying the wind shear exponent α for the next set of test cases.

Figure 11 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases B which are given in Table 4 that are compared to the test cases A2 given in Table 3.

4.3 Effect of the wind shear

A few test cases are synthesized similar to A2 but using different wind shear exponents while keeping all other parameters the same. The test cases are shown in Table 5.

From Figure 12, it can be seen that the influence of the wind shear on the observed SPL is small and the audible difference is subtle. The cases close to the crosswind (θ = 80°, C1,2-2) show nearly the same response for the SPL as well as the AM for all frequencies, which is expected as the effect of refraction due to the wind shear is low for this orientation. For the downwind cases C1-1 and C2-1, the relatively equal AM is explained by the fact that the maximum difference in the level of the refracted noise observed by each segment is almost the same for the different wind shear exponents and also there is no shadow zone for any source height observed. The effect of the different wind shear exponents in the upwind condition C1,2-3 is significant in the AM above 1 kHz (Fig. 12b). The high AM generated for a particular system upwind, as explained earlier, is the fluctuation of the received noise from inside to outside the limits of the shadow zone. This difference can also be heard clearly in the respective signals of this case. As it is known that the wind shear effectively modifies the limits of the shadow zone, the difference between the AM upwind for cases C1-3 and C2-3 is explained. The parameter that is seen to influence the level of the noise is the scattering of the atmospheric turbulence as seen in Figures 3c and 5. We investigate the influence of this parameter in the next set of test cases.

Figure 12 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases C which are given in Table 5.

4.4 Effect of the scattering due to turbulence

Here, we synthesize a few cases similar to the case A2-3 in the upwind direction, but with two other levels of turbulence while the other parameters remain the same. To study the influence of the scattering due to turbulence, as before we apply the turbulence to the source (TIN) as well as the scattering (ΔL_scat) in the cases D1 (Source and scattering) and compare it to the cases where we now apply the turbulence only to the source in D2 (Source only). The test cases are summarized in Table 6. The TDR (ϵ) is modified and corresponds to the low and high levels of turbulence in Table 1.

It can be seen directly from Figure 13 that the influence of the scattering due to the turbulence is almost negligible for the SPL as well as the AM. This inaudible difference difference is also observed in the signals provided. For a high turbulence it can be seen that SPL is higher due to TIN source model but provides almost the same AM as compared to the case with low turbulence. An explanation of this low scattering level can be attributed to the fact that not all the segments of the rotating blade are inside the computed shadow zone and thus experience low scattering due to turbulence, as can be seen in Figure 5. The noise from the segments that are outside the shadow zone dominates the generated noise and thus the influence of scattering due to turbulence is low. As observed until now, the limits of the shadow zone are crucial for the level of AM upwind. We finally investigate this influence by extending the distance x_R away from the wind turbine in the next set of test cases.

Figure 13 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases D which are given in Table 6.

Figure 14 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases E which are given in Table 7.

4.5 Shadow zone

We synthesize again two cases similar to case A2-3 with a moderate turbulence but at different distances x_R and the other parameters remaining the same as before. The test cases are given in Table 7.

The influence of the distance on the SPL is clearly seen in Figure 14a. The level of the SPL is seen to decrease as the distance from the wind turbine is increased, as is expected. However, the level of the AM as in Figure 14b evolves in a different way with respect to distance. The level of the AM for the receiver with the shortest distance (A2-3) is seen to be intermediate between the cases E1-1,2. This is because for the case A2-3, most of the sources above the height of 80 m are observed to be outside the shadow zone as can be seen in Figure 5. On the other hand it can be seen that for the case E1-1 (Fig. 6), the source at height 80 m is already in the shadow zone and so the moving blade produces a greater difference in the amplitude while moving in and out of the shadow zone, an outcome of which produces AM even in the higher frequencies. For the case E1-2, all of the high frequencies above 300 Hz are seen to remain in the shadow zone (Fig. 7), thus AM is produced only at lower frequencies. For these cases, the decrease in the level can be heard along with the corresponding change in the AM as mentioned above. It should be noted that the sharp transition of the amplitude heard in the audio signal is due of the number of source heights N_h = 7. For a larger number of source heights, the transition through the shadow zones would be smoother.

5 Conclusions

In this article, we have presented a physics-based synthesis model of wind turbine noise in an inhomogeneous moving medium where the influence of the ground impedance, atmospheric refraction and turbulence scattering has been included. The source part of the model based on Amiet’s theory is used to predict the noise from the trailing edge and the leading edge of the airfoil. The prediction of the trailing edge noise uses the wall pressure spectra model of Lee and Goody while the prediction of the turbulent inflow noise is based on the Kolmogorov spectrum. The propagation modelling uses the parabolic equation in moving medium in order to accurately capture the refraction effect due to the wind speed gradient. To consider the effect of scattering due to turbulence the Harmonoise model is implemented using the Kolmogorov spectrum. The synthesis of the signals are done based on the moving monopole model and using the method of overlapping grains with an adaptive window.

From the comparison of the individual test cases, the influence of the ground effect, wind shear coefficient and scattering due to the turbulence on the wind turbine noise have been studied. Even with the inclusion of the propagation effect, the characteristic large amplitude modulation of the wind turbine noise crosswind is observed throughout all the studied cases. The amplitude modulation produced in the crosswind direction originates from the Doppler effect and convective amplification caused by the rotating blades. Among the test cases studied, the effect of the implemented ground parameters and the wind shear seems to produce a significant difference mainly in the upwind direction. This is attributed to the influence of the parameters on the limits of the shadow zones that are observed due to the different source heights during the synthesis. The limits of the shadow zone as observed by the receiver upwind for each source height is seen to greatly influence the amount of the amplitude modulation. It is clear that the receiver distance influences the level of amplitude modulation in the upwind direction as it relates to the shadow zones observed from each source height. The effect of the scattering due to turbulence is seen to have a small effect on the synthesized noise. With the inclusion of prerecorded background noise at the correct level, the scattering due to turbulence may possibly be negligible.

This physics-based synthesis model can be easily used to obtain wind turbine noise signals for a receiver in a desired scenario. This makes it a useful tool for further studies that can avail of this model. One application of this tool is its use for psychoacoustic analysis of wind turbine noise. To have a better understanding of the perception of wind turbine noise and its psychoacoustic impact, it is beneficial to obtain signals with controlled parameters. This is achievable with the help of the presented model. With this model we have a straightforward way to synthesize the noise observed from one wind turbine. We can proceed in this same course to synthesize an entire wind farm with this tool and thus study the noise of a complete wind farm in the design phase itself. Along with obtaining the spectral levels of the expected noise from the wind farm, the tool also produces audible signals that can be further used for the desired purpose. Thus making it a useful tool that bridges the gap between citizens, local authorities and acousticians.

Conflict of interest

Author declared no conflict of interests.

Data availability statement

This article includes selected audio files embedded in the article. The sound files associated with this article are also available in https://zenodo.org/, under the reference doi: 10.5281/zenodo.7044407 [47].

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 812719.

Appendix

Numerical solution of the wide-angle parabolic equation in moving medium

Here, we explain how the parabolic equation presented in Section 2.2.1 can be numerically solved. First, the solution of equation (10) is advanced along x using the Crank-Nicholson (CN) scheme:

$$\left[{\mathrm{\Psi }}_1-\frac{i{k}_0\Delta x}{2}{\mathrm{\Psi }}_2\right]\overline{\varphi }\left(x+\Delta x\right)=\left[{\mathrm{\Psi }}_1+\frac{i{k}_0\Delta x}{2}{\mathrm{\Psi }}_2\right]\overline{\varphi }\left(x\right),$$(A1)where the terms Ψ₁ and Ψ₂ can be written:

$${\mathrm{\Psi }}_1=1+\frac{\eta }{4}+\frac{1}{4k_0^2{\zeta }_x^2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2},$$(A2)

$${\mathrm{\Psi }}_2=\frac{{\zeta }_x^2\eta }{2}-\left(1+\frac{\eta }{4}\right)\tilde{\tau }+\frac{\left(2{\zeta }_x^2-\tilde{\tau }\right)}{4k_0^2{\zeta }_x^2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}.$$(A3)

The domain is now discretized with mesh sizes Δx and Δz: ${\varphi }_m^n=\overline{\varphi }(m\Delta x,n\Delta z)$, with m = 1…M and n = 1…N. The second derivative with respect to z is estimated using a second order finite difference scheme:

$$\left(\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}\right){\varphi }_m^n=\frac{{\varphi }_m^{n+1}-2{\varphi }_m^n+{\varphi }_m^{n-1}}{k_0^2\Delta z^2}.$$(A4)

The numerical scheme associated with the CN for the WAPE method is thus:

$$M_1{\varphi }_{m+1}^n=M_2{\varphi }_m^n,$$(A5)where the matrices M₁ and M₂ are given by:

$$M_1{\varphi }_m^n=\left[1+\frac{{\eta }_m^n}{4}-\frac{i{k}_0\Delta x}{2}\left(\frac{({\zeta }_x^2{)}_m^n{\eta }_m^n}{2}-(1+\frac{{\eta }_m^n}{4}){\tilde{\tau }}_m^n\right)\right]{\varphi }_m^n+\left[\frac{2-i{k}_0\Delta x\left(2\right({\zeta }_x^2{)}_m^n-{\tilde{\tau }}_m^n)}{8k_0^2({\zeta }_x^2{)}_m^n}\right]\frac{{\varphi }_m^{n+1}-2{\varphi }_m^n+{\varphi }_m^{n-1}}{\Delta z^2},$$(A6)

$$M_2{\varphi }_m^n=[1+\frac{{\eta }_m^n}{4}-\frac{i{k}_0\Delta x}{2}\left(\frac{({\zeta }_x^2{)}_m^n{\eta }_m^n}{2}-(1+\frac{{\eta }_m^n}{4}){\tilde{\tau }}_m^n\right)-\frac{2-i{k}_0\Delta x\left(2\right({\zeta }_x^2{)}_m^n-{\tilde{\tau }}_m^n)}{4k_0^2({\zeta }_x^2{)}_m^n\Delta z^2}]{\varphi }_m^n +\left[\frac{2-i{k}_0\Delta x\left(2\right({\zeta }_x^2{)}_m^n-{\tilde{\tau }}_m^n)}{8k_0^2({\zeta }_x^2{)}_m^n\Delta z^2}\right]\left({\varphi }_m^{n+1}+{\varphi }_m^{n-1}\right).$$(A7)

The matrix M₁ in equation (A6) is tridiagonal with diagonal elements

$$b_n=\left[1+\frac{\eta }{4}-\frac{i{k}_0\Delta x}{2}\left(\frac{{\zeta }_x^2\eta }{2}-\left(1+\frac{\eta }{4}\right)\tilde{\tau }\right)-\frac{2-i{k}_0\Delta x\left(2{\zeta }_x^2-\tilde{\tau }\right)}{4k_0^2{\zeta }_x^2\Delta z^2}\right],$$(A8)and off-diagonal elements

$$a_n=c_n=\left[\frac{2-i{k}_0\Delta x\left(2{\zeta }_x^2-\tilde{\tau }\right)}{8k_0^2{\zeta }_x^2\Delta z^2}\right].$$(A9)

Similarly, the matrix M₂ in equation (A7) is tridiagonal with diagonal elements

$$e_n=\left[1+\frac{\eta }{4}+\frac{i{k}_0\Delta x}{2}\left(\frac{{\zeta }_x^2\eta }{2}-\left(1+\frac{\eta }{4}\right)\tilde{\tau }\right)-\frac{2+i{k}_0\Delta x\left(2{\zeta }_x^2-\tilde{\tau }\right)}{4k_0^2{\zeta }_x^2\Delta z^2}\right],$$(A10)and off-diagonal elements

$$d_n=f_n=\left[\frac{2+i{k}_0\Delta x\left(2{\zeta }_x^2-\tilde{\tau }\right)}{8k_0^2{\zeta }_x^2\Delta z^2}\right].$$(A11)

The boundary condition at z = 0 (n = 1) written with respect to the normalized admittance β = 1/Z can be obtained by using the centered second order scheme at the fictitious point ${\varphi }_m^0$ with z = −Δz:

$$\frac{{\varphi }_m^2-{\varphi }_m^0}{2\Delta z}+i{k}_0\beta {\varphi }_m^1=0.$$(A12)

The first lines of the matrices M₁ and M₂ are changed accordingly, with modified coefficients:

$$\begin{array}{c}{c}_{1g}=2c_1,\hspace{1em}{b}_{1_g}=b_1+2i{k}_0\Delta \mathrm{z\beta }{c}_1,\\ f_{1g}=2f_1,\hspace{1em}{e}_{1_g}=e_1+2i{k}_0\Delta \mathrm{z\beta }{f}_1.\end{array}$$

References

All Tables

All Figures

	Figure 1 Schematics for (a) the modelled wind turbine blades with the receiver represented as a sphere, and (b) the segmented blade in motion represented as a set of point sources whose heights are interpolated to the nearest neighbor in the propagation calculation based on a set of Nh source heights (Nh = 5 here).
In the text

	Figure 2 The window function W[k] of length N with the overlap functions f[χ] and g[χ]. The white area indicates the grain in consideration and the grey shaded area indicates the previous and next grains. The blue dashed lines represent the overlap functions of the adjacent grains.
In the text

Figure 3 Relative SPL ΔL calculated at f = 500 Hz for the receiver at height zR = 2 m and xR upwind (−x) and downwind (+x) with Uref = 8 m/s at zref = 80 m for (a) different wind shear exponents α with the source at zS = 80 m and a hard ground with no turbulence in the atmosphere, (b) α = 0.3 and with the grass ground in the summer and different source heights with no turbulence in the atmosphere, and (c) α = 0.3 with the source at zS = 80 m and different ground conditions and including turbulence in the atmosphere. The shaded area shows where the implemented WAPE-MM is not validated (solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).

In the text

	Figure 4 Narrowband spectrum of the relative SPL ΔL calculated for the receiver downwind at xR = 500 m with (a) a grass ground in the summer and (b) grass ground in the winter (α = 0.3 used with Uref = 8 m/s at zref = 80 m. Solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).
In the text

	Figure 5 Narrowband spectrum of the relative SPL ΔL calculated for the receiver upwind at xR = 500 m with (a) a grass ground in the summer and (b) grass ground in the winter (α = 0.3 used with Uref = 8 m/s at zref = 80 m. Solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).
In the text

	Figure 6 Narrowband spectrum of the relative SPL ΔL calculated for the receiver upwind at xR = 800 m with a grass ground in the summer (solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).
In the text

	Figure 7 Narrowband spectrum of the relative SPL ΔL calculated for the receiver upwind at xR = 1000 m with a grass ground in the summer (solid line with no turbulence, dashed line with low turbulence and dotted line with strong turbulence).
In the text

	Figure 8 Spectrogram for the cases A in free field: (a) A1-1 and (b) A1-2.
In the text

	Figure 9 Spectrogram for the cases A with the inclusion of the propagation effects: (a) A2-1 and (b) A2-2.
In the text

	Figure 10 Narrowband spectra of the (a) time-averaged SPL and (b) AM of the test cases A which are given in Table 3.
In the text

	Figure 11 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases B which are given in Table 4 that are compared to the test cases A2 given in Table 3.
In the text

	Figure 12 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases C which are given in Table 5.
In the text

	Figure 13 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases D which are given in Table 6.
In the text

	Figure 14 Narrowband spectra of the (a) Time-averaged SPL and (b) AM of the test cases E which are given in Table 7.
In the text