Development of an efficient numerical method for wind turbine flow, sound generation, and propagation under multi-wake conditions

: The propagation of aerodynamic noise from multi-wind turbines is studied. An efﬁcient hybrid method is developed to jointly predict the aerodynamic and aeroacoustics performances of wind turbines, such as blade loading, rotor power, rotor aerodynamic noise sources, and propagation of noise. This numerical method combined the simulations of wind turbine ﬂow, noise source and its propagation which is solved for long propagation path and under complex ﬂow environment. The results from computational ﬂuid dynamics (CFD) calculations not only provide wind turbine power and thrust information, but also provide detailed wake ﬂow. The wake ﬂow is computed with a 2D actuator disc (AD) method that is based on the axisymmetric ﬂow assumption. The relative inﬂow velocity and angle of attack (AOA) of each blade element form input data to the noise source model. The noise source is also the initial condition for the wave equation that solves long distance noise propagation in frequency domain. Simulations were conducted under different atmospheric conditions which showed that wake ﬂow is an important part that has to be included in wind turbine noise propagation. modeling of ﬂow, ANSource and ANPropagation, and the detailed phenomenon of ANPropagation through multi-wind turbine wakes are studied.


Introduction
The generation and propagation of wind turbine aerodynamic noise exhibit several special characteristics as compared to some other industrial noise problems. Considering the generation of wind turbine aerodynamic noise, the rotor aerodynamic noise level depends on airfoil profiles at cross sections, blade shape (the spanwise distribution of chord, twist, and airfoils), rotor size, rotational speed, yaw angle, tilt angle, pitch setting, the angular position of blades, wind speed profile, inflow turbulence level, density, and viscosity of air. From the perspective of aerodynamic noise propagation (ANPropagation), aerodynamic noise source (ANSource) of large wind turbines is often located at high altitude which leads to more significant problem for ANPropagation over long range. Intrinsically, the study of wind turbine noise is a multi-disciplinary subject which should consider wind farm aerodynamic noise generation (ANGeneration), ANPropagation, and energy production. However, ANGeneration from a single airfoil has a similar basis in the nature of aerodynamics and aeroacoustics as that from a large wind farm. Sophisticated numerical methods are needed to understand the physical mechanisms of flow induced noise. Computational aero-acoustics (CAA) methods were widely applied favored by the high-performance computing technology. However,

RANS/AD Method
To simulate the flow passes wind turbine rotor, the general equations are modified by adding the aerodynamic loadings of the rotor as external volume force terms in the momentum equations of the NS equations where i x is the displacement in the inertial coordinate system and i u is the velocity components.
t, ρ , and P represent time, air density and pressure respectively. ν is the kinematic viscosity term, ext f is the added source term. With the assumption based on the axisymmetric property of the HAWT, 2D-RANS/AD method is utilized such that the computational effort can be greatly decreased.
The rotor is assumed as permeable disc such that the mass equation remains unchanged. To avoid numerical instability the volume force in Equation (2) should be smeared along the flow direction using 1D Gaussian function. For the elements which are away from the disc with a normal distance of d, the smeared force f ' is computed by convolution computation

RANS/AD Method
To simulate the flow passes wind turbine rotor, the general equations are modified by adding the aerodynamic loadings of the rotor as external volume force terms in the momentum equations of the NS equations ∂u i ∂x i = 0 (i = 1 , 2 ) (1) where x i is the displacement in the inertial coordinate system and u i is the velocity components. t, ρ, and P represent time, air density and pressure respectively. ν is the kinematic viscosity term, f ext is the added source term. With the assumption based on the axisymmetric property of the HAWT, 2D-RANS/AD method is utilized such that the computational effort can be greatly decreased. The rotor is assumed as permeable disc such that the mass equation remains unchanged. To avoid numerical instability the volume force in Equation (2) should be smeared along the flow direction using 1D Gaussian function. For the elements which are away from the disc with a normal distance of d, the smeared force f ' is computed by convolution computation where the free parameter ε = 3∆z is three times larger than the axial size of reference grid ∆z. What is more, the body forces of rotor can be computed using blade element theory with the aerodynamic data of airfoil and the simulated information of flowfield. The present 2D-AD method solves the ideal axisymmetric flowfield. The axial and tangential induced velocities are naturally combined with the NS equations, such that the inflow angle and the relative inflow velocity are determined from The local AOA is easily computed by α = φ − γ, where γ is the real twist angle (added by pitch) at local blade element. The lift and drag forces per unit span are computed using Equation (7) with the aerodynamic data of airfoil.
The blade forces are updated after each iteration using newly computed velocity field, then the external volume force terms of the NS equations apply renewed blade forces in the next iteration.

BPM Model for Noise Source or Generation
The original BPM noise model was developed by Brooks et al. [14] which is able to predict airfoil self-noise. The model was based on scaling the wind tunnel experimental data using the NACA 0012 airfoil with chord lengths varies from 2.5 cm to 61 cm. Here, the BPM model was integrated with the BEM method [15] in order to deal with the complicated flow for wind turbine. The required boundary layer thickness data in the BPM model is obtained from the Xfoil computation so that this data varies with the airfoils rather than fixed with NACA 0012. A general expression of the BPM model is shown in Equation (8) SPL = 10log 10 δM n lD r 2 In the equation, it is shown that rotor aerodynamic noise is mainly a function of parameter δ which represents boundary layer thickness. The boundary layer thickness depends on airfoil shape, Reynolds number, AOA, turbulence level etc. l is the spanwise length of an airfoil segment. G 1 , G 2 , and G 3 are related to the Strouhal number St, Mach number M, Reynolds number Re and boundary layer thickness parameter δ. Detailed descriptions of noise mechanisms can be found in references [14,15]. ANSource level of rotor generally depends on the blade shape and operational conditions, such as rotational speed, wind speed, yaw, tilt, etc. Instead of combining the engineering BEM model, the BPM noise model is directly integrated into the NS equations with the AD model. The BPM model is a passive subroutine that can be activated anytime with the inflow velocity and AOA provided from AD simulations. The trailing edge boundary layer thickness is interpolated from the existing database prepared with Xfoil through a wide range of Reynolds numbers and AOAs. To take clean and rough blade surface conditions into account, the boundary layer thickness database is created under different surface roughness cases, such that wind turbine noise under free transition and fully turbulent flows can be modelled. A fully-turbulent boundary layer thickness database is utilized to calculate the noise of rotor whose blade surface is rough, and free transition database is applied for a rotor with a clean surface.

