Comparison of aerodynamic models for horizontal axis wind turbine blades accounting for curved tip shapes

Curved tip extensions are among the rotor innovation concepts that can contribute to the higher performance and lower cost of horizontal axis wind turbines. One of the key drivers to exploit their advantages is the use of accurate and efficient computational aerodynamic models during the design stage. The present work gives an overview of the performance of different state-of-the-art models. The following tools were employed, in descending order of complexity: (i) a blade-resolved Navier Stokes solver, (ii) a lifting line model, (iii) a vortex-based method coupling a near-wake model with a far-wake model, and (iv) two implementations of the widely used blade element momentum method (BEM), with and without radial induction. The predictions of the codes were compared when simulating the baseline geometry of a reference wind turbine and different tip extension designs with relatively large sweep angle and/or dihedral angle. Four load cases were selected for this comparison, to cover several aspects of the aerodynamic modeling: steady power curve, pitch step, extreme operating gust impact, and standstill in deep stall. The present study highlighted the limitations of the BEM-based formulations to capture the trends attributed to the introduction of curvature at the tip. This was true even when using the radial induction submodel. The rest of the computational methods showed relatively good agreement in most of the studied load cases. An exception to this was the standstill configuration, as the blade-resolved Navier-Stokes solver was the only code able to capture the highly unsteady effects of deep stall.


| INTRODUCTION
The design of modern horizontal axis wind turbines poses several challenges, due to their significant rotor size, that could be overcome through innovative design concepts. 1When focusing on the blades, tip extensions are seen as one efficient solution to optimize the performance of the machine, while keeping the loads under the design requirements.Of particular relevance to the present work are the add-on modifications of the blade geometry through the installation of site-specific tip extensions.This represents a modular way of altering existing blade geometries, in a retrofit fashion, and several studies have anticipated significant benefits after the installation of the tips.Some of these studies simply consider straight tip extensions in order to, for instance, extend the lifetime of a rotor. 2 However, more interesting seem to be the improvements attributed to curved tips, such as winglets 3,4 or complex shapes accounting for both sweep and dihedral angular deflections. 5,6These geometries can be designed to improve the performance of the machine, through for instance the steering of the aeroelastic tailoring.However, the consideration of these highly curved tip shapes during the design stage of wind turbine blades, which often relies on fast engineering aerodynamic models, may result in suboptimal designs.
In other words, if the curvature effects are not properly taken into account in the simulations, the resulting geometry could underperform with respect to the baseline design.Along these lines, recent work has extended engineering aerodynamic models applied for wind energy, that traditionally assume straight blades and a flat rotor plane, to account for both dihedral angles 7,8 and sweep angles. 9,10The extended models showed improved agreement with results from higher fidelity numerical methods, namely, a blade-resolved Navier-Stokes code and a lifting line free vortex wake solver, when comparing the sectional loading variations attributed to the curvature of the geometry.For the particular case of swept tips, a validation campaign of a wide range of numerical models was also recently conducted, where a representative tip shape was studied in a wind tunnel. 11e present work aims at characterizing the performance of several state-of-the-art solvers for the prediction of aerodynamic loads of wind turbine blades.In order of descending complexity, these solvers can be listed as follows: 1.A finite volume blade-resolved Navier-Stokes solver, run for different rotor-resolved meshes.This method is labeled in the present document as CFD, which is the acronym for computational fluid dynamics 2. A lifting line aerodynamic model coupled to a hybrid filament-particle-mesh solve to model the flow, labeled as LL.
3. A computationally efficient vortex-based method coupling a near-wake model with a far-wake model, referred to as NW.

