Development and assessment of a method to determine the angle of attack on an operating wind turbine by matching onboard pressure measurements with panel method simulations

Wind Energy is substantially growing in recent years and is now one of the most competitive renewable energy sources on the market. To further foster the growth of this energy source, increasing effort is put into building accurate numerical models. Most models compute the loads acting on the turbine as a dependence of some sort to the angle of attack (AoA). Accurate AoA measurements would allow for comparison with experiments and would be of great benefit for the improvement of numerical models and the investigations of aerodynamic phenomena such as stall delay. However, the determination of the angle of attack during operation is troublesome to the present day. In addition to what was already mentioned, the AoA is key to evaluate loads acting on the wind turbine and assessing experiments, computational models, and aeroelastic models. This paper proposes a simple comparative method to estimate the AoA based on pressure distributions. The proposed method is tested using data from different numbers of pressure taps placed on the Berlin Research Turbine (BeRT) at the Hermann Föttinger Institut of the Technische Universität Berlin. The predicted results are in line with those from other methods while the operating conditions to which the model can be applied are improved.


Introduction
Turbine upscaling, retrofit on installed rotors and development of new projects both onshore and offshore have made wind the most competitive and fast-growing renewable energy source in the market. To improve wind turbines further, a better understanding of the actual Angle of Attack (AoA) on a wind turbine blade is most important, since it is key to correctly evaluate the loads acting on it, which are in turn pivotal for proper aeroelastic development. Operating in a non-optimal AoA in fact would come with the drawback of decreasing the aerodynamic performance and also affecting the lifetime of the turbine.
In the case of a wing, the AoA can be evaluated by measuring the angle of the incoming and undisturbed flow with respect to the wing itself. Conversely, in the case of a wind turbine, determining the AoA is not straightforward: as an effect of the energy extraction, the local inflow velocity changes in the proximity of the rotor, and therefore, the geometrical angle between the blade section and the undisturbed flow is not representative of the AoA. This phenomenon is commonly referenced as induction and its magnitude depends on the blade loads and thus, on the AoA itself, which brings inherent problems to any iteration model between loads and AoAs such as inverse Blade Element Momentum (BEM).
In other terms one would need to include the effects of the wake on the local inflow but not the flow deformations caused by the blade itself in order to correctly evaluate AoA; a non-trivial task. The latter, also known as the upwash effect, has been included in corrections [1] and also, neglected [2] leading in both cases to good agreements between experiments and computational methods, as highlighted and discussed in detail in Sheppers et al. [3] Several authors have proposed methods to estimate the AoA of a rotating blade. Starting from the pressure distribution the inverse BEM method is one of the most common ways to determine the AoA. From pressure measurements, normal and tangential forces are defined. Assuming them uniform over an annular ring containing the blade section, the wake-induced velocities are calculated with the momentum theory. The free stream wind speed and the rotational velocity are known, and thus it is possible to obtain the e 1 ffective velocity vector and the desired angle of attack.
Sant et al. [4] suggested an algorithm to estimate the AoA from blade pressure distribution based on the iteration of an AoA guess, lift calculation, bound circulations and induced velocities. Soto-Valle et al. [5] proposed an application of the Gaunaa's method [6] on the rotating blade, using a reduced amount of pressure sensors positioned chordwise on the surface of a wind turbine model.
Both methods showed positive results. However, the former implies a complex implementation, while the latter exhibits problems with the estimation of the dynamic pressure, reducing the effectiveness in a narrow range. This study provides a simpler prediction of the AoA, based on the matching up of the pressure distribution on a rotating blade with airfoil simulations made with the industry-standard panel methods, Drela et al. [7].

Experimental setups
The analyses shown in the following chapters rely on two existing datasets. The first one was obtained from measurements on an isolated wing in a wind tunnel. The second was obtained from wind tunnel measurements on a rotating wind turbine model. Both datasets were used to derive assumptions and to define the parameters for the methodological development.
The 2D wing dataset was used for validation in a case where the AoA is known. The rotating framework dataset was used to test the methodology, while it is compared with other AoA prediction methods, such as Gaunaa's method. The following sections provide the context of the experiments but for a more detailed description the original references are provided the setups [8,9,10].