Noise Propagation Method
Sound speed in an ambient flow is written as c e f f = c + v x where c is the constant sound speed at a given air temperature, v x is the wind speed component through the propagation path. The wind speed components may come from field measurements, empirical expressions or numerical simulations. For an incompressible wind flow, the air density is assumed as a constant, the wave equation in a moving atmosphere reads [12] where k = ω/c 0 , ω is the radian frequency of sound wave, c 0 is the reference sound speed, P (r) denotes the monochromatic field of sound pressure and ε = (c 0 /c) 2 − 1. If the effect from ambient flow is neglected, this equation can be further reduced to Helmholtz equation. The ANSource of a wind turbine is placed on the left boundary of the computational domain depicted in Figure 2. A classical approach to represent a monopole ANSource is to distribute the source along the vertical direction which is also called starting field. In this study, the rotor noise source is located at the hub height as sketched at the left side of the boundary. Based on the initial source strength, the equation or the sound pressure amplitude is computed along the x-direction, such that solution is renewed from P(x) to P(x + ∆x). The boundary condition applied on the ground surface is calculated from the impedance of the ground. The empirical impedance model for absorbing materials developed by Delany and Bazely [21] is used in this study, where Z is the normalized acoustic impedance, σ represents the flow resistivity and f is the frequency. The similar approach can be applied on the top surface. A constant of the normalized acoustic impedance Z = 1 is implemented on the top grid, such that the vertical waves through the top surface is vanished without any reflection. However, to sufficiently eliminate wave reflections, an absorbing layer is defined on top of the grid. As shown in the figure, the grey region sketched on top of the domain is the absorbing layer. In the study, the width of such an absorbing layer is defined as 50λ. This implies that for lower frequency sound propagation, a larger absorbing layer is needed. The sound pressure along the vertical line is updated from step n to n + 1 until the receiver position is reached. The solution procedure starts with the separation of the forward and backward propagation waves.
Here, only the forward wave propagates from source towards receiver is considered (in the positive x-axis), which is in the downwind direction. Detailed mathematical manipulations of the equation are found in [12,13]. According to the methodologies described in Sections 2.1 and 2.2, the noise spectrum of a given wind turbine is the logarithmic sum of sound pressure level (SPL) from all the blade elements such that (11) where n is the number of the blade elements and SPL i total is the SPL of the ith blade element, which accounts for all the noise mechanisms. The SPL of wind turbine ANSource L p ( f ) is the initial condition for calculating long-path ANPropagation. The sound power level L w ( f ), which is a representation of wind turbine total noise strength, it does not change with propagation range. The SPL at the receiver is where the first term on the right side is the sound power level, the second term −10 log 10 4πD 2 is the geometric attenuation term that corresponds to the sound intensity loss when the sound wave sphere is expanding. The third term shows the attenuation due to air absorption, which is proportional to the absorption coefficient α and the propagation length D. ∆L is the relative SPL differences caused by factors such as ground reflection/absorption (flat/irregular terrain), atmospheric refraction, wind and turbulence and sound barriers, etc. It does not include the effects of geometric attenuation and air absorption. The ∆L varies along the propagation path with these factors, which is beneficial to clearly manifest the effects of propagation. Most of the discussions and results in Section 3 are dealing with this relative pressure level which is abbreviated as RSPL. The values of α depend on sound frequency, air temperature, pressure and humidity. In the following studies, unless specified separately, a relative humidity of 70% and temperature of 10 • C are assumed.
where the first term on the right side is the sound power level, the second term is the geometric attenuation term that corresponds to the sound intensity loss when the sound wave sphere is expanding. The third term shows the attenuation due to air absorption, which is proportional to the absorption coefficient α and the propagation length D. L Δ is the relative SPL differences caused by factors such as ground reflection/absorption (flat/irregular terrain), atmospheric refraction, wind and turbulence and sound barriers, etc. It does not include the effects of geometric attenuation and air absorption. The L Δ varies along the propagation path with these factors, which is beneficial to clearly manifest the effects of propagation. Most of the discussions and results in Section 3 are dealing with this relative pressure level which is abbreviated as RSPL. The values of α depend on sound frequency, air temperature, pressure and humidity. In the following studies, unless specified separately, a relative humidity of 70% and temperature of 10 °C are assumed. The wave equation is solved based on each frequency. The ANSource and wake information of given wind turbines from BPM/AD computations are combined into the ANPropagation solver. For each frequency, the corresponding ANSource, known as the starting field, is individually applied as the initial condition to the computational domain. Detailed descriptions about starting field formulation can be found in publication from Salomons [11]. The horizontal size of the simulation range is determined by the location of the receiver. However, the size of the vertical domain is flexible such that the height of domain can vary with the simulated frequency. The maximum vertical size of the grid should be less than / 8 x z λ Δ = Δ = so that a sound wave can be resolved with at least eight grid points. As a result, the domain height is inversely proportional to frequencies for a fixed number of vertical grids. Besides, a larger CFD domain should be used in the simulation such that the flow field can fully cover the noise propagation domain.