Several implementations of the commonly used blade element momentum model (BEM).
More in particular, the focus is put on the implications that curved tip shapes may have on the blade and rotor loads.Indeed, while engineering models have been found to give accurate loading predictions for straight blades, 12 the implications that sweep angle and dihedral angle may have in a blade design context are still a research subject.Additionally, it is also known that the accuracy of all models relying on airfoil data is compromised in deep stall conditions. 13The present study aims at performing a first assessment of the implications of these factors in the aerodynamic modeling of wind turbine blades accounting for complex tip shapes.For that, a selection of representative curved tip shapes was used, and a range of load cases was studied in detail: • Steady power curve: operational conditions, assuming uniform axial flow.Used to analyze the influence of the curved tip shapes on the static loads.
• Pitch step: perturbation of a converged state through a fast pitch variation.Used to study the time evolution of the aerodynamic loads.
• Gust impact: modeling of an extreme operating gust wind profile.Used to study the influence of inflow dynamics on the loading.
• Standstill: considering the blades to be pitched, and under a large yaw misalignment, so that the flow impinges the blade from the pressure side.
Used to highlight the implications of deep stall on the modeling of aerodynamics.
The IEA 10 MW offshore reference wind turbine (RWT) 14 was taken as a baseline, and its tip was subsequently modified to account for representative sweep and dihedral angle variations.To ensure a purely aerodynamic comparison, all the simulations considered the structure to be stiff.It should be noted that neglecting the blade elasticity and considering a reduced number of load cases do limit the possibilities of concluding on the performance of the selected blade tips.However, this objective was not pursued in the present work, as the focus was put strictly on the aerodynamic modeling.

| METHODOLOGY
In this section, the methodology followed in the present work is described.The considered curved tip shapes are first introduced in Section 2.1.
The different load cases are then described in Section 2.2, both in terms of inflow conditions and operational settings.Finally, Section 2.3 presents the computational methods employed, together with the numerical models that were built.

| Blade and tip geometries
In the present study, four tip extensions of the IEA 10 MW RWT blade, here labeled as BASELINE, were considered.All the tip extensions had the same projected length on the pitch axis (i.e., 5% of the BASELINE blade projected length starting at 97.5% and ranging until 102.5% projected length).The blade tip geometries were designed through a parametrization based on the Akima cubic interpolation, through three spanwise points. 15The aforementioned points were located at (i) 2.5% projected length on the pitch axis inboard of the baseline tip, (ii) at the baseline tip, and (iii) at the extended tip.In order to dimension the tip extensions, a reference chord c ref was used.c ref corresponded to the chord at 99% of the total radius of the baseline geometry.In this way, the chord of all the tip extensions, at the location of the baseline tip, was set to be c ref .The chord at the tip of the extension was defined as 0:5c ref , resulting in a more slender tip geometry compared to the baseline.The twist distribution was the same as the BASELINE in the overlapping part and constant at 0 in the extended part.The relative thickness of the tip extension profiles remained the same as the baseline tip profile.The four considered tip extensions represent the parametric space expected for curved tip shape design optimization.These tip extensions were characterized by smooth shapes, representing specific geometrical design concepts: • EXT: a straight tip extension of the blade.It is used as a basis to understand the effects of introducing curved tip shapes, as its projected length on the pitch axis is the same as the other designs.
• SWEEP: a tip extension only accounting for backwards sweep angle.It is used to study the effects of in-plane geometrical modifications on the blade performance.
• DIHE: a tip extension only accounting for a dihedral angle, in the opposite direction of the blade prebend.It is used to study the effects of outof-plane geometrical modifications on the blade performance.
• COMBINED: accounting for the geometrical modifications of both SWEEP and DIHE.It is used to understand their combined effects and implications on the numerical modeling.
It should be noted that the resulting blade geometries were not the outcome of any detailed design process, so that their optimality from an operational point of view was not ensured.A summary of the blades considered in the present work is included in Table 1, and a visualization of the tip geometries is depicted in Figure 1.The reader is referred to Appendix A for more details about the geometries, where plots of the distributions of the following quantities are included: (i) blade centerline coordinates, (ii) twist, (iii) chord, and (iv) relative thickness.

| Load cases
The different load cases considered for the present study, which were run for every geometry and every computational method, are listed below.
Note that, in all cases, the inflow was assumed to be laminar (i.e., with a turbulence intensity level of 0 %.).

| Steady power curve
Simulations assuming a uniform wind speed, a constant rotational speed and a fixed pitch angle were performed.An overview of the selected operating points for this load case is shown in Figure 2.These values correspond to the steady state operating conditions included in the IEA 10 MW RWT definition document 14 and were used for every numerical method of the present study.

| Pitch step
Two pitch motions were taken into account, a nose-up and a nose-down, which, respectively, accounted for an increase and decrease of the angle of attack and therefore the loading.In particular, a change of ±2 was imposed over 0.1 s, starting at a time where the solution of the different models had reached steady conditions.Three particular wind speeds were considered for this load case (i.e., 6, 10, and 14 m/s), with the rotational speeds and initial pitch angles as compiled in Figure 2.