2D Wing
The experimental measurements were carried out on an isolated wing with the Clark Y airfoil profile, Gorrell et al. [11]. The wing has a chord length of 0.305 and a span length of 0.800 . The test campaign was performed in the 2D active grid wind tunnel of the University of Oldenburg. The test section has a length of 2.6 m and the outlet has a cross section of 1.0 x 0.8 m 2 . The experiments were conducted with inflow velocity U∞≅ 15 m/s. This value of velocity resulted in a Reynolds number of 2.9×10 5 using as reference length, the chord of the wind turbine model blade described below at 75% . The wing was positioned vertically on the top of a turntable in the middle of the test section 1.1 m downstream the active grid. The AoA was controlled by a stepper motor placed on the top of the wing.
The wing was equipped with two force gauges, a torsion sensor, a custom made three-hole probe, a trailing edge flap (TEF), and 32 surface pressure taps. Due to manufacturing reasons, not all the pressure ports were in the same spanwise position. The chord positions of the pressure taps are shown in Figure 9 when pressure distributions are compared. The reference static pressure was obtained from two wall-mounted pressure taps. More information about the tunnel and the wing model can be found in Bartholomay et al. [8].

Berlin Research Turbine -BeRT
The experimental measurements were carried out on the Berlin Research Wind Turbine (BeRT) test rig, placed in the wind tunnel facility at the Hermann Föttinger Institut of the Technische Universität Berlin. BeRT is a three bladed wind turbine model with a 1.5 m diameter, designed with a single profile airfoil along the entire span and without cylindrical root section, The model was placed in the settling chamber which has a cross-section of 4.2 × 4.2 m 2 . All the experiments were conducted with inflow velocity U∞ ≅ 6.5 m/s with a level of turbulence up to 1.5%.
The static pressure was measured with two Prandtl tubes and used as a reference for the pressure measurements on the blade; they were located on the wind tunnel walls at 2.7 m height and 0.43 rotor diameters upstream of the turbine. In the experimental campaign, one of the blades was equipped with, 3 three-hole probes mounted at 65%R, 75%R, 85%R, TEFs located between 60%R and 90%R and pressure taps.
Experimental data used in this work were collected using the pressure taps mounted on the equipped blade. The taps were located chordwise at 45%R. Pressure taps positions on the airfoil chord are shown in Figure 2a. Pressure transducers (HCL0025E)) were located in a pressure box inside the hub.
In order to compare the performance of BeRT at several operational points, three alignments (yaw angle) and eleven pitch angles were analyzed ( Table 1). The turbine was fixed by means of the yaw and pitch angles in each case. The measurement duration was 100 rotations at each operational point to have a significant statistical average in the data. The data is phase averaged with respect to the azimuth angle with steps of Δϕ =1°. More information about the tunnel and the turbine model can be found in Pechlivanoglou et al. [9] and Vey et al. [10].

Methods
In this section, the methodology of this research is described. Two methods are compared throughout the study and applied to the pressure distributions from the rotating blade. On one side, the method proposed by Gaunaa and Andersen [5], which in the interest of brevity, only a brief description is presented and the details of its derivation and implementation can be found in Soto-Valle et al. [5] and Gaunaa [12,13]. On the other hand, the here proposed method is described at the end of this chapter. Gaunaa [13] presented an analytical expression for the unsteady 2D force distribution on a variable geometry airfoil undergoing arbitrary motion under the assumptions of incompressible, irrotational and inviscid flow. The local pressure difference could be reduced to Eq. 1 when the the high order terms and the added mass from the pitching are neglected by using only the pressure information at the 12.5% of the chord:

Gaunaa's method
where Δ is the pressure difference between pressure and suction side at 12.5% chord, the density, the inflow velocity; and 1 and 2 fit parameters obtained from static measurements. Therefore, the effective AoA can be derived by solving when the constants are known.

Proposed method: matching-up 2D and 3D pressure distributions
The method is based on matching pressure distributions measured on the wind turbine blade to the database from the computational tool XFOIL, where the AoA is known. The validation of the database is presented in Sect. 4.1. After the simulations, a linear interpolation of XFOIL data was done in order to obtain data at the same chord positions of the pressure taps during the measurements on BeRT. Given the high number of panels, linear interpolation should not introduce additional modeling errors.
The predicted AoA was selected as the one for which the Fit Parameter (FP) was minimized between the matching pressure distributions. The parameter FP is based on the point-to-point comparison of measured and simulated data and is defined as: where the subscript i refers to the chord position. Additionally, the coefficient of determination R 2 , between the two pressure distribution curves is calculated. The aim of R 2 introduction is to verify if the AoA predicted by using FP was the same as the one obtained by choosing the corresponding to the maximum value of R 2 . A flowchart of the proposed method is shown in Figure 3. The experimental pressure distributions for a given yaw angle () and pitch angle (), for wich the AoA has to be determined are compared to the XFOIL database using the FP parameter. The XFOIL pressure distribution that minimizes FP is selected as best match, and thus the AoA is determined since this parameter is known in the numerical simulations. More detail about the steps to follow to perform such comparison are shown in Figure 4.

Figures 3 (left) and 4 (right): Flowcharts of the proposed method
In comparison to Gaunaa's method, the proposed method, while requiring more complex experimental data to compare since a full Cp distribution is required, is conceptually simple to apply and, as will be discussed in the following of this paper, is able to predict the AoA for a wider range of operating pitch and yaw angles. In fact, when applying Gaunaa's method, the relative velocity (Eq. 1) needs to be estimated. This can be challenging in a rotating framework with high yaw and pitch angles ultimately leading to highly uncertain predictions when applying Eq. 1.

Results
The results are presented in this section. The general tuning of the numerical panel method parameters that was done as a preliminary step is discussed. The database is then initially validated using the 2D wing experimental dataset. Finally, the method is tested in the rotating framework (BeRT) on the aligned-flow case (ψ =0°) and the misaligned-flow cases.

