Comparative Study of the Aerodynamic Performance of the New MEXICO Rotor under Yaw Conditions

The influence of yaw misalignment on the aerodynamic performance of the New MEXICO rotor is investigated using blade-resolved Computational Fluid Dynamics (CFD) approaches with three wind speeds considered at a fixed yaw angle of 30 degrees. The air-load predictions and near wake characteristics from the numerical results are compared and discussed against the most recent wind tunnel test data. The nacelle impact, dynamic stall phenomenon and wake characteristics are analyzed, demonstrating the yaw effects and numerical issues raised from Reynolds-Averaged Navier Stokes (RANS) and Detached Eddy Simulation (DES) computations.


Introduction
A yaw situation occurs when a rotor is not perpendicular to the wind direction, which could happen to any wind turbine on every wind farm at any time. Even with a controlling system implemented, the rotor is not capable of following the continuous changing of the wind direction and hence spends much of its operating time misaligned with the wind. The air-loads on each blade vary in time when the turbine is yawed due to large fluctuations in the relative velocity at the leading edge. Results from both experimental and computational studies have shown that unsteady aerodynamics exist during all operating conditions and that dynamic stall can occur for high yaw angle operation while stall hysteresis occurs for even small yaw angles [1][2][3][4]. These unsteady aerodynamic phenomena are coupled with a long and flexible blade together with a tall tower for modern wind turbine, resulting severe structural dynamic response and hence wind turbine failures, reduced machine life and increased operating maintenance. Understanding the unsteady aerodynamic behaviour of the wind turbine in complicated yawed flow is essential to reliable and efficient designs.
Many experiments have been carried out on wind turbine aerodynamics. One is the well-known Unsteady Aerodynamics Experiment (UAE) Phase VI conducted by National Renewable Energy Laboratory (NREL) in the NASA Ames 24.4 m × 36.6 m wind tunnel [5,6], where comprehensive wind tunnel data of a two-bladed, stall-regulated horizontal axis wind turbine were obtained under a wide range of wind inflow speeds and yaw angles. Another is the ongoing Model Experiment in Controlled Conditions (MEXICO) project [7][8][9][10] carried out in Large Scale Low Speed Facility (LLF) of the German-Dutch Wind Tunnel (DNW), in which pressure and load data of a three-bladed model were measured under various inflow conditions together with detailed flow field data taken by Particle Image Velocimetry (PIV) technique. Several follow-up experiments, namely the Mexnext-II and Mexnext-III New Mexico project [11] have been performed to enlarge the yaw flow wind tunnel database. In addition, experiments have been carried out by Schreck [12,13] to explore the 3D dynamic stall processes on the NREL UAE Phase VI rotor through the analysis of blade surface pressure data and the local inflow angles when operating under yaw conditions. The wind tunnel experiments of these two wind turbines have accumulated detailed data, which provides the modeling basis and comparison reference for wind turbine aerodynamic computation and modeling.
The present numerical study is built on the recent yawed inflow wind tunnel experiments in Mexnext-III New MEXICO project [11]. Unsteady aerodynamic performance of yawed wind turbine under three wind speeds of 10, 15, and 24 m/s, is investigated using computational fluid dynamic (CFD) method. The yaw angle between the rotor axis and wind direction is fixed at 30 • . Full rotor geometry is resolved in the computations. Computational results of the rotor power and load variations are compared against the wind tunnel measurements. Sensitivity analysis of the influence of the grid resolution and nacelle effect is presented. Reynolds Averaged Navier-Stokes (RANS) and Detached Eddy Simulations (DES) computations are both employed, and evaluation regarding the adaptability and accuracy of the two methods has been made. The advancing and retreating effect and the flow physics of dynamic stall phenomenon on the blade inboard section are discussed.