| Gust impact
Starting from a stabilized simulation of the 10 m/s operating point included in Figure 2, an extreme operating gust was introduced.Its wind speed profile was defined according to the IEC 61400-1 standard, 16 and it was aligned with the rotor axis.An amplitude of 62.1 m/s and a duration of 10.5 s were considered, as depicted in Figure 3.For all the numerical models included in the present work, any time shifts corresponding to the convection of the gust were completely disregarded, so that the changes of the inflow were instantaneously seen by all the blade sections.

| Standstill
For this case, the rotor was braked and the blades were pitched by 90 .Only the blade pointing directly upwards was studied.In addition to that, a large yaw misalignment of 90 was considered, so that flow impinged the targeted blade from the pressure side.A single wind speed of 25 m/s was simulated.

| Numerical models
This section introduces the different numerical methods and describes the turbine and blade models that were built in order to simulate each of the load cases listed in Section 2.2.It should be noted that, to simplify the comparison, both the rotor tilt and the blade coning were neglected in all the performed simulations.Tower shadow effects were not accounted for in any of the models.
All the described numerical methods, apart from the blade-resolved CFD solver, required a set of airfoil polar data as an input.The polar data defined in the IEA 10 MW RWT report 14 was used for this purpose.The FFA-W3 airfoil data were computed using the two-dimensional finite- Parameters characterizing the different simulation points included for the steady power curve.Includes the resulting tip speed ratio for the BASELINE configuration volume CFD solver EllipSys2D v16.0. 17Fully turbulent flow was assumed, and the data were 360 extrapolated using the pre-processor AirfoilPreppy, 18 with no three-dimensional corrections.A single Reynolds number of 10 million was assumed.

| Blade element momentum BEM
The results labeled as BEM in the present study correspond to the model described in Madsen et al, 19 which is implemented in the in-house finite element multibody aeroelastic code HAWC2. 20This aerodynamic solver includes, among others, submodels for dynamic inflow, yawed inflow, and unsteady 2-D airfoil aerodynamics.The casting of the BEM problem into a polar grid allows to model turbulent inflow conditions.
As described in Madsen et al, 19 in the present study care is taken to ensure that the thrust forces for the induction computations are correctly computed as thrust per unit radius, instead of per unit span.In addition, a radial induction model based on the lateral induction for a 2-D actuator disk was included as an option. 19This last feature was found to be particularly relevant for the present work, as it is developed to take into account part of the effects attributed to the out-of-plane geometries. 7That is the reason why the simulations of the Steady power curve were also run without any radial induction and labeled as BEM0 in this study.For the particular case of Standstill, where the rotor is braked, the induction model was completely deactivated.Therefore, these simulations corresponded to the application of the blade element model with the geometric inflow (i.e., disregarding the coupling with the momentum theory), and they are referred to as BE.The unsteady airfoil aerodynamic model 21,22 (usually referred to as the dynamic stall model) was used for all load cases, including the Steady power curve.As shown by Li et al, 8 directly using quasi-steady 2-D airfoil data will introduce errors for blades with out-of-plane geometries, even for steady state calculations.

| Near wake and vortex cylinder NW
The coupled near-and far-wake model [23][24][25][26] is a computationally efficient vortex-based method.In the near wake, which is defined as the first quarter revolution, the trailed vortices from each blade are modeled using helical trailed vortices.The indicial function approach is used, and the computational effort for each iteration remains constant throughout the simulation.The far wake is modeled using a far-wake BEM method that is without the tip-loss correction.The near wake model and the far wake model are coupled together using a coupling factor, in order to get the total induction. 25,27The coupling factor is calculated so that the rotor thrust is comparable to the one computed with the BEM method.The near wake model was recently modified to model the blade sweep effects, 10 by also accounting for the influence of the curved bound vortex. 28In this work, an earlier version of the modification of Li et al 10 is used, as its final implementation was not available during the generation of the results.In particular, another set of influence coefficient tensors is employed, which was pre-calculated based on a smaller design space of the blade sweep geometry.While this resulted in a slightly lower accuracy when compared to Li et al, 10 the regeneration of a few selected results indicated that the use of the different set of coefficients does not compromise the conclusions of the present work.A vortex cylinder model is also employed 7 to account for the influence of out-of-plane geometries, here implemented as a correction of the far-wake induction model.The system closure described in Branlard and Gaunaa 29 was adopted.As for the BEM method, the unsteady airfoil aerodynamic model 21,22 is included for all load cases.The described near wake and vortex cylinder coupled model is labeled as NW in the followings of this work.
F I G U R E 3 Wind speed time series for the studied gust profile

