**Experimental Study on Aerodynamic Characteristics of a Gurney Flap on a Wind Turbine Airfoil under High Turbulent Flow Condition**

**Junwei Yang 1,2,3, Hua Yang 1,2,\*, Weijun Zhu 1,2, Nailu Li 1,2 and Yiping Yuan <sup>4</sup>**


Received: 20 September 2020; Accepted: 14 October 2020; Published: 16 October 2020

**Abstract:** The objective of the current work is to experimentally investigate the effect of turbulent flow on an airfoil with a Gurney flap. The wind tunnel experiments were performed for the DTU-LN221 airfoil under different turbulence level (T.I. of 0.2%, 10.5% and 19.0%) and various flap configurations. The height of the Gurney flaps varies from 1% to 2% of the chord length; the thickness of the Gurney flaps varies from 0.25% to 0.75% of the chord length. The Gurney flap was vertical fixed on the pressure side of the airfoil at nearly 100% measured from the leading edge. By replacing the turbulence grille in the wind tunnel, measured data indicated a stall delay phenomenon while increasing the inflow turbulence level. By further changing the height and the thickness of the Gurney flap, it was found that the height of the Gurney flap is a very important parameter whereas the thickness parameter has little influence. Besides, velocity in the near wake zone was measured by hot-wire anemometry, showing the mechanisms of lift enhancement. The results demonstrate that under low turbulent inflow condition, the maximum lift coefficient of the airfoil with flaps increased by 8.47% to 13.50% (i.e., thickness of 0.75%), and the Gurney flap became less effective after stall angle. The Gurney flap with different heights increased the lift-to-drag ratio from 2.74% to 14.35% under 10.5% of turbulence intensity (i.e., thickness of 0.75%). However, under much a larger turbulence environment (19.0%), the benefit to the aerodynamic performance was negligible.

**Keywords:** wind tunnel experiment; wind turbine airfoil; turbulence; Gurney flap; aerodynamic characteristics

#### **1. Introduction**

Wind power generation technology has been maturely developed in the past decades. Airfoil is a basic element of a wind turbine blade, and its aerodynamic characteristics have a major influence on the wind energy conversion efficiency. Among the conventional rotor aerodynamic design strategy, the blade add-ons were of particular interest to further improve wind energy efficiency. Therefore, mounting flap to the airfoil trailing edge was one of the most feasible methods to improve the aerodynamic performance of wind turbines. In addition to power production, such a technique can also effectively reduce the aerodynamic loads both of wind turbine blades and tower. If a sophisticated controller was implemented, the flap can further reduce turbulence-induced fatigue loads, so that longer lifetime was guaranteed. After the pioneering work of Liebeck [1], a large number of studies have been conducted to explain the phenomena induced by the presence of this device. More recently,

research objects were most focused on airfoil attached various shapes of flaps, such as Gurney flaps [2,3], triangular flaps [4], separate trailing edge flaps [5,6] and deformable trailing edge flaps [7,8]. For experiment tests, Zhang et al. [9] and Amini et al. [10] studied the aerodynamic effect of a Gurney flap based on the airfoil through wind tunnel experiments. T Lee et al. [11] carried out a wind tunnel test of the lift force and pitching moment coefficients of both trailing edge flaps and Gurney flaps with different shape parameters. Based on a 5 MW reference wind turbine, Chen et al. [12] designed and optimized a trailing edge flaps such that the blade mass can be further reduced but still maintain the desired power performance. Medina et al. [13] explored the flow mechanisms of a flap at a high deflection angle. When the flap works in a separated flow region, he provided some ideas for realizing instantaneous action or alleviating extra aerodynamic loads on wind turbines. Elsayed et al. [14] studied the flap tip vortexes and characterized the flow structures behind a flap in a low-speed wind tunnel by using particle image velocimetry. Little et al. [15] designed a trailing edge flap by using a single medium plasma driver which resulted in a higher lift force. Bergami et al. [16] designed an active controller of a trailing edge flap on a 5 MW reference wind turbine and proved that the flap could effectively control aerodynamic loads. Edward et al. [17] also conducted experiments on a flaps noise drop. According to the wind tunnel tests based on a full-size rotor, Straub et al. [18] found that flaps could also be used to control noise and vibration such that the noise generated by blade vortex interaction could be reduced by 6 dB. For numerical simulations and theoretical analyses, Traub et al. [19] fitted a semi-empirical equation by summarizing the performance of a large number of flaps. Lario et al. [20] numerically solved the unsteady flow field of Gurney flaps at high Reynolds numbers through the discontinuous Galerkin method. By analyzing dynamic characteristics of an airfoil with Gurney flaps through the numerical simulation, Li et al. [21] found that flaps were capable to reduce unsteady aerodynamic loads of wind turbines. Zhu et al. [22] numerically simulated the airfoil with trailing edge flaps using the immersed boundary method and found that flaps could be combined with the paddle movement to adjust the aerodynamic loads of a wind turbine airfoil. Ng et al. [23] investigated the trailing edge flap together with an aeroelastic analysis. It was noted that the trailing edge flap could be a smart device to control aeroelastic deformation.

All the above researches on the flaps were built on the uniform inflow of low turbulence intensity. In the past, there were few studies carried out by flow over flaps under high turbulence intensity [24,25], however, such a flow condition often occurs on wind turbines operating in a wind farm. Considering wind turbines operate in a turbulence environment, the conclusions obtained from the previous studies might not be accurate. Therefore, to simulate wind turbines under turbulence environment, active and passive wind tunnel turbulence generation methods can be used. The active technology includes a vibrating grille and multi-fan wind tunnel; the maximum turbulence intensity can reach more than 20% [26]. The passive control structure was relatively simple, which can be divided into grille, wedge, and rough square types among the passive turbulence generators, it was more convenient to construct grilles, which have gain very popular use. In this experiment, specific turbulence levels were passively controlled by a grille with proper grid size.

The investigations presented in this paper were focused on the coupled effects of Gurney flap and turbulence inflow. The Gurney flaps were experimentally investigated under various turbulence intensities and flap configurations. The desired turbulent flow passes the airfoil was achieved by changing the grille size as well as the distance between the grille and the airfoil model. The hot-wire anemometer was used to record the wind speed and turbulence intensity in a flow cross-section. The quantitative information obtained during the experiment includes: (1) turbulent field descriptions, (2) airfoil pressure coefficients, (3) lift-to-drag coefficients, (4) wake measurements. On that basis, the aerodynamic performance of Gurney flaps with different heights and thicknesses was researched. The values of lift and pitching-moment coefficients were obtained through the integration of surface pressures. The wake rake array was also used to determine the values of drag coefficient. Furthermore, the flow fields near the trailing edge of airfoil were tested which could further

verify the reason for lift improvement. Finally, concluding remarks accompany the discussion of the experimental investigations.

#### **2. Experimental Setup**

The experiments were carried out in a wind tunnel located at Yangzhou University. It was a close-loop type wind tunnel which contains two experimental sections. The measurements were conducted in the smaller section with the cross-section parameters of 3 m × 1.5 m, and the length is 3 m. The operational wind speed range was 0~50 m/s and the calibrated maximum turbulence intensity in the free stream was 0.2%. The DTU-LN221 airfoil model [27], as shown in Figure 1, was adopted which has a chord length of 0.6 m and a span length of 1.5 m. The airfoil model was vertically placed in the test section. The bottom part of the airfoil section was connected by the rotating shaft, which was fixed on a rotational plate, as demonstrated in Figure 2. Therefore, the angle of attack can be remotely controlled via a shaft connection with a motor below the wind tunnel. The airfoil model was made of aluminum alloy where small taps were drilled on the surface of the middle section, with the location sketched in Figure 1. The arc length between the taps was approximately spaced every 2.5% of chord length starting from the leading edge to 90% of the trailing edge. The hollow plastic hoses were connected on the reverse side of the airfoil surface which has an outer diameter of 1.2 mm. Thus, a total of 77 pressure taps were arranged for the pressure measurement on the surface of the airfoil. The pressures on the airfoil surface were sensed by the electronic pressure acquisition system of PSI (Pressure Systems Inc., Hampton, VA, USA). The pressure system used in the experiment has a sampling frequency of 333.3 Hz, a measuring range of ±2.5 kPa and a measuring accuracy of ±0.05%. Figure 1 also shows the chord-wise position, so the lift and pitching-moment of the test airfoil were achieved by integrating the pressure and the position of the taps on the surface. The lift coefficient and pitching-moment coefficients were given by

$$\begin{aligned} \mathbf{C}\_{l} &= \mathbf{C}\_{n} \cos \alpha - \mathbf{C}\_{l} \sin \alpha\\ \mathbf{C}\_{\text{m}} &= \int\_{0}^{1} (\mathbf{C}\_{pl} - \mathbf{C}\_{\text{pu}}) (0.25 - \hat{\mathbf{x}}) d\hat{\mathbf{x}} \end{aligned} \tag{1}$$

where

$$\begin{aligned} \mathsf{C}\_{n} &= \underset{\hat{0}}{\overset{\circ}{\underset{0}{\mathbb{C}}}} (\mathsf{C}\_{pl} - \mathsf{C}\_{pu})d\hat{x} \\ \mathsf{C}\_{t} &= \underset{\underset{\hat{0}}{\underset{0}{\text{max}}}}{\overset{\circ}{\underset{0}{\underset{0}{\text{max}}}}} (\mathsf{C}\_{p,he} - \mathsf{C}\_{p,af})d\hat{y} \end{aligned} \tag{2}$$

where *Cl* and *Cm* are the lift coefficient and pitching-moment coefficient, respectively; αrepresents the angle of attack; *x*ˆ is the relative chord length; *y*ˆ is the thickness value relative to the chord length; *Cpu* and *Cpl* are the pressure coefficient on the suction and pressure sides of the airfoil, respectively; *Cn* and *Ct* are the normal force coefficient and tangential force coefficient, respectively. *Cp*,*be* and *Cp*,*af* are the pressure coefficient before and after the maximum thickness of airfoil; *y*ˆ*u max* and *y*ˆ*l max* represent the maximum thickness values of the suction and pressure sides relative to chord length.

**Figure 1.** Schematic diagram of measuring position.

**Figure 2.** Experimental test device.

At a small angle of attack, the aerodynamic drag largely consists of frictional drag, while the data measured by surface pressure taps cannot accurately represent the drag force, so the drag can be measured by the momentum method more precisely. As shown in Figure 2, a wake rake array was placed at 0.7 chord length behind the trailing edge of the airfoil and at the same vertical level as the pressure taps. The measurement range of the wake probes was 80.8 cm, 102 total pressure pipes (with an outer diameter of 1.2 mm) and 4 static pressure pipes (with an outer diameter of 2 mm) were averagely arranged, with an 20 cm apart for the static pressure pipes. To prevent air leakage, all pressure tubes were connected by plastic hoses to the pressure measurement device. Besides, two Pitot tubes were installed at 1.6 m measured from the downstream of grilles where the free stream velocity was recorded. The drag coefficients were given as follows