Validation
The values of the XFOIL parameters chosen are based on the experimental setup. In this way, the range of AoAs between -10° and 20° has been selected based on previous measurements made on BeRT blades and chosen wide enough to cover all the possible AoAs experienced by the blade during the test cases. The step Δα =0.1° used for the computations was chosen as a better compromise between accuracy and computational time, without compromising on the generality of the prediction method.
As the BeRT blades and the wing were completely clean and without any trip strip tape during the test campaign, only free-transition based on Ncr was used and no forced transition was imposed.
To evaluate the Reynolds number, density was calculated considering the air temperature (T = 18±1.5°C), pressure (P∞ = 100kPa) and relative humidity (H = 40 ± 5%). The viscosity value has been chosen according to Sutherland's formula, Chapman et al. [14]. The calculation of the chord length at the spanwise position of the pressure taps (45%R) has been done through a linear interpolation of the BeRT blade chord distribution, Pechlivanoglou et al. [9]. The definition of local relative velocity, Urel, was based on three different criteria of evaluation: geometry, pressure taps measurements and three 3-hole probes measurements, Bartholomay et al. [15] and Klein et al. [16]. Table 2 summarizes the XFOIL parameters while Table 3 the methods to estimate the relative velocity with their resulting Reynolds number.
It must be noted that the estimation from the 3-hole probes was made through a linear extrapolation, since the positions of the three three-hole probes were different from the position of pressure taps.
The influence of the Reynolds number is shown in Fig. 5, using the range of Re numbers from Table 3. The lift curve shows no significant differences in the Re-range between the AoA chosen to simulate. Therefore, considering that the measurements in the wind tunnel of the wing were conducted at Re = 2.9 × 10 5 the same value was chosen for this analysis.
The effect of the Mach number is also negligible in this case, with no significant changes in the polar curves ( Figure 6), thereby Ma = 0 has been assumed for the following part of the work. The free transition on the airfoil surface is described in XFOIL, using the Ncr method. Thus, the parameter Ncr, required for the computation, is the exponent of the transition method. In the present case, according to the previous knowledge on previous experiments made on the BeRT, the values of Ncr evaluated were from 6 to 9 which correspond to dirty and average clean test cases, respectively.   Furthermore, as a preliminary evaluation of the matching method, a comparison between Cp distributions from XFOIL and from the 2D wing tested in the Oldenburg university wind tunnel, as discussed in section 2.1 was done. The AoAs were known for both cases. In Figure 9, the comparison of the Cp distributions from XFOIL at α = 5° and from the wing at α = 4°, 5°, 6° and 7° are shown. The experimental and numerical Cp distributions at 5° of AoA show very good agreement, and thus the method can be considered effective in determining the AoA in non-rotating conditions. While more extensive testing could be carried out on this testcase, as previously discussed, this work focuses on AoA determination in a rotating framework and thus this result is considered sufficient. From the results, it can be noticed that the airfoil employed in the turbine posses a large suction peak located close to the Leading Edge (LE). For this reason, several approaches were tested to perform the matching-up procedure. In particular, referring to the side of the blade either the Suction Side (SS) or the Pressure Side (PS), or the number of pressure tap data used. The matching approaches listed below will be considered in the following part of the work: • PS -SS with x/c ≥ 0 or x/c ≥ 0.01 that will be indicated respectively as PS -SS 0% and PS -SS 1% • SS with x/c ≥ 0 or x/c ≥ 0.01 that will be indicated respectively as SS 0% and SS 1% The previous approaches are chosen after some consideration referred to the baseline case: • The AoAs predicted using just the PS data significantly deviate with respect to the other predictions ( Figure 10). The differences in the prediction are caused by the great difference of Cp values corresponding to the chord position x/c = 0.005 on the PS (see Figure 11). This difference could possibly be justified by small manufacturing errors of the pressure taps, which is likely at these small dimensions. To overcome this problem a strategy could be to ignore the first two points on PS (x/c=0.005, 0.01) but, in this case, the matching would be based on very few data points because on PS there are just ten pressure ports after that position, and it would be less reliable. Therefore, it was decided to test the method using all the PS pressure ports available for matching (all the chord positions).
• The decision to exclude the first 1% of the blade from the matching is mainly due to the presence of the pressure peak close to the LE on the SS. The shape of the peak is difficult to describe with the linear interpolation as it is shown in Figure 12. For example, in 12 a) XFOIL Cp distribution (red) is close to the experimental data, but with the interpolation the interpolated point (blue) in x/c = 0 goes well below the measured one. In Figure 12 c) the shape of XFOIL pressure peak is lost due to the interpolation, even if the interpolated point remains over the measured one in x/c=0. From Figure 12 b) instead, even if the XFOIL pressure peak (red) is greater than the experimental one (black), the XFOIL interpolated point (blue) goes below the experimental one. Therefore, due to the poor resolution of the experimental pressure taps, the position of the interpolated point in x/c = 0 moves around the measured point without a real dependence on the XFOIL Cp distribution shape.
• The method is tested also just considering SS data, to reduce the amount of data and thereby to reduce computational time for future uses of the tool. In the baseline case, differently from the PS, the SS shows the same behaviour of the PS-SS from 0% and 1% of the blade chord ( Figure 10)