Results and Discussions
A few validations about ANPropagation cases are shown first. Under various flow conditions, the flow passing wind turbine rotor, ANGeneration and ANPropagation problems are discussed in detail. The RSPL L Δ over a long distance is computed for each frequency. After the validations, the wind turbine ANPropagation across multi-wakes is shown. In all the cases, the temperature at the ground level is 15 °C, such that the equivalent sound speed is 340 m/s. The wave equation is solved based on each frequency. The ANSource and wake information of given wind turbines from BPM/AD computations are combined into the ANPropagation solver. For each frequency, the corresponding ANSource, known as the starting field, is individually applied as the initial condition to the computational domain. Detailed descriptions about starting field formulation can be found in publication from Salomons [11]. The horizontal size of the simulation range is determined by the location of the receiver. However, the size of the vertical domain is flexible such that the height of domain can vary with the simulated frequency. The maximum vertical size of the grid should be less than ∆x = ∆z = λ/8 so that a sound wave can be resolved with at least eight grid points. As a result, the domain height is inversely proportional to frequencies for a fixed number of vertical grids. Besides, a larger CFD domain should be used in the simulation such that the flow field can fully cover the noise propagation domain.

Results and Discussions
A few validations about ANPropagation cases are shown first. Under various flow conditions, the flow passing wind turbine rotor, ANGeneration and ANPropagation problems are discussed in detail. The RSPL ∆L over a long distance is computed for each frequency. After the validations, the wind turbine ANPropagation across multi-wakes is shown. In all the cases, the temperature at the ground level is 15 • C, such that the equivalent sound speed is 340 m/s.

Validation of Long-Range Sound Propagations without Ambient Flow
This validation case considers the idealized long-range sound propagation without any moving medium or temperature gradient effects. The terrain contains a flat and acoustically soft surface with a total length of 200 m. The flow resistivity of a grassland are usually in range from 100 kPa · s · m −2 to 300 kPa · s · m −2 . The absorbing ground surfaces in the current simulation has a flow resistivity of 250 kPa · s · m −2 . The ANSource locates at x = 0 m, z = 0.5 m and the receiver at x = 200 m, z = 1.5 m. The current results are compared with the data presented in reference [22] where either measured data or data obtained from other models are available. The sound transmission losses over 200 m are shown in Figure 3 where the selected frequencies in the 1/3-octave band are gathered. At very low frequencies, i.e., f = 20 Hz and f = 40 Hz, the RSPL ∆L is increased with the propagation distance. At higher frequencies, ∆L gradually decreases along the propagation distance. The ∆L values at 200 m are compared in Figure 4 which shows that there is a large propagation loss between 400 Hz and 1000 Hz. The present simulation agrees with the reference data and the simulation obtained from the Nord2000 commercial model [23]. medium or temperature gradient effects. The terrain contains a flat and acoustically soft surface with a total length of 200 m. The flow resistivity of a grassland are usually in range from 100 2 kPa s m − ⋅ ⋅ to 300 2 kPa s m − ⋅ ⋅ . The absorbing ground surfaces in the current simulation has a flow resistivity of 250 2 kPa s m − ⋅ ⋅ . The ANSource locates at x = 0 m, z = 0.5 m and the receiver at x = 200 m, z = 1.5 m. The current results are compared with the data presented in reference [22] where either measured data or data obtained from other models are available. The sound transmission losses over 200 m are shown in Figure 3 where the selected frequencies in the 1/3-octave band are gathered. At very low frequencies, i.e., f = 20 Hz and f = 40 Hz, the RSPL L Δ is increased with the propagation distance. At higher frequencies, L Δ gradually decreases along the propagation distance. The L Δ values at 200 m are compared in Figure 4 which shows that there is a large propagation loss between 400 Hz and 1000 Hz. The present simulation agrees with the reference data and the simulation obtained from the Nord2000 commercial model [23].   The current results are compared with the data presented in reference [22] where either measured data or data obtained from other models are available. The sound transmission losses over 200 m are shown in Figure 3 where the selected frequencies in the 1/3-octave band are gathered. At very low frequencies, i.e., f = 20 Hz and f = 40 Hz, the RSPL L Δ is increased with the propagation distance. At higher frequencies, L Δ gradually decreases along the propagation distance. The L Δ values at 200 m are compared in Figure 4 which shows that there is a large propagation loss between 400 Hz and 1000 Hz. The present simulation agrees with the reference data and the simulation obtained from the Nord2000 commercial model [23].

Validation of Long-Range Sound Propagations with a Positive Temperature Gradient
In this case, the terrain distance, source, and receiver locations remain the same, but the ground is an acoustically hard surface, meaning that the flow resistivity is close to infinity. The temperature gradient was measured as 0.0846 K/m. If a linear temperature variation is assumed, the temperature and sound speed profiles are easily computed using the following relations Appl. Sci. 2019, 9, 100 8 of 22 Such a temperature gradient leads to a downward sound propagation. As a result, the RSPL along the propagation direction will probably increase. In Figure 5, it is seen that there is a general increase of ∆L along the propagation path except near the source locations (x < 40 m). The ∆L values are gathered at x = 200 m distance and are depicted in Figure 6. It is found that the ∆L values are positive between 20 Hz and 5 kHz. The current prediction fits well with the reference data. gradient was measured as 0.0846 K/m. If a linear temperature variation is assumed, the temperature and sound speed profiles are easily computed using the following relations Such a temperature gradient leads to a downward sound propagation. As a result, the RSPL along the propagation direction will probably increase. In Figure 5, it is seen that there is a general increase of L Δ along the propagation path except near the source locations (x < 40 m). The L Δ values are gathered at x = 200 m distance and are depicted in Figure 6. It is found that the L Δ values are positive between 20 Hz and 5 kHz. The current prediction fits well with the reference data.
Such a temperature gradient leads to a downward sound propagation. As a result, the RSPL along the propagation direction will probably increase. In Figure 5, it is seen that there is a general increase of L Δ along the propagation path except near the source locations (x < 40 m). The L Δ values are gathered at x = 200 m distance and are depicted in Figure 6. It is found that the L Δ values are positive between 20 Hz and 5 kHz. The current prediction fits well with the reference data.