$$\mathcal{C}\_{d} = \frac{2}{c} \int\_{w} \sqrt{\frac{P\_{01} - P}{P\_{0} - P\_{\infty}}} \Big| 1 - \sqrt{\frac{P\_{01} - P\_{\infty}}{P\_{0} - P\_{\infty}}} ds \tag{3}$$

where *Cd* is the drag coefficient; *c* is the chord length, *w* is the range of integration; *s* is the coordinate along the thickness direction of the wake rake array; *P*<sup>∞</sup> and *P*<sup>0</sup> are the static pressure and total pressure measured by the Pitot tubes; *P* and *P*<sup>01</sup> are the static pressure and total pressure measured by the wake rake array. It should be noted that there is a total pressure loss along the Pitot tubes to the wake rake array, therefore the total pressure loss should be added to each total pressure measuring point of the wake rake array.

Figure 3 presents the schematic diagram of the grille geometry. The grilles assembled with many squared alloys with geometry specified by four parameters a, b, c and d. These small squares were

bolted together and on top of the grille the rubber pads were attached. During the measurements, two types of grilles were implemented. The dimensions of the two grilles were provided in Table 1. The width of the longitudinal grille and the transverse grille were denoted by a and c, respectively; the width and the height of the spacing were denoted by b and d, respectively, the thickness along the flow direction was e.

**Figure 3.** Schematic diagram of experimental grille. (**a**) Front view; (**b**) side view.


**Table 1.** Dimensions of the experiment grilles.

Behind the grille, there were eight measurement sections designed for velocity recording. At each of the section, 15 measurement points were uniformly spaced. The center of the topmost measurement point corresponds to the geometrical center of the wind tunnel section. The horizontal and vertical spacing of measuring points were 10 cm, as displayed by the black spots in Figure 3a. In addition, we set eight measurement sections one to two meters behind the grille, as shown in Figure 3b. The wind speed was measured through a Dantec hot-wire probe (type 55P61). The sampling frequency was 5 kHz, and the test wind speeds were consistent with that in the airfoil experiment. Figure 4 displays the device of turbulence measurement. Besides, wake measurements in this paper were carried out on an electronically controlled traverser consisting of a two-dimensional probe, as displayed in Figure 5.

**Figure 4.** Grille turbulence generated device.

**Figure 5.** Wake measurements.

The influence of the flap height and thickness were two key parameters which influence the overall aerodynamic performance of the airfoil. For this reason, a total of nine types of flaps were designed in the experiments, namely, the flap with heights h of 1%, 1.5% and 2% chord length (6 mm, 9 mm and 12 mm, respectively) and with relative thicknesses d of 0.25%, 0.5% and 0.75% chord length (1.5 mm, 3 mm and 4.5 mm, respectively). The diagram of the flap used in the experiment was given in Figure 6. As shown, the flap has an L shape; the bottom part has a width of 7 mm which was attached to the pressure side of the airfoil trailing edge. During the experiments, the flap was vertically attached near the position of 100% of the chord length, see Figure 2 with the specific experimental setup. The coordinate system behind the trailing edge was defined for subsequent wake measurements, also the x and y axis were set in the vertical and the forward directions, respectively. The black spots in the Figure 6 represented the hot-wire measuring positions set at downstream of the trailing edge.

**Figure 6.** Flap geometry.

#### **3. Experimental Results and Analysis**

In the experiment, the free stream flow velocity was 20 m/s, and the corresponding Reynolds number based on the airfoil chord length was 0.8 <sup>×</sup> 106, although wind turbines often operate at wind speed below 20 m/s, but the magnitude of the Reynolds number was up to an order of 106. Therefore, the experimental value of the Reynolds number was chosen to approach an order of magnitude corresponding to those obtained from full-scale wind turbines. Two types of grilles generate different turbulence levels in the wind tunnel, such that the experiments were mainly separated into low and high turbulence cases. To be specific, a grille was placed at the upstream of the airfoil test section which results in different aerodynamic characteristics of flow over the airfoil with and without a Gurney flap.

#### *3.1. Experimental Results of Grille Turbulence Generator*

The measurements were first conducted for the two grilles. Figure 7 illustrates the time series of the instantaneous velocity at the central measuring point of the 6th section. As shown in Figure 7, the incoming wind became highly disturbed, for this reason, a statistical analysis was done for the flow field. The decay of turbulence intensity from the grilles towards the airfoil test section was reported in Figure 8a. Within the scope of the test, the averaged turbulence intensity was recorded along 8 downstream sections (each measured with 15 points). Besides, the margin of error in Figure 8a was expressed as the standard deviation of each section data. According to Figure 8a, the average turbulence intensity of downstream in scheme one dropped from 13.5% at the position of 1.0 m behind the grille to 9.4% at the position of 2.0 m; the average turbulence intensity of the downstream in scheme two was higher than scheme 1, which decreased from 27.7% at the position of 1.0 m behind the grille to 15.9% at the position of 2.0 m. The standard deviations of turbulence intensity of the two schemes drops quickly towards the test section. The change in turbulence level indicates that there were very high turbulence fluctuations just behind the grille which yielded very unstable flow field. The flow became gradually stabilized owing to the mutual dissipation over a propagation distance. To ensure a homogeneous and isotropic turbulent field, so the distance 1.6 m (from the leading edge of the airfoil to the grille) was selected for placing the leading edge of the airfoil. The turbulence intensity measured at the 6th section of the two schemes was 10.5% and 19.0%, respectively. Figure 8b illustrates the power spectral density (PSD) of the wind speed at the two grille schemes measured at the central measuring point of the 6th section (1.6 m measured from grille to the hot-wire probes). Obviously, the larger the turbulence intensity was generated, the larger the amplitude of the power spectral density was observed. Moreover, the power spectra amplitude in the downstream direction was consistent with the corresponding Karman spectra, and it dropped down rapidly after the frequency over 30 Hz, which suggests that the grille distance was more proportional to the turbulent kinetic energy.

**Figure 7.** The time series of the instantaneous streamwise velocity (the 6th section, central measuring point).

**Figure 8.** Description of the flow field behind the grilles. (**a**) Turbulence intensity and average velocity (Re: 0.8 <sup>×</sup> 106); (**b**) Comparison of the PSD (1.6 m behind the grille).

#### *3.2. Experimental Results of the Baseline Airfoil*

To verify the accuracy of the experiment, the experimental data from the wind tunnel of LM Wind Power were selected as a reference [28]. The aerodynamic data of the DTU-LN221 airfoil under uniform inflow were measured in the LM wind tunnel experiment at the approximate Reynolds number of 1.5 <sup>×</sup> 106. The lift and drag data were selected as a cross-validation case. In consideration of the interference effect of wind tunnel wall, the maximum blockage ratio in this experiment was about 8.4%, so the experiment measured data of the airfoil were corrected by after the reference [29], and the experimental data analyzed below were all corrected. In order to perform a comparison, Figure 9 shows the comparison of experimental data of the DTU-LN221 airfoil under uniform inflow condition. To minimize the uncertainty, measurements were performed several times for comparison test at the approximate Reynolds number of 1.5 <sup>×</sup> 106, YZU1, YZU2, YZU3 represents the results of the repeated experiments. As shown in Figure 9, the lift and drag coefficients were very similar before the stall angle, while there were some discrepancies in the stall state. Consider the cause of flow separation, the deviation of results of drag is quite considerable, but the stall angles tested in the two wind tunnels were nearly the same. The reason for this may be that although the airfoils were of the same type, they were individually manufactured, for example the roughness of the airfoil surface may affect the experimental results. Besides, there were some differences in drag measurements, LM adopted the model surface pressure distribution in a certain range of the angle of attack, while YZU used the wake probes at all angles. Similar phenomena were seen in different wind tunnel laboratories [30–32]. In the following, for reason of comparisons, the measured data are averaged.

**Figure 9.** Comparison of experimental data of the DTU-LN221 airfoil (uniform inflow, Re: 1.5 <sup>×</sup> 106). (**a**) Lift coefficient comparison; (**b**) drag coefficient comparison.

#### *3.3. Experimental Results of the Gurney Flap under Uniform Inflow*

This subsection starts to investigate the aerodynamic characteristics of the effect of Gurney flap in low turbulent inflow condition, the Reynolds number was decreased to 0.8 <sup>×</sup> 106. Seen from the Figure 10, three lift and drag curves were compared with the original at the angle of attack ranged from −9.6◦ to 14.4◦. As shown, the lift coefficients increased with the increase of the height of the Gurney flap before the stall angle, while it was also realized that the effect of the Gurney flap was limited at a high angle of attack. The lift coefficient slopes and the stall angles of the flapped airfoil remain essentially unchanged compared to the baseline airfoil, but the flapped airfoil caused a leftward shift compared to the baseline airfoil. In order to show the changes of aerodynamic characteristics clearly, the following lift and drag differences were defined:

$$
\Delta cl = \frac{\left| cl\_{\rm bs} - cl\_{\rm ff} \right|}{cl\_{\rm bs}} \times 100\% \tag{4}
$$

$$
\Delta cd = \frac{\left| cd\_{\rm bs} - cd\_{\rm gf} \right|}{cd\_{\rm bs}} \times 100\% \tag{5}
$$

$$
\Delta cld = \frac{|cld\_{\rm bs} - cld\_{\rm gf}|}{cld\_{\rm bs}} \times 100\% \tag{6}
$$

where *cl*bs and *cd*bs are the lift and drag coefficients obtained from the baseline airfoil, respectively. *cl*gf and *cd*gf are the lift and drag coefficients obtained from the flapped airfoil. *cld*bs and *cld*bs are the lift-to-drag ratio obtained from the baseline airfoil and flapped airfoil.

**Figure 10.** Variation of the aerodynamic characteristics with different heights of Gurney flaps (uniform flow, Re: 0.8 <sup>×</sup> 106). (**a**) Lift coefficient comparison; (**b**) drag coefficient comparison; (**c**) pitching-moment coefficient comparison; (**d**) lift-to-drag ratio comparison.

When the flap height changes to 6 mm, 9 mm and 12 mm, the maximum lift coefficients were increased by 8.47%, 9.56% and 13.50% at 9.4 degrees, respectively. Looking into the drag coefficients shown in Figure 10b, within the range of −9.6~10.4◦, the drag coefficients show regular change where the frictional drag dominates as expected. The drag coefficients rise rapidly in the stall state, at this time the change in drag coefficient was dominated by pressure drag. For example, at an angle of attack of 12.4◦, with the flap heights of 6 mm, 9 mm, 12 mm the drag coefficients were increased by 39.67%, 59.92% and 113.43%, respectively. According to Figure 10c, the Gurney flap generated a prominent increase in the pitching-moment coefficient compared to the baseline airfoil, the value increased with the heights of the Gurney flaps and reached negative peak values at approximately 4.4◦. Based on Figure 10d, the lift-to-drag ratios were compared. When the angle of attack was less than 8.4◦, the lift-to-drag ratios were all larger than the baseline airfoil. However, the maximum lift-to-drag ratios of the Gurney flap at three heights were smaller than the baseline airfoil. When the flap heights were 6 mm, 9 mm and 12 mm, the corresponding maximum lift-to-drag ratios were decreased by 17.16%, 22.79% and 24.47%, respectively. In the range after the stall angle, the presence of the Gurney flap reduces the lift-to-drag ratios and the higher the flap height, the more the lift efficiency decreases. When the angle of attack was 12.4◦, for example, the lift-to-drag ratios of the flap with heights of 6 mm, 9 mm and 12 mm were decreased by 26.55%, 38.11% and 52.20%, respectively. Meanwhile, we observed the same trend from the other thicknesses of the Gurney flaps (i.e., of 1.5 mm), When the flap height changes to 6 mm, 9 mm and 12 mm, the maximum lift coefficients were increased by 9.27%, 9.78% and 14.08% at 9.4 degrees.