Axial flow
The behaviour of the AoA as a function of the azimuthal angle is depicted in Figure 13 for the baseline case. The proposed method appears to be able to capture the effect of the tower. The tower influence appears as a reduction of the AoA of Δα = 0.5° around ϕ = 180° (tower location) and where the AoA takes its minimum value. The AoA decrement is in line with experimental results exhibited by Soto-Valle et al. [5], where a drop of 0.7° +/-0.3° is shown. Results are within the uncertainty band of AoA predictions from Gaunaa's method (Figure 13) that show an AoA drop of 0.7° +/-0.2° of uncertainty [5].
The approaches PS-SS 1% and SS 1% give very similar results for all azimuthal angles. In both cases a difference of 0.75° ± 0.05° with respect to the approach PS -SS 0% is shown. This difference can be justified, with the effect of the pressure peak on the SS close to the LE. In Figure 14, the pressure coefficient distribution of the predictions that best matched the experimental data, obtained using the different approaches of the proposed method for the baseline case at ϕ = 315° are shown. This azimuthal angle is chosen as it is the farthest from the wind tunnel walls and tower, and thus it is assumed is the lesser influenced by blockage and tower shadow. Firstly, as also observed in Figure 14, the AoA predicted by the method, using the PS-SS 1% and SS 1% approaches, is the same, as can be seen by the perfect overlapping of the pressure distributions corresponding to these approaches.
Secondly, the main differences between the approaches 1% and 0% are located near the LE, where the pressure peak on the SS is located. The pressure peak shown by the pressure distributions corresponding to the 0% approaches is higher than that shown by those corresponding to the 1% approaches. Moving forward to the TE, the pressure distributions show more similar behaviours between each other. An exception is made by the point where the transition is located on the SS and for which the Cp values differ the most.
Both FP and R 2 obtained for several azimuthal angles are shwn in Table 3. It can be noted that the lower the FP, the higher the R 2 . The coefficient of determination magnitude is over 0.9 in all cases. Furthermore, considering the azimuthal angles shown in Table 3 the approaches PS-SS 1% and SS 1%, are those for which the coefficient of determination assumes its maximum values (R 2 ≥ 0.99), and the fitting parameter its minimum values. By observing of the Cp distributions at = 315°, when the pitch angle increases the pressure peak on the SS corresponding to Cp distribution of SS 0% is lower than that of the Cp distribution PS-SS 0%. The former is also closer to the peak of BeRT data. This behaviour is shown until θ = 6°; when θ ≥ 8° this consideration is no longer valid. Indeed, the trend is confirmed by the results obtained θ = -2°, 2°, 4° and 6° when shown in Figures 15 and 16 (blue lines).
It is worth noting that in Figure 16, there are no curves corresponding to the results obtained using Gaunaa's method since it was not able to predict the AoAs outside the range -6° ≤ θ ≤ 2°. Therefore, the proposed method can be used to predict the AoA in a wider range of pitch angle conditions. Furthermore, even if the results obtained with the PS-SS 0% approach agree well with those obtained using Gaunaa's method (see Figures 13 and 15) the best predictions seem to be those obtained with the 1% approaches.
It is worth noting that when comparing results, in terms of AoA, obtained with different methods, the agreement does not necessarily mean that the prediction is accurate or reliable because the actual AoA is never known since it is not directly measurable. Therefore, the results obtained with Gaunaa's method can be used as reference, but they have not to be considered as rigorous valid since they are influenced by the conditions and limits of the methodology development.
In Figure 17 the evolution of the azimuthal average AoA over the pitch angles is shown. A linear trend is shown. Point to point it is found that the slope is not constant, changing as moving to the pitch angle. These changes can be justified by the fact that for each pitch angle the induction factor is different since the position of the blade with respect to the inflow is different. Thereby it is different how the blade affects the streamlines, and consequently the velocities. The overall relation between these angles is fitted as y = m * x+q, for each approach. The slopes m found within the range ∆m = 0.76 ± 0.04. First, it must be noted that the mutual behaviour of the evaluated approaches for the proposed method shows some common aspects with the baseline case previously discussed. Therefore, the results achieved with SS 0% are more similar to PS-SS 0% results than those obtained with PS-SS 1% and SS 1%. However, the small constant difference between the AoAs predicted by the SS 0% approach with respect to the PS-SS 0% approach, shown in the baseline case, is no longer present in this case. The biggest differences between the SS 0% and PS-SS 0% are shown for 160° ≤ ϕ ≤ 240°, and they are approximately of 0.35°±0.05°.
Outside that range the difference between the AoAs predicted becomes smaller and at some points the AoA obtained with the two approaches are the same. Furthermore, even with yawed inflow, a difference between the AoAs predicted by the PS-SS 1% and the SS 1% approaches regarding the PS-SS 0% approach, can be found. However, even this difference is no longer constant over the azimuthal angle. In detail, the difference between the AoAs predicted by the approaches at 1% and 0% increases from ϕ = 0° until more than ϕ = 180°, it remains constant and when ϕ = 180°, it decreases again.
With the introduction of a yaw misalignment, no approach can capture the influence of the tower. Its effect on the AoA predictions, consisting in a reduction of the AoA when ϕ ≈ 180° is not visible since it is covered by the crossflow. However, the effect of the crossflow is less marked than the one shown from Gaunaa's method and it results in a smoother peak of the AoA for the proposed method respect to Gaunaa's method (see Figures 18 and 19).
If the yaw misalignment is greater, as shown in Figures 20 and 21, where ψ = −30°, the results obtained with the proposed method are similar to those shown for, and ψ=−15°. In general, it can be affirmed that they present the same behaviour but in a more pronounced way, and the same considerations are valid.  Even with increased yaw angle the difference between the curves corresponding to the 0% approaches and 1% approaches is not constant. The ∆α between the curves is about 1.3° in the range of azimuth positions ϕ ≈ 180° and much lower outside that range.
In Figures 22 and 23 the evolution of the mean AoA over the azimuth, with the pitch angle is shown. The results obtained in both cases are similar to those shown with axial inflow. To assess the relation between the angle of attack and the pitch angle a linear fit, y = m * x + q, has been found. The slopes are in the range ∆m = 0.61 ± 0.01 when ψ= -15° and ∆m = 0.66 ± 0.03 when ψ=-30°. In both cases, R 2 ≥ 0.99, indicating a linear relationship between pitch angle and AoA.

