Combining Unsteady Blade Pressure Measurements and a Free-Wake Vortex Model to Investigate the Cycle-to-Cycle Variations in Wind Turbine Aerodynamic Blade Loads in Yaw

Abstract

Prediction of the unsteady aerodynamic flow phenomenon on wind turbines is challenging and still subject to considerable uncertainty. Under yawed rotor conditions, the wind turbine blades are subjected to unsteady flow conditions as a result of the blade advancing and retreating effect and the development of a skewed vortical wake created downstream of the rotor plane. Blade surface pressure measurements conducted on the NREL Phase VI rotor in yawed conditions have shown that dynamic stall causes the wind turbine blades to experience significant cycle-to-cycle variations in aerodynamic loading. These effects were observed even though the rotor was subjected to a fixed speed and a uniform and steady wind flow. This phenomenon is not normally predicted by existing dynamic stall models integrated in wind turbine design codes. This paper couples blade pressure measurements from the NREL Phase VI rotor to a free-wake vortex model to derive the angle of attack time series at the different blade sections over multiple rotor rotations and three different yaw angles. Through the adopted approach it was possible to investigate how the rotor self-induced aerodynamic load fluctuations influence the unsteady variations in the blade angles of attack and induced velocities. The hysteresis loops for the normal and tangential load coefficients plotted against the angle of attack were plotted over multiple rotor revolutions. Although cycle-to-cycle variations in the angles of attack at the different blade radial locations and azimuth positions are found to be relatively small, the corresponding variations in the normal and tangential load coefficients may be significant. Following a statistical analysis, it was concluded that the load coefficients follow a normal distribution at the majority of blade azimuth angles and radial locations. The results of this study provide further insight on how existing engineering models for dynamic stall may be improved through the integration of stochastic models to be able to account for the cycle-to-cycle variability in the unsteady wind turbine blade loads under yawed conditions.

1. Introduction

Horizontal axis wind turbines (HAWTs) are bound to operate under a 3D complex unsteady aerodynamic environment caused by various unsteady sources, such as variations of wind speed, atmospheric turbulence, wind shear, yawed flow, and the operational states of the wind turbine. The cumulative effect of these phenomena leads to significant variation of the aerodynamic loads on the blades, with azimuth angle and instantaneous change in the blade angle of attack. Unsteady effects appear when the reduced frequency of the effective pitching of the turbine blade becomes large. The effects of unsteadiness lead to dynamic stall, which impacts on a significant region of the rotor disc, especially on the inboard sections of the wind turbine blades [1]. Since these effects play an important role in predictions of wind turbine performance and fatigue lifetime of the structure, it is therefore essential to have an in-depth understanding of the aerodynamic load phenomena before the structural response can be simulated.

The National Renewable Energy Laboratory (NREL) in the USA [2] conducted a detailed experiment in the NASA Ames wind tunnel to produce very detailed measurements using high quality instrumentation under controlled conditions to avoid disturbances caused by the real environment. These measurements formed a basis for further understanding 3D flow effects and dynamic stall in a rotating environment, which is crucial for improving engineering modelling tools. Schreck et al. [3] analyzed the experimental surface pressures and normal force histories to characterize the dynamic stall vortex kinematics and normal force amplification for a rotating blade observed in the NREL Phase IV Unsteady Aerodynamics Experiment. Stall vortex convection was found to vary along the blade, leading to rapid deformation of the vortex and amplification of normal forces. When the static stall angle of attack threshold is exceeded, these forces were observed to occur frequently at all span locations except those at the blade tip. The results showed that exceeding the static stall angle leads to a two dimensional dynamic stall vortex close to the leading edge. Subsequently, it undergoes rapid three-dimensional deformation as it convects rapidly downstream. In a further work, Schreck and Robinson [4] also confirmed that, under yawed rotor conditions, dynamic stall is present when the blade angle of attack exceeds the static stall threshold line. The growth and convection of this vortex aft along the blade significantly amplifies blade aerodynamic loads. In addition, under yawed flow operation, pseudo-sinusoidal inflow angle oscillations generate dynamic stall which is highly unsteady, leaving a significant impact on the aerodynamic loads. Similarly, Simms et al. [5] showed that there were measurement uncertainties associated with the wind tunnel test data. However, the measurements were both accurate and repeatable. The aerodynamic forces, including the lift, drag, and pitching moment, were analyzed against angle of attack at five span locations for 36 consecutive blade rotation cycles. It was shown that, when angle of attack is below the static stall angle (<15°), the flow is attached and cycle-to-cycle variation in the aerodynamic loads is not significant. When angle of attack exceeds the static stall angle of attack, however, unsteady separation dominates the flow field leading to dynamic stall. As a result, significant cycle-to-cycle variations and deviations about the mean values are seen.

Several semi-empirical models have been developed in the past few decades to simulate unsteady aerodynamics loads on wind turbine and helicopter rotors under the influences of dynamic stall. These are 2D correction models that are usually integrated with simplified 3D potential flow models to predict the rotor aerodynamics. Examples include the Leishman and Beddoes model [1,6], the Onera model [7], Øye model [8], and the Risø fgh-model [9]. However, these models still depend on 2D aerofoil wind tunnel data and various semi-empirical parameters that usually cater for 3D aerodynamic phenomena, particularly those under yawed flow conditions. These models are not able to simulate the cycle-to-cycle variations in the aerodynamic loads experienced by the rotor under yawed conditions. Unlike helicopter rotors, modelling such variations is more relevant for wind turbines, given that the latter usually operate at higher angles of attack in order to maximize power extraction from the wind.