The following study in this section will be focused on the effect of the thickness of the Gurney flap. Figure 11 illustrates the changes in the aerodynamic characteristics of the baseline airfoil and airfoil with a height of 9 mm of Gurney flaps under uniform flow. As seen from the figure, by changing the thickness of the Gurney flap, the lift-to-drag curves of the flapped airfoils with different thicknesses almost overlaps, it seems that the flap thickness does not have obvious effect in the aerodynamic performance. The results imply that the small increase of the flap thickness does not change the degree of downward turning of the mean flow and recirculation flow around the trailing edge of the airfoil.

**Figure 11.** *Cont*.

**Figure 11.** Variation of the aerodynamic characteristics with different thicknesses of Gurney flaps (uniform flow, Re: 0.8 <sup>×</sup> 106). (**a**) Lift coefficient comparison; (**b**) drag coefficient comparison; (**c**) pitching-moment coefficient comparison; (**d**) lift-to-drag ratio comparison.

#### *3.4. Experimental Results of Gurney Flap under 10.5% Turbulence Intensity*

This subsection starts to investigate the aerodynamic characteristics of the Gurney flap under a turbulent inflow with an intensity of 10.5%. Figure 12 shows the changes in the aerodynamic characteristics of the baseline airfoil and the Gurney flaps with a thickness of 4.5 mm when the turbulence intensity of the incoming flow was 10.5%. The measured angle of attack ranged from −9.6◦ to 20.4◦. As shown, under the present turbulence intensity of 10.5%, all cases induced a stall delay. The stall angle was nearly delayed from the original 9.4◦ to about 16.4◦ under turbulent flow. In Figure 12a, with the increase of turbulence intensity, the maximum lift coefficient of the baseline airfoil reaches 1.611 and increases by 17.59% as compared with uniform inflow condition. An interesting observation was that in such a turbulence environment, the lift coefficients were not significantly different from each other under the current three flapped configurations. For the flap heights of 6 mm, 9 mm and 12 mm, the maximum lift coefficients were increased by 15.21%, 17.19% and 17.26% at 18.4 degrees, respectively. The increase rate was found larger than that under the uniform turbulence condition. The margin of error in Figure 12a was expressed as the variance of the time-averaged lift coefficient. As can be seen from the Figure 12a, the pressure uncertainty increased with the increase of the angle of attack, and the lift coefficient fluctuations caused by the pressure fluctuations were the most obvious pronounced when the attack angle reaches the stall angle. According to Figure 12b, the increase in the lift of the flapped airfoil was at a cost of increasing drag. The drag coefficient of the flapped airfoils was also non-linearly raised in the stall state. For example, when the angle of attack was 20.4◦, the drag coefficients were increased by 12.50%, 37.50% and 44.53%, respectively. As can be seen from Figure 12c, compared with the baseline airfoil, the Gurney flap also generated a prominent increase in the pitching-moment coefficient, the negative peak values delay to about 6.4◦compared with the uniform inflow condition. As demonstrated in Figure 12d, the peak value of the lift-to-drag ratio becomes less as compared with the uniform flow condition in all cases. Unlike the low turbulence condition, the maximum lift-to-drag of Gurney flap with different heights slightly increased relative to the baseline airfoil. The maximum lift-to-drag raised 14.35%, 9.15% and 2.74% respectively when the flap heights were 6 mm, 9 mm, and 12 mm, respectively, with the maximum lift-to-drag ratio occurred at the smallest flap height. Besides, in contrast to the situations of 1.5 mm thickness of the Gurney flaps, the Gurney flaps also achieved a lift increase under 10.5% turbulence intensity. The maximum lift coefficient was increased by 14.77%, 14.83%, and 16.44% for height equals to 6 mm, 9 mm and 12 mm, respectively. The two thickness configurations both had maximum lift-to-drag efficiency at 6 mm height, for the three heights, respectively.

**Figure 12.** Variation of the aerodynamic characteristics with different heights of Gurney flaps (T.I. of 10.5%, Re: 0.8 <sup>×</sup> 106). (**a**) Lift coefficient comparison; (**b**) drag coefficient comparison; (**c**) pitching-moment coefficient comparison; (**d**) lift-to-drag ratio comparison.

In addition, the effects of various flap thicknesses were investigated. Figure 13 indicates the change of the aerodynamic characteristics of the baseline airfoil and airfoil with a height of 4.5 mm of the Gurney flaps when the turbulence intensity was 10.5%. As Figure 13 shows, like low turbulence conditions, after changing the thickness of the Gurney flap under turbulence condition of 10.5%, the aerodynamic characteristics of flapped airfoil almost show no change. The main difference in the Figure 13d was concentrated in the large lift-to-drag ratio zone, obviously because the pressure fluctuations caused by the turbulence and the flow separation at the airfoil surface result in the deviation of the measurement.

**Figure 13.** Variation of the aerodynamic characteristics with different thicknesses of Gurney flaps (T.I. of 10.5%, Re: 0.8 <sup>×</sup> 106). (**a**) Lift coefficient comparison; (**b**) drag coefficient comparison; (**c**) pitching-moment coefficient comparison; (**d**) lift-to-drag ratio comparison.

#### *3.5. Experimental Results of Gurney Flap under the 19.0% Turbulence Intensity*

To study the aerodynamic characteristics of the Gurney flap under high turbulence conditions, Figure 14 shows the changes of aerodynamic characteristics under the turbulence condition of 19.0%. The measured angle of attack ranged from −9.6◦ to 24.4◦. As seen in Figure 14, the presence of turbulence enhanced the lift coefficient and drag coefficient for all cases. In Figure 14a, with the increase of turbulence intensity, the maximum lift coefficient continually increased, and the maximum lift coefficient of the baseline airfoil reaches 1.781 at 18.5 degree, with a 30.00% increase relative to that under low turbulence condition. The maximum lift coefficients increased by 2.13%, 3.37%, 5.39% at 18.5 degrees, respectively relative to the baseline airfoil, when the heights of the flaps were 6 mm, 9 mm, and 12 mm. Furthermore, since the pressure fluctuations were determined by the fluctuation's energy of the turbulent flow, the fluctuations of lift coefficient were increased relative to that under the turbulence condition of 10.5%. According to Figure 14b, the drag coefficient of the baseline airfoil changes more obvious than those under the low turbulence condition. The drag coefficients were no longer gentle within the range before stall, and it also rises sharply in the stall state. Such a trend was more noticeable after adding the Gurney flap. For example, at the angle of attack was 22.4◦, the drag coefficients increased by 32.11%, 37.22% and 37.96%, when the flap heights were 6 mm, 9 mm, and 12 mm respectively. As can be seen from Figure 14c, the Gurney flap also produced a significant increase in the pitching-moment coefficient and such situation still existed in turbulent flow condition. Based on Figure 14d, unlike that under the condition of 10.5%, the lift gain obtained with Gurney flap is, unfortunately, coupled to a larger increase in drag. Adding Gurney flap with different heights at most of the angles reduces the lift-to-drag relative to the baseline airfoil, and the Gurney flap doesn't

show obvious effect under the turbulence condition of 19.0%. At the same time the same situation occurred at other thicknesses of the Gurney flaps.

**Figure 14.** Variation of the aerodynamic characteristics with different heights of Gurney flaps (T.I. of 19%, Re: 0.8 <sup>×</sup> 106). (**a**) Lift coefficient comparison; (**b**) drag coefficient comparison; (**c**) pitching-moment coefficient comparison; (**d**) lift-to-drag ratio comparison.

Figure 15 exhibits the changes in the aerodynamic characteristics of the baseline airfoil and airfoil with Gurney flaps of height 4.5mm under the turbulence condition of 19.0%. As can be seen from it, with the increase of turbulence intensity, obviously, there were more pressure fluctuations. The curves of flapped airfoils with different thickness no longer overlap totally, in the majority of angles, the difference from 2% to 5%, though it may sometimes reach 20%.

**Figure 15.** *Cont*.

**Figure 15.** Variation of the aerodynamic characteristics with different thicknesses of Gurney flaps (T.I. of 19%, Re: 0.8 <sup>×</sup> <sup>10</sup>6). (**a**) Lift coefficient comparison; (**b**)drag coefficient comparison; (**c**) pitching-moment coefficient comparison; (**d**) lift-to-drag ratio comparison.

#### *3.6. Surface Pressure Characteristics*

To analyze the mechanism of turbulent inflow coupled with Gurney flap, Figure 16 illustrates the pressure coefficient distribution of the baseline airfoil and airfoil attached with the Gurney flap under different turbulent inflow when Reynolds number Re = 0.8 <sup>×</sup> 106. The x-axis shows the chord-wise position x/c and the y-axis shows the pressure coefficients *Cp*. The pressure coefficients *Cp* on the airfoil surface were given

$$\mathcal{C}\_{\mathbb{P}} = \frac{p\_{\mathrm{i}} - p\_{0}}{0.5 \rho \mathcal{U}\_{0}^{2}} \tag{7}$$

where *p*<sup>i</sup> is the pressure at the pressure tap of i (i = 1:77), *p*<sup>0</sup> is the free stream static pressure at the airfoil tested by the Pitot tube, ρ is the air density and *U*<sup>0</sup> is the free stream velocity.

(**b**)

**Figure 16.** *Cont*.

**Figure 16.** The surface pressure coefficients under different angle of attack (Re: 0.8 <sup>×</sup> <sup>10</sup>6). (**a**) <sup>α</sup> <sup>=</sup> <sup>−</sup>9.6◦; (**b**) α = −3.6◦; (**c**) α = 2.4◦; (**d**) α = 8.4◦; (**e**) α = 12.4◦; (f) α = 18.4◦.

Figure 16 illustrates that the pressure coefficients of the baseline airfoil under different turbulence conditions with small angles were similar to each other. However, with the increase of angle of attack, the differences become more pronounced at the leading edge. In contrast to the Gurney flaps, Figure 16a shows that at the angle = −9.6◦, the differences between the flapped airfoils were small. With the increase of the angle of attack, as the angle of −3.6◦, compared with the baseline airfoils, the intersection of the airfoil with Gurney flap moves to the leading edge under different conditions. Besides, the pressure coefficients (with the angle = −3.6◦ and 2.4◦), the pressure distribution indicate a larger the pressure difference in the range of 0.8 < x/c < 0.9. When the angle of attack was 2.4◦ and 8.4◦, both the pressure coefficient absolute values of the baseline airfoil and the Gurney flap tendency to increase with the increase of turbulence intensity, and the flapped airfoil becomes more obvious than the baseline airfoil. At a large angle of attack, due to the condition of stalling (Figure 16e,f), there was significant flow separation both on the suction side of flapped airfoil and baseline airfoil under uniform