Validation of Long-Range Sound Propagations with Upward Refraction under Wind Condition
Similar as shown in the previous test case, the ground acoustic impedance is also infinity. However, an upward refraction case is considered which is caused by a moving atmosphere. The wind speed recorded at 10 m height above the ground is −4.165 m/s, and the surface roughness length is 0.1 m. In a moving medium, the sound speed can be computed using the equation below where c e f f is the profile of effective sound speed along the vertical direction, c 0 = 340 m/s is the sound speed and z 0 = 0.1 m is the roughness height. Knowing the wind speed at 10 m height and the roughness length, it is easy to show that b ≈ −1. Such a sound speed profile represents a slightly upward atmospheric refraction. In Figure 7, except at very low frequencies, the RSPL decays very fast along the propagation path. Figure 8 shows a very good agreement with the reference data which indicates the capability of the current method to handle noise propagations through moving media.
is 0.1 m. In a moving medium, the sound speed can be computed using the equation below where eff c is the profile of effective sound speed along the vertical direction, 0 c = 340 m/s is the sound speed and 0 z = 0.1 m is the roughness height. Knowing the wind speed at 10 m height and the roughness length, it is easy to show that 1 b ≈ − . Such a sound speed profile represents a slightly upward atmospheric refraction. In Figure 7, except at very low frequencies, the RSPL decays very fast along the propagation path. Figure 8 shows a very good agreement with the reference data which indicates the capability of the current method to handle noise propagations through moving media.
where eff c is the profile of effective sound speed along the vertical direction, 0 c = 340 m/s is the sound speed and 0 z = 0.1 m is the roughness height. Knowing the wind speed at 10 m height and the roughness length, it is easy to show that 1 b ≈ − . Such a sound speed profile represents a slightly upward atmospheric refraction. In Figure 7, except at very low frequencies, the RSPL decays very fast along the propagation path. Figure 8 shows a very good agreement with the reference data which indicates the capability of the current method to handle noise propagations through moving media.

Wake Modeling of a Wind Turbine Cluster
Wind turbine flow are simulated for a Vestas wind turbine called NM80. The AOA and the relative velocity at all blade elements from RANS/AD simulations are inputted to the BPM model. The AD solver is much more efficient because it represents the rotor by volume force term rather than the dense body-fitted-mesh around blade surface. The flow computations are carried out with the in-house developed flow solver [24][25][26] which is an incompressible flow solver.
The grid for AD simulation and a typical flowfield with multi-wake structures are depicted in Figure 9. In Figure 9a, two wind turbines, depicted with thick white lines, are used for demonstration. The inflow boundary is applied on the left and top domain edge (numbered as 1 in Figure 9), the axisymmetric boundary on bottom edge (numbered as 2) and outflow boundary on the right edge (numbered as 3). The computational domain consists of the several blocks whose number is the same as that of wind turbines to be investigated. Each block contains grid points of 64 × 64. Using this approach, multi-wake flowfield can be efficiently modelled. A typical horizontal velocity contour plot is shown in Figure 9b.
The AD solver is much more efficient because it represents the rotor by volume force term rather than the dense body-fitted-mesh around blade surface. The flow computations are carried out with the inhouse developed flow solver [24][25][26] which is an incompressible flow solver.
The grid for AD simulation and a typical flowfield with multi-wake structures are depicted in Figure 9. In Figure 9a, two wind turbines, depicted with thick white lines, are used for demonstration. The inflow boundary is applied on the left and top domain edge (numbered as 1 in Figure 9), the axisymmetric boundary on bottom edge (numbered as 2) and outflow boundary on the right edge (numbered as 3). The computational domain consists of the several blocks whose number is the same as that of wind turbines to be investigated. Each block contains grid points of 64 × 64. Using this approach, multi-wake flowfield can be efficiently modelled. A typical horizontal velocity contour plot is shown in Figure 9b. To show the model accuracy, in Figure 10, some results from the current simulations using RANS turbulence model are compared with the LES results under same flow conditions [27,28]. On the modeling of wind turbine wake, many studies can be found [29][30][31]. In the present study, it is not of our primary interest to go into more details on this topic. Some more details of the AD modeling approach and cross comparisons can be found from our previous works [32][33][34]. To show the model accuracy, in Figure 10, some results from the current simulations using RANS turbulence model are compared with the LES results under same flow conditions [27,28]. On the modeling of wind turbine wake, many studies can be found [29][30][31]. In the present study, it is not of our primary interest to go into more details on this topic. Some more details of the AD modeling approach and cross comparisons can be found from our previous works [32][33][34].

Modeling of the Wind Turbine ANSource
As is shown in the above AD simulation, every wind turbine has different inflow conditions because some wind turbines are in the wake of another, which leads to different ANSources. The simulation were conducted for a free-stream wind velocity of 8 m/s. For a receiver at 120 m far away, the detailed noise spectra solely from the first upstream wind turbine are shown in Figure 11. The noise spectra for other downstream wind turbines will be shown in the later sections. As shown in the figure, the total SPL spectrum is name as SPL-TOTAL which is logarithmically summed from several noise mechanisms. These noise mechanisms include: noise from blade tip (SPL-TIP); from

Modeling of the Wind Turbine ANSource
As is shown in the above AD simulation, every wind turbine has different inflow conditions because some wind turbines are in the wake of another, which leads to different ANSources.
The simulation were conducted for a free-stream wind velocity of 8 m/s. For a receiver at 120 m far away, the detailed noise spectra solely from the first upstream wind turbine are shown in Figure 11. The noise spectra for other downstream wind turbines will be shown in the later sections. As shown in the figure, the total SPL spectrum is name as SPL-TOTAL which is logarithmically summed from several noise mechanisms. These noise mechanisms include: noise from blade tip (SPL-TIP); from trailing edge bluntness (SPL-TEBLUNT); from the shedding of laminar boundary layer vortex (SPL-LBLVS); from flow separation (SPL-SEPARATION); from turbulent boundary layer trailing edge (SPL-TBLTE); from turbulent inflow (SPL-INFLOW) [14][15][16]. Later, the total noise spectrum will be applied into the noise propagation computations. Figure 10. Comparisons of the normalized wake velocity from LES and RANS simulations, for downstream locations of 5D and 7D after a single wind turbine.

