{'fraction_non_alphanumeric': 0.023832945805759256, 'fraction_numerical': 0.008138079055625112, 'mean_word_length': 4.623585617299472, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Investigation of the current yaw engineering models for simulation of wind turbines in BEM and comparison with CFD and experiment

The aim of this work is to investigate the capabilities of current engineering tools based on Blade Element Momentum (BEM) and free vortex wake codes for the prediction of key aerodynamic parameters of wind turbines in yawed flow. Axial induction factor and aerodynamic loads of three wind turbines (NREL VI, AVATAR and INNWIND.EU) were investigated using wind tunnel measurements and numerical simulations for 0 and 30 degrees of yaw. Results indicated that for axial conditions there is a good agreement between all codes in terms of mean values of aerodynamic parameters, however in yawed flow significant deviations were observed. This was due to unsteady phenomena such as advancing & retreating and skewed wake effect. These deviations were more visible in aerodynamic parameters in comparison to the rotor azimuthal angle for the sections at the root and tip where the skewed wake effect plays a major role.


Introduction
Wind turbines in the atmospheric boundary layer (ABL) are practically operating in fluctuating wind, which leads to a misalignment of the velocity vector with respect to the rotor plane. This effect is commonly known as yawed flow effect. Modelling of wind turbines in yaw is still one of the major design challenges of wind turbines [1].
One of the most common approaches to simulate the aerodynamics of wind turbines is the Blade Element Momentum (BEM) theory. Low computational costs for simulations, which is a result of this steady and two dimensional theory, make the simulation affordable even for more than thousand load cases. Nevertheless, complex three dimensional flows can not be captured accurately by the basic BEM theory. Several studies show that codes based on BEM with different correction models for yawed inflow are often not sufficiently accurate and reliable for predicting yawed flow acting on wind turbine blades [2,3].
Yawed flow reduces the effective projected area with respect to the wind direction and therefore in most cases reduces the energy extracted from the wind. Additionally, yawed flow causes two unsteady aerodynamic phenomena. First, it causes the advancing and retreating blade effect: In case of positive yaw, the blade will advance in the lower half of the rotor plane and retreat in upper side. This leads to 1p variation of the angle of attack and the effective inflow velocity [4]. In this case, the horizontal component of the wind velocity is no longer negligible and should be considered for computing the angle of attack. This effect has a cosine type and its symmetric around the 180 degree azimuth (blade pointing down in vertical position). Therefore, the maximum thrust will occur at this point due to the higher effective velocity. Moreover, advancing and retreating mainly occurs at high wind speeds i.e. low tip speed ratio and it's more significant at the inboard area of the blade since for the tip the rotational speed is dominating the horizontal component of the wind velocity.
Second, there is the skewed wake effect: Due to asymmetric position of the wake (root and tip free vortex) relative to the rotor plane the blade in downwind side will be more into the wake and therefore the induction over the rotor plane will vary. This will produce a lower load for the downstream blade and hence a yawing moment [5]. There were various researchers who developed models for this effect. The first model was proposed by Glauert [6]. He found an equation to correct the induction factor at the blades based on a yaw angle. This equation uses the wake skew angle defined as the angle of the actual flow leaving the turbine and the rotor axis [7] a = a 1 (1 + k r R sinφ).
In equ. (1), k is a function of the radial station r, R rotor radius and φ the azimuthal position of the rotor blade. Many attempts were made in the past to improve the Glauert model by more accurate modelling of the wake expansion and deflection, see [8,9,10,11]. However, dynamics of these models show discrepancies between each other [3]. They mainly assume the variation in induced velocity to be purely sinusoidal, which is not always the case. The root vortex also induces axial velocities and this causes deviation at the root to be different and not purely sinusoidal [2]. Schepers [4] used experimental results to develop a new model which predicts the dependency of the induction variation on the azimuthal position. The main difference in comparison with the previous models lies in inclusion of the effects of the root vortex [7]. An oversimplification of the aerodynamic phenomena governing the power extraction at yawed conditions can lead to large errors in the design process and reduces the efficiency of wind turbines. Therefore, it is essential to conduct more realistic simulations in order to capture all the important details, which play a role in modelling of the wind turbines correctly. The current work has the following contributions to the correct modelling of yawed inflow in wind turbines: • The limits of current correction models will be shown by comparing BEM simulations with other numerical methods and also experimental data. • A detailed investigation for different unsteady phenomena such as advancing and retreating blade effect and skewed wake effect will be presented and explained. • A proposal for a new correction model will be presented.

Test Cases and Methods
In this paper three different turbines namely, NREL VI [12], the AVATAR [13] and INNWIND.EU [14] were investigated under uniform inflow conditions. The computational grids were created by the BladeBlockMesher code of ForWind and Fraunhofer IWES [15]. The meshes are shaped spherically and fully structured based on hexahedral cells. For the NREL VI turbine, the total grid size is 22 million cells. It has 300 cells around the airfoils, 250 cells in the span-wise direction for each blade, and 150 cells in the wall normal direction. The y + values at the surface are kept below 1 everywhere on the blade surface. For the AVATAR and INNWIND turbine, a total grid size of 30 million cells is used, having 240 cells around the airfoils, 380 cells in the span-wise direction for each blade, and 120 cells in the wall normal direction. For all the cases the effects of tower and nacelle are neglected and computations are performed in a uniform flow condition. To examine the grid-independent behaviour of the solution Grid Convergence Index (GCI) study were used for all the grids presented here [15].
All the CFD calculations are performed using the open-source CFD software OpenFOAM [16] with the Spalart-Allmaras-DDES model [17].
In order to solve the pressure-velocity coupling, the PIMPLE algorithm which is a combination of the loop structures of SIMPLE [18] and PISO [19] is used. The second-order linearUpwind scheme is used for discretizing the convective terms. To perform the numerical simulations, the facility for large-scale computations in wind energy research (FLOW) of the University of Oldenburg was used. In total 264 cores are used, resulting in approximately 96 hours for NREL VI turbine and 120 hours for the AVATAR and INNWIND turbine. The physical time step for NREL VI is ∆t = 2 × 10 −4 and for the AVATAR and INNWIND ∆t = 2 × 10 −3 . The converged numerical result is achieved after 5 rotation.
The AVATAR and INNWIND turbine are simulated using ECN's BEM code and free vortex wake code AWSM [22]. For the ECN BEM code, the skewed wake effects are modeled using the Glauert [6] and Schepers model [4], in AWSM the skewed wake effect is modeled intrinsically by the wake definition.

Results
In this work, first NREL VI turbine is investigated for inflow velocity of 5 m/s and 30 degrees of yaw. The numerical results in terms of normal and tangential force coefficients from BEM using different models will be compared with CFD for two radial sections of r/R = 30 and 95 %.
Next, the numerical results for AVATAR and INNWIND turbine using a BEM code, free vortex wake code and CFD are compared for both axial and yawed flow in terms of mean value of axial induction factor, thrust, and tangential force, and also the axial induction factor, thrust and tangential force vs azimuthal position for several sections along the span.

NREL VI
In Figure 1 the normal and tangential force coefficient vs azimuthal angle are presented. The CFD results show a reasonable concurrence with the experimental data for both normal and tangential force coefficients, for both sections presented here. Though, the aerodynamic properties of larger rotor could be diverse from small sized rotors, still CFD can be used as a reference to compare the BEM results for larger rotors where experimental data is not available.
In terms of the normal force coefficient at 30% span, all three correction models predict the maximum forces at the same azimuthal position. The maximum is in downwind side of the rotor (plane). It should be noted that NREL VI turbine is rotating counter clockwise, therefore the upwind side is at 90 degree of azimuthal position. The Pitt & Peters model shows a smaller amplitude for the maximum value and predicts a pure sinusoidal behaviour. This can be due  to the fact that Pitt & Peters model does not take the root free vortex effect into account. The GDW model is underestimating the forces at the upper half of the rotor blade significantly. The model from Schepers has a closer agreement to the experimental results in terms of amplitude and phase. The result for the tangential force coefficient using these three correction models show more deviation from the experimental results. All the correction model results are shifted towards the downwind side of the rotor (plane). The model from Schepers has a closer agreement to the experiments since this model considers the root free vortex effect as a function of blade radius.
Comparing advancing and retreating with skewed wake effect, it can be concluded that, in this section of the blade, both effects contribute to the load variation. This can be seen from the fact that the maximum forces are shifted towards the downwind side of the rotor (plane).
For the section at 95% of the span, for both normal and tangential force coefficients, experimental and CFD data show sinusoidal behavior which comes from tip vorticity. The Pitt & Peters model predicts different phase and amplitude, and the GDW model predicts the correct phase but not the amplitude. The Schepers model also miss-predicts the position of maximum force and underpredicts the forces on average. This could be due to the fact that, for this model, the effect of root vortex is still present at the tip.
Comparing advancing and retreating with skewed wake effect for the tip section shows that the contribution of advancing and retreating effect to the force coefficient is less than that of the skewed wake effect. This can be seen from the fact that the maximum forces are shifted to the upwind side of rotor (plane). This leads to stabilizing yaw moment. The other possible source of error for this section is the tip correction model. The capability of these correction models in prediction of the loads in the yawed flow needs to be investigated.

AVATAR and INNWIND
For the AVATAR and INNWIND turbines the CFD calculation will be compared with the results from AWSM and BEM code of ECN. BEM calculations are performed with two different skewed wake correction models, namely Glauert [6] and Schepers [4]. For the purpose of verification and validation, aerodynamic key parameters, such as induction factor and local forces at different span-wise positions along the rotor span, are compared.

Axial Induction Factor
Several methods have been introduced so far for calculating the induced velocity at the rotor blade in axial condition; these range from theoretical models, such as inverse BEM [23], to the numerical models, such as Average Azimuthal Technique (AAT) [24,25] and the model from Shen [26]. In the following section, 4 different methods are implemented and investigated: • Inverse BEM Technique: The computed load from CFD is used to estimate the induction and AoA using general BEM theory formulation [23]. In Figure 2 the azimuth averaged axial induction factor, a, is plotted as a function of radial position for both the AVATAR and INNWIND turbine for 0 degrees of yaw. The first conclusion is that BEM and AWSM are in good agreement almost everywhere along the span. However, the CFD methods used for calculating the axial induction factor show deviation from each other. The 3-point technique and inverse BEM are in good agreement with each other and closer to the results from BEM and AWSM. The method of Shen and AAT show 15% less induction from all other methods. Moreover, the trend of CFD results at the tip is different from the results of BEM and AWSM. BEM and AWSM show an increase in induction, however CFD methods show a decrease. This indicates that more research is still needed for the calculation of the CFD induction factor. One possible solution would be to conduct more experimental data in this topic.
In Figure 3 the azimuth averaged axial induction factor is plotted as a function of radial position for both the AVATAR and INNWIND, for 30 degrees of yaw. For this case, the 3-point technique predicts 5% more induction than all other calculations. At the very tip, BEM and AWSM results still deviate from the CFD results. This brings up the question of accuracy of tip correction models for BEM models. In Figure 4 the axial induction factor vs azimuth positions for section 60 and 90 % of the span are shown. In both sections, AWSM present the closest results in terms of the location of the maximum and also phase to the CFD calculations. For the mid-span section, the skewed wake model of BEM-Schepers and BEM-Glauert deviate from each other. BEM-Schepers predict a maximum in the upwind side of rotor (plane), which is in contradiction with the results from CFD and AWSM. BEM-Glauert demonstrates purely sinusoidal behaviour, which is not the case for other methods. In the section close to tip, the BEM-Glauert and BEM-Schepers predict less induction than AWSM and CFD. In the downWind side of rotor (plane) (around 90 degree), the trailed tip vortices of the skewed wake are closer to the rotor plane and therefore increase the induced velocites. In this section, the skewed wake effect plays a major role. This effect is much more visible for the INNWIND turbine, since it operates with a much higher induction factor. Figure 5 the azimuth averaged normal and tangential force per unit length for axial flow is shown. In terms of thrust force along the rotor span, all the codes are coincide for both the AVATAR and INNWIND turbine. However, there is still deviation in the codes at the very tip (95% of span). This deviation was also observed when comparing the averaged axial induction factor in Figure 2, where CFD results where predicted much lower induction at the very tip. This brigs up the question of accuracy of the tip correction model in BEM, as well as the 3D effects on the polars used by BEM and AWSM. In terms of tangential force for the AVATAR and INNWIND turbines, deviation is observed between CFD and the results of AWSM and BEM code, for up to 30% of the span. This might be due to the 3D effects. For the INNWIND turbine, there is also a deviation at the tip. At mid span, the CFD results predict 10 to 15% more force than the other two methods. This deviation might be related to the input polars of BEM and AWSM. It is worth mentioning that, since the INNWIND turbine is a high induction turbine, the forces are much higher than the AVATAR turbine.

Normal and tangential force In
In Figure 6, the azimuth averaged normal and tangential force per unit length for 30 degrees of yaw is shown. As the case of axial flow, the codes seems to concur in terms of thrust force except at the tip, where CFD predicts around 10% lower. In case of tangential force for the AVATAR turbine much lower forces are predicted up to 40% of span. For the INNWIND turbine, CFD and BEM codes coincide along the span except the sections close to the tip, where CFD predicts much lower forces than BEM. The predication of forces from AWSM showed much lower results than CFD and BEM, which still requires investigation.
In Figure 7, the thrust force is shown over a rotor revolution for 3 different sections (R=40,60,90%) along the rotor span. For section R=40% the skewed wake model of Schepers predicts a very similar phase and amplitude as AWSM for the downwind side of the rotor and also the location of maximum force. In contradiction, the skewed correction of Glauert predicts the maximum force on the upwind side of rotor (plane). This observation is also valid for the INNWIND turbine, where the maximum for the BEM-Schepers model and AWSM are on the downwind side, and the Glauert model on the upwind side. The CFD results, in both case, predict the maximum close to the results of AWSM and BEM-Schepers. This could be due to the fact that the BEM Schepers and AWSM models include the root vortex, while BEM-Glauert does not. For the R=60% the maximum force is shifted for CFD, AWSM, and BEM-Glauert towards the upwind side, however the BEM-Schepers still predicts the maximum on the downwind side. In terms of amplitude and phase, the AWSM has closer agreement with CFD, while BEM-Glauert predicts much higher loads at the upwind side of rotor (plane). For the section at the very tip (R=90%), the prediction of BEM-Schepers and BEM-Glauert are very similar to each other; they both predict much higher force than AWSM and CFD. As the section before, the results of AWSM are very close to CFD in terms of amplitude and phase.
In Figure 8, the tangential force is shown over a rotor revolution for 3 different sections (R=40,60,90%) along the rotor span. The first conclusion is that the AWSM results in terms of phase and location of maximum force are very close to the results of CFD. As is the case of thrust in the 40 % section the phase and the location of maximum force for BEM-Schepers and BEM-Glauert are different. Due to the fact that BEM-Schepers models the presences of the root vortex, the results are very similar to CFD and AWSM. At the mid span BEM-Schepers mispredicts the location of maximum force, and also the phase. At the very tip the results, of BEM-Schepers and BEM-Glauert are very close to each other. The qualitative behaviour of the induced velocity for BEM-Schepers is slightly poorer since the effect of the root vortex is still present for this section.

Conclusion and future work
In this paper, the performance of different engineering BEM based code on predicting the aerodynamic parameters of wind turbine in yawed flow were investigated. The aerodynamic key parameters, such as axial induction factor and sectional loads, were compared with CFD and experimental data (when available). From the presented results, the following conclusion can be made: • In general concurrence was observed when comparing the mean value of aerodynamic quantities, such as axial induction factor and blade loads, for axial and yawed flow by all the codes. However at the very tip, there was still deviation between the results of BEM and AWSM, compared to the CFD method, with regard to axial induction factor. This brings up the question of accuracy for tip correction models for BEM models. • From four different methods used to calculate the axial induction factor, the 3-point method has an excellent agreement to the results from AWSM for axial and yawed condition. Underprediction with respect to other results were observed using the AAT and Shen methods. Although the authors didn't find more than 10% dependency of the results to the position of monitoring point for Shen and 3 point methods, this still needs to be investigated in more detail. • Sectional results over a rotor revolution showed a significant difference between different codes. In the root area, the skewed wake model from Schepers coincides more with the results from CFD and AWSM, since this model takes the root vortex into account. Hence, the azimuthal variation of the induced velocities at the root predicted by this model deviates from the sinusoidal behaviour, induced by tip vortices. Glauert skewed wake correction models are demonstrate the axial induction factor and loads with a sinusoidal behaviour which is in contrast with the results of CFD and AWSM. It should be noted that the deviation from sinusoidal behaviour is due to presence of the root vortex. By CFD, AWSM and BEM Schepers maximum velocity predicted is induced on the upwind side and the minimum velocity is induced at the downwind side. However, since by the Glauert model velocity is induced by a tip vortex wake, only the maximum induced velocity is reached at the downwind side of the rotor plane. • In the mid span where the advancing and retreating effect is dominant reasonably similar agreement was observed between different codes. ECN-Schepers is slightly poorer since the effect of the root vortex is still present for this section. The advancing and retreating effect leads to destabilizing yawing moment. • Close to the tip the behaviour is closer to the sinusoidal behaviour however some deviation is still observed, due the effect of advancing and retreating. The qualitative behaviour of the induced velocity for ECN-Schepers is slightly poorer since the effect of the root vortex is still present for this section. This offers room for the improvement of ECN's yaw model in a way that the radial dependency of the parameters should be changed, so that the root vortex effects are damped out at the tip [27]. • In the future, these CFD calculations will be used for extracting the axial induction factor for different TSRs, yaw angles, sections along the span and azimuthal positions. Afterwards, a fourier series analysis will be performed to find an empirical model which covers the variations in the wake induced axial velocity in the rotor plane as a function of rotor (plane) azimuth angle, radial position, tip speed ratio and yaw angle. The new model will cover the root and tip vortex, and is therefore expected to be more accurate than the previous ones.