| Lifting line LL
Mid-fidelity simulations were performed with the in-house multi-fidelity vortex solver MIRAS. 30,31The investigation is carried out with a lifting line (LL) aerodynamic model coupled to a hybrid filament-particle-mesh solver to model the flow. 32The rotor blades are therefore modeled as discrete vortex segments that account for the bound vortex strength and act releasing vorticity into the flow in the form of straight vortex filaments.
In the hybrid method, the filaments are transformed into vortex particles, which vorticity is later on interpolated into an auxiliary Cartesian mesh.
The motion of the vortex elements is determined by the velocity of their markers, calculated by a superposition of the free-stream velocity, the velocity induced by the blade bound vortex, the filament wake, and the particle wake.The flow is governed by the vorticity equation, obtained by taking the curl of the Navier-Stokes equation.It describes the evolution of the vorticity of a fluid particle as it moves with the flow.Moreover, MIRAS has been recently modified to accurately account for blade curvature effects. 9Øye's dynamic stall model 33 is used to account for the dynamics of separated flows subjected to oscillatory motions.A 10R Â 4R Â 4R Cartesian mesh is employed in all cases, being R the blade radius of the BASELINE geometry.The mesh has a constant spacing of 5 m in the three directions, which adds up to a total of 1.28 million cells.The bound vortex is discretized with 80 straight segments with a constant spacing.The leading segments of the bound vortex rings are placed along the blade quarter chord line, with the calculation point defined at the center of such segments.The length of the first row of trailing filaments is defined by the velocity computed at the filament end points.For all the cases simulated for the present study, a time step size of 0.04 s was used.
A free boundary condition is used in all directions.Moreover, an eight-order stencil is used to spatially discretize the vorticity equation; a particle re-meshing is forced every time step to maintain a smooth field in combination with a periodic re-projection of the vorticity field, which is applied to maintain it divergence free.

| Blade-resolved Navier-Stokes solver CFD
For the simulations labeled as CFD, the incompressible Navier-Stokes pressure-based solver EllipSys3D is used. 17,34,35EllipSys3D is a threedimensional finite volume code that solves the equations on multiblock structured grids.For the present work, an automatized process is followed to generate a series of body-fitted meshes for each geometry.This ensures comparable grid qualities and was achieved by first creating a surface mesh based on the blade center line and the airfoil coordinates with the PGL library. 36In a second step, a hyperbolic volume grid is generated with Hypgrid. 37Depending on the targeted load case, a different meshing strategy is followed, and the assumptions taken for the solution of the Navier-Stokes equations are also adapted: • Steady power curve: These simulations are based on three-bladed rotor-resolved meshes.Every grid accounts for 14.2 million cells, and the reader is referred to Li et al 7 for more details about their topology.A Reynolds-averaged Navier-Stokes (RANS) turbulence model is used, based on k-ω SST. 38Pitch step and Gust impact: The same rotor meshes used for the Steady power curve are employed for the simulations of these load cases.However, due to the unsteadiness of the problem, the Navier-Stokes equations are solved following an unsteady RANS (URANS) dual time-stepping approach.A time step size of 4.9 ms is used for the Pitch step load case and of 3.45 ms for the Gust impact.
• Standstill: An improved delayed detached eddy simulation (IDDES) turbulence model is employed 39,40 for these simulations, with a time step size of 1.5 ms.For the region involving the RANS method, the k-ω SST turbulence model was also used.The body-fitted grids generated for the analysis of the Standstill case only consider a single blade, and account for a total of 12.6 million cells.More details about the mesh discretization can be found in Horcas et al. 41 All the aforementioned CFD models assume fully turbulent flow and follow a wall-resolved approach.Prior to the generation of the CFD results based on the described setups, sensitivity studies concerning mesh refinement and time step size were conducted.