inflow. The larger the turbulent inflow, the larger the peak value on the leading edge of the suction side, and also the larger pressure coefficient on the pressure side.

As seen from the above, the influence of turbulent inflow coupled with Gurney flap on the pressure coefficients of airfoil was complicated. Turbulent inflow changes the flow around the airfoil surface and restrains the flow separation, after attaching the Gurney flap, the pressure distribution on both suction and pressure sides changed more significantly. Especially in the turbulence condition, the trailing edge pressure distributions were farther apart for the flapped airfoil. In these situations, the increased peak values of pressure coefficients at the leading edge increased the integral area of the pressure coefficients, so the lift force was therefore increased.

#### *3.7. Wake Profile Characteristics*

Figure 17 presents the distribution of wake velocity measured by the wake rake array under uniform inflow. As can be observed, when the angle of attack was 8.4◦, which was before the stall angle of attack, the wake of airfoil with the Gurney flap was inclined to the pressure side of the airfoil. When the angle of attack was 11.4◦, which was after the stall angle, the Gurney flap weaken the ability of the wake position deflection. Moreover, the wake velocity deficit was significantly increased, and the Gurney flap began to have side effects on the aerodynamic characteristics of the airfoil. The higher the Gurney flap is, the larger the wake velocity deficit is, so it indicates a higher mean drag.

**Figure 17.** The wake velocity measured by the wake rake array (uniform inflow, Re: 0.8 <sup>×</sup> 106). (**a**) α = 8.4◦; (**b**) α = 11.4◦.

By comparing the flows in the near trailing edge region of the airfoil with Gurney flap. The resultant velocity distribution of the baseline airfoil and the airfoil attached with Gurney flap (9 mm in height and 4.5 mm in thickness) at the angle of attack 8.4◦ (before the stall angle) and 11.4◦ (after the stall angle) under uniform inflow were depicted in Figure 18. The vertical direction Y denoted the forward distance from the hot-wire probe to the trailing edge, and the horizontal direction X denoted the lateral distance from the hot-wire probe to the trailing edge. As shown in Figure 18, the black solid spots were the actual measurement positions by using the electronically controlled traverser. As shown, the wake flows shift to the pressure side in both the two situations. When the angle of attack was 8.4◦, the airfoil with Gurney flap had a larger velocity deficit than the baseline airfoil, which suggests an increased drag. When the angle of attack was 11.4◦, the velocity deficit of the baseline airfoil was also less than the Gurney flap, but there was no obvious difference between the offset direction caused by the flapped airfoil and the baseline airfoil which were being similar to the far field results measured by the wake rake array.

**Figure 18.** Velocity distribution of the region behind trailing edge (uniform inflow, Re: 0.8 <sup>×</sup> 106). (**a**) α = 8.4◦, baseline airfoil; (**b**) α = 8.4◦, Gurney flap 9 mm × 4.5 mm; (**c**) α = 11.4◦, baseline airfoil; (**d**) α = 11.4◦, Gurney flap 9 mm × 4.5 mm.

Figure 19 presents the power density spectrum measured by the region 1 cm directly behind the trailing edge under uniform inflow. As can be observed, at the angle of attack 8.4◦, both along-wind and across-wind display an obvious peak in the power density spectrum, meanwhile the across-wind energy increased in the flapped airfoil. It indicates there are counter rotating vortices generated by the Gurney flap, while the wake of baseline airfoil does not show the existence of vortex shedding [33]. Gurney flaps also have a wake structure similar to the Karman vortex street behind the cylinder at bigger Reynolds number (i.e., of 0.8 <sup>×</sup> 106), the flow oscillations were caused by the separating of the shear layer on the pressure side and suction side of the trailing edge. The Karman vortex street both can be found in uniform flow and turbulent flow, so the increased lift produced by the Gurney flap seems can be explained in turbulent flow condition, either. Moreover, from the lift coefficients shown above, the size of structures was affected by turbulence level. Besides, as compared with the reference data under similar condition [33], with the increase of the Reynolds number, from Re = 1 <sup>×</sup> 105 to 0.8 <sup>×</sup> 106, the Strouhal number decreased from 0.151 (Re: 1 <sup>×</sup> <sup>10</sup>5, 2% chord length) to 0.128 (Re: 2 <sup>×</sup> <sup>10</sup>5, 2% chord length) and to 0.088 (Re: 0.8 <sup>×</sup> 106, 1.5% chord length).

**Figure 19.** Comparison of the power density spectrum (uniform inflow, Re: 0.8 <sup>×</sup> 106). (**a**) <sup>α</sup> = 8.4◦, baseline airfoil; (**b**) α = 8.4◦, Gurney flap 9 mm × 4.5 mm.

#### **4. Conclusions**

In summary, this paper presented a study of the aerodynamic performance of airfoils with/without the Gurney flap under different turbulence conditions by means of wind tunnel experiments. The experimental observations were as follows:

The Gurney flap deflects the wake position from the pressure side of the airfoil, increasing the vertical distance between the airfoil chord and the middle arc. The surface pressure characteristics revealed that the pressure difference of the trailing edge was larger for the flapped airfoil in the turbulence condition. Moreover, due to the special wake structure, the Gurney flap also indicates the vortex generated at the trailing edge, which results in a greater lift of the airfoil. The lift increment obtained from Gurney flap was influenced by the turbulence intensity, meanwhile a large increase in drag was observed as well. Under a medium range of turbulence intensity (i.e., of 10.5%), the maximum lift-to-drag ratios increased from 14.35% (flap 6 mm × 4.5 mm, a = 12.4◦) to 14.47% (flap 6 mm × 1.5 mm, a = 12.4◦), the 1.0% Gurney flap case has the best performance, that's probably because the decrease of boundary layer thickness under turbulent conditions, so the height of 1.0% was optimal. However, under a much larger turbulence inflow (i.e., of 19.0%), all experimental data have shown that the impact was negligibly small. The drag increment was found to surpass the benefit of the increase in lift, resulting in a worse lift-to-drag ratio than the baseline airfoil.

It was of importance to note that there are two kinds of measurement uncertainties could affect the acquired results. One was stochastic uncertainties caused by the discrete pressure data which can be eliminated by multiple measurements. The other one was systematic uncertainties caused by measuring instruments which can be reduced by calibration of the instruments before the experiment.

Based on the observations from the present measurements, it was found that under very high turbulence level, e.g., with 19% of turbulence intensity, the Gurney flap has very small influence on the airfoil aerodynamic performance. The measurement results show the Gurney flap may be used in areas with slightly lower turbulence. Nevertheless, wind turbines in real life were running in the turbulent atmospheric boundary layer, as shown above, the effects of Gurney flap are highly dependent on turbulent inflow, it seems that the Gurney flap would be installed as a scalable form and preferred to have smaller size, such that it would get more potential benefits at all periods of rotation. In addition, it should be noted that the incoming flow is often a combination of rotational wake from the previous wind turbines. Research on the lift increment of Gurney flap was not isolated, it always interacted with such highly unsteady aerodynamic issues. Therefore, to analyze more precisely the effects of combining such turbulence levels with Gurney flap, further experiments should be carried out for the rotor blades in highly turbulent flows. This was in progress on a wind turbine blade composed of several DTU-LN221 wing sections.

**Author Contributions:** Conceptualization, J.Y.; conducted the data collection, J.Y. and Y.Y.; funding acquisition, H.Y.; supervision, N.L., H.Y. and W.Z.; formal analysis, J.Y., H.Y., W.Z.; writing—original draft preparation, J.Y.; writing—review and editing, H.Y., W.Z.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Postgraduate Research & Practice Innovation Program of Jiangsu Province under grant number KYCX19\_2104 and Youth Fund Project of Jiangsu Natural Science Foundation under grant number BK20170510.