Modeling of the Wind Turbine ANSource
As is shown in the above AD simulation, every wind turbine has different inflow conditions because some wind turbines are in the wake of another, which leads to different ANSources. The simulation were conducted for a free-stream wind velocity of 8 m/s. For a receiver at 120 m far away, the detailed noise spectra solely from the first upstream wind turbine are shown in Figure 11. The noise spectra for other downstream wind turbines will be shown in the later sections. As shown in the figure, the total SPL spectrum is name as SPL-TOTAL which is logarithmically summed from several noise mechanisms. These noise mechanisms include: noise from blade tip (SPL-TIP); from trailing edge bluntness (SPL-TEBLUNT); from the shedding of laminar boundary layer vortex (SPL-LBLVS); from flow separation (SPL-SEPARATION); from turbulent boundary layer trailing edge (SPL-TBLTE); from turbulent inflow (SPL-INFLOW) [14][15][16]. Later, the total noise spectrum will be applied into the noise propagation computations. Figure 11. Aerodynamic noise source spectra of a single wind turbine.

Sound Profiles Inside wind Turbine Wakes
Although the flow and the noise source are modelled without wind and temperature gradients, in reality the effective sound speed profile will be affected by wind and temperature profiles. It was shown that the combinations of wind speed and temperature gradient forms different types of sound speed profiles [35]. Three idealized initial sound profiles are considered in the following sound propagation simulations. In Figure 12a, a uniform sound speed at temperature of 15 °C is assumed Figure 11. Aerodynamic noise source spectra of a single wind turbine.

Sound Profiles Inside wind Turbine Wakes
Although the flow and the noise source are modelled without wind and temperature gradients, in reality the effective sound speed profile will be affected by wind and temperature profiles. It was shown that the combinations of wind speed and temperature gradient forms different types of sound speed profiles [35]. Three idealized initial sound profiles are considered in the following sound propagation simulations. In Figure 12a, a uniform sound speed at temperature of 15 • C is assumed at the left boundary. The wind turbines are separated with a normalized distance of seven rotor diameters. For example, the second wind turbine is located at x = 560 m and the third wind turbine is located at x = 1120 m. Inside the computational domain, the sound speed is superimposed with the wind turbine wake field. Similarly, Figure 12b,c represent downward and upward sound speed contours. These sound profiles merge with the wind turbine multi-wakes and show distinct characteristics. The sound profiles have significant effects on wind turbine noise propagation as shown in the next subsection.
diameters. For example, the second wind turbine is located at x = 560 m and the third wind turbine is located at x = 1120 m. Inside the computational domain, the sound speed is superimposed with the wind turbine wake field. Similarly, Figure 12b,c represent downward and upward sound speed contours. These sound profiles merge with the wind turbine multi-wakes and show distinct characteristics. The sound profiles have significant effects on wind turbine noise propagation as shown in the next subsection.

Wind Turbine ANPropagation across Multi-Wake
As shown in Figure 13, the RSPL or ∆L is simulated along 1500 m propagation distance, at f = 300 Hz and with the noise spectrum in Figure 11. The first wind turbine locates at x = 0 m, the second at 560 m and the third at 1120 m. The distance between each wind turbine are relative large. The wind turbine ANSource are placed at a height of 70 m. The ANSources from the second and third wind turbines, depicted as white lines, are not considered so as to focus on the propagation phenomenon through the wakes. In Figure 13, RSPL results are shown under several conditions: (1) Case 1, ANPropagation under homogeneous atmospheric condition, no wake effect; (2) Case 2, ANPropagation under homogeneous atmospheric condition, considering multi-wakes; (3) Case 3,