| RESULTS
This section summarizes the results of the numerical campaign performed in the present study, comparing both the integrated rotor loads and the sectional blade loading.While the former quantity describes the overall performance of the wind turbine, the latter gives additional information on the structural implications for design and allows to better characterize the differences between the numerical methods.The sectional loading per meter span is presented as function of the normalized projected length on the pitch axis, in percentile.This quantity is referred to in this document as Norm.radius [%].
The same methodology was used for the post-processing of the results from all aerodynamic solvers.They were launched with the in-house Python-based framework DTU coupling. 42This framework coupled the sectional loading per span (normal to the center line of the blade) of the different aerodynamic models, with a stiff model built in HAWC2.This allowed to establish a common set of sensors, defining the projection of the loads in a consistent manner, thus limiting potential errors from data manipulation.

| Steady power curve
Figure 4 depicts the CFD predicted aerodynamic thrust, torque, and power of the BASELINE rotor operating under the different conditions presented in Figure 2. Results for the other fidelity models are presented in Figure 5.These are expressed as percentile relative differences compared to CFD, which is the reference model for the benchmark.
The LL method overestimates the rotor thrust by up to 3.5% when the wind speed is above 15 m/s and underestimates the rotor power by 6% at the wind speed of 4 m/s.For the other conditions, the LL method general has good agreement with the CFD.Regarding the BEM and the NW predictions, both solvers overestimate the aerodynamic thrust and power in below-rated operational conditions, i.e., wind speeds lower than 10 m/s.For example at 4 m/s, both methods overestimate the rotor thrust and power by approximately 4.5% and 17%, respectively.At this low wind speed, the turbine is running far above the design tip speed ratio at very high thrust coefficient above 0.9.This is a challenge for both BEM and NW, because the empirical high thrust correction of the C T ðaÞ curve (i.e., thrust coefficient versus induction), which is described in section 2.2 of Madsen et al, 19 becomes active.In above rated conditions, BEM has a general good agreement with the CFD predictions.NW overestimates the rotor thrust, especially at a wind speed of 25 m/s, where the rotor thrust is overestimated by 5.5%.The agreement of NW with CFD regarding the rotor power is slightly better than BEM below rated wind speed, and slightly worse above.
The capabilities of the different codes to model curved geometries (i.e., SWEEP, DIHE, and COMBINED designs) are addressed in what follows.
Figure 6 shows the difference in % in thrust and power predictions, relatively to the results obtained for the EXT configuration by the same method.To isolate the effects of the blade tip geometry, only the influence of the aerodynamic loading in the outer half of the blade has been considered.Regarding the thrust, all the methods show similar tendencies for the three geometries.There is a general thrust increase in the first half of the wind speed range (4 to ≈15 m/s) and a decrease at high wind speeds.At below rated conditions, the SWEEP geometry has the smallest Integrated rotor loads found with the CFD aerodynamic model, as a function of the wind speed F I G U R E 5 BASELINE configuration.Relative differences found for the integrated loads of every computational method, with respect to the rotor loads obtained for CFD thrust increase while the COMBINED geometry has the highest one.At higher wind speeds, the codes predict a thrust decrease.This is because the loading close to the tip at high wind speed is negative due to the large blade pitch angle.The same changes in tip geometry that cause an increase in thrust at low wind speeds will therefore lead to a decrease in thrust at high wind speeds.The trends are not as clear in terms of power prediction.In below rated conditions, there are large differences between the power computed by each method, and the results are changing abruptly with the wind speed.For above-rated wind conditions, the three lower fidelity methods generally predict lower power differences, when compared the CFD.
In what follows, a more detailed comparison of the solvers is carried out, focusing on the sectional loads of the outer half of the blade and starting from the BASELINE geometry.Again, the predictions for the different models are presented as relative differences compared to the CFD results, which are depicted in Figure 7. Three operational conditions with the wind speed of 10, 12, and 14 m/s are chosen for the comparison.As shown previously in Figure 5, the integrated rotor loads from different methods are in generally good agreement for these operational points.
Therefore, the chosen cases are expected to better assess the capability of the different solvers to capture the correct blade loading.The relative differences of the BASELINE sectional loads on the other models, compared to CFD, are visualized in Figure 8.The LL method is in good agreement with CFD, up to 90 % of the radius.That is also the case for the NW code, even if some noticeable discrepancies are also found at lower radii.
The BEM method slightly overpredicts both the in plane and out of plane loads.Closer to the tip, there is a clear overestimation of the loads by the three methods, especially for the NW.The aforementioned discrepancies are likely related to the 3D flow effects near the tip, which were not correctly modeled, as 2D airfoil data were employed.Additionally, the relative comparison is also negatively affected by the fact that close to the tip the blade loading tends to 0 (as previously shown in Figure 7).
In what follows, the relative percentile variation predicted between each of the curved shapes and the EXT straight extension is presented, for each of the considered aerodynamic models.Besides the already introduced methods, i.e., BEM, NW, LL, and CFD, results for the BEM0 model (i.e., without the radial induction correction) are included.This is to benchmark the performance of this correction for different curved blade geometries.The predicted blade loading for the SWEEP configuration is presented in Figure 9, for wind speeds of 10, 12, and 14 m/s.For all three operational conditions, BEM and BEM0 predict almost identical results.For the out of plane loading, there is an excellent agreement between the NW, LL, and CFD up to the last few percent of the radius.On the other hand, the lower fidelity BEM is not able to capture the change on the blade loading.Regarding the in plane loads, LL and NW follow the trends established by CFD relatively well, for the particular wind speed of 10 m/s.With the increase of wind speed, the two methods still predict similar trends, but those deviate from the CFD results.The trends computed by BEM deviate from the other methods.Moreover, the BEM trends in the relative differences for the in-plane loads and the out-of-plane loads align, because only the effect of the geometrical projections of the sectional velocity and loading was captured.
F I G U R E 6 Difference on the integrated loads for every tip geometry (outer half of the blade), when passing from the EXT geometry to every other tip variant.The results are normalized by the loads found for the BASELINE configuration, when using the same computational method The blade loading for the dihedral design DIHE is presented in Figure 10.All codes show similar trends in terms of wind speed evolution for both out of plane loads and in plane loads.However, both BEM and NW underestimate the increase of the loads, when compared to CFD.For 10 m/s, LL predicts almost identical results as BEM and NW.For the higher wind speeds, LL predictions appear mid-way between the lower and higher fidelity codes, capturing similar trends as the CFD solution.BEM0 overestimates the increase of the load for all operational conditions compared to CFD.Especially for the wind speed of 10 m/s, the trends of the BEM0 in plane loads differ from the results of the other methods.
The blade loading for the COMBINED case is presented in Figure 11.For the out of plane load, the predictions of LL, NW, BEM, and BEM0 have similar trends.BEM predicts a lower value of the relative difference near the tip, especially for 12 and 14 m/s.NW predicts similar results to LL, except for the tip-most region at 14 m/s.The relative difference predicted by BEM0 is overestimated, when compared to LL.The change in out of plane loads predicted by CFD is larger than any other lower fidelity model, with the exception of BEM0.For the in plane loads, LL predicts similar trends to CFD but tends to underestimate the relative difference at the blade tip, especially for 14 m/s.NW is able to predict some of the trends, at Sectional loads computed with the CFD numerical method for the BASELINE geometry.Legend entries correspond to selected wind speeds, in [m/s] Relative difference of sectional loads predicted by every method with respect to the CFD solution.Wind speeds of 10, 12, and 14 m/s are included.Only results for the outer half of the blade are shown least for 10 m/s, but it further underestimates the relative increase of loads at the tip.On the other hand, BEM predicts different trends for all conditions, when compared to the rest of the methods.The same comment applies to BEM0 for the studied conditions, with the particular exception of 14 m/s.
F I G U R E 9 SWEEP configuration.Percentile differences of the sectional loading at the vicinity of the tip, with respect to the results obtained for the EXT geometry by the same method.BEM0 shows the effect of disregarding the radial induction of the BEM method F I G U R E 1 0 DIHE configuration.Percentile differences of the sectional loading at the vicinity of the tip, with respect to the results obtained for the EXT geometry by the same method.BEM0 shows the effect of disregarding the radial induction of the BEM method This section focuses on the different model performance under a step change in the blade pitch angle, both positive and negative.First, a detailed investigation of the transient phenomena is shown for the BASELINE case.Later, results for all five blade designs at three different wind speeds are presented, in order to assess the influence of the blade geometry and wind conditions in the maximum attained loading.Note that in the fig- ures, the start of the pitch step was set to be at 0 s and the simulations were carried out for 30 s more.
Figure 12 shows the time variation of the difference in rotor thrust and mechanical power, with respect to the steady values before applying the pitch variation, and for the BASELINE blade design.At first sight the results of all the involved methods are consistent.There is an excellent agreement between BEM and CFD loads after 30 s, and a similarly good comparison was found for the peak loading of LL and CFD.Regarding the time evolution of the quantities, it is observed that LL and CFD show very similar trends in all cases, especially in the power output.On the other hand, the BEM and NW computations move more quickly towards the new equilibrium after the initial overshoot.The most obvious difference between the models is that only the LL and CFD simulations capture the rapid oscillations in the loading, attributed to the effect of vortex passing.
This is an expected outcome since only these two solvers can model the downstream development of the vorticity released by the blades and their mutual interaction.This effect was previously observed in IEA Wind Task 29 Phase III, 43 where the CFD and vortex codes predicted a clear staircase shape of the power and thrust time series, after the initial over-or under-shoot of a stiff turbine.The overshoots caused by blade pitching can have a role in the estimation of the fatigue loading during the lifetime of a wind turbine.For this reason, the range of the rotor thrust time series was computed, for each of the performed simulations (Figure 13).Generally, the higher the model fidelity the larger the predicted range value, with the models aligning fairly consistently for both nose-up and nose-down pitching motions at wind speeds of 6 and 10 m/s.At a wind speed of 14 m/s, there is a trend change with CFD predicting the lowest range and followed by NW, BEM, and LL.Regarding the influence of the blade geometry, it is observed generally that the BASELINE configuration shows the lowest range values, while a similar plateau is seen for all four tip extensions.CFD predicts the largest differences among the considered tip configurations.Still, the different tip configurations behave very similarly, which indicates that changing the tip shape only has a very limited effect on the dynamic inflow behavior.