Wind Turbine Geometrical and Operational Descriptions
The recent wind tunnel experiments in the MEXICO project, which was carried out in Large Scale Low Speed Facility (LLF) of the German-Dutch Wind Tunnel (DNW) with an open test section of 9.5 m × 9.5 m, involved extensive measuring of air loads, surface pressure and near wake flow field data. The wind tunnel blockage ratio was 18% and breathing slots were placed behind the collector to decrease the tunnel effect. The experimental model was a three-bladed wind turbine with a diameter of 4.5 m. The geometry profile of the blade is given in Table 1. The rotor was designed for an optimum tip speed ratio of 6.7 with the pitch angle of −2.3 • at the tunnel velocity of 15 m/s. The tip speed ratio is defined as λ = ωR/U ∞ , where ω is the rotational speed, R stands for rotor radius and U ∞ is the wind speed. The model was placed on a six component balance, measuring three forces and three moments at the foot of the tower. A total of 148 pressure sensors were instrumented and distributed at five blade spanwise positions 25%R, 35%R, 60%R, 82%R and 92%R, where R stands for the rotor radius. A zig-zag tape at 5% of the chord was used both on the blade upper and lower surfaces to trigger the laminar-to-turbulent transition. Phase locked stereo PIV measurements were performed at the 9 o'clock plane of the rotor for a variety of configurations and locations to characterize near wake flow. A complete description of the experiment setups is documented in the report of Boorsma and Schepers [11]. The nacelle geometry is simplified in this study. It is modelled as a cylinder with a diameter of 0.54 m and a length of 4 m. The nacelle nose and tail take the form of a tapered cylinder ending with a spherical cap.
In the most recent yawed condition MEXICO experiments, on which this study is based, the rotating speed of the rotor was fixed at 425.1 rpm and the blades were pitched with −2.3 • . Measurements were performed at the yaw angle of 30 • under the tunnel speeds of 10, 15 and 24 m/s, corresponding to tip speed ratios of 10, 6.7 and 4.2, respectively.

Grid Generation
Numerical investigations in this work are built on the yawed inflow wind tunnel experiments as mentioned above. Full computational domain is required in the transient yawed inflow simulation. Body-fitted mesh is applied in the present computations. The wind turbine is modeled without consideration of either zigzag or the tower in this numerical study. The sliding mesh technique is utilized to enable the rotation movement of the rotor. A pair of non-conforming axisymmetric interfaces is constructed between the static region and the rotating region. The structured multi-blocked computational meshes are generated using Pointwise. Mesh parameters are detailed in the following.

•
Rotational region: The rotational zone has a shape of a cylinder. The rotor rotates around the +x axis with the origin located at point (0, 0, 0). The radius of the rotating region is 1D, and the region extends from −D to 3D in the x direction, where D is the rotor diameter. The boundary mesh topology around the blade surface is O-type. The height of the first boundary layer cells is kept at 5 × 10 −6 m to ensure y + ≈ 1. • Stationary region: The outer stationary zone has a cuboid shape in order to ensure the inlet and outlet boundaries easy to specify. The computational domain ranges from [−5D, −5D, −5D] to [10D, 5D, 5D] in the x, y, z directions, respectively. The snapshots of the computational mesh are shown in Figure 1.

Numerical Methods
The tip Mach number of the rotor blade does not exceed 0.3, and consequently compressibility is neglected. The governing equations are the incompressible form of Navier-Stokes equations: where → U is the fluid velocity vector, p the density-normalized pressure, ν the fluid kinematic viscosity, and ν T the fluid turbulent viscosity.
With the inflow wind speed ranging from 10 m/s to 24 m/s, unsteady flow structures, such as massive flow separations and dynamic stall phenomenon, can be found near the blades. Therefore, the SST k-ω DES model [31] is employed in our study in addition to the commonly used Stress Transport (SST) k-ω RANS model [31]. The governing equations for the SST k-ω RANS model are described as follows: where: The value of the coefficients α and β can be found in [31].
Modification to the dissipation term in the k equation is made to switch the model from RANS to DES: where ∆ is the maximum local grid spacing (∆ = max(∆x, ∆y, ∆z)) and C DES = 0.61 is a calibration constant of the DES formulation. Numerical simulations are conducted based on OpenFOAM-v2.3.x, in which the finite volume method is used to solve the Navier-Stokes equations. The convective terms are discretized using a third-order Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme, and the viscous terms are discretized using central difference scheme. Temporal discretization of unsteady solution is treated through a second-order iterative time-stepping method. A uniform inflow velocity boundary condition is set at the inlet. A convective outflow boundary condition is applied at the outlet. Rotor blades and nacelle surface are treated as no-slip walls. The interpolation of the flux and variables between the non-conforming mesh interfaces is treated via an arbitrary mesh interface (AMI) projection technique [32]. The time step is set at 1.96 × 10 −4 s, equivalent to 720 steps per revolution. All the numerical results are sampled after 20 revolutions with the scaled residual of all equations below 1 × 10 −5 and time-averaged over the last five revolutions to ensure the convergence of the simulation. All our numerical simulations have been carried out using both RANS and DES models.