**Acknowledgments:** The authors wish to express acknowledgement to the Postgraduate Research &Practice Innovation Program of Jiangsu Province under grant number KYCX19\_2104 and Youth Fund Project of Jiangsu Natural Science Foundation under grant number BK20170510.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Evaluation of Tip Loss Corrections to AD**/**NS Simulations of Wind Turbine Aerodynamic Performance**

#### **Wei Zhong 1, Tong Guang Wang 1, Wei Jun Zhu 2,\* and Wen Zhong Shen <sup>3</sup>**


Received: 17 October 2019; Accepted: 12 November 2019; Published: 15 November 2019

**Abstract:** The Actuator Disc/Navier-Stokes (AD/NS) method has played a significant role in wind farm simulations. It is based on the assumption that the flow is azimuthally uniform in the rotor plane, and thus, requires a tip loss correction to take into account the effect of a finite number of blades. All existing tip loss corrections were originally proposed for the Blade-Element Momentum Theory (BEMT), and their implementations have to be changed when transplanted into the AD/NS method. The special focus of the present study is to investigate the performance of tip loss corrections combined in the AD/NS method. The study is conducted by using an axisymmetric AD/NS solver to simulate the flow past the experimental *NREL Phase VI* wind turbine and the virtual *NREL 5MW* wind turbine. Three different implementations of the widely used Glauert tip loss function *F* are discussed and evaluated. In addition, a newly developed tip loss correction is applied and compared with the above implementations. For both the small and large rotors under investigation, the three different implementations show a certain degree of difference to each other, although the relative difference in blade loads is generally no more than 4%. Their performance is roughly consistent with the standard Glauert correction employed in the BEMT, but they all tend to make the blade tip loads over-predicted. As an alternative method, the new tip loss correction shows superior performance in various flow conditions. A further investigation into the flow around and behind the rotors indicates that tip loss correction has a significant influence on the velocity development in the wake.

**Keywords:** wind turbine aerodynamics; actuator disc; AD/NS; tip loss correction; blade element momentum

#### **1. Introduction**

Wind energy is nowadays an important and increasing source of electric power. It has been the biggest contributor of renewable electricity except for hydropower, sharing about 5.5% of the global electricity production in 2018 [1]. As a primary subject of wind turbine technology, aerodynamics [2] largely determines the efficiency of wind energy extraction of an individual wind turbine or a wind farm. Along with the extensive development of high quality wind resources onshore, there are several trends in wind power industry: A large number of newly installed wind turbines have to be installed in areas with lower wind speeds and complex terrain; the layout optimization of wind turbine array becomes very important for wind farm design; offshore wind power development is accelerating, and the rotor size is continuously increasing. These trends need to be supported by more advanced and refined aerodynamic tools. More accurate aerodynamic load prediction is required for the design of a new generation of wind turbines with high efficiency and relatively low weight. Furthermore, a wind farm with dozens of wind turbines needs to be studied as a whole in order to take into account the complex interference between wind turbines.

Blade Element Momentum Theory (BEMT) [3,4] and Computational Fluid Dynamics (CFD) [5,6] are two essential methods of wind turbine aerodynamic computation. BEMT is no doubt the key method for rotor design [7], while CFD gives the most refined data of aerodynamic loads and flow parameters. A full CFD simulation with resolved rotor geometry is usually employed for an individual wind turbine [8,9]. However, full CFD is not suitable for a wind turbine array, due to the huge computational cost. The Actuator Disc/Navier-Stokes (AD/NS) method [10–13] was developed for this situation, in which the rotor geometry is not resolved, and thus, the number of mesh cells is greatly reduced. AD/NS is based on a combination of the blade-element theory and CFD. The flow field is still solved by CFD, while the rotor entity is replaced by a virtual actuator disc on which an external body force is acted. If the blade entities are represented by virtual actuator lines, it is called the Actuator Lines/Navier-Stokes (AL/NS) method [14–17]. As compared to AD/NS, the flow around the rotor solved by AL/NS is closer to the reality, since the rotor vortices are simulated. However, AL/NS requires more grid cells to describe the actuator lines, and thus, is seldom employed in wind farm simulations.

The AD/NS method has played a significant role in wind farm simulations involving wake interaction [18,19], complex terrain [20,21], atmospheric boundary layer [22,23], noise propagation [24,25], etc. Nevertheless, the AD model assumes that the number of blades is infinite, and thus, no tip loss is simulated, causing an over-prediction of the blade tip loads and power output. In order to take tip loss into account, a reliable engineering model has to be embedded into the numerical solver, which is usually called tip loss correction. However, all tip loss corrections were originally proposed for BEMT, and there is no tip loss correction specially developed for AD/NS. Most of the literature about AD/NS simulations either did not mention tip loss or declared that the Glauert tip loss correction [26] was employed. It is worth mentioning that AD/NS and BEMT solve the momentum of the flow past the rotor by using two completely different approaches. The former employs CFD, while the latter applies the momentum theory, though they share the actuator disc assumption and the blade-element theory. That leads to different implementations of tip loss correction for the two methods. In the BEMT, tip loss correction is realized by applying a correction factor *F* into the induced velocity through the rotor. In the AD/NS, the velocity is naturally found by the NS solver, and factor *F* can only be used to modify the external body force. The question is whether a tip loss correction has a good global performance when it is transplanted into AD/NS. Even for the BEMT itself, evaluations indicate that accurate tip loss correction is still not well achieved [27], and the development of new correction models is still going on [28–33]. In contrast with the massive study in the BEMT, tip loss corrections applied to AD/NS lack a comprehensive evaluation.

In the present work, the 2-Dimensional (2D) axisymmetric AD model is employed, and steady-state simulations are performed. The code solves the incompressible axisymmetric NS equations for the experimental *NREL Phase VI* wind turbine [34] and the virtual *NREL 5MW* [35] wind turbine under various axial inflow conditions. The main purpose of the numerical study is to evaluate the performance of the tip loss corrections applied to AD/NS. Three different implementations of the Glauert tip loss correction are discussed. In addition, a new tip loss correction recently proposed by Zhong et al. [33] is introduced and compared. The normal and tangential forces of blade cross-sections are chosen as the key parameters for the present evaluation. Because of the long arm of force, computational errors of the forces acted on the tip region are most likely to cause non-ignorable errors of the blade bending moment and the rotor torque (power generation). A BEMT study on the *NREL Phase VI* wind turbine by Branlard [27] showed that the power generation would be overestimated by 15% if no tip loss correction was made and a deviation of about 5% exists when various existing tip loss corrections were applied. That highlights the significance of the present study on tip loss correction.

The innovation of the present study lies in the following items. (a) Different implementations of the Glauert tip loss factor *F* are gathered, discussed and compared, and finally, one of them is recommended according to its best performance. Such kind of work has never been reported in the existing literature. (b) The Glauert correction is found to have a similar performance when it is transplanted from BEMT to AD/NS. (c) The new tip loss correction of Zhong et al. [33] is for the first time employed in AD/NS simulations. It is found to be generally superior to the Glauert-type corrections. That provides an alternative choice for a more accurate prediction of the blade tip loads. (d) Tip loss correction is found to have a significant influence on the velocity field, which highlights the importance of an accurate tip loss correction for not only the blade loads, but also the wake development.

The paper is organized as follows: In Section 2, the axisymmetric AD/NS method is introduced, including the governing equations of the NS approach and the formulae of the AD model; In Section 3, the Glauert and the new tip loss corrections are introduced; In Section 4, the implementations of the tip loss corrections used in AD/NS are described; In Section 5, the numerical setup, as well as the involved simulation cases, are introduced; In Section 6, all simulation results are presented and discussed; Final conclusions are made in the last section.

#### **2. Axisymmetric AD**/**NS Method**

#### *2.1. Governing Equations*

The governing equations of the present AD/NS simulation are the incompressible axisymmetric NS equations. In a cylindrical coordinate system as shown in Figure 1, the axial direction is defined along the *z*-axis, the radial direction is represented by *r*, and the tangential direction is represented by θ, the continuity equation is written as

$$
\nabla \cdot \overrightarrow{u} = \frac{\partial u\_z}{\partial z} + \frac{\partial u\_r}{\partial r} + \frac{u\_r}{r} = 0,\tag{1}
$$

the axial and radial momentum equations read

$$\frac{\partial u\_z}{\partial t} + \frac{1}{r} \frac{\partial}{\partial z} (r u\_z^2) + \frac{1}{r} \frac{\partial}{\partial r} (r u\_l u\_z) = -\frac{1}{\rho} \frac{\partial p}{\partial z} + \frac{1}{r} \frac{\partial}{\partial z} \left[ r \nu \left( 2 \frac{\partial u\_z}{\partial z} - \frac{2}{3} (\nabla \cdot \overrightarrow{u}) \right) \right] + \frac{1}{r} \frac{\partial}{\partial r} \left[ r \nu \left( \frac{\partial u\_z}{\partial r} + \frac{\partial u\_r}{\partial z} \right) \right] + f\_z,\tag{2}$$