Conclusions
In this work a simple method to determine the angle of attack on a rotating blade is proposed. The method is based on the comparison between measured and simulated pressure distributions, for which the AoA is known. A 2D panel method is used to compute the numerical pressure distributions. The method is tested with respect to wind tunnel data, first comparing the predicted AoA to the one set on a 2D static wing, and then comparing pressure distributions measured on a rotating wind turbine model. Firstly, a calibration of the numerical parameters of the panel method is performed. Particular focus is posed on Reynolds, Mach and Ncr, which determines where the boundary layer free transition will be located along the suction-side. The present test case was found to be insensible to the range of Reynolds and mach numbers tested, and the same can be said for the Ncr parameter, which was found to not influence results significantly. Nonetheless it is recommended that a similar calibration procedure is performed if this method shall be replicated.
When testing the proposed method against 2D wind tunnel measurements, a satisfactory match between pressure distributions predicted in the linear AoA-range was found, and thus, comparing pressure distributions in this case would be an effective way of determining AoA.
When testing the method on the dataset derived from the BeRT model wind turbine, predictions were compared to the ones derived with Gaunaa's method. Several yaw and pitch angles were used in the comparison, and several comparison methods were tested. It was found that the proposed method is able to predict AoA for a wider range of pitch angles with respect to Gaunaa's method, while predicting AoAs in line with this method in its range of applicability. The best fit between measured and computed Cp distributions was achieved when excluding the first 1% of the chord, due to the strong gradients in this area, which can easily be afected by errors.
In conclusion, the proposed method can predict AoA on rotating blades with results comparable to Gaunaa's method, while being applicable in a broader range of inflow conditions and simple to apply.