Wind Turbine ANPropagation across Multi-Wake
As shown in Figure 13, the RSPL or ∆L is simulated along 1500 m propagation distance, at f = 300 Hz and with the noise spectrum in Figure 11. The first wind turbine locates at x = 0 m, the second at 560 m and the third at 1120 m. The distance between each wind turbine are relative large. The wind turbine ANSource are placed at a height of 70 m. The ANSources from the second and third wind turbines, depicted as white lines, are not considered so as to focus on the propagation phenomenon through the wakes. In Figure 13, RSPL results are shown under several conditions: (1) Case 1, ANPropagation under homogeneous atmospheric condition, no wake effect; (2) Case 2, ANPropagation under homogeneous atmospheric condition, considering multi-wakes; (3) Case 3, ANPropagation under atmospheric condition of down-ward refraction, considering multi-wakes; (4) Case 4, ANPropagation under atmospheric condition of upward refraction, considering multi-wakes. It is shown that the wind turbine ANPropagation is clearly affected by the complicated flow environment. For the case without any wake effect shown in Figure 13a, the ANPropagation is mainly influenced by the reflection and absorption of ground. As is shown in Figure 13b, the contours of RSPL is more complicated than that of Figure 13a which distinctly shows the influences of wake. For atmospheric condition of downward refracting depicted in Figure 13c, the 'noise pollution' band also shows a downward reflecting behavior. Additionally, a low-noise zone, shown in blue, can be found near x = 1000 m. In Figure 13d, a much quieter situation is observed under an upward refraction condition. In such a case, the quiet zone is located between wind turbine #2 and #3. It should be emphasized that uniform inflow boundary condition is applied in the flowfield simulations of Figure 13c,d. However, the un-uniform or down/up-ward atmospheric conditions, which means atmospheric altitude-dependent variation of the sound speed, are used in the ANPropagation computations.
ANPropagation under atmospheric condition of down-ward refraction, considering multi-wakes; (4) Case 4, ANPropagation under atmospheric condition of upward refraction, considering multi-wakes. It is shown that the wind turbine ANPropagation is clearly affected by the complicated flow environment. For the case without any wake effect shown in Figure 13a, the ANPropagation is mainly influenced by the reflection and absorption of ground. As is shown in Figure 13b, the contours of RSPL is more complicated than that of Figure 13a which distinctly shows the influences of wake. For atmospheric condition of downward refracting depicted in Figure 13c, the 'noise pollution' band also shows a downward reflecting behavior. Additionally, a low-noise zone, shown in blue, can be found near x = 1000 m. In Figure 13d, a much quieter situation is observed under an upward refraction condition. In such a case, the quiet zone is located between wind turbine #2 and #3. It should be emphasized that uniform inflow boundary condition is applied in the flowfield simulations of Figure  13c,d. However, the un-uniform or down/up-ward atmospheric conditions, which means atmospheric altitude-dependent variation of the sound speed, are used in the ANPropagation computations. Manifesting the noise received by humans, the RSPL or ∆L at h = 2 m is extracted from Figure 13 and shown in Figure 14.The sound attenuation ∆L of Case 1 starts to vary linearly with the propagation distance at around x = 400 m, because the reflection effect becomes weaker after that. In contrast, the ∆L of Case 2, 3, and 4 vary obviously even after a long distance. In Cases 2 and 3, ∆L fluctuates over distance due to the wake effect and the strong downward refraction respectively. Whereas in Case 4, a more significant RSPL loss is shown which is caused by the upward refraction. Manifesting the noise received by humans, the RSPL or ∆L at h = 2 m is extracted from Figure 13 and shown in Figure 14.The sound attenuation ∆L of Case 1 starts to vary linearly with the propagation distance at around x = 400 m, because the reflection effect becomes weaker after that. In contrast, the ∆L of Case 2, 3, and 4 vary obviously even after a long distance. In Cases 2 and 3, ∆L fluctuates over distance due to the wake effect and the strong downward refraction respectively. Whereas in Case 4, a more significant RSPL loss is shown which is caused by the upward refraction. It should be mentioned that the ∆L values might be quite different when the observer height changes. From this comparison, it can be concluded that ANPropagation can be greatly influenced by the environmental conditions. What is more, the change of SPL in the near downstream is much smaller than that at far downstream (x > 400 m), which shows the significance of long distance wind turbine noise predictions under complex atmospheric conditions. The near field (x < 400 m) simulations and measurements only belong to ANSource study which is not identical to the long-path ANPropagation. Manifesting the noise received by humans, the RSPL or ∆L at h = 2 m is extracted from Figure 13 and shown in Figure 14.The sound attenuation ∆L of Case 1 starts to vary linearly with the propagation distance at around x = 400 m, because the reflection effect becomes weaker after that. In contrast, the ∆L of Case 2, 3, and 4 vary obviously even after a long distance. In Cases 2 and 3, ∆L fluctuates over distance due to the wake effect and the strong downward refraction respectively. Whereas in Case 4, a more significant RSPL loss is shown which is caused by the upward refraction. It should be mentioned that the ∆L values might be quite different when the observer height changes. From this comparison, it can be concluded that ANPropagation can be greatly influenced by the environmental conditions. What is more, the change of SPL in the near downstream is much smaller than that at far downstream (x > 400 m), which shows the significance of long distance wind turbine noise predictions under complex atmospheric conditions. The near field (x < 400 m) simulations and measurements only belong to ANSource study which is not identical to the long-path ANPropagation.  Figure 14 represents the ∆L solely due to propagation which is computed at a fixed frequency of f = 300 Hz. Repeating the calculation for many frequencies and extracting the data at a receiver 1500 m away and 2 m height, a full spectrum of ∆L will be obtained. Using Equation (12) and considering all the effects (including ∆L), the SPL at receiver's location under Case 1 and Case 3 is gathered and shown in Figure 15. Besides, the spectrum of ANSource and at receiver's location (without ∆L) are also depicted in Figure 15. In this work, the ANSource is considered in an 1/3-octave band from 20 Hz to 5000 Hz. Since a minimal eight grid points are needed to resolve a sound wave, the grid number for high-frequency simulations increases sharply if the same domain size is adopted.  Figure 14 represents the ∆L solely due to propagation which is computed at a fixed frequency of f = 300 Hz. Repeating the calculation for many frequencies and extracting the data at a receiver 1500 m away and 2 m height, a full spectrum of ∆L will be obtained. Using Equation (12) and considering all the effects (including ∆L), the SPL at receiver's location under Case 1 and Case 3 is gathered and shown in Figure 15. Besides, the spectrum of ANSource and at receiver's location (without ∆L) are also depicted in Figure 15. In this work, the ANSource is considered in an 1/3-octave band from 20 Hz to 5000 Hz. Since a minimal eight grid points are needed to resolve a sound wave, the grid number for high-frequency simulations increases sharply if the same domain size is adopted. For instance, there is a 100-fold increase in the grid number for frequency of 3000 Hz compared with that of 300 Hz. To reduce the computational effort, frequencies above 5000 Hz are not considered here. Another reason is that the noise reduction due to air absorption, −αD in Equation (12), is much larger than the other factors at large frequency (f > 5000 Hz). The SPL at the source location is shown with black line with asterisk and labeled as 'Source: Lp' in Figure 15. The black lines with stars represent the SPL at receiver (2 m height at x = 1500 m) without considering the propagation loss ∆L in Equation (12). The blue line (Case 1) and red line (Case 3) consider all the factors, including geometric attenuation, air absorption and ∆L, in the long-path ANPropagation simulations. The big deviation between the two black lines demonstrates that the geometrical attenuation (10 log 10 4πD 2 ) contributes most to the ANPropagation loss. The even larger reduction of noise, when f > 1000 Hz, means that air absorption has a larger effect at high frequencies. The differences between blue/red lines and the black line with stars implies the value of ∆L. The value of blue and red lines are larger than that of black line with stars, which means positive ∆L near low frequency (f = 30 Hz). However, Case 1 and Case 3 show opposite ∆L around f = 200 Hz which implies the influences of wake and downward refraction. geometric attenuation, air absorption and ∆L, in the long-path ANPropagation simulations. The big deviation between the two black lines demonstrates that the geometrical attenuation ( contributes most to the ANPropagation loss. The even larger reduction of noise, when f > 1000 Hz, means that air absorption has a larger effect at high frequencies. The differences between blue/red lines and the black line with stars implies the value of ∆L. The value of blue and red lines are larger than that of black line with stars, which means positive ∆L near low frequency (f = 30 Hz). However, Case 1 and Case 3 show opposite ∆L around f = 200 Hz which implies the influences of wake and downward refraction. Similar as the previous simulations, the following study focuses on the noise propagation area between 1120 m and 2500 m, and a higher frequency f = 1000 Hz is selected. In this region, the three individual wind turbine noise sources propagate and superimpose on each other. In Figure 16a, the results cover the entire domain behind the last wind turbine. The solution differs from Figure 13a in the sense that it is the RSPL from all the three wind turbines. In Figure 16b, the uniform sound profile at the inlet boundary is affected by the multi-wake (see Figure 12a), such that the resulted RSPL is enhanced around wake center. With the prescribed sound profile of Figure 12b, the RSPL is refracted downwards, as shown in Figure 16c. Similarly, using the sound profile given in Figure 12c, the resultant sound propagation is refracted upwards, see Figure 16d. Similar as the previous simulations, the following study focuses on the noise propagation area between 1120 m and 2500 m, and a higher frequency f = 1000 Hz is selected. In this region, the three individual wind turbine noise sources propagate and superimpose on each other. In Figure 16a, the results cover the entire domain behind the last wind turbine. The solution differs from Figure 13a in the sense that it is the RSPL from all the three wind turbines. In Figure 16b, the uniform sound profile at the inlet boundary is affected by the multi-wake (see Figure 12a), such that the resulted RSPL is enhanced around wake center. With the prescribed sound profile of Figure 12b, the RSPL is refracted downwards, as shown in Figure 16c. Similarly, using the sound profile given in Figure 12c, the resultant sound propagation is refracted upwards, see Figure 16d. Figure 17 better illustrates the superposition of the sound pressure propagates from each wind turbine, which is a decomposition of Figure 16d. In the first figure, noise is generated from wind turbine #1 and propagates through the wake flow. The other two wind turbines are depicted with white lines, which imply that they are currently excluded for noise propagation. In the next simulation, shown in Figure 17b, ANPropagation from wind turbine #2 is individually considered, which excluded noise propagation from wind turbine #1 and #3. Similarly, wind turbine #3 is then considered for ANPropagation, which does not take account of the two upstream wind turbine noise propagations. The total RSPL is the summation of sound pressure following the logarithmic sumation rule, which is shown in Figure 17d. Figure 18 shows the RSPL or ∆L extracted from a line at 2 m high. The overall sound pressure variation for different cases is smaller as compared to the previous results shown in Figure 14. This is due to the mixing of three individually simulated propagation field. Among all the test cases, Case 3 can be regarded as the worst-case scenario. For example, around the distance of x = 1800 m, there is an obvious increase of SPL due the strong downward refraction at this region, which is also observed in Figure 16c at x = 1800 m.