Most aerodynamic models implemented in wind turbine design codes still rely on the definition of the angle of attack to be able to estimate blade loads. Such models include the Blade-Element-Momentum and free-wake vortex methods [1,10]. Analysis of measurements to accurately derive the angle of attack in a rotating environment is, however, a major challenge due to difficulty associated with the impacts of the bound circulation at the blades and the complex 3D wake vorticity distribution. These would lead to significant differences in between the local angle of attack and the local flow angle measured by flow direction probes normally installed on wind turbine blade sections at a fixed distance upstream from the leading edge. A correction factor obtained from simple 2D wind tunnel calibration procedures is usually used to estimate the angle of attack from the flow angle, as discussed in the work of Rooij et al. [11]. However, this correction is invalid for yawed flow conditions due to the significant impacts of the unsteady shed vorticity and the effects of the skewed wake. Tangler et al. [12,13] established an iterative process to derive the angle of attack distributions by a prescribed-wake vortex model and the experimental aerodynamic loads on the NREL rotor under axial flow conditions. The process mainly starts with an initial radial distribution for the angle of attack and, subsequently, the unsteady normal and tangential forces coefficients and the angle of attack are used to determine the lift force coefficient. The bound circulation is then calculated using the Kutta–Joukowski theorem and is further used to determine the induction on the blades, and a new angle of attack distribution is finally obtained. A successive repetition for this process is made until convergence in the angle of attack is met for all time steps and radial positions. Sant et al. [14,15] applied a similar approach to estimate the unsteady angle of attack using a combination of a free-wake vortex model and the blade measurements of the NREL wind turbine rotor Phase VI under axial and yawed flows. Shen et al. [16] established two approaches in which the Biot–Savart integral is used to calculate the influence of the bound vorticity on the velocity field. In the first approach it is assumed that the load and the velocity distributions at a monitoring point in the vicinity of a rotating blade are known from measurements or CFD. The experimental load distribution is used to determine the bound circulation and, hence, the induced velocities are obtained. These were then subtracted from the velocity at the monitor point in order to find the angle of attack along the blade. The second approach involves calculating the local circulations on the surface contour of the blade from the pressure distribution. Subsequently, the influence of the bound vorticity was evaluated in order to enable the velocity monitor points to be located near the blade and, therefore, the angle of attack is determined. Although the above-mentioned approaches led to a better insight into the aerodynamics of rotors in axial and yawed flows, the cycle-to-cycle variation in the aerodynamic loads and the resulting change in the blade angle of attack and the induction on the blades were not presented. This is because only the experimental aerodynamic loads and the dynamic pressure distributions averaged over a number of rotor rotations (cycles) were used in the analyses.

The main objective of this study is to investigate the cycle-to-cycle variations in the aerodynamic loads, angles of attack, and the axial-induced velocities experienced by yawed wind turbine rotors operating at constant speed with a fixed yaw angle and a steady and uniformly distributed wind flow. This work is an extension to the work of Sant et al. [14,15] who investigated the yaw aerodynamics of the NREL Phase VI rotor by coupling the blade normal and tangential coefficients (C_N and C_T) derived from the measurements to a lifting-line free-wake vortex model through an iterative algorithm. Sant et al. only based their analysis on the averaged values for C_N and C_T at each blade azimuth position obtained from multiple revolutions. The present work adapts this iterative algorithm using a different vortex model to account for cycle-to-cycle differences in these coefficients by basing the free-wake computations on the entire time history of the experimentally derived C_N and C_T values obtained over 30 successive rotor revolutions. The revised method made it possible to quantify the cycle-to-cycle variability in the angle of attack and wake-induced velocities corresponding to the variability in C_N and C_T values over different blade locations and operating conditions. This analysis provides insight of how existing engineering models for dynamic stall currently embedded in wind turbine design software packages may be upgraded to be able to model the stochastic behavior on the rotor blades under yawed conditions.

2. NASA Ames UAE Wind Tunnel Data Used

The NREL Phase VI rotor is a two-bladed HAWT that was tested comprehensively in the NASA Ames 24.38 m × 36.57 m wind tunnel (The National Renewable Energy Laboratory, Golden, CO, USA) under a wide range of operating conditions. The turbine had a rotor diameter equal to 10.1 m, with the blade sections having the geometric profile of the S809 aerofoil. The tower has three sections: 0.6096 m outer diameter at the base (from 0 to 3.04 m), transition section (from 3.04 to 3.54 m), and 0.4064 m outer diameter at the top (from 3.54 to 11.14 m). The rotor overhang distance (from yaw-axis to blade-axis) is 1.401 m. More details about the wind turbine geometry can be found in [17]. The surface pressure distributions were measured on one blade using 22 pressure sensors distributed on the upper and lower surfaces of the S809 aerofoil sections at 0.30R, 0.47R, 0.63R, 0.80R, and 0.95R. The aerodynamic normal and tangential loads were subsequently determined by integrating the surface pressure sensors at these radial locations. The local inflow angle (LFA) was acquired on one blade using five-hole pressure probes at five radial locations: 0.34R, 0.51R, 0.67R, 0.84R, and 0.91R. Further details about the measurements are available in [2,18]. The measurements were conducted for axial and yawed flows at wind speeds of 5–25 m/s and yaw angles of 0°–180°. For each wind speed a data record of 30 s duration was taken, which is equivalent to 36 rotor rotations (cycles) [2]. Tests considered in this paper were conducted on an upwind, rigid turbine with a 0° cone angle, a blade tip pitch angle of 3°, a constant rotor rotation of 72 rpm and yawed flow conditions at 5 and 15 m/s, and yaw angles 10°, 30° and 45°. These two wind speeds are considered, as the flow field at 5 m/s (tip speed ratio = 7.58) is usually attached, while at 15 m/s (tip speed ratio = 2.53) the flow is highly unsteady and the effect of the dynamic stall is evident. Figure 1 illustrates typical cycle-to-cycle variations of the measured unsteady normal and tangential force coefficients against the blade azimuth angle for a total of 36 rotor rotations at three spanwise locations. In Figure 1, zero azimuth angle denotes position of the blade when it is located at 12 o’clock. At the inboard section of the blade, r/R = 0.3, significant cycle-to-cycle variation was observed between 0°–120° and 240°–360°. Cycle-to-cycle variations in blade azimuth angle ranging from 120° to 240° are marginal and remain close to the mean, suggesting fully attached flow, as previously shown in the work of Simms et al. [5] and Schreck and Robinson [19]. It will be indicated later in the discussion section that the local angle of attack in the separated region is above the static stall angle of the S809 aerofoil of 15.3°, while those in the attached region are below 15.3°. The cycle loads become smaller when approaching the tip of the blade. The cycle-to-cycle variations observed for fixed and uniform wind speeds were predominantly due to the presence of dynamic stall separation that dominates the flow field [5,19]. Other factors, including wind tunnel turbulence, measurement errors in the pressure measurement techniques, and blade elasticity effects could have also affected such cycle-to-cycle variations. However, their influence was found to be marginal [2]. The turbulence intensities for the tunnel wind speeds of 5 and 15 m/s were found to be equal to only 3.28% and 0.41%, respectively. The errors in the measurements due to the installation of the five-hole probes and the pressure sensors on the upper and lower surfaces were evaluated by repeating the tests with and without the five-hole probes to ensure such errors were negligible. Furthermore, the turbine was purposely designed to be as rigid as possible, with the hub supported to minimize the blade deflection.

3. Free-Wake Vortex Model

