Eeffect of Blade Thickness on Noise Pollution of H-type Darrieus Wind Turbines: A Numerical study

Noise pollution is one of the biggest problems of wind turbines, especially when these turbines are located near residential areas. In this article, the effect of blade thickness is numerically investigated on the noise pollution of an H-type Darrieus wind turbine. The flow is first simulated using the unsteady Reynolds averaged Navier-Stokes equations and the SST-kω model at the tip speed ratio of 2.64. Then, the noise is calculated using Ffowcs Williams-Hawkings equations. Blade thickness is changed using NACA airfoils from NACA 0008 up to NACA 0024. It is concluded that noise calculation at only one point, known as a routine method in noise investigation of wind turbines, is insufficient to investigate the noise of this turbine. Here, maximum noise in directivity is defined as the criterion of noise pollution. The results show that changing the blade profile of the benchmark turbine from NACA 0021 to NACA 0015 increases the power coefficient from 0.318 to 0.371 and reduces the maximum noise from 95.67 dB (76.35 dB) to 90.19 dB (71.01 dB) at R = 2 m (8m). For NACA 0018, the power coefficient is 0.353, and the maximum noise is 89.78 dB (70.47 dB) at R = 2 m (8m). Overall, the highest output power is for NACA 0015, and the lowest noise pollution is for NACA 0018.