$$\begin{split} \frac{\partial u\_r}{\partial t} &+ \frac{1}{r} \frac{\partial}{\partial z} (r u\_z u\_r) + \frac{1}{r} \frac{\partial}{\partial r} (r u\_r^2) = \\ &- \frac{1}{\rho} \frac{\partial p}{\partial r} + \frac{1}{r} \frac{\partial}{\partial z} \Big[ r \nu \left( \frac{\partial u\_r}{\partial z} + \frac{\partial u\_z}{\partial r} \right) \Big] + \frac{1}{r} \frac{\partial}{\partial r} \Big[ r \nu \left( 2 \frac{\partial u\_r}{\partial r} - \frac{2}{3} (\nabla \cdot \vec{u'}) \right) \Big] - 2 \nu \frac{u\_r}{r^2} + \frac{2}{3} \frac{\nu}{r} \Big( \nabla \cdot \vec{u'} \Big) + \frac{u\_\theta^2}{r} + f\_r \end{split} \tag{3}$$

and the tangential momentum equation is

$$\frac{\partial \boldsymbol{u}\_{\boldsymbol{\theta}}}{\partial t} + \frac{1}{r} \frac{\partial}{\partial z} (r \boldsymbol{u}\_{z} \boldsymbol{u}\_{\boldsymbol{\theta}}) + \frac{1}{r} \frac{\partial}{\partial r} (r \boldsymbol{u}\_{r} \boldsymbol{u}\_{\boldsymbol{\theta}}) = \frac{1}{r} \frac{\partial}{\partial z} [r \boldsymbol{v} \frac{\partial \boldsymbol{u}\_{\boldsymbol{\theta}}}{\partial z}] + \frac{1}{r^{2}} \frac{\partial}{\partial r} [r^{3} \boldsymbol{v} \frac{\partial}{\partial r} (\frac{\boldsymbol{u}\_{\boldsymbol{\theta}}}{r})] - \frac{\boldsymbol{u}\_{r} \boldsymbol{u}\_{\boldsymbol{\theta}}}{r} + f\_{\boldsymbol{\theta}}.\tag{4}$$

**Figure 1.** Definition of coordinates and the stream tube through an actuator disc.

In the above equations, *uz*/*ur*/*u*<sup>θ</sup> is the axial/radial/tangential velocity, *p* is the static pressure, ρ is the air density, ν is the kinematic viscosity coefficient of air, *fz*, *fr*, *f*<sup>θ</sup> are the source term of the axial, radial, tangential external body forces, respectively.

#### *2.2. Force on Actuator Disc*

The conceptual idea of the AD/NS method is to solve the aerodynamic force of rotor blades by using the blade-element theory and then applying its counterforce as an external body force into the momentum equations.

At each time step of solving the NS equations, the velocity passing through the actuator disc is detected and used to determine the axial induced velocity *Wa* and the tangential induced velocity *Wt*,

$$\mathcal{W}\_a = V\_0 - u\_{z\prime} \tag{5}$$

$$\mathcal{W}\_t = -\boldsymbol{\mu}\_0. \tag{6}$$

The axial interference factor *a* and tangential interference factor *a* are then determined by,

$$a = \frac{W\_a}{V\_0} \,\prime \tag{7}$$

$$a' = \frac{\mathcal{W}\_{\text{f}}}{\Omega r'} \tag{8}$$

where *V*<sup>0</sup> is the wind speed, and Ω is the rotating speed of the wind turbine rotor.

The above interference factors are azimuthally unchanged according to the axisymmetric condition, which implies an infinite number of blades and zero tip loss. In order to estimate the tip loss of the realistic rotor with a finite number of blades, the interference factors need to be corrected as

$$d = f\_{corr}(a),$$

$$
\vec{a}' = f'\_{\
avr}(a'),
\tag{10}
$$

where *a*˜ and *a*˜ denote the corrected axial and tangential interference factors, *fcorr* and *f corr* represent correction functions.

The axial and tangential induced velocities after the correction are written as

$$
\tilde{\mathcal{W}}\_{\mathfrak{a}} = V\_0 \tilde{a}\_{\mathfrak{a}} \tag{11}
$$

$$
\vec{W}\_t = \Omega \,\hat{n}'.\tag{12}
$$

According to the velocity triangle shown in Figure 2, the flow angle φ, angle of attack α, and relative velocity *Vrel* are then determined by

$$\phi = \tan^{-1} \left( \frac{V\_0 - \bar{W}\_a}{\Omega r + \bar{W}\_t} \right) \tag{13}$$

$$
\alpha = \phi - \beta,\tag{14}
$$

$$V\_{rel}^2 = \left(V\_0 - \vec{W}\_a\right)^2 + \left(\Omega r + \vec{W}\_t\right)^2,\tag{15}$$

where β is the local pitch angle of a blade cross-section.

**Figure 2.** Illustration of the velocity and the aerodynamic force for a blade cross-section.

With the determined angle of attack and relative velocity, the lift coefficient *Cl* and drag coefficient *Cd* of each cross-section of the blade can be obtained from tabulated airfoil data. According to the relationship of force projection, the normal and tangential force coefficients are determined by

$$\mathbf{C}\_{\text{H}} = \mathbf{C}\_{\text{I}} \cos \phi + \mathbf{C}\_{\text{d}} \sin \phi,\tag{16}$$

$$C\_t = C\_l \sin \phi - C\_d \cos \phi\_\prime \tag{17}$$

The axial force *Fz* and tangential force *F*<sup>θ</sup> per radial length of the rotor are then given by

$$F\_z = \frac{1}{2}\rho V\_{\text{ref}}^2 cB \mathbb{C}\_n = \frac{1}{2}\rho V\_{\text{ref}}^2 cB (\mathbb{C}\_l \cos\phi + \mathbb{C}\_d \sin\phi),\tag{18}$$

$$F\_0 = \frac{1}{2}\rho V\_{rel}^2 c B \mathbb{C}\_l = \frac{1}{2}\rho V\_{rel}^2 c B (\mathbb{C}\_l \sin\phi - \mathbb{C}\_d \cos\phi),\tag{19}$$

where *c* is the local chord length, and *B* is the number of blades. Ignoring the effect of the radial flow, the vector form of the counterforce acting on the air is

$$
\overrightarrow{F} = (-F\_{z\prime}0\_{\prime} - F\_{\partial}).\tag{20}
$$

Rather than distributing the force only on the disc, the above force is regularized by the following Gaussian distribution along the axial direction in order to avoid a numerical singularity. The Gaussian distribution has been widely used and proven to be proper for AD/NS and AL/NS simulations [10,14–17].

$$\stackrel{\rightarrow}{F}\_{\varepsilon} = \eta\_{\varepsilon} \stackrel{\rightarrow}{F}\_{\prime} \eta\_{\varepsilon} = \frac{1}{\varepsilon \sqrt{\pi}} \exp\left[-\left(\frac{z - z\_{0}}{\varepsilon}\right)^{2}\right] \tag{21}$$

where *z*<sup>0</sup> is the axial position of the disc. The parameter ε serves to adjust the concentration of the regularized force, which is in the present study set to ε = 0.02*R* where *R* is the rotor radius.

The source terms of the external body force in the momentum Equations (2)–(4) need to be replaced with the above force per unit volume. In the axisymmetric coordinate system, a 2D grid cell with an axial length of Δ*z* and a radial height of Δ*r* represents a volume of 2π*r*Δ*r*Δ*z*, and thus, the force per unit volume is

$$
\overrightarrow{f} = \frac{\overrightarrow{F}\_{\varepsilon} \Delta r \Delta z}{2\pi r \Delta r \Delta z} = \frac{\overrightarrow{F}\_{\varepsilon}}{2\pi r} \tag{22}
$$

Using Equations (20) and (21), it reads

$$\overrightarrow{f} = \frac{\eta\_{\epsilon}\overrightarrow{F}}{2\pi r} = \left(-\frac{\eta\_{\epsilon}F\_{z}}{2\pi r}, 0, -\frac{\eta\_{\epsilon}F\_{\theta}}{2\pi r}\right) \tag{23}$$

i.e.,

$$\begin{cases} \ f\_z = -\frac{\eta\_e F\_z}{2\pi r} \\ \ f\_r = 0 \\ \ f\_\partial = -\frac{\eta\_e F\_\partial}{2\pi r} \end{cases} \tag{24}$$

#### **3. Tip Loss Corrections**

#### *3.1. Glauert Correction*

The tip loss correction of Glauert [26] was developed from the study of Prandtl [36]. A function *F*, which was later recognized as the first tip loss factor, was derived by Prandtl making Betz's optimal circulation [37] go to zero at the blade tip,

$$F = \frac{2}{\pi} \arccos\{ \exp\left[ -\frac{B}{2} \left( 1 - \frac{r}{R} \right) \sqrt{1 + \lambda^2} \right] \}\tag{25}$$

where *B* is the number of blades, *r* is the local radial location, *R* is the rotor radius, and λ is the tip speed ratio. The original derivation of Prandtl was written very briefly and was elaborated more in [4,27,38] as for a more detailed derivation.

Later on, Glauert made further contributions. First, he interpreted the physical meaning of factor *F* as the ratio between the azimuthally averaged induced velocity and the induced velocity at the blade position, leading to

$$F = \frac{\overline{a}}{a\_B} = \frac{\overline{a'}}{a'\_B},\tag{26}$$

where *a* = <sup>1</sup> 2π <sup>2</sup><sup>π</sup> <sup>0</sup> *a*dθ and *a* = <sup>1</sup> 2π <sup>2</sup><sup>π</sup> <sup>0</sup> *a* dθ are the azimuthally averaged axial and tangential interference factors, and *aB* and *a <sup>B</sup>* are the interference factors at the blade position. Second, the local inflow angle φ was introduced into the function in order to make it consistent with the local treatment of the BEMT, leading to the following new formula of factor *F*,

$$F = \frac{2}{\pi} \arccos\left\{ \exp\left[ -\frac{B(R-r)}{2r\sin\phi} \right] \right\}.\tag{27}$$

The factor *F* was then applied into the momentum theory through a straightforward way of multiplying the axial velocity at the rotor plane with factor *F*, resulting in the following thrust and torque for an annular element,

$$\mathrm{d}T = 4\pi r \rho V\_0^2 a (1 - a) F \mathrm{d}r,\tag{28}$$

$$\text{rdd} = 4\pi r^3 \rho V\_0 \Omega a' (1 - a) F \text{d}r. \tag{29}$$

Using another two equations of d*T* and d*M* derived from the blade-element theory,

$$\mathbf{d}T = \frac{1}{2}\rho V\_{rel}^2 \mathbf{c} \mathbf{C}\_n \mathcal{B} \mathbf{d}r = \frac{1}{2}\rho V\_{rel}^2 \mathbf{c} (\mathbf{C}\_I \cos\phi + \mathbf{C}\_d \sin\phi)\mathcal{B}\mathbf{d}r,\tag{30}$$

$$\text{d}M = \frac{1}{2}\rho V\_{\text{rel}}^2 \text{c} \mathbf{C}\_l \text{Br} \mathbf{d}r = \frac{1}{2}\rho V\_{\text{rel}}^2 \mathbf{c} (\mathbf{C}\_l \sin\phi - \mathbf{C}\_d \cos\phi) \text{Brd}r,\tag{31}$$

and the velocity triangle at the rotor plane, the equations for the interference factors in the BEMT approach was derived to be:

$$\frac{a}{1-a} = \frac{\sigma(\mathcal{C}\_l \cos \phi + \mathcal{C}\_d \sin \phi)}{4F \sin^2 \phi},\tag{32}$$

$$\frac{a'}{1+a'} = \frac{\sigma \left(\mathbb{C}\_l \sin \phi - \mathbb{C}\_d \cos \phi\right)}{4F \sin \phi \cos \phi}. \tag{33}$$

The above two equations lead to the final iterative formulae of the BEMT approach with the Glauert tip loss correction. Their difference from the baseline BEMT approach is only the appearance of factor *F* in the equations. Obviously, the application of the Glauert correction in BEMT is very simple, which is also a great advantage. There are several variations of the Glauert correction in which the way of applying factor *F* is different [39,40], but the original Glauert tip loss correction is the most commonly used form till today.

#### *3.2. A Newly Developed Correction*

A new tip loss correction was recently proposed for BEMT by Zhong et al. [33], based on a novel insight into tip loss. In contrast with the Prandtl/Glauert series corrections that begin with an actuator disc and estimate the effect of the finite number of blades, the new correction begins with a non-rotating blade and estimates the effect of rotation on tip loss. It has been validated in BEMT computations and showed superior performances in the cases involving flow separation or high axial interference factor.

The correction was realized by using two factors of *FR* and *FS* that treat the rotational effect and the 3D effect, respectively.

$$F\_R = 2 - \frac{2}{\pi} \arccos\left\{ \exp\left[ -2B(1 - r/R)\sqrt{1 + \lambda^2} \right] \right\},\tag{34}$$

$$F\_S = \frac{2}{\pi} \arccos\left\{ \exp\left[ -\left(\frac{1 - r/R}{\overline{c}/R}\right)^{3/4} \right] \right\}, \ \overline{c} = \frac{S\_\mathbf{t}}{R - r},\tag{35}$$

where *c* denotes a geometric mean chord length, and *S*t is the projected area of the blade between the present cross-section and the tip. The purpose of introducing *c* is to deal with tapered blades and those with sharp tips.

The factor *FR* was applied to the BEMT approach by using:

$$\frac{aF\_R(1-aF\_R)}{(1-a)} = \frac{\sigma \Big(\mathcal{C}\_l \cos\phi + \mathcal{C}\_d \sin\phi\Big)}{4\sin^2\phi},\tag{36}$$

$$\frac{a'F\_R(1-aF\_R)}{(1+a')(1-a)} = \frac{\sigma(\mathbb{C}\_l \sin \phi - \mathbb{C}\_d \cos \phi)}{4 \sin \phi \cos \phi},\tag{37}$$

in which the employed lift and drag coefficients were corrected by factor *FS*, rather than the direct use of the airfoil data.

$$\bar{\mathbf{C}}\_{l} = \frac{1}{2} \left[ \mathbf{C}\_{l}(a) F\_{S} + \mathbf{C}\_{l}^{\*} \right] \tag{38}$$

$$\bar{\mathbf{C}}\_{d} = \frac{1}{\cos^{2}a\_{l}} \left[ \mathbf{C}\_{d}(a\_{\ell}) \cos a\_{l} + \left( \bar{\mathbf{C}}\_{l} \cos a\_{l} + \mathbf{C}\_{d}(a\_{\ell}) \tan a\_{l} \right) \sin a\_{l} \right] \tag{39}$$

where

$$\mathbf{C}\_{l}^{\*} = \frac{1}{\cos^{2}\alpha\_{l}} [\mathbf{C}\_{l}(\alpha\_{\ell}) \cos \alpha\_{l} - \mathbf{C}\_{d}(\alpha\_{\ell}) \sin \alpha\_{l}],\tag{40}$$

$$a\_l = \frac{C\_l(\alpha)}{m}(1 - F\_S)\_l \tag{41}$$

where *m* is the slope of the lift-curve of the airfoil before flow separation, α*<sup>i</sup>* is called the downwash angle, and α*<sup>e</sup>* is called the effective angle of attack.

By comparing with the Glauert tip loss correction, the new tip loss correction appears more complicated in application. It is simplified in the present study: First, Equations (36) and (37) will not be employed naturally in the AD/NS method; Second, Equations (38)–(41) are further simplified to

$$\mathcal{C}\_{l} = \frac{1}{2} [\mathbb{C}\_{l}(\alpha) F\_{\mathbb{S}} + \mathbb{C}\_{l}(\alpha\_{\mathfrak{e}})], \ a\_{\mathfrak{e}} = \alpha - a\_{i\mathfrak{e}} \tag{42}$$

$$
\bar{\mathbf{C}}\_d = \mathbf{C}\_d(\alpha\_\ell) + \bar{\mathbf{C}}\_l \tan \alpha\_{i\nu} \ a\_i = \frac{\mathbf{C}\_l(\alpha)}{m} (1 - F\_\mathcal{S}). \tag{43}
$$

The simplification is derived by using the fact that α*<sup>i</sup>* and *Cd* are relatively small. More detailed applications are shown in the next section.

#### **4. Applying Corrections to AD**/**NS Simulation**

#### *4.1. Application of Glauert Tip Loss Factor F*

Equations (32) and (33) where the Glauert tip loss factor *F* is applied are not employed in the AD/NS method because the flow is simulated by the NS solver instead of the momentum theory. As a result, an alternative way has to be employed correctly to apply the tip loss factor *F*. Nevertheless, among the literature studies, little literature mentions the detail of how the tip loss factor is applied to the AD/NS simulations. After an extensive literature review, we have found three representative documents in which Sørensen et al. [41], Mikkelsen [42], and Shen et al. [43] described their implementations. In order to facilitate the distinction, the implementations adopted by the three studies are denoted as Glauert-A, Glauert-B and Glauert-C in the present paper, respectively.

#### 4.1.1. Glauert-A Correction

Sørensen et al. [41] performed a tip loss correction by applying factor *F* to modify the aerodynamic force that determines the external body force in the NS equations. They replaced the lift coefficient *Cl* by *Cl*/*F* (there was no need to modify the drag in their study as the drag was assumed not to produce the external body force). In the present Glauert-A correction, *Cn* and *Ct*, in Equations (18) and (19), are replaced by:

$$
\tilde{\mathbb{C}}\_n = \mathbb{C}\_n / \mathbb{F}\_\prime \tag{44}
$$

$$
\vec{\mathsf{C}}\_{l} = \mathsf{C}\_{l} / \mathsf{F}.\tag{45}
$$

That is equivalent to replacing *Cl* by *Cl*/*F* and *Cd* by *Cd*/*F*, according to Equations (16) and (17).

Sørensen et al. [41] did not explain the reason why factor *F* could be directly used to modify the force. They might be inspired by Equations (32) and (33) in which the existence of *F* can be looked as corrections to *Cn* and *Ct*, although from a physical point of view it is a correction to the interference factors. It is important to note that the corrected force is only used for determining the external body force, so as the flowfield, while the force acting on the blade should be regarded as the original one. In addition, the interference factors in Equations (9) and (10) should no longer be corrected, leading to *a*˜ = *a* and *a*˜ = *a* .

This implementation involves a division operation with denominator *F*, which causes unreasonable big values of *<sup>C</sup>*˜*<sup>n</sup>* and *<sup>C</sup>*˜*<sup>t</sup>* at the extreme tip where *<sup>F</sup>* <sup>→</sup> 0. However, there is no exact criterion for defining what a big value is unreasonable because the *C*˜*<sup>n</sup>* and *C*˜*<sup>t</sup>* themselves are introduced as an engineering correction rather than a physical concept. Sørensen et al. [41] did not mention this problem, and no obvious numerical fluctuation is observed in his result. That is possible because the problem is limited to a very small area at the tip, and thus, the integral effect of the resulted unreasonable body force is also small.

#### 4.1.2. Glauert-B Correction

Mikkelsen [42] compared Equations (32) and (33) to the corresponding baseline equations without factor *F*. The baseline equations are

$$\frac{a}{1-a} = \frac{\sigma(\mathcal{C}\_I \cos \phi + \mathcal{C}\_d \sin \phi)}{4 \sin^2 \phi},\tag{46}$$

$$\frac{a'}{1+a'} = \frac{\sigma \left(\mathcal{C}\_l \sin \phi - \mathcal{C}\_d \cos \phi\right)}{4 \sin \phi \cos \phi}. \tag{47}$$

In order to distinguish from the parameters in the baseline equations, we here rewrite Equations (32) and (33) to

$$\frac{F\vec{a}}{(1-\vec{a})} = \frac{\sigma \left(\mathbb{C}\_{l}\cos\vec{\phi} + \mathbb{C}\_{d}\sin\vec{\phi}\right)}{4\sin^{2}\vec{\phi}},\tag{48}$$

$$\frac{F\vec{a}^{\prime}}{(1+\vec{a}^{\prime})} = \frac{\sigma \big(\mathbb{C}\_{l}\sin\vec{\phi} - \mathbb{C}\_{d}\cos\vec{\phi}\big)}{4\sin\vec{\phi}\cos\vec{\phi}}.\tag{49}$$

The implementation of Mikkelsen [42] adopts the following equations which implies an assumption of φ = φ˜ which in fact does not accurately hold,

$$\frac{a}{1-a} = \frac{\sigma(\mathcal{C}\_l \cos \phi + \mathcal{C}\_d \sin \phi)}{4 \sin^2 \phi} = \frac{F\overline{a}}{(1-\overline{a})'} \tag{50}$$

$$\frac{a'}{1+a'} = \frac{\sigma \left(\mathbb{C}\_l \sin \phi - \mathbb{C}\_d \cos \phi\right)}{4 \sin \phi \cos \phi} = \frac{F \vec{a}'}{(1+\vec{a}')}.\tag{51}$$

The corrected interference factors were, thus, determined by

$$\vec{a} = \frac{a}{F(1-a) + a'} \tag{52}$$

$$\overline{a}' = \frac{a'}{F(1+a') - a'}.\tag{53}$$

The tip loss correction was completed as the above functions for *a*˜ and *a*˜ were used to replace Equations (9) and (10).

#### 4.1.3. Glauert-C Correction

Shen et al. [43] performed their correction by directly using Equation (26) which represents the physical meaning of factor *F*. Using the azimuthally uniform condition (*a* = *a*), the correction was completed by replacing Equations (9) and (10) with

$$
\vec{a} = a\_B = a/F,\tag{54}
$$

$$
\overline{a}' = a'\_B = a'/\mathcal{F}.\tag{55}
$$

Similar to the Glauert-A correction, this implementation involves a division operation with denominator *F*. That causes an unphysical big value of *a*˜ at the extreme tip where *F* → 0. A limiter is set in the present study to force the result to be *a*˜ = 1 when it is greater than 1. The value of *a*˜ is usually much less than *a*˜ and is not limited here.

#### *4.2. Application of New Correction*

The application of the new correction consists of two steps: The first is using factor *FR* determined by Equation (34) and the second is using factor *FS* determined by Equation (35).

In the first step, an implementation similar to the Glauert-C correction is adopted, which is performed by replacing Equations (9) and (10) with

$$
\vec{a} = a/F\_{R\nu} \tag{56}
$$

$$
\bar{a}' = a'/F\_{\mathbb{R}}.\tag{57}
$$

A limiter of *a*˜ = 1 is used when *a* is greater than 1. The angle of attack α can then be determined by calculations using Equations (11)–(14).

The second step is to correct the lift and drag coefficients by using Equations (42) and (43), the result of which is used to replace the *Cl* and *Cd* in Equations (18) and (19).

#### **5. Computational Setup**

#### *5.1. Flow Solver and Mesh Configuration*

The 2D axisymmetric NS equations are solved by using an in-house code EllipSys2D [44,45] developed at Technical University of Denmark (DTU). The code is a general incompressible flow solver with multi-block and multi-grid strategy. The equations are discretized with a second-order finite volume method. In the spatial discretization, a central difference scheme is applied to the diffusive terms and the QUICK (Quadratic Upstream Interpolation for Convective Kinematics) upwind scheme is applied to the convective terms. The SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations) scheme is used for the velocity-pressure decoupling. The turbulence flow is simulated using the method of Reynolds-averaged Navier-Stokes equations (RANS) in which the *k-*ω turbulence model of Menter [46] is employed with a modification for a better simulation of the turbulence quantities in the free-stream flow [47,48].

The coordinate for the AD/NS simulations is defined in the *z*-*y* plane. A computational mesh is generated, as shown in Figure 3 where the inflow, outflow and axisymmetric boundaries are indicated, the length of the computational domain is 30, which is nondimensionalized with the rotor radius. The blade is positioned at *z* = 0 and the grid cells are clustered around *z* = 0 and *y* = 1 to ensure a better resolution near the blade tip. The mesh is composed of four blocks where 64 × 64 grid points are used for each block. There are 64 grid points on the AD along the radial direction, and 20 points in the range of [-0.02, 0.02] in the axial direction where the Gaussian distribution of the body force plays a significant role. Since the axisymmetric boundary condition is applied, the 2D flow solution is regarded as an azimuthal slice of a full 3D field.

**Figure 3.** Mesh configuration for the Actuator Disc/Navier-Stokes (AD/NS) simulations.

The above computational setup has been proven to perform well in AD/NS simulations, as shown in the previous work of Cao et al. [49] where both blade loads and wake flows were validated against experiments.

#### *5.2. Simulation Cases*

Simulations are performed for two different wind turbines of the *NREL Phase VI* [34] and the *NREL 5MW* [35] that represent a small and a large rotor size, respectively. Additionally, the *NREL Phase VI* rotor has a very blunt blade tip shape as compared with the *NREL 5MW* rotor.

The operational conditions of the two rotors are listed in Table 1. Cases 1, 2 and 3 are set to be consistent with the *NREL UAE* experiment in axial inflow conditions [34]. Various wind speeds of 7 m/s, 10 m/s and 13 m/s are considered, while the rotating speed remains unchanged. The flow is fully attached on the blade surface at a wind speed of 7 m/s, but is partly separated at 10 m/s and 13 m/s (Higher wind speed leads to heavier flow separation) [8]. Measured pressure distributions (from which

the force can be obtained by pressure integral) at five blade sections (*r*/*R* = 0.30, 0.47, 0.63, 0.80, 0.95) are available for these cases [50]. Cases 4 and 5 are two typical points on the designed power curve of the *NREL 5MW* wind turbine [35]. The wind speed and rotating speed of Case 5 are the rated parameters of this wind turbine. Case 4 represents a condition with a lower wind speed of 8 m/s at which the rotating speed is reduced to 9.22 rpm for tracking the optimum tip speed ratio.


**Table 1.** Simulation cases for two different wind turbines.

These cases cover the following multiple situations: Wind turbines with remarkably different rotor sizes and tip shapes; flow with fully attached and separated conditions; operating conditions with various wind speeds and rotating speeds. It is, therefore, interesting to see the performance of the tip loss corrections in these cases.

Experimental data are preferred as the reference for the computational results of the *NREL Phase VI* rotor, as well as the full CFD data are also displayed in some results. For the cases of the *NREL 5MW* rotor, there is no experimental data available such that the full CFD data are the only reference. The related full CFD simulations were previously conducted by the authors, see Zhong et al. [8,33].

#### **6. Simulation Results**

#### *6.1. Results of Glauert-Type Corrections*

#### 6.1.1. Glauert-A/B/C Corrections

The blade loads for Cases 1, 2 and 3 are gathered in Figure 4, which represent results obtained from the *NREL Phase VI* rotor. The *NREL Phase VI* blade has a linear change of chord distribution starting from *r* = 1.3 m with a constant slope of 0.1. The chord length is 0.358 m at the tip, which is about 48% of the largest chord length in the blade inboard part. Therefore, the loads do not converge to zero without a tip loss correction.

**Figure 4.** *Cont*.

**Figure 4.** Normal force *Fn* and tangential force *Ft* along the *NREL Phase VI* blade.

It is noticed that there is a certain degree of difference between the results of the Glauert-A/B/C corrections in the tip region (e.g., at *r*/*R* > 0.8). In general, Glauert-C leads to lower loads which are closer to the experimental data, Glauert-A gives results close to, but slightly higher than, those of Glauert-C, while Glauert-B results in relatively higher loads near the tip. Taking the *Fn* in Figure 4a as an example of quantitative comparison, the relative difference between the results of Glauert-B and Glauert-C is about 4% at *r*/*R* = 0.95. The difference in other cases is not larger than this value.

At a wind speed of 7 m/s, all the curves of *Fn* generally agree with the five experimental data except that at *r*/*R* = 0.95 where an overestimation is observed. This overestimation becomes much more remarkable as the wind speed increases to 10 m/s and 13 m/s. At 10 m/s, there is no much difference around *r*/*R* = 0.9 between the results with no correction and with the Glauert-type corrections, indicating that almost no effective correction is made here. At 13 m/s, the curves of the Glauert-type corrections even exceed the uncorrected curve in a *r*/*R* range of about [0.66, 0.88]. It is clear that the corrections perform much worse at the higher wind speeds as compared with 7 m/s. The overestimation at *r*/*R* = 0.95 is also observed in the results for *Ft*, although it looks better to some extent, especially at 7 m/s.

Considering the fact that the Glauert correction was derived based on the potential hypothesis, it is reasonable to believe that the unusually poor performance at the higher wind speed is caused by the flow separation. According to the results of the *NREL UAE Phase VI* experiments [50], the angle of attack at most cross-sections of the blade exceeds the linear range of their aerodynamic polar when the wind speed is increased to be higher than 10m/s. Taking the cross-section of *r*/*R* = 0.63 as an example, the angles of attack are 5.9◦, 11.8◦, 16.9◦ at a wind speed of 7 m/s, 10 m/s, 13 m/s, respectively. The exact operating points for three cases on the lift polar of the S809 airfoil (the airfoil of all cross-sections

of the *NREL Phase VI* blade) are depicted in Figure 5. It clearly shows that Case 2 and Case 3 are in the nonlinear lift region corresponding to flow separation, while Case 1 is in the linear lift region corresponding to attached flow, which implies that flow separation occurs at 10 m/s and 13 m/s.

**Figure 5.** Operating points of the cross-section of *r*/*R* = 0.63 on the lift polar of the S809 airfoil. (The Reynolds number of the airfoil is 1 <sup>×</sup> <sup>10</sup><sup>6</sup> which is close to that of the cross-section).

The normal and tangential forces along a blade of the *NREL 5MW* rotor are plotted in Figure 6. As a widely studied virtual rotor, full CFD solutions [33] are often used as the reference data. Although the typical rotating speed of the *NREL 5MW* rotor is a few times less than the *NREL Phase VI* rotor, the resulted tip speed ratio is higher, which is closer to the optimal design point. The blade tip of the *NREL 5MW* rotor is smoothly sharpened where the slope at the last 4 m is 0.46, and the chord length at the tip is only 4.3% of the largest chord length in the blade inboard part. It can be seen in Figure 6 that the loads at the tip naturally approach to zero even without a tip loss correction.

<sup>ȑ</sup> **Figure 6.** *Cont*.

**Figure 6.** Normal force *Fn* and tangential force *Ft* along the *NREL 5MW* blade.

As a large wind turbine with pitch control, there is no flow separation on the most area of the blade, including the tip. That means the flow pattern ideally keeps the same on the blade as the wind speed increases, which is significantly different from the situation of the *NREL Phase VI* rotor. As seen in Figure 6, by increasing the wind speed from 8 m/s to 11.4 m/s, the forces are greatly increased, whereas, their distributions near the tip region remain similar. All the Glauert-A/B/C corrections over-predict the normal forces near the blade tip, whereas, better agreement with the reference data is found for the tangential forces. Such a result is to some extent similar to that of the *NREL Phase VI* rotor at 7 m/s. Slight discrepancies are also noticed between the three Glauert-type corrections, and the Glauert-C performs better among them.

#### 6.1.2. Comparison with BEMT Results

The Glauert tip loss correction was originally proposed for BEMT. As applied to AD/NS simulations, the correction factor *F* works in a different way. It is worth to mention that the specific difference in the final results would be caused by the transplantation from BEMT to AD/NS. In order to perform a comparison, BEMT computations under the same conditions of the present study are conducted.

The comparisons were shown by percentages, which represent the relative change of blade loads before and after using a tip loss correction,

$$
\delta\_n = \frac{F\_{n, \text{No}} - F\_{n, \text{Gl}}}{F\_{n, \text{No}}} \times 100\%\_{\text{s}} \tag{58}
$$

ȑ

$$\delta\_{t} = \frac{F\_{t, \text{No}} - F\_{t, \text{GI}}}{F\_{t, \text{No}}} \times 100\%\_{t} \tag{59}$$

where *Fn*,*No* and *Fn*,*Gl* are the normal forces obtained from computations with no tip loss correction and with the Glauert correction (the standard Glauert correction in BEMT and the Glauert-type corrections in AD/NS), respectively, *Ft*,*No* and *Ft*,*Gl* are the corresponding tangential forces.

The above parameters are calculated independently using BEMT and AD/NS. The results are shown in Figure 7 where the Glauert-BEMT denotes the results for the standard Glauert correction used in BEMT, and the Glauert-A/B/C denotes the results for the Glauert-type corrections used in AD/NS. A certain degree of difference between the Glauert-BEMT and the Glauert-A/B/C results is observed. The Glauert-C again shows the best performance, since it is generally most consistent with the Glauert-BEMT, which means the transplantation of the tip loss correction from BEMT to AD/NS does not notably influence its performance. The results for Cases 3 and 4 also agree with the above conclusion, which is not displayed here for simplicity.

**Figure 7.** Distributions of δ*n* and δ*t* along a rotor blade.

#### *6.2. Results of New Tip Loss Correction*

The simplified form (see Section 3.2) of the new correction developed by Zhong et al. [33] is used in the present simulations. The results of blade loads for Cases 1, 2 and 5 are displayed in Figure 8a–c, respectively. Case 1 represents a typical condition of the *NREL Phase VI* wind turbine with the fully attached flow on its blades, while Case 2 represents another condition of flow separation (stall). Case 5 represents the rated operating condition of the *NREL 5MW* wind turbine. The results for Cases 3 and 4 are not displayed here for simplicity. (The performance of the correction in Case 3 is similar to that in Case 2, while the performance in Case 4 is similar to that in Case 5.) The Glauert-C correction, which performs best in Glauert-A/B/C, is compared with the new correction. Full CFD results [33] are employed as the reference data for the *NREL 5MW* wind turbine. As compensation for the sparse points of the experimental data, full CFD data [8,33] are also employed for the *NREL Phase VI* wind turbine.

(**b**) Case 2: *NREL Phase* Ď, *V*0 = 10 m/s, ȍ= 72 rpm

<sup>ȑ</sup> **Figure 8.** *Cont*.

(**c**) Case 5: *NREL 5MW*, *V*0 = 11.4 m/s, ȍ= 12.06 rpm

**Figure 8.** Normal force *Fn* and tangential force *Ft* along the *NREL 5MW* blade.

It is seen in the result of Case 1, as shown in Figure 8a, that the full CFD result achieves the best agreement with the experimental data and the new correction is superior to the Glauert-C correction. As compared to the Glauert-C correction, the overestimation of *Fn* in the tip region is properly corrected towards the reference data by the new correction. A similar tendency is seen from the tangential force distribution. In Case 2, as shown in Figure 8b, the new correction maintains its accuracy, while the Glauert-C correction produces larger positive error, due to the occurrence of flow separation. The performance of the two corrections is similar to those in BEMT computations [33].

For the *NREL 5MW* rotor, as shown in Figure 8c, the new correction again makes a better correction to the blade loads in the tip region. The corrected curve of *Fn* for the new correction is closer to the reference data compared to that for the Glauert-C correction. The difference between the two corrections in tangential force (*Ft*) is small, and both the corrections can be considered to be doing well.

#### *6.3. Comparison of Velocity Field*

The blade loads are affected by the local flow around the blade, and the forces consequently modify the flow around the blade and further in the wake. The flowfield simulated by the AD/NS demonstrates how the tip loss corrections affect the velocity field. In order to clearly show the influence of using different tip loss corrections, instead of directly showing the velocity field, the following velocity difference is defined,

$$
\delta\omega\_z = \frac{\left|\mu\_{z,\text{new}} - \mu\_{z,\text{Gl}}\right|}{\mu\_{z,\text{new}}} \times 100\%,\tag{60}
$$

ȑ

ȑ

where *uz*,*new* and *uz*,*Gl* are the axial velocity obtained from the simulations using the new tip loss correction and the Glauert-C tip loss correction, respectively.

Figure 9 displays the contours of δ*uz* in the flow field of the studied two rotors, which highlights the influence of the tip loss corrections in terms of changing the velocity field not only around the blade tip, but also in the wake. The influence of δ*uz* is more concentrated near the blade tip in the rotor plane, whereas, the inner part is hardly affected. In the wake, the influence range of δ*uz* tends to cover a wider radial space downstream. The magnitude is even increased in the wake of the *NREL 5MW* rotor.

**Figure 9.** Contours of the axial velocity difference between the results of the new and the Glauert-C tip loss corrections.

#### **7. Conclusions**

This work presented AD/NS simulations for the *NREL Phase VI* and the *NREL 5MW* rotors, which represent the typical small size and modern size rotors. The primary purpose of the numerical study is to evaluate the performance of tip loss corrections employed in the AD/NS simulations. Three different implementations of the widely used Glauert tip loss factor *F*, denoted as Glauert-A/B/C, are discussed and evaluated in the study. In addition, a newly developed tip loss correction is also evaluated. As an extension, the influence of employing different tip loss corrections to the velocity field is studied as well. The following conclusions are drawn from the present study:


The above observations help to understand the performance of different tip loss corrections employed in AD/NS simulations. It is indicated that the blade loads, as well as the wake velocity, can be simulated better by choosing the best implementation of the Glauert tip loss factor or employing the new tip loss correction. That is the major contribution of this study.

**Author Contributions:** Conceptualization, W.Z.; Formal analysis, W.Z. and W.J.Z.; Funding acquisition, T.G.W.; Investigation, W.Z. and W.J.Z.; Methodology, W.Z.; Supervision, T.G.W. and W.Z.S.; Writing—Original draft, W.Z.; Writing—Review and editing, W.J.Z. and W.Z.S.

**Funding:** This work was funded jointly by the National Natural Science Foundation of China (No. 51506088 and No. 11672261), the National Basic Research Program of China ("973" Program) under Grant No. 2014CB046200, the Priority Academic Program Development of Jiangsu Higher Education Institutions, and the Key Laboratory of Wind Energy Utilization, Chinese Academy of Sciences (No. KLWEU- 2016-0102).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