The free-wake vortex model WInDS (Wake Induced Dynamics Simulator) was used in this study. WInDS is an open-source code developed and validated by Sebastian and Lackner [20,21]. This model was further improved by Farrugia et al. [22] to reduce the large computational times when modeling a large number of rotor rotations (cycles). A wake cut-off parameter was introduced such that the number of wake filaments in the near wake is fixed, whereas those after the reference point are discarded when applying the Biot–Savart law. This process showed lower computational times, enabling the efficient simulation of a large number of rotor rotations while maintaining the same level of accuracy [22]. WInDS uses a lifting line method with discrete line filaments at c/4 to model the blade bound circulation distribution, Г_B, determined by the Kutta–Joukowski theorem. A time-marching algorithm based on an Euler predictor-corrector method is also incorporated to generate a vortex sheet for each blade. Each sheet consists of a number of straight-line filament segments to model the wake vorticity as a trailing circulation, Г_T, and a shed circulation, Г_S, as illustrated in Equations (1) and (2), respectively. The trailing circulation presents the variation in spanwise bound circulation, while the shed circulation presents the variation in bound circulation strength with time. The filaments are connected by nodes which are updated continuously after each rotor time step.

$${\mathsf{\Gamma }}_T=\frac{\partial {\mathsf{\Gamma }}_B}{\partial r}\mathrm{d}r$$

(1)

$${\mathsf{\Gamma }}_S=\frac{\partial {\mathsf{\Gamma }}_B}{\partial t}\mathrm{d}t$$

(2)

A vortex model is implemented in WInDS, the Lamb–Oseen vortex model [23,24], and is used to cater for the viscous effects in the evolution of the near wake. This model is a solution to the one-dimensional laminar Navier–Stokes equations in which an axisymmetric solution for the swirl velocity was used, assuming that the axial and radial velocities are zero. The form of the Lamb–Oseen vortex model for the swirl velocity is:

$$V_{\mathsf{\theta }}(r)=\frac{\mathsf{\Gamma }}{2\mathsf{\pi }{r}_v}\left[1-\text{exp}(-\frac{r_v{}^2}{4\mathsf{\nu }t})\right]$$

(3)

The viscous core radius (the radial location at which the swirl velocity reaches the maximum value) can be determined by differentiating Equation (3) with respect to r_v (radial location from vortex axis) and substituting the derivative equal to zero. The viscous vortex core radius as a function of a gradual growth with time is defined as:

$$r_c=\sqrt{4{\mathsf{\alpha }}_L\mathsf{\nu }t}$$

(4)

where α_L = Oseen empirical constant equal to 1.25643. The Lamb–Oseen model was further developed by Squire [25], where an average apparent or “eddy” viscosity coefficient δ_v was incorporated into the Lamb–Oseen model to cater for effects of turbulence development on the vorticity diffusion rate. A time-offset parameter, S_c, was also proposed to give a finite vortex core at the origin of the trailing vortex with a more gradual increase with time as follows:

$$r_{\mathrm{c}}=\sqrt{4{\mathsf{\alpha }}_L{\mathsf{\delta }}_v\mathsf{\nu }(t+S_c)}$$

(5)

where S_c is a time-offset parameter defined as:

$$S_c=\frac{r_0^2}{4{\mathsf{\alpha }}_L\mathsf{\nu }{\mathsf{\delta }}_v}$$

(6)

where r₀ is the initial core radius, which can be assumed to be approximately equal to half the blade tip thickness.

In the free-wake code, the 3D-induced velocities are determined by the Biot–Savart law as given by:

$$\overrightarrow{du}=K_v.\frac{\Gamma }{4\pi }.\frac{\overrightarrow{ds}\times \overrightarrow{r}}{\dot{{\left|\overrightarrow{r}\right|}^3}}$$

(7)

where K_v is a viscous parameter to account for viscous effects in rotor wake.

4. Estimation of the Unsteady Angle of Attack Using the NREL Experimental Data and the Free-Wake Vortex Model

For the present study, the main algorithm implemented in the free-wake model in WInDS, (see Section 3), was modified to estimate the angle of attack variation, with time from the known unsteady aerodynamic loads derived experimentally by integrating blade surface measurements. As already highlighted in Section 1, the approach established by Sant et al. [14,15] was used. However, a major difference in this paper is that the procedure is successively applied for each individual rotor rotation (cycle) over 30 revolutions. The time-series data for C_N and C_T over such rotations are used by the free-wake code. This involves a more complex analysis than that of Sant et al., who simply used the average values for such coefficients for each blade azimuth position obtained from a multitude of rotor rotations. In this way, the cycle-to-cycle differences in the blade aerodynamics loads observed in the measurements, especially in the presence of dynamic stall, are taken into account, thereby capturing more realistically the unsteady wake phenomena occurring in the close proximity of the rotor. Indeed, the effects of the cycle-to-cycle variation in the aerodynamic loads are considered within the wake development downstream of the wind turbine over 30 rotor rotations. Figure 2 shows the main methodology adopted. It involves the following steps:

• The unsteady measured aerodynamic loads, C_N and C_T, and the unsteady measured dynamic pressure, Q_NORM, at the five stations (0.30R, 0.47R, 0.63R, 0.80R, and 0.95R), and the unsteady measured local flow angle at (0.34R, 0.51R, 0.67R, 0.84R, and 0.91R) are first collected for blade one. Given that measurements were only carried out on one blade, a phase shift of 180° is applied to the time-series data to obtain corresponding estimates for the second blade. This assumption is reasonable given that the measurements were conducted in a wind tunnel at a fixed speed and a uniform and steady wind flow.
• The relative flow velocity is computed from the unsteady measured dynamic pressure, Q_NORM, at the five given stations (0.30R, 0.47R, 0.63R, 0.80R, and 0.95R) as expressed by:

  $$V_{rel}=\sqrt{2Q_{NORM}/\mathsf{\rho }}$$

  (8)