Result shown in
individual wind turbine noise sources propagate and superimpose on each other. In Figure 16a, the results cover the entire domain behind the last wind turbine. The solution differs from Figure 13a in the sense that it is the RSPL from all the three wind turbines. In Figure 16b, the uniform sound profile at the inlet boundary is affected by the multi-wake (see Figure 12a), such that the resulted RSPL is e  Figure 17 better illustrates the superposition of the sound pressure propagates from each wind turbine, which is a decomposition of Figure 16d. In the first figure, noise is generated from wind turbine #1 and propagates through the wake flow. The other two wind turbines are depicted with white lines, which imply that they are currently excluded for noise propagation. In the next simulation, shown in Figure 17b, ANPropagation from wind turbine #2 is individually considered, which excluded noise propagation from wind turbine #1 and #3. Similarly, wind turbine #3 is then considered for ANPropagation, which does not take account of the two upstream wind turbine noise propagations. The total RSPL is the summation of sound pressure following the logarithmic sumation rule, which is shown in Figure 17d.  Figure 17 better illustrates the superposition of the sound pressure propagates from each wind turbine, which is a decomposition of Figure 16d. In the first figure, noise is generated from wind turbine #1 and propagates through the wake flow. The other two wind turbines are depicted with white lines, which imply that they are currently excluded for noise propagation. In the next sim (a)  Figure 18 shows the RSPL or ∆L extracted from a line at 2 m high. The overall sound pressure variation for different cases is smaller as compared to the previous results shown in Figure 14. This is due to the mixing of three individually simulated propagation field. Among all the test cases, Case 3 can be regarded as the worst-case scenario. For example, around the distance of x = 1800 m, there is an obvious increase of SPL due the strong downward refraction at this region, which is also observed in Figure 16c at x = 1800 m. Similar to the procedures previously mentioned for Figure 15, Figure 19a-c show the noise spectra from wind turbines #1, 2, and 3, respectively. The difference is that the three wind turbines are individually considered for its noise generation and propagation at the given frequencies (20 Hz-5000 Hz). Due to the presence of velocity deficit in the wake, the inflow conditions of the two downstream wind turbines are slightly different, and thus the ANSources of the three wind turbines are slightly different from each wind turbine. In each figure, the original ANSource is shown on top of the other noise spectra. Similar to the situations in Figure 15, geometrical attenuation plays an primary role in the large noise reduction for ANPropagation at f < 1000 Hz, and the noise absorption by air has a strong effect at f > 1000 Hz. As expected, Case 3 shows more distinct characteristics for Similar to the procedures previously mentioned for Figure 15, Figure 19a-c show the noise spectra from wind turbines #1, 2, and 3, respectively. The difference is that the three wind turbines are individually considered for its noise generation and propagation at the given frequencies (20 Hz-5000 Hz). Due to the presence of velocity deficit in the wake, the inflow conditions of the two downstream wind turbines are slightly different, and thus the ANSources of the three wind turbines are slightly different from each wind turbine. In each figure, the original ANSource is shown on top of the other noise spectra. Similar to the situations in Figure 15, geometrical attenuation plays an primary role in the large noise reduction for ANPropagation at f < 1000 Hz, and the noise absorption by air has a strong effect at f > 1000 Hz. As expected, Case 3 shows more distinct characteristics for the frequencies. In Figure 19, the receiver is fixed at x = 2500 m, z = 2 m, therefore the sound pressure is largely attenuated for wind turbine #1 due to the longest propagation distance. A similar behavior is seen for wind turbines #2 and #3, where the noise level is the highest for wind turbine #3 due to the shortest relative distance from source to receiver. Figure 19a-c also show that at different relative propagation distances 2500 m (#1), 1940 m (#2) and 1380 m (#3), the distribution of SPL maintains some similarity. The differences of these sound spectra come from the combined effect of wake flow and the propagation distance. Finally, the total noise received from the three wind turbines can be got through logarithmically adding the sound spectra from individual wind turbines. In Figure 19d, the sources Lp1, Lp2, and Lp3 represent the ANSource level of the three wind turbines. The major difference is seen at the middle frequency range, which is due to the different input velocity data obtained from flow simulation. Case 1 and Case 3 all show that the most significant influence comes from wind turbine #3, which is not a surprise. To highlight more detailed difference in the noise spectra, a few typical values are compared in Table 1. At f = 20 Hz, the summated SPL from three turbines (Case 3) is 37.92 dB which is 3 dB higher than that from wind turbine #3 (Case 3). This difference becomes smaller when frequency increases. This implies that the low frequency noise generated from the far upstream wind turbines still plays an important role. On the other hand, if a simple sound propagation algorithm is applied, such as shown in the last column of Table 1, a large error is observed in a frequency range below 1000 Hz. It indicates that the high frequency noise will be heavily attenuated over long propagation distance except in special cases such as sound propagation over water (hard surface), stronger downward refraction under special climate conditions.   Finally, it is interesting to predict the noise spectra at different receivers. The receiver height is again fixed at z = 2 m, and the horizontal distances measured from wind turbine #1 are: x = [1620, 1820, 2020, 2220, 2420] m. The noise spectra are shown in Figure 20 and marked with different colors. The general trend is that the noise level decreases as observer distance increases. However, it does not always hold true at some distances and frequencies. For example, at x = 1820 m (22.75D), except at low frequencies, increase of noise level is observed as compared to the x = 1620 m case. The simulations show that the low frequency noise might be the major wind turbine noise source in the farfield. Since the spectra are obtained in farfield, the overall SPLs are quite low. According to the wind turbine noise regulations in many countries (such as 45 dB), the calculated noise evidently fulfill the limit. From the results, it is seen that SPLs over 1000 Hz is below 20 dB. For future studies, it is indicated that more focus should be paid on long-path ANPropagation at low frequencies.   Finally, it is interesting to predict the noise spectra at different receivers. The receiver height is again fixed at z = 2 m, and the horizontal distances measured from wind turbine #1 are: x = [1620, 1820, 2020, 2220, 2420] m. The noise spectra are shown in Figure 20 and marked with different colors. The general trend is that the noise level decreases as observer distance increases. However, it does not always hold true at some distances and frequencies. For example, at x = 1820 m (22.75D), except at low frequencies, increase of noise level is observed as compared to the x = 1620 m case. The simulations show that the low frequency noise might be the major wind turbine noise source in the farfield. Since the spectra are obtained in farfield, the overall SPLs are quite low. According to the wind turbine noise regulations in many countries (such as 45 dB), the calculated noise evidently fulfill the limit. From the results, it is seen that SPLs over 1000 Hz is below 20 dB. For future studies, it is indicated that more focus should be paid on long-path ANPropagation at low frequencies.