A B S T R A C T
Noise pollution is one of the biggest problems of wind turbines, especially when these turbines are located near residential areas. In this article, the effect of blade thickness is numerically investigated on the noise pollution of an H-type Darrieus wind turbine. The flow is first simulated using the unsteady Reynolds averaged Navier-Stokes equations and the SST-kω model at the tip speed ratio of 2.64. Then, the noise is calculated using Ffowcs Williams-Hawkings equations. Blade thickness is changed using NACA airfoils from NACA 0008 up to NACA 0024. It is concluded that noise calculation at only one point, known as a routine method in noise investigation of wind turbines, is insufficient to investigate the noise of this turbine. Here, maximum noise in directivity is defined as the criterion of noise pollution. The results show that changing the blade profile of the benchmark turbine from NACA 0021 to NACA 0015 increases the power coefficient from 0.318 to 0.371 and reduces the maximum noise from 95.67 dB (76.35 dB) to 90.19 dB (71.01 dB) at R = 2 m (8m). For NACA 0018, the power coefficient is 0.353, and the maximum noise is 89.78 dB (70.47 dB) at R = 2 m (8m). Overall, the highest output power is for NACA 0015, and the lowest noise pollution is for NACA 0018. doi: 10 Wind energy utilization has grown 233%, equal to 513 GW, in 2010s [1]. It is also predicted that the global wind energy market will annually have a 4% increase [2], which means the use of wind turbines will considerably increase in the future. Depending on the rotor axis position, wind turbines are divided into Vertical-Axis *Corresponding Author Email: bozorgi@arakut.ac.ir (A. Bozorgi) Wind Turbines (VAWTs) and Horizontal-Axis Wind Turbines (HAWTs). VAWTs are used in smaller sizes and are becoming increasingly popular in residential areas [2]. However, several studies [3][4][5] have shown that wind turbines adversely affect human health. Wind turbines are known as the cause of distress and anxiety for people living around wind farms [6,7]. Long-term exposure to the noise of wind turbines could cause several psychological and physiological disorders, such as stress, sleep disturbance, irritability, difficulty concentrating, fatigue and blood pressure [8]. The noise of VAWTs is even more annoying than HAWTs since VAWTs are usually installed in residential areas and have a higher rotational speed, generating humming, swooshing, and whistling sounds [9][10][11]. With the expanding use of these turbines in urban environments, their noise can constitute an important component of urban noise [12]. According to environmental standards used in many countries, the A-weighted Sound Pressure Level (SPL) in residential areas must be lower than 40 dBA for nights and lower than 50 dBA for days. In addition, some countries, such as Denmark, define a limit of 20 dBA for low-frequency noises . According to Iran's national noise pollution standard, the SPL is limited to 45 dBA for nights and 55 dBA for days in residential areas. However, there is no limitation for low-frequency noises. On the other hand, there is no effective deterrence for preventing lawbreakers, not regarding the standard of noise pollution.
H-type Darrieus Wind Turbine (HDWT) is one of the widely used VAWTs installed in low-capacity potentials. In recent years, several studies have been done on the aerodynamic analysis of HWDTs. For example, Hashem et al. [13] and Parakal et al. [14] investigated the effect of blade shape on the aerodynamic performance using Computational Fluid Dynamics (CFD). Hashem and Mohamed [13] numerically simulated the aerodynamic performance of a three-bladed HDWT by using the Unsteady Reynolds Averaged Navier-Stokes (URANS) equations. This two-dimensional simulation exhibited that the HDWT with S1046 airfoil has the highest output power. Parakkal et al. [14] also simulated a three-bladed HDWT in two dimensions and investigated the effect of blade geometry on the output power of the turbine. In their study [14], a Joukowski airfoil was compared against NACA 0012 and NACA 4312 airfoils. The results revealed that the Joukowski airfoil has higher torque and power coefficient almost over the entire operable Tip Speed Ratio (TSR), but negatively affects the self-starting ability of the turbine. Huang et al.'s study [15] showed that reducing the radius distance for part of the blades increases the self-starting ability by improving the output power in low TSRs. It is a significant advantage for using these turbines in residential areas. Nemati [16] and Abid et al. [17] showed that the combination of Savonius and Darrieus turbines improves the self-starting ability. Singh et al. [18] also investigated self-starting performance for a three-bladed HDWT. The results revealed that it is a self-staring turbine for solidities of 0.8 to 1.2. A numerical simulation performed by Celik et al. [19] about the solidity showed that the increase in blade number improves the self-starting performance of the turbine; however, this may reduce the output power. Moreover, numerical studies of Tong et al. [20] and Maalouly et al. [21] confirm that solidity significantly affects the aerodynamic performance of HDWTs. Maalouly et al. [21] revealed that solidity significantly impacts the aerodynamic performance for transient and steady operations.
The effect of blade profile on the noise pollution of an HDWT has been evaluated by Mohamed [22]. The simulation has been performed for NACA 63418, S1046, FXLV152, and NACA 0018 airfoils in two dimensions using the URANS and Ffowcs Williams-Hawkings (FW-H) equations [23]. The results showed that blade profile significantly affects noise generation, and the lowest noise belongs to S1046 airfoil. Ghasemian and Nejat [24] calculated noise pollution of an HDWT using Large Eddy Simulation (LES) and FW-H equations and showed that the noise is logarithmically reduced by increasing observer distance from the turbine. Su et al. [25] simulated in three dimensions an HDWT and showed that the noise increases with increasing the TSR. The effect of pitch angle on the aerodynamic performance and noise pollution of an HDWT was investigated by Karimian et al. [26]. The results showed that positive angles reduce the output power, but the noise reduces for angles of 1 o to 3 o .
In this paper, the effect of blade thickness on the output power and noise pollution of an HDWT is investigated using NACA 0008, NACA 0012, NACA 0015, NACA 0018, NACA 0021 and NACA 0024 airfoils. For this purpose, the flow is first simulated in two dimensions using the continuity and URANS equations and the SST-kω turbulence model. Then the noise is calculated by using the FW-H equations. In this work, several observers are defined on full circles around the turbine to investigate directivity. In most past studies, observers have been located in only one angle position (in wind direction and downstream of a turbine) [22,[24][25][26] that cannot comprehensively present the noise behavior of wind turbines.

Aerodynamic simulation
Here, the aim of the aerodynamic simulation is to calculate the power coefficient Cp at TSR = 2.64 and also to prepare parameters that the FW-H equations need for noise calculation. The Cp and TSR are defined by: where, R is rotor radius, ω rotational speed, U wind speed, P output power, ρ air density, and L rotor height [27]. Geometric characteristics and operating conditions of the benchmark HDWT are exhibited in Table 1. The blade profile corresponds to NACA 0021 airfoil. The geometry was modeled using ANSYS Design Modeler, and numerical simulation was performed using ANSYS Fluent. The two-dimensional continuity and URANS equations [27] were applied for flow simulation. The equations are written for incompressible flow as: where u i and p denote resolved velocity components and pressure, respectively. Here, the Reynolds stress term (ρu i ′ u j ′ ) is solved using the SST-kω turbulence model. This model was developed by Menter [28] to effectively blend the robust and accurate formulation of the k-ω model in the near-wall region with the free-stream independence of the k-ɛ model in the far field. The SSTkω model is more accurate and reliable than the standard k-ω model for a wider class of flows, such as adverse pressure gradient flows and airfoils. Flow domain and boundary conditions are exhibited in Figure 1a. Mohamed et al. [29] investigated the convenient size of the domain and explained that the domain must spread more than ten times the rotor diameter in all directions. It has been obeyed in the current simulation. The domain consisted of rotary and stationary domains. A boundary layer mesh was defined on blade walls (Figure 1b) such that the mean of Y + on the wall of the blades was 1.9, which is agreeable for the SST-kω model. In Figure 2, Y + on the wall of the blades has been shown. Velocity inlet boundary condition was used for the upstream and pressure outlet for the downstream. Mesh independency was checked using grid sizes from 119.597 to 779.687 cells. The study revealed that over 232,658 cells, the relative deviation of Cp calculation was less than 1% ( Figure 3). Here, the grid size of 381,324 cells was selected for all simulations to reduce the computation time. In this grid, the average orthogonal quality on the rotating and stationary zones was 0.943 and 0.995, respectively, and the equiangle skewness on the rotating and stationary zones was 0.162 and 0.014, respectively.
The SIMPLE algorithm was used for pressurevelocity coupling, and the Green-Gauss cell-based   method was used for gradient calculation. The secondorder upwind method was selected to discretize convection terms and turbulence equations. The secondorder method was used for pressure terms and the secondorder implicit method for transient terms.
Each period was divided into 128 time-steps. with 20 internal iterations in each time-step. It means a whole rotation of the three-bladed turbine was divided into 384 time-steps (less than one-degree rotation per timestep).
The simulations were carried out on 12 processors, 4.7 GHz PC, which requires a total CPU time of about 24 min for every revolution. All simulations were carried out for at least ten revolutions, while the results converged after five revolutions.
Comparing the CFD results of the current work and other numerical simulations [13,27,[30][31][32] with experimental data [30] (Figure 4) exhibits that the CFD results of the current work are the nearest curve to the experiment. The maximum difference in Cp is at the TSR of 2.33, equal to 0.42, while the minimum difference is at TSR of 2.64, equal to 0.006. The results show that the current numerical simulation using the SST-kω model is acceptable for predicting the aerodynamic performance of this HDWT.

Noise calculation
The noise of the HDWT was calculated in far-field using the FW-H equations. The equation is a rearrangement of the Navier-Stokes equations widely used for calculating the noise of moving bodies [23], such is wind turbines. The equation by defining blade walls as data surface (where f = 0) is given by: Here, Q, Li and Tij are monopole, dipole and quadrupole sources. c0 is the speed of sound, vn the normal component of data surface velocity, Pij compressive stress tensor, σij compressive stress tensor, and nj the normal unit vector of data surface. A "0" subscript represents mean quantities. Because of the Dirac delta function δ(f), the Q and Li terms are only defined on the data surface, while the Tij term must be defined out of the data surface because of the step function H(f). It should be noted that the quadrupole sources are ignored in noise calculations [22,[24][25][26] since calculating them requires surface integration, which is expensive. Moreover, Ansys Fluent uses an integral solution of the FW-H equations presented by Farassat [33] (formulation 1A) that applies only Q and Li terms. In this work, the far-field density was defined as equal to the density of the incoming flow, i.e., 1.225 kg/m 3 . The speed of sound was 340 m/s; the reference pressure was 2e-5 Pa, and the source correlation length was 1 m. The convective effect option was activated to consider the effect of mean flow (wind speed) on the speed of sound.
Finally, 24 observers were defined on two circles with radiuses of 2 m and 8 m ( Figure 5).

RESULTS AND DISCUSSION
The results of Cp at TSR = 2.64 are shown in Figure 6 for different blade thicknesses. It is observed that the blade profile has a significant effect on Cp. The Cp increases by changing the blade profile of the validated case from NACA 0021 to NACA 0012, NACA 0015, NACA 0018 or NACA 0024. In contrast, NACA 0008 decreases the Cp. The results for NACA airfoils thinner Than NACA 0008 are not presented because they had negative Cp, and therefore, cannot be used. The numerical results show that the change in blade profile from NACA 0021 to NACA 0015 increases the Cp from 0.318 to 0.371, a considerable improvement in the output power (16.66% increase). The numerical study of Mohamed et al. [27] on several VAWTs showed that the maximum Cp in VAWTs is in the range of 0.3-0.4. As shown in Figure 4, the CFD results of Mohamed et al. [27] are more overpredicted than the current CFD results, and therefore, the Cp of 0.371 obtained in the current CFD work can be considered as a high value.   According to Betz's law, the limit of Cp in HAWTs is 0.593, but practically the Cp of existing HAWTs is lower than it. Moreover, the Cp of VAWTs is lower than the Cp of HAWTs. For example, Fadil and Ashari [34] compared the performance of a VAWT with a HAWT that both of them had a swept area of 3.14 m 2 , three blades with NACA 4412 profile, and showed that the maximum Cp of the HAWT is 0.54, while the maximum Cp of the VAWT is 0.34.
The pressure contour is shown in Figure 7 for all blade thicknesses at t = 0.680943 s. The results show that the pressure distribution in the flow domains depends on blade thickness. Moreover, the pressure distribution around the blades of a rotor is different from each other. As shown in Equation 6, dipole sources Li define noise sources generated by pressure distribution around the wall of the blades. Therefore, it is concluded that the dipole sources generated by the blades of a rotor are different form each other at identical times. The results of directivity are shown in Figure 8. It is concluded that there is no similarity in the directivity profiles for different distances. It means that the rotor is not similar to a point noise source. Moreover, it is shown that the noise reduces with increasing the observer distance, except where observers located in ϴ = 0 o for NACA 0015 or NACA 0018 blade profiles. For these cases, the time history of sound pressure received by the observers during one period (T) is shown in Figures 9 and  10. The rotor's sound pressure is equal to the sum of all blades' sound pressure. Comparing Figure 9a ( Figure  10a) with Figure 9b (Figure 10b) shown that the sound pressure of corresponding blades reduces with increasing the distance from 2 m to 8 m (at ϴ = 0 o ). However, the interaction of the blades is such that the sound pressure of the rotor increases. The directivity for identical distances is shown in Figure 11. The results show that noise pollution comparison of the different blade thicknesses cannot be performed by calculating the noise in only one observer point. For example, a comparison of the noise of NACA 0008 and NACA 0021 shows that at ϴ = 0 o , the noise of NACA 0008 is higher, while at ϴ = 90 o , the noise of NACA 0021 is higher (Figures 11a and 11b).  The IEC 61400-11 standard defines only one necessary point with ϴ = 0 o for noise calculation of such wind turbines, and according to it, in most past studies, the noise calculation has been performed only at ϴ = 0 o [22,[24][25][26]. However, Figure 11 shows that the maximum noise is not located at ϴ = 0 o , and even in some cases, the minimum noise is at this ϴ, e.g. NACA 0024 for R = 2 m (Figure 11a). For more detailed investigations, the maximum noises for all thicknesses are shown in Table 2.
It is observed that the ϴ of maximum noise is not identical for all blade thicknesses, and moreover, in most blade thicknesses, the ϴ of maximum noise changes when R changes. Therefore, defining a specific ϴ for the angle position of maximum noise is impossible, and directivity must be calculated for this purpose. HDWTs are usually used in areas with variable wind direction. Since the directivity profile rotates equal to the change of wind direction, and wind direction is variable, maximum noise in directivity profile can be located in all observer points. Here, maximum noise in the directivity profile was defined as the criterion of noise pollution. The results show that by changing the blade profile from NACA 0021 to NACA 0012, NACA 0015, NACA 0018 or NACA 0024, the noise pollution reduces. In contrast, NACA 0008 increases noise pollution. Overall, it is concluded that NACA 0008 has the most noise pollution and NACA 0018 has the least noise pollution. Considering the results of Figure 6 (Cp) and Table 2 (Maximum noise), it can be concluded that NACA 0008 has the lowest output power and the most noise pollution, while NACA 0015 has the highest output power and NACA 0018 has the lowest noise pollution. The performance of NACA 0015 and NACA 0018 are close to each other. They have less than 0.02 difference in Cp (the Cp of NACA 0015 is more) and less than 1.0 dB difference in SPL for all observer points (the SPL of NACA 0018 is lower). The Tij term in Equation (6), ignored in the integral solution of the FW-H equations (formulation 1A [33]), represents the effect of vortices and wakes on noise generation [35,36]. The vorticity contours in Figure 12 show that the vorticity near the blades is stronger in all rotors. In VAWTs, vortex shedding from one of the blades can affect the other blades. As shown here, vortex shedding from the 1 st blade affects the 3 rd blade in all rotors (blade numbering is according to Figure 1). Moreover, a qualitative comparison between these contours depicts that the strongest and widest vortices are for NACA 0008. Therefore, it is expected that the Tij term is the highest for NACA 0008 and causes the highest increase in the noise of this blade profile. The frequency distribution for some observer points is depicted in Figures 13 and 14. The frequency distribution is calculated by using Fast Fourier Transform (FFT). The results show that the first peak in all points is located at three times the rotor frequency, which is Blade Passing Frequency (BPF) and is equal to: BPF = blade numbers × ω(rad/s) 2π = 3 × 46.14 2π = 22.03 Hz The results show that the first peak is the highest in all cases, but the shape of frequency distribution is not similar for different blade thicknesses. For example, the directivity profiles of NACA 0015 and NACA 0018 are very close (Figure 11), but the shape of frequency distribution for identical observers is not similar (Figures  13 and 14). On the other hand, the shape of frequency distribution also depends on observer position. Therefore, in addition to directivity, frequency distribution is needed to have a comprehensive study on the noise behavior of this wind turbine.

CONCLUSION
In the present work, the effect of blade thickness on the output power and noise pollution of HDWTs was numerically investigated using CFD. For this purpose, a benchmark HDWT having three blades with NACA 0021 profile was selected. The present CFD results were in very good agreement with the experiment compared to other numerical works. In the following, the effect of blade thickness was investigated using NACA 0008, NACA 0012, NACA 0015, NACA 0018, NACA 0021 and NACA 0024. The flow was first simulated using the continuity and URANS equations and the SST-kω turbulence model at TSR = 2.64. Then, the FW-H equations were used to calculate SPL in 24 observer points defined on circles with 2 m and 8 m radiuses. The results were used to obtain directivity at R = 2 m and 8 m.
The results showed that blade thickness considerably affects the output power and noise pollution. Moreover, it was concluded that noise calculation at only one point (a routine method in noise calculation of wind turbines) is insufficient for investigating the noise pollution of this HDWT. For this purpose, maximum noise in directivity profile was defined as the criterion of noise pollution.
Changing the blade profile of the benchmark HDWT from NACA 0021 to NACA 0015 increased the Cp from 0.318 to 0.371 and reduced the maximum noise from 95.67 dB (76.35 dB) to 90.19 dB (71.01 dB) at R = 2 m (8m). For NACA 0018, the Cp was 0.353, and the maximum noise was 89.78 dB (70.47 dB) at R = 2 m (8m). Overall, the highest output power was for NACA 0015, and the lowest noise pollution was for NACA 0018.
Finally, frequency distribution was investigated using FFT, and the results showed that the highest peak was located in the first BPF in all cases. It was observed that the shape of frequency distribution depends on blade thickness and observer position. It means that in addition to directivity, frequency distribution should be investigated to have a comprehensive study on the noise behavior of this HDWT.