• A linear interpolation scheme is then applied to the data in Steps 1 and 2 to derive the various parameter values at different radial locations and blade azimuth angles, as required by the discretization scheme adopted in the free-wake vortex model.
• An unsteady angle of attack distribution, α_initial, is assumed along all radial locations and blade azimuth angles for a pre-defined number of rotor revolutions. In this step, the unsteady measured local flow angle is recommended to be used as an initial assumption because it can lead to a faster convergence in the angle of attack.
• For the first time step (t = 0), initial values for the free wake vortex model are calculated, including: the position of blade nodes determined using the kinematic chain, aerodynamic loads determined using the blade element momentum (BEM) theory, bound circulation strength, initial shed and trailing circulation, Г_S and Г_T, vortex core size, and velocities induced by shed and trailing filaments on all nodes in the wake calculated by the Biot–Savart law.
• For each new individual time step (or azimuth angle), the following are calculated: (i) shed and trailing circulation, Г_S and Г_T; (ii) vortex core size; (iii) update vortex core size due to filament strain; (iv) the 3D induction at all wake nodes via the Biot–Savart law; and (v) convect wake nodes using numerical integration.
• The collected data in Steps 3 and 4 are then prescribed to the free-wake vortex model. For each individual time step (or azimuth angle), the following process is applied until successive convergence in the angle of attack is achieved for all radial locations:

(a)

Compute the unsteady lift and drag coefficients, C_L and C_D, using the initial angle of attack, α_initial, and Equation (9).

$$\begin{array}{l}{C}_L=C_N\cos ({\mathsf{\alpha }}_{\text{initial}})+C_T\sin ({\mathsf{\alpha }}_{\text{initial}})\\ C_D=C_N\sin ({\mathsf{\alpha }}_{initial})-C_T\cos ({\mathsf{\alpha }}_{initial})\end{array}$$

(9)

(b)

Calculate the bound circulations, Г_B, using the Kutta–Joukowksi theorem with zero values at the blade root and tip as given by:

$${\mathsf{\Gamma }}_B=\frac{1}{2}.C_L.V_{rel}.c$$

(10)

(c)

Generate the free wake downstream for the current time step and determine the 3D unsteady induced velocities on the lifting line. Subsequently, estimate a new angle of attack distribution along the radial locations, as illustrated by:

$${\mathsf{\alpha }}_{new}={\tan }^{-1}(V_n/V_t)$$

(11)

(d)

Check convergence between the initial and the new angle of attack distributions using the fixed point iteration method according to

$$\mathrm{Error}=|{\mathsf{\alpha }}_{initial}-{\mathsf{\alpha }}_{new}|$$

(12)

(e)

If the convergence is met by Equation (12) within a specified tolerance, the results are saved and the Steps (a)–(d) are repeated for the next time step. In case the convergence is not achieved for any time step, the initial angle of attack is replaced by the new one (α_initial = α_new) and the Steps (a)–(d) are repeated until convergence in the local angle of attack.

When a converged solution in the unsteady angles of attack is obtained for all radial locations and the blade azimuth angles, the unsteady aerodynamic lift and drag coefficients are determined by Equation (9). The time-averaged results obtained by this approach (see Subsection 5.2) agree well with results of Sant et al. [14,15]. Slight differences can be seen as the latter work only considered an average of 36 cycles for C_N, C_T, and Q_NORM as input values.

5. Results and Discussion

This section presents the results of the iterative procedure (Figure 2) for deriving the unsteady angle of attack on the NREL rotor at wind speeds of 5 and 15 m/s at yaw angles of 10°, 30°, and 45°. These wind speeds were chosen as they result in attached and separated flow conditions. In the numerical computations with the lifting-line free-wake vortex model, each blade was discretized into 41 radial elements (N_b) and an azimuth step (N_θ) of 10° was used. The relative numerical error for the verification tests (N_θ = 5° – N_θ = 10°) and (N_b = 41 – N_b = 37) are less than 1.9% in the axial induction factor (a_x = u_ax/U_∞) and less than 0.064% in the bound circulation. In total, 33 revolutions were modeled for each operating condition, with the first three rotor rotations always used as initial values for the free wake model to make sure that the wake extended downstream by more than 2R. Only the results from cycle 4 to cycle 33 were mainly used for the analysis.

In the absence of detailed near-wake measurements, it is very difficult to accurately establish the optimum values for the viscous parameters (δ_v, S_c). These parameters determine the initial core radius and core growth with time for each wake vortex filament, leaving a considerable impact on the development of the modeled helical vortex wake structure. A detailed sensitivity analysis of how parameters (δ_v, S_c) affect the induced velocities at the blades may be found in Sant et al. [15]. In this study, these were set to (50, 10) as these yield a viscous core radius that is within the order of the boundary layer thickness at the trailing edge of the blades.

5.1. Limitations in the Developed Approach for Deriving the Unsteady Angle of Attack

5.1.1. The Classic Kutta–Joukowski Theorem

The free-wake vortex model, WInDS, in Section 3, applies the classic Kutta–Joukowski theorem to model the blade bound circulation distribution under the assumption of steady potential flow. This theorem relates the lift of an airfoil to circulation. The unsteady term in the classic Kutta–Joukowski theorem (Equation (10)) is not included in WInDS and this may be a major limitation within the iterative procedure (Figure 3) for deriving the unsteady angle of attack. On the other hand, it has been observed that excluding the unsteady term yielded a small error when it was applied for the unsteady flow conditions. A case study using the parameters (N_θ = 10°, N_b = 41, and N_rev = 5) showed that a maximum relative error of less than 0.2% at (U_∞ = 5 m/s and Φ = 30°) and less than 2% at (U_∞ = 15 m/s and Φ = 30°) exists only at the blade azimuth angle 0°/360° when the steady and unsteady terms in the classic Kutta–Joukowski theorem were both involved, as depicted in Figure 3. It should be understood that the error is only small in the current simulations, as the unsteady experimental pressure distributions/aerodynamic loads were already affected by the unsteady flow, which may eliminate the effect of the unsteady term used in the classic Kutta–Joukowski theorem. Furthermore, it was shown by Sant et al. [15] and Micallef et al. [26] that the contribution of shed vorticity to the induction at the rotor blades is small as compared to that resulting from trailing circulation.

5.1.2. Effect of Tower Disturbance

The presence of the tower causes transient changes in angle of attack on the wind turbine blade sections. The effect is significant in downwind turbines, while this has less impact on upwind turbines, especially if the distance between the center of the tower and the rotor is large. This section briefly addresses the unsteady effects of the tower disturbance on the experimental data of the upwind NREL rotor Phase VI under yawed flow conditions used in this paper. The dimensions of the tower of the upwind wind turbine Phase VI are presented in Section 2. The 5.029 m blade length is affected by the disturbance of the top part of the tower (0.4064 m outer diameter) in the height from 3.54 to 11.14 m. The distance from the yaw-axis to the blade-axis is 1.401 m. Using these dimensions and the potential flow theory [27], the flow velocity upwind of the tower can be calculated by:

$$\text{U}={\text{U}}_{\infty }(1-\frac{{(\text{D}/2)}^2({\text{x}}^2-{\text{y}}^2)}{{({\text{x}}^2+{\text{y}}^2)}^2})$$

(13)

where D is the tower diameter, and x and y are the longitudinal and lateral distances to the tower center. Using Equation (13) and the tower dimensions, it was found that approximately 2.1% relative error results from neglecting the effect of the tower in the numerical approach presented in Section 4. It should be noted that although the impacts of the tower were not physically modeled, its effects were to a certain extent still accounted for since their influence was captured by the blade pressure measurements.

5.2. Variation of the Unsteady Angle of Attack with the Blade Azimuth Angle

Figure 4 depicts the variation of the local angle of attack at different radial locations as a function of the blade azimuth angle at wind speeds 5 and 15 m/s, respectively. The variation is presented over the entire 30 rotor rotations. For the lower wind speed (Figure 4a–c), all the predicted local angles of attack are below the static stall angle of 15.3° [28], indicating attached flow at all radial locations. The sinusoidal variation in the angle of attack has the highest values near the blade root at the blade azimuth angle close to 0°/360°. At the higher wind speed (Figure 4d–f), the static stall angle is exceeded, even at the outboard sections. This is clearly shown for blade azimuth positions between around 0°–120° and 240°–360°, leading to dynamic stall with massive unsteady flow separation that dominates the flow field. On the other hand, for blade azimuth angles between approximately 120°–240°, the flow remains fully attached, as indicated in [5,19].

The cycle-to-cycle variation in the angle of attack at 5 and 15 m/s is very small. This may be observed from

Table 1. The maximum normalized standard deviation in the unsteady angle of attack at 5 and 15 m/s.

U∞	Maximum Normalized Standard Deviation (σα/μα) (%)
	Radial Location	Φ = 10°	Φ = 30°	Φ = 45°
U∞ = 5 m/s	r/R = 0.30	3.92	2.06	1.51
	r/R = 0.47	2.29	2.70	2.82
	r/R = 0.63	1.93	1.51	2.26
	r/R = 0.80	1.28	1.26	1.80
	r/R = 0.95	1.38	1.46	1.75
U∞ = 15 m/s	r/R = 0.30	3.81	2.87	1.41
	r/R = 0.47	1.69	1.59	1.56
	r/R = 0.63	3.15	2.55	1.84
	r/R = 0.80	2.25	1.65	1.16
	r/R = 0.95	1.13	1.12	1.11

which shows the maximum normalized standard deviation in the unsteady angle of attack (σ_α/μ_α) at the given five radial locations. For U_∞ = 5 m/s, this increases marginally at the inboard blade sections and for the larger yaw angles. However, the maximum value reached is estimated to be lower than 4%. At a wind speed of 15 m/s, the cycle-to-cycle variations in the angle of attack are also small, with the highest normalized standard deviation reaching 3.81% at r/R = 0.30 at Φ = 10°. Although turbulence in the undistributed wind approaching the rotor may lead to a cycle-to-cycle variation in the angle of attack, such influences are likely to be very small given the fact that the turbulence levels in the NASA Ames wind tunnel were very small (see Section 2). The variation between the different rotor rotations is primarily induced by the unsteady flow phenomena at the blades and is thus self-induced by the rotor when operating in the yawed state. The resulting unsteadiness in the bound circulation along each blade results in unstable vorticity in the near wake (Equations (1) and (2)). Consequently, variations in the induced velocities, (hence also in the angle of attack) are experienced at the blades between successive rotations.

The variation in the reduced frequency across the entire rotor disc for the different operating conditions is presented in Figure 5. It may be noted that the cycle-to-cycle variation in the angle of attack shown in Figure 4 is well correlated to the reduced frequency (k = Ωc/(2V_r)). Higher reduced frequencies imply higher levels of unsteadiness in the flow regime over the blade aerofoil sections, which lead to dynamic stall when the local angle of attack exceeds the static stall angle (in this case at U = 15 m/s). Under yawed rotor conditions, their occurrence shifts to the blade inboard regions over the upper half of the disc, at azimuth positions (0° ≤ Ψ≤ 40°) and (320° ≤ Ψ ≤ 360°). These coincide with those experiencing large cycle-to-cycle variations on the angle of attack (compare Figure 5d–f with Figure 4d–f).

5.3. Variation of C_N and C_T with α (Hysteresis Loops) and the Blade Azimuth Angle

Although the cycle-to-cycle variations in α were found to be small, the corresponding variations in the normal and tangential force coefficients were significant for certain conditions. Using the approach presented in this study, it was possible to derive directly from the measurements the hysteresis loops, C_N-α and C_T-α, for 30 rotor revolutions. These are presented in Figure 6 and Figure 7 for Φ = 30°. The 2D steady wind tunnel aerofoil data provided by Ohio State University (OSU) are also included in these plots for comparison [28]. At lower wind speeds (U_∞ = 5 m/s and Φ = 30°), the boundary layer is attached on the suction surface of the aerofoil during the cycles, as the angle of attack is below the static stall angle. The C_N and C_T behavior is essentially quasi-static, with elliptical hysteresis loops, in the sense that departures from the static behavior are not particularly apparent throughout the cycle. At the inboard sections the flow is generally unsteady and becomes quasi-steady when approaching the blade tip. The loops of C_N and C_T remain close to the static value, except at r/R = 0.95 where the slope of the normal force coefficient is smaller than that of the 2D aerofoil. This is related to the strong 3D effects resulting from the tip vortex formations in the near wake on the flow, suggesting that special treatment should be given to the tip sections when modeling dynamic stall phenomenon by possibly using a smaller slope for the aerofoil data curve than that of the 2D static data curve. Under 3D flow conditions, a rotating blade experiences unsteady axial, tangential, and radial velocity components. The relative wind speed on a rotating blade increases with the radial locations and the blade exhibits two external forces during rotation, called centrifugal and Coriolis forces. The blade also exhibits an increase in dynamic pressure from the root to the blade tip [29]. For this wind speed (5 m/s), the flow is attached and the angle of attack is below the static stall angle and, hence, the chord-wise flow is more effective compared to the radial flow [30]. The lift and drag force coefficients obtained from the 2D wind tunnel can be, therefore, used to predict the aerodynamic forces on rotating blade sections with a reasonable level of accuracy.