Conclusions
Combing the flowfield aerodynamic simulation and the aeroacoustics computation of multiwind turbines is the focus of this paper. To achieve this objective, ANSource modeling, flow simulation and ANPropagation are individually described and then integrated. The numerical computations are carried out using the in-house developed tool EllipSys. The ANSource modeling and ANPropagation modeling are implemented in the platform of EllipSys. The steady CFD

Conclusions
Combing the flowfield aerodynamic simulation and the aeroacoustics computation of multi-wind turbines is the focus of this paper. To achieve this objective, ANSource modeling, flow simulation and ANPropagation are individually described and then integrated. The numerical computations are carried out using the in-house developed tool EllipSys. The ANSource modeling and ANPropagation modeling are implemented in the platform of EllipSys. The steady CFD simulations are performed to obtain the flowfield with a relatively efficient 2D-AD solver. The ANSource of wind turbine is calculated using the flowfield information from CFD computations. Then the ANPropagation is simulated which needs the data from both CFD simulations and ANSource computations. Results from different case studies show that a long-path ANPropagation are greatly influenced by different flow conditions, atmosphere conditions and wake condition. An overall trend is that the ANSource from wind turbines is broadband, but the received noise level at long travelling distances from multi-wind turbines is complicated. The noise spectra show that the low frequency noise from a wind turbine cluster might be an important issue to be further investigated. It has been shown in the results that at a distance over 1 km (about 12 rotor diameter), the high frequency components are largely attenuated, say at f > 1000 Hz. In a large wind farm with superimposition of noise sources from multi-wind turbines, the computational load will be heavy if high frequency noise is simulated. The ray theory, as compared in Sections 3.1-3.2, is a high-frequency approach therefore the low frequency solutions are not very accurate if strong refraction and diffraction are represented. The present numerical approach has same accuracy over all frequency range. However, to solve high frequency ANPropagation, higher grid density is required, thus increasing computational time. Fortunately, results presented in Figures 18 and 19 indicate that it is the low frequency noise components that contribute to the long range wind turbine ANPropagation. The receiver distances shown in Figure 19 ranges from 20 to 30 rotor diameters. It is clear that at such a large distance, attenuation of low frequency noise is still not significant.