| Gust impact
This section summarizes the results of the different numerical models when simulating an idealized gust, which is introduced in the simulations once a converged state is reached.Figure 14 depicts the time series of the integrated loads, which were made relative to the corresponding Percentile differences of the sectional loading at the vicinity of the tip, with respect to the results obtained for the EXT geometry by the same method.BEM0 shows the effect of disregarding the radial induction of the BEM method converged values, for the computations based on the BASELINE geometry.Only the period of time when the wind turbine loads are affected by the impact of the gust is shown.For clarity, the time is set to 0 at the start of this time window.Small differences (i.e., below 2%) are identified for the maximum and minimum loads.A fraction of these differences can be attributed to the discrepancies found for the different solvers in steady conditions (see Section 3.1).It should be highlighted that all the methods predict the same time evolution pattern.As commented when introducing this particular case, this good matching is affected by the fact that an idealized gust was modeled, as it was understood as an instantaneous change on the incoming wind speeds for all the numerical methods.
Similar trends are found for the rest of the considered geometries.This is summarized in Figure 15, which includes the range of the integrated loads time series for all the performed simulations.It is worth mentioning that all the geometries accounting for a tip extension predict similar values, irrespective of the numerical method employed.These values are also relatively small, when compared to the BASELINE configuration.An exception to these observations is the particular case of CFD, where noticeable deviations for every curved tip shape are predicted.However, the differences introduced by these geometries remain relatively low, even for this higher fidelity solver.Consistently to the study of the pitch step (Section 3.2), these results indicate that the instantaneous dynamic variation of the inflow can be captured in a similar way by all the involved numerical methods.Moreover, at a certain extent, these variations are independent from the dihedral and sweep angle modifications considered in the present study.