At the higher wind speed (U_∞ = 15 m/s and Φ = 30°), the mean angle of attack decreases from root to tip, down to below 15.3° at r/R = 0.95. Compared to the hysteresis loops for wind speed 5 m/s, the shape and trend of the hysteresis loops for 15 m/s are wider and cycle-to-cycle differences are larger. Stall is delayed to a higher angle of attack than the static stall angle and, hence, dynamic effects became pronounced with significant C_N undershoots and overshoots. At the inboard sections, the maximum angle of attack greatly exceeds the static stall angle and the deviation from the static characteristic is evident. The impact of a dynamic stall vortex (DSV) on the upper surface of the aerofoil is also evident, which leads to a significant increase in the normal force coefficient hysteresis loops. When the DSV sheds into the wake, the aerofoil experiences a sudden drop in lift. The presence of a secondary peak in the normal force coefficient hysteresis loops at high angles of attack is commonly known to be a result of the shedding of the secondary vortex, as was shown in Choudhry et al. [31]. The flow reattachment process at the inboard sections is delayed to a higher angle of attack than the static stall angle, with generally a high level of unsteadiness. The kinematics of the dynamic stall vortex exhibits a different behavior from the cycle-to-cycle, leading to a significant variation in the pressure distribution on the blade surfaces. It is common that the blade at high wind speed experiences high angles of attack, especially at the inboard section of the blade. In addition, 3D stall delay due to rotation contributes to radial flow, due to the centrifugal force which forces the air to move along the radial locations of the blade. Under separated flow conditions, the radial movement of the separated flow is further impacted by Coriolis forces, which act as a favorable pressure gradient in the chord-wise direction. The latter process causes a reduction in the thickness of the boundary layer, a delay in the onset of stall, and an increase in the lift force coefficient to a high angle of attack [32]. Therefore, the lift and drag force coefficients obtained from the 2D wind tunnel cannot properly model the aerodynamic loads on a rotating blade, as the 3D effects are not included in these data.

5.4. Variation of Axial-Induced Velocities at Rotorplane

Figure 8 present the unsteady axial induction factor, a_x = u_ax/ U_∞, at the blade lifting line estimated from the measurements and the free-wake model at U_∞ = 5 m/s and U_∞ = 15 m/s. In each figure, the mean and the corresponding standard deviation estimated from the 30 rotor revolutions at each blade radial location and azimuth position are plotted. The maximum normalized standard deviations in a_x (σ_ax/μ_ax) at five radial locations are tabulated in

Table 2. The maximum normalized standard deviation in the axial induction factor at 5 and 15 m/s.

U∞	Maximum Normalized Standard Deviation (σax/μax) (%)
	Radial Location	Φ = 10°	Φ = 30°	Φ = 45°
U∞ = 5 m/s	r/R = 0.30	4.71	5.58	3.33
	r/R = 0.47	3.07	2.94	2.55
	r/R = 0.63	2.66	1.87	2.25
	r/R = 0.80	2.27	2.01	2.32
	r/R = 0.95	6.34	5.75	5.46
U∞ = 15 m/s	r/R = 0.30	14.92	15.99	15.50
	r/R = 0.47	28.88	13.05	10.06
	r/R = 0.63	39.11	25.63	10.71
	r/R = 0.80	31.08	14.32	7.33
	r/R = 0.95	22.67	17.61	13.60

. In the free-wake model, the induced velocity components are computed from the near wake trailing and shed circulation. Since the measured bound circulation at each time step and radial location is used as input value for the free-wake model, as previously explained in Section 4, the cycle-to-cycle variation in the normal forces and the bound circulations is also reflected in the axial-induced velocities on the blade lifting line.

5.4.1. Variation of the Mean Axial-Induced Velocities

At U_∞ = 5 m/s the mean axial-induced velocity variation at the blade lifting line follows a quasi-sinusoidal trend with azimuth angle. At the outboard sections of the blade, the induced velocities at 90° are higher than at 270°, given that, under the influence of the skewed wake formation, the outboard sections are, on average, closer to the strong tip vorticity at the 90° blade position. Trailing vorticity emerging from the blade root is also known to influence the flow field at the blades under yawed conditions [10]. The inboard blade sections are closer to the wake root vorticity at a blade angle of around 270°, hence causing a higher axial induction to be experienced by these sections at this blade angle than at 90°. At U_∞ = 15 m/s, the sinusoidal trend in the mean values of a_x is no longer present and is replaced by an irregular one under the influence of unsteady flow separation, resulting from large blade angles of attack. In fact, the inboard sections, which are more prone to dynamic stall, exhibit higher levels of unsteadiness in the circulation shed by the blades in the near wake region.

5.4.2. Cycle-to-Cycle Variations of the Axial-Induced Velocities

It is observed from Figure 8a–c that the cycle-to-cycle variation in the axial induction factor, a_x, along the blade lifting line is generally small at U_∞ = 5 m/s. This is noted from the error bars, which represent the standard deviations, which are small, mainly due to the fact that the flow around the blades is generally attached. The cycle-to-cycle variation is largest at the blade tip and root regions, where the largest radial gradients in the bound circulation (dГ_B/dr) are known to occur (Table 2). Consequently, any variations in loading due the turbulence in the wind tunnel speed, despite being very small (Section 2), would be amplified in these regions. Generally, higher values of the standard deviation in bound circulation imply higher cycle-to-cycle variations in the axial induction factor. Figure 9 plots the standard deviation in the bound circulation (${\sigma }_{{\mathsf{\Gamma }}_B}$) as a function of the radial location and blade azimuth position at yaw angle of 10°, 30°, and 45° and wind speeds of 5 and 15 m/s. In these plots it should be noted that, for the blade azimuth positions between around 50°–100°, the boom installed at the hub during the tunnel measurements resulted in a spike inboard of the blade, especially at high yaw angles, which may contribute to the larger standard deviations in a_x, presented in Figure 8. Excluding the latter effects, the relation between the variation in a_x (Figure 8a–c) and ${\sigma }_{{\mathsf{\Gamma }}_B}$ (Figure 9a–c) at 5 m/s is very clear, as it can be observed that the higher variation in ${\sigma }_{{\mathsf{\Gamma }}_B}$ at 10° leads to a higher variation in a_x. When the yaw angle increases, ${\sigma }_{{\mathsf{\Gamma }}_B}$ is reduced and, hence, a_x follows the same trend (see Table 2). At U_∞ = 5 m/s, ${\sigma }_{{\mathsf{\Gamma }}_B}$ is small, reaching a maximum value of not more than 0.12 at the blade tip and root regions. This explains why the largest standard deviations in a_x at this wind speed were obtained at these regions. At the higher wind speed (U_∞ = 15 m/s), however, the cycle-to-cycle variation in a_x is found to be significantly larger (Figure 8d–f and Table 2). Such variations are largest at the inboard sections, which coincide with the regions where significant cycle-to-cycle variations in the normal and tangential coefficients measurements were recorded.

The unsteadiness in a_x is found to be related to the cycle-to-cycle variations in the bound circulation experienced by the rotating blades (Figure 9d–f). Figure 9d–f shows that ${\sigma }_{{\mathsf{\Gamma }}_B}$ was small for blade positions between 100° and 240° at yaw angle 30° and 45° (Figure 9e,f) resulting in small standard deviations for a_x over this range. Given that the trailing circulation is directly related to the spanwise gradients in the bound circulation (Γ_T = dΓ_B/dr), it also follows that ${\sigma }_{{\mathsf{\Gamma }}_B}$ at a particular radial position (r/R) also influences the neighboring blade sections. Consequently, the high levels of ${\sigma }_{{\mathsf{\Gamma }}_B}$ observed at the mid-span locations at Φ = 30° (Figure 9e) also influence the cycle-to-cycle variations in a_x at the adjacent radial locations. A trend showing a decrease in the cycle-to-cycle variations in a_x with yaw angle could be also noted in Figure 8d–f and Table 2. In fact, ${\sigma }_{{\mathsf{\Gamma }}_B}$ was found to be significantly lower, on average, for Φ = 45° as compared to Φ = 10° and 30°. This is due to the decreasing of the angle of attack and the effective wind speed components when the projected rotor area is decreased with yaw angle.

5.5. Normality Test of the Randomly Fluctuating Aerodynamic Loads

All lifting line based models for representing blade loads rely on the use of engineering aerofoil data models which basically relate the load coefficients C_N and C_T (or C_L and C_D) to an angle of attack. It should be appreciated that developing improved engineering models to cater for cycle-to-cycle variability due to dynamic stall would be more complex if both the cycle-to-cycle variability in the load coefficients (C_N and C_T) and the angle of attack are large. The iterative approach adopted in this study (Section 4) enabled a detailed quantitative evaluation of the variability in the angle of attack over successive rotor revolutions induced by rotor yaw. It was evident that the cycle-to-cycle variability in the derived angles of attack was found to be very small, even for conditions when the corresponding variability in C_N and C_T is large. This may be observed from the plots in Figure 10, which show the variation in the normalized standard deviations for C_N, C_T and α with the blade azimuth position. This behavior was observed at all blade radial locations and operating conditions of the NREL Phase VI rotor considered in this study. Consequently, the probability distribution of the aerodynamic loads should not be significantly dependent on the probability distributions of the other aerodynamic quantities, like those for the angle of attack. This would considerably simplify the use of engineering modeling for simulation of the cycle-to-cycle variation in the aerodynamic loadings.

In this section, a statistical analysis is presented for the aerodynamics load coefficients based on 30 rotor revolutions. This made it possible to examine the statistical distribution expected for a given angle of attack and radial/azimuthal position. Normality tests were also carried out at 5 and 15 m/s and yaw angles 10°, 30°, and 45°, in order to establish the extent to which the cycle-to-cycle variability in C_N and C_T measurements at each rotor blade azimuth angle and radial location exhibits a normal probability distribution. This would further simplify the integration of stochastic modeling for wind turbine design tools. This analysis was based on different criteria and approaches: (1) the skewness and kurtosis z-values (−1.96 < z-value < 1.96) [33]; (2) the Shapiro–Wilk p-value (p-value > 0.05) [34]; and (3) the Normal Quantile–Quantile plots and Histograms. The data throughout these tests were binned per one degree azimuth angle. Hence, each bin includes a total of 30 data points.

Table 3. A summary of descriptive statics and Shapiro–Wilk’s test for normal force coefficient at U_∞ = 5 m/s, Φ = 30°, and r/R = 0.30.

Bin	Skewness		Kurtosis		Shapiro–Wilk (p-value)
	Statistic	z-value	Statistic	z-value
Bin 1 (10°)	−0.1889	−0.5425	−0.0913	0.3584	0.8515
Bin 2 (135°)	0.1516	0.4270	−1.0290	1.3758	0.2191
Bin 3 (210°)	−0.3923	−1.1200	−0.1894	0.4955	0.5583
Bin 4 (300°)	−0.5576	−1.6060	0.0089	0.2356	0.2191

and

Table 4. A summary of descriptive statics and Shapiro–Wilk’s test for normal force coefficient at U_∞ = 15 m/s, Φ = 30°, and r/R = 0.30.

Bin	Skewness		Kurtosis		Shapiro–Wilk (p-value)
	Statistic	z-value	Statistic	z-value
Bin 1 (10°)	−0.7162	−1.9380	1.1735	1.5162	0.0674
Bin 2 (135°)	0.0553	0.1545	−0.8802	1.2355	0.5298
Bin 3 (210°)	−0.2033	−0.5772	−0.6635	1.0425	0.1879
Bin 4 (300°)	0.0938	0.2705	−0.7593	1.1457	0.5298

show some selected samples at the azimuthal bins shown in Figure 11 at r/R = 0.30.

Figure 12 and Figure 13 illustrate the corresponding histograms and the normal Quantile–Quantile plots for each bin. When constructing histograms, a free important parameter is the number of bins. The number of bins used in a histogram is sensitive to the width of the bins. If the bin width is too large, the resulting probability distribution will then tend to be uniform, giving the wrong distribution, while using too small a bin width may lead to gaps in the resulting probability distribution, which is probably not representative of the correct distribution. For this reason, the optimal number of bins and the bin width should be selected carefully. The optimal bin width ∆* was calculated using the approach presented by Shimazaki and Shinomoto [34,35] as a minimizer of the formula [(2${\mu }_{C_N}$– ${\sigma }_{{}_{C_N}}^2$)/∆²], where ${\mu }_{C_N}$ and ${\sigma }_{{}_{C_N}}^2$ are mean and variance of C_N sample counts across bins with width ∆. The optimal number of bins is then calculated as [(max(X) − min(X))/∆*], where X is C_N sample data vector. The same also holds for C_T data.

It can be seen that, from Table 3 and Table 4, the Shapiro–Wilk’s test (p-value > 0.05) and the Skewness and kurtosis z-values (−1.96 < z-value < 1.96) agree well, identifying the same normal distributions. For example, it can be observed from Table 3 that all the skewness and kurtosis z-values are within the range (−1.96 < z-value < 1.96), while the results of Shapiro–Wilk’s test show that all p-values are greater than the significance level of 0.05. This hypothesis indicates that the data come from a population that can be assumed to be normally distributed. The same observation can be also seen in Table 4 for the higher wind speed of 15 m/s. The histograms and the normal Q–Q plots presented in Figure 12 and Figure 13 also corroborate with results of the Shapiro–Wilk’s test and the Skewness and kurtosis z-values as regards normality. It should also be understood that a management decision is issued regarding normality when more than two statistical tests agree with each other. The example in Figure 13 includes results of the Normal Quantile–Quantile plot of bin 3. For this bin, the skewness and kurtosis z-values, Shapiro–Wilk’s test p-value in Table 4 and the corresponding histogram indicate that the data are normally distributed. The normal Q–Q plots, on the other hand, show that the data do not seem to be normally distributed and, hence, such distribution is considered as a normal distribution.

For a complete assessment, the Shapiro–Wilk’s test was carried out for the full azimuth angle range and all radial locations at both U_∞ = 5 and 15 m/s and yaw angles 10°, 30°, and 45°. The results for C_N and C_T are presented in Figure 14 and Figure 15. The different colors shown in the plots indicate the probability distribution (p-value) as computed by the Shapiro–Wilk’s test. The percentage number of data points across the entire rotor disc failing the Shapiro–Wilk’s test is also shown in

Table 5. The percentage number of data points across the entire rotor disc failing the Shapiro–Wilk’s test.

U∞	Force Coefficient	Φ = 10°	Φ = 30°	Φ = 45°	Average
U∞ = 5 m/s	CN	10.47%	6.41%	9.25%	8.71%
	CT	7.02%	6.89%	11.28%	8.39%
U∞ = 15 m/s	CN	11.14%	10.74%	10.67%	10.85%
	CT	17.90%	13.51%	15.87%	15.76%

. This was computed at all radial locations using the formula:

$$\%\text{\hspace{0.17em}number\hspace{0.17em}of\hspace{0.17em}data\hspace{0.17em}bins}=\frac{\text{total\hspace{0.17em}of\hspace{0.17em}bins\hspace{0.17em}with\hspace{0.17em}}(p-\text{values}\le 0.05)}{\text{total\hspace{0.17em}number\hspace{0.17em}of\hspace{0.17em}bins\hspace{0.17em}}}\times 100$$

For wind speed 5 m/s, the average percentage of data bins (C_N data points) that are not normally distributed is less than 9%, while for wind speed 15 m/s it is less than 11% (see Table 5). In the case of the C_T data, the average percentage of data bins that are not normally distributed is less than 9% for wind speed 5 m/s and less than 16% for wind speed 15 m/s (see Table 5). It should be noted, however, that for data points failing the Shapiro–Wilk’s test, normality is still usually observed within ± two standard deviations. From Table 5, it can be also seen that the average of the percentage number of data points showing non-normality in C_N and C_T for the three azimuth angles at 5 m/s is approximately the same (8.71% for C_N and 8.39% for C_T). This trend is, however, not seen at 15 m/s, where it is noted that the average of the percentage number of data points for C_T having non-normal distribution is higher than this for C_N (10.85% for C_N and 15.76% for C_T).

As observed from the results presented in this subsection, the majority of the data points of C_N and C_T at various blade azimuth angles and radial locations (r/R) are normally distributed. This matter has so far not been addressed sufficiently in literature related to yawed wind turbine aerodynamics (see Section 1), where the analysis for stalled conditions has mainly concentrated on the mean values. The current paper demonstrates the following: (i) aerodynamic modeling of yawed turbines based solely on the average estimates at the different blade azimuth positions would neglect some important information about unsteady flow phenomena; (ii) the simulated wake downstream of the yawed turbine over a number of rotor rotations is considerably affected by the cycle-to-cycle aerodynamic loads on the blades, hence, it is important to understand the unsteady effects on the blades under stalled flow conditions which are self-induced by the rotor; and (iii) the characterization of the type of distribution for the cycle-to-cycle aerodynamic loads at a given mean angle of attack and radial location provides insight for determining the type of stochastic models (Monte Carlo methods) required to account for cycle-to-cycle differences in blade loads. The assumption that both C_N and C_T are normally distributed for a given angle of attack simplifies the stochastic modeling process significantly. Yet, the standard deviations under several operating conditions may still be required. Reasonable estimates may be derived from experiments or from a Computational Fluid Dynamics (CFD) code. In the present paper, the standard deviations of C_N and C_T for the different angles of attack and radial locations have been collected from the developed approach (Figure 2) for all radial locations, blade azimuth angles, wind speeds, and yaw angles, and arranged as shown in Figure 16.

6. Conclusions

This paper presents a numerical approach in which a free-wake vortex model and the blade pressure measurements are used to derive the unsteady angle of attack distributions over successive rotor rotations under yawed flow conditions. The new approach is able to capture the cycle-to-cycle differences in various aerodynamic parameters, which become significant at high angles of attack as a result of the complex 3D dynamic stall kinematics, which differ between one rotor revolution and the next. It is shown that, while cycle-to-cycle differences in the normal and tangential coefficients may be significant, the corresponding differences in the angle of attack are marginal. The maximum turbulence intensities for the tunnel wind speeds were found to be less than 4% and the effects of the tower disturbance are negligible, indicating that the main cause of the cycle-to-cycle variation in the aerodynamic loads is the local unsteady blade aerodynamics and the resulting unstable mechanism from one cycle to another. Larger yaw angles are found to reduce the cycle-to-cycle variations in the axial induction factors at high wind speeds. Normality tests have shown that the cycle-to-cycle variations in the normal and tangential force coefficients exhibit almost the same degree of normality at low wind speeds. Less than 9% of the analyzed data bins across the entire rotor disc was found to fail the Shapiro–Wilk test. For high wind speeds, the degree of normality in the tangential force coefficient was found to be lower than that for the normal coefficient. However, most of the non-normally distributed data, especially the normal force coefficients, follow a near-normal variation within ± two standard deviations, which may indicate that it can be reliably assumed that all data are normally distributed for engineering analysis. The present study was conducted on a model wind turbine subjected to a fixed speed and a uniform and steady wind flow in a wind tunnel. However, analysis of full-size turbines exposed to an atmospheric inflow with shear and turbulence/gust may result in different unsteady aerodynamic phenomena, due to variation of the wind speed and its direction, continuously with space and time. The results should provide motivation for developing improved aerofoil performance prediction models to include probabilistic models to account for the stochastic behavior under unsteady flow conditions. It is expected that such work would potentially reduce the uncertainty in aerodynamic load prediction in wind turbine design.