Results and Discussion
To facilitate the interpretation of the results, definitions of the two coordinate systems in yawed flow are specified in Figure 2a,b. The rotor rotates in clockwise direction when looking from the upstream. The blade azimuthal angle is defined to be zero when pointing to +z direction. The normal force on a blade section is normal to the local chord, positive pointing in downwind direction. The tangential force is parallel to the local chord, positive pointing from the trailing edge to the leading edge, see Figure 2c. The experimental data presented in the following discussion come from the MEXICO Experiment wind tunnel database from ECN [9,10] and the wind tunnel experiments described in the report of Boorsma and Schepers [11].

Grid Dependency
Four different resolution meshes as listed in Table 2 have been applied to investigate the influence on the load predictions. Computations are performed using the method described in Section 2.3 with turbulence modeled by SST k-ω RANS model. Wind inflow of 15 m/s is chosen to conduct this round of grid dependency study. Comparisons of the rotor torque and axial force between computations and experiments are shown in Figure 3. The air-loads on the blade experience a periodic change as the blade incidence flow angle varies with azimuthal position. The overall shape and fluctuation amplitude of the rotor torque and axial force are well predicted. Over-estimation can be observed comparing to the measured data, which is also been found in similar studies by other research groups using CFD methods. The reason for the discrepancy could be tunnel effects, which cause the velocity accelerated. It has been pointed out in the numerical study of Réthoré [24], Bechmann [25] and Shen [26] that normalized velocity one rotor diameter upstream is different for simulations and measurements in the aligned flow condition. In addition to that, recent study [33] has also show that zig-zag tape is one of the probable factor resulting large discrepancy between fully turbulent CFD results and measured ones.
In Figure 3, the largest deviation can be found between the coarse mesh case and the measured data. The mean levels of the axial force and torque are over-predicted by 10%. With the increase in rotor resolution, the deviations gradually decrease as more detail flow structures can be captured in the simulation. The average difference is reduced down to around 6% for medium grid resolution case, and further decreased to less than 5% for the fine mesh case. The reduction of the mean value is less than 1% when the grid resolution changes from the fine to finer one, indicating that, with respect to the integrated loads, the results can be considered to be grid independent.

Rotor Performance Validation at Different Tip Speed Ratios
Rotor performance under inflow wind speeds of 10, 15 and 24 m/s, corresponding to tip speed ratios of 10, 6.7 and 4.2, is investigated. Comparisons of torque and axial force variations with blade azimuth are presented in Figures 4-6. Large separations may occur when high wind speed is considered, hence DES is employed in this round of numerical computations. Air-loads of the rotor at the optimal inflow 15 m/s from RANS computations have already been presented in Section 3.1, and so we pay more attention to the other two "off-design" conditions. For inflow wind speed of 15 m/s, corresponding to tip speed ratio of 6.7, the results of torque and thrust predictions using the two turbulence models are very close. The over-estimations of axial thrust and rotor torque are both within 6% compared to the measured data, which is better than in previous studies. Even though the differences are small, the amplitudes of the force and torque fluctuations from DES results fit the experimental data a little better than the RANS results. For instance, the fluctuation amplitude of the axial force with azimuthal angle calculated from RANS is 11 N, while this result from DES is 17 N closer to the experimental result of 15 N.  The scenario changes for the "off-design" state of a lower wind speed of 10 m/s, corresponding to high tip speed ratio of 10.0. In Figure 5, the axial thrust predictions are in good agreement with the measurements. The maximum over-estimation in RANS result is within 3%. DES result is slightly improved over RANS with 0.5 percentage reduction in maximum error. Strong discrepancy can be observed in torque prediction compared to measurement, although the relative error of the thrust prediction remains the same level as that at the wind speed of 15 m/s. The predicted torque is about 20% higher than the experimental result. At this low wind speed the blade inflow angle is rather low, and therefore the drag plays a more important negative role in torque production. It is reasonable that drag prediction may not be accurate enough using either RANS or DES technique, resulting in notable relative error in torque prediction, despite the fact that its absolute error is decreased compared to that in the design state. In this case, DES results were not significantly improved compared with RANS results.
For the wind speed of 24 m/s, flow separations are expected to occur as the blade operates at high incidence angles. The difference between RANS and DES results is significant, as expected. In Figure 6, predictions of total force and torque from RANS computations have a maximum over-estimation of 30% compared to the measurements. On the other hand, DES gives a better prediction of the error level below 20% compared to the wind tunnel data, although this is still unsatisfactory, at least to some extent. The dynamic stall phenomenon, the three-dimensional rotational augmentation and advancing and retreating effect all together pose more difficulties in rotor load predictions.