| Standstill
This section presents the results for the standstill case, which assumed the studied blade to be pitched and under a large yaw misalignment, so that the flow came from the pressure side.This particular load case was only simulated with the higher fidelity blade-resolved Navier-Stokes solver (CFD) and with a version of the BEM method where the induction was deactivated (here labeled as BE).The momentum part had to be neglected, because the disc assumption that is the basis for the BEM formulation is not valid in standstill conditions.Vortex-based models have not been employed in this study since, as pointed out by Chattot, 44 their use in stall conditions can cause numerical instabilities when the local angle of attack is larger than the maximum lift angle.As shown by the cited author, this could be avoided by the use of an artificial viscosity term.
Figure 16 depicts the sectional loads.
For brevity, only the direction aligned with the inflow (i.e., the in-plane) is shown.In the same line of thought, the BASELINE and COMBINED configurations were selected, but similar observations could be made for the other geometries involved in the present study.It is important to highlight that only the CFD method was able to model any type of unsteadiness for this particular load case.Indeed, the only source for F I G U R E 1 2 BASELINE configuration at a wind speed of 10 m/s.Relative variation in rotor loads with respect to the stabilized values before pitching, for every computational method.The abscissa values were shifted so that the pitch step starts at t = 0 s unsteadiness for this standstill condition case was the vortex shedding, which is illustrated in Figure 17 through the visualization of the BASELINE wake, and that is not modeled in the BE engineering model.However, the relative differences of the mean loading predicted by the blade-resolved solver CFD and the engineering model BE remain relatively low.
A detailed study of the vortex shedding and its potential effects on, e.g., standstill vibration is outside of the scope of the present study, and the reader is referred to a published work from the same authors. 41However, it is interesting to illustrate the implications of such phenomenon on the sectional loading of CFD.This is done in Figure 18 by means of the time series at the tip vicinity.trends would be captured no matter the complexity of the numerical method.However, it should be stated that the selected load case disregarded any effect attributed to the convection of the flow, which could eventually lead to larger discrepancies.
• Standstill: This case considered a blade pointing upwards, under deep stall.The CFD method was able to capture the unsteadiness of the flow related to vortex shedding.This effect, which could potentially lead to vibrations, was completely disregarded when applying an engineering method based on airfoil data.However, such a simple approach was found to give relatively accurate predictions of the mean loading.
Overall, the presented results highlighted the importance of selecting the right level of fidelity of the aerodynamic model, when simulating wind turbine blades with curved tip shapes.In that regard, the usage of BEM-based methods, even when including the radial induction method, can fail to model the effect of curvature.That seemed to apply to both static and dynamic load cases.The presented LL and NW implementations were found to better capture the trends imposed by the dihedral and the sweep angles, when compared to the results of the blade-resolved CFD.
While a rigorous study of the computational effort of each of the numerical methods was outside of the scope of this project, it should be noted that both LL and especially NW required much less resources than CFD.That makes this type of approaches very attractive to, for instance, their integration in an optimization design loop.
Future work could be devoted to expand the considered geometrical space and to assess the performance of the methods for more abrupt curved shapes.Along these lines, it could be interesting to carry out a similar comparison around winglet shapes, or around tips accounting for very large backward or forward sweep angles.
APP E NDIX A: BLADE DESCRIPTIONS

1
Visualization of the tips included in this study.The freestream velocity vector in normal operating conditions U, is included for reference.Left: isometric view; center: side view (seen from the left cheek of the nacelle); right: front view (seen from the depicted inflow)

F I G U R E 1 4
BASELINE configuration.Relative variation of integrated loads with respect to the stabilized values before introducing the gust, for every computational method.The depicted time window has the maximum gust speed at the center of the plot F I G U R E 1 3 Range of the rotor thrust due to the pitch step.Values normalized based on the stabilized value of every simulation (i.e., before pitching) and shown in percentile.The considered geometries are included in the abscissa, where B refers to the BASELINE geometry, E to the straight tip extension, S to the tip accounting for sweep, D to the tip accounting for dihedral, and C to the tip extension combining both dihedral and sweep anglesF I G U R E 1 5Range of the integrated loads due to the extreme operating gust.Values normalized based on the stabilized simulations (i.e., prior to the impact of the gust) and shown in percentile.The considered geometries are included in the abscissa, where B refers to the BASELINE geometry, E to the straight tip extension, S to the tip accounting for sweep, D to the tip accounting for dihedral, and C to the tip extension combining both dihedral and sweep angles F I G U R E 1 6 Standstill load case.Sectional in-plane loads distribution for CFD and the engineering model labeled as BE.For the former numerical method, a range of results is shown, which limits are the mean value AE one standard deviation, for the last 6 s of simulations.The results of the BE method are simply plotted by a line, as no unsteadiness was observed in the loads F I G U R E 1 7 Wake visualization of the BASELINE simulation.CFD method for the standstill load case.Iso-surfaces of Q-criterion,45 for a value of 0.5, at the end of the simulation.Colored by the sign of the vertical vorticity (positive values in red, negative values in blue).The direction of the freestream velocity U was included for reference

Figure
Figure A1 depicts the blade centerline and the airfoil stacking parameters, for each of the geometries employed in the present work.