Nacelle Effect
The influence of nacelle on the rotor aerodynamics and near wake flow tends to be neglected in axial flow because the flow is axisymmetric and the nacelle frontal area is relatively small compared to the rotor swept area. According to the numerical study carried out by Sørensen [34], the presence of the nacelle would increase the power and the axial load by approximately 1-2%. He also pointed out that the scales of the blade and nacelle in the wind tunnel experiment were not in proportion to a full-size wind turbine, and therefore the effect is much larger than it should be expected for a full-size turbine [34]. However, when yaw misalignment is present, the flow state changes from axisymmetric to non-asymmetric. The impact area of nacelle is not confined to the frontal area, and it is possible that the nacelle interacts with the turbine wake. Under this circumstance, the presence of nacelle is more likely to bring more impact on the near wake flow than the axial inflow condition.
In order to evaluate exactly how much influence nacelle would bring on the rotor aerodynamic performance and near wake flow field, a trial of computations is conducted with and without the nacelle. The wind speed of 15 m/s is chosen in this case. The exact same resolution grids and numerical scheme set-ups are applied. Comparisons of the loads from the two cases are shown in Figure 7. The mean level of total torque gets a 3% increase when nacelle geometry is considered. Meanwhile, the axial force is less sensitive and the presence of nacelle geometry only causes a less than 2% increase.  To elaborate on the nacelle influence on the near wake flow field statistically, radial and axial profiles of velocity components are compared against the wind tunnel data obtained using PIV technique [8,11] in Figures 9-11. The radial profiles of averaged velocity components extracted parallel to the rotor plane at x = −0.3 m (upstream) and x = 0.3 m (downstream) in the model coordinate system are presented in Figures 9 and 10. Results in both cases are very close and agree well with the measured data in the outboard section. Discontinuous peaks in the center area of no nacelle case occur due to the root vortices. The radial profiles of velocity components do not vary much whether the nacelle geometry is considered or not.   The velocity distribution in the axial direction is sampled at y = +1.5 m with respect to the wind tunnel system when the blade is pointing straight up. The influence of nacelle geometry on the axial profile of the near wake can be easily observed in the range between x = 2 m and x = 4 m in Figure 11. The results of the velocity components when nacelle is considered are in good agreement with the experimental data. In the meantime, large deviations can be observed in the "nacelle impact zone" between simulated results without nacelle and the measured data, especially the over-estimation of velocity component u. Therefore, for the future study of rotor flow field under yawed inflow, nacelle presence needs to be included.

Advancing and Retreating Phenomenon
In contrast the axial inflow, there is a contribution of the wind in both x and y directions when the rotor is positioned in yawed inflow. Contour of pressure on the suction side of rotor blades at different azimuthal angles is shown in Figure 12. The influence of the extra y direction contribution on the rotor blade is changing due to the rotational motion, and consequently increases or decreases the velocity in chordwise direction periodically. This advancing and retreating phenomenon is characterized by a minimum and maximum at 0 • and 180 • azimuthal angles, respectively. At these two positions, the flow is either fully aligned with or opposite to the rotational speed. The local incidence angle reaches the highest and lowest at 0 • and 180 • azimuthal angle. From Figure 12, we can observe that the flow is mainly attached at 180 • and almost fully separated at 0 • . The change of incidence angles further induces dynamic stall. On the other hand, the extra y-contribution also has a periodic impact on the flow in the radial direction. The radial contribution either weakens or enhances the centrifugal pumping caused by rotational augmentation. This effect is important for blade inboard sections, where the inflow speed is relatively large compared to the rotational speed at low TSR. This effect reaches the maximum and minimum at the azimuthal angle of 90 • and 270 • , respectively. In addition to that, the axial induction on the downwind side becomes larger than the upwind rotor side resulting in an additional deflection beyond the wind direction to the skewed wake (marked by the two dashed lines) due to yawed flow, see Figure 13.

Dynamic Stall
Dynamic stall is expected to occur in the blade inboard section when the turbine is yawed with respect to the wind [8]. Large load fluctuations are known to be imposed on the airfoil section when compared to the static loading characteristics. The wind speed of 24 m/s is chosen for the dynamic stall investigation. Three typical cross sections with two airfoil types are discussed here, namely the DU91-W2-250 airfoil at the 25% and 35% span and the RISØ A1-21 airfoil at the 60% span.
Comparisons of normal and tangential forces at the three spanwise positions are presented in Figures 14-16. Strong discrepancies can be easily detected between RANS and DES results especially in the downwind side azimuth range from 0 • to 180 • . The same pattern can also be found in tangential force predictions between 0 • and 180 • . At the 25% span, despite the good agreement in the range of 0 • to 90 • , RANS result fails to reproduce the drop in normal force prediction compared to the measurements, and strong discrepancy can be observed after the drop as the RANS prediction rises to the maximum, which is almost 20% higher than the measured value. On the other hand, the general shape and mean level from the DES prediction are closer to the measured data. The fluctuation-related difference between the measurements and DES results may be related to the data filtering approach, although the reason for this difference still remains open. As for the tangential force comparison of the 25% span, it is obvious that DES results are closer to the measured data than URANS results. RANS results fail to reproduce either the overall shapes or the mean level. At the 35% span positions, the fluctuations of the normal forces when the incidence angle is retreating from the maximum value are not predicted in RANS computations. DES method, by contrast, is able to give a better estimation of both mean value and fluctuations. However, when it comes to the 60% span location, the agreements of both RANS and DES results with the measured data are deteriorated. The reason for this still remains open, but it may relate to the three-dimensional vortex shedding phenomenon at this particular mid-span section due to the blade design [35].   The lift and drag characteristics of the three sections extracted from the simulations have been compared with the 2D case provided by Mexnext MEXICO project [7][8][9][10]. The exact angle of attack (AOA) of the blade cross section is known to be difficult to determine, and various approaches have been proposed, such as free wake vortex [36], inverse BEM [37,38], 3-Point method [39], Shen model [40]. Evaluation of these methods can be found in [41]. The sectional lift and drag coefficients are derived using the inverse BEM theory through the normal and tangential coefficients obtained from CFD simulation. The amplitude of the changes in AOA is high at the inboard sections and gradually decreases toward the outboard region.
The DU91-W2-250 airfoil is used at the 25% and 35% span positions. It can be observed from the 2D CL plots (see Figure 17) that the airfoil stalls at the critical angle of attack of 25 • , where the lift has a sudden drop and the drag increases significantly. The changes in angle of attack at the 25% and 35% span positions are very large during one rotational cycle in yawed inflow. At the 25% span, rotational augmentation causes the lift coefficient at 10 • to be higher than 2D case. The lift coefficient increases linearly up to the maximum value of 2.2 at a 3D stall angle of attack about 32 • , and then it drops with the increase of AOA. The three dimensional drag coefficient is much larger than the 2D result with no sudden increase at the critical angle of attack. Similar patterns can be found in the lift and drag characteristics at the 35% span location in Figure 18. The increase of lift coefficient is also linear and the hysteresis loop at the 35% span is wider than that at the 25% span. The RISØ A1-21 airfoil is used at the 60% span. The measured 3D lift demonstrates a much different scenario from both the computations and 2D static lift line, see Figure 19. There is no sudden drop in 2D lift when the stall angle is reached. Neither the RANS nor the DES results agree with the measured dynamic data, they both give much larger predictions and wider loops. The reason of these differences is still not well determined, but they may be related to the blade design and a special aerodynamic behaviour at this particular spanwise station, where a mid-span vortex exists [35]. This requires further investigations.

Near Wake Characteristics Analysis
The predicted flow field velocity profiles in the near wake region are discussed in the following sections. The velocities and coordinates discussed here are given with respect to the wind tunnel coordinate system. Instantaneous velocity components u, v, w in the wind tunnel coordinates at different positions along the wind direction are obtained via wind tunnel PIV data. Axial profiles of velocity components from numerical computations at two positions are compared against the measured data, one is at y = 1.5 m when blade azimuth is 0 • and the other one is at y = −1.5 m when blade azimuth is 60 • [8,11]. All the velocity components are extracted on the horizontal plane at the height of the rotor center. Unfortunately, this round of wind tunnel measurements did not provide much data for the two off-design state (10 m/s and 24 m/s) and the wake deficit after 1D. Hence, the analyses of the axial velocity distribution are primarily concentrated on the very near wake region.
For the axial velocity distributions at the windward side y = 1.5 m in Figures 20-22, good agreements can be observed between the computational results and measured data before the wake structure interacts with nacelle geometry. The obstruction of the nacelle between x = 2 m and x = 4 m can be clearly observed in the numerical results. After that, strong discrepancy can be found between RANS and DES results in all the three wind speeds because of the unsteadiness in the interacting wake.   While in the leeward side of y = −1.5 m, the influence of the nacelle obstruction disappears. For optimal wind speed of 15 m/s, other than the overestimation of w-velocity component in RANS computation, the axial profiles of the computational velocity components agree well with the measurements, see Figure 23. At the wind speed of 10 m/s corresponding to the tip speed ratio of 10, a strong discrepancy can be found in Figure 24 between the RANS and DES results, probably due to the quick dissipation of RANS model leading to a quick breakdown of the tip vortices. At higher wind speed of 24 m/s, larger deviations can be found among the experimental results and different numerical results in Figure 25. It seems from the velocity components comparison that DES result fits the experimental data better than the RANS one. For example at x = 1 m, a very near wake position, the error of the u component prediction using DES is reduced almost by a half compared to that using RANS. The similar improvement can be found in the comparison of v and w components. The DES model shows advantage in capturing the evolution of tip vortex structure and the breakdown.

Conclusions
A detailed investigation of the effects from the yaw misalignment has been presented based on the New MEXICO rotor using blade-resolved CFD approach. Three inflow wind speeds are considered at a fixed yaw angle of 30 • . Two types of numerical computation approaches, i.e., RANS and DES, are used in the study. Sensitivity studies on the grid dependency and nacelle effect have been performed. The simulation results indicate that yaw misalignment has significant influence on aerodynamic torque, thrust behaviour. Complex flow phenomena have been observed in the wind turbine flow field. The following conclusions are obtained:

1.
Good agreements can be found in the total force and thrust prediction in optimal wind speed of 15 m/s with the maximum error below 6%. Improvements have been made in off-design state of 10 m/s and 24 m/s.

2.
The presence of nacelle causes a 3% increase in rotor torque and a 2% increase in thrust.
Interactions between the blade tip/root vortices and nacelle are demonstrated. The axial profile shows less wake deficit of the mainstream velocity component when the nacelle obstruction is not taken into consideration.

3.
Results of the normal and tangential forces at the 25%, 35% and 60% spanwise stations have been validated against the experimental data, and the three dimensional lift and drag characteristics have been analyzed to demonstrate the influence of dynamic stall on the rotor load predictions at the wind speed of 24 m/s in yawed flow. 4.
The near wake velocity components show a good agreement with the measurement. The DES model shows its advantage in both blade load predictions and wake characteristics.