Vortex induced vibrations of wind turbine blades: Influence of the tip geometry

The present investigation used numerical simulations to study the vortex induced vibrations (VIVs) of a 96 m long wind turbine blade. The results of this baseline shape were compared with four additional geometry variants featuring different tip extensions. The geometry of the tip extensions was generated through the variation of two design parameters: the dihedral angle bending the blade out of the rotor plane and the sweep angle bending the blade in the rotor plane. The applied numerical methods relied on a fluid structure interaction (FSI) approach, coupling a computational fluid dynamics solver with a multi-body structural solver. The methodology followed for locating VIV regions was based on the variation of the inclination angle. This variable was defined as the angle between the freestream velocity and the blade axis, being 0° when these vectors were normal and positive when a velocity component from tip to root was introduced. For the baseline geometry, the FSI simulations predicted significant blade vibrations for inclination angles between 47.5° and 60° with a maximum peak-to-peak amplitude of 2.3 m. The installation of the different tip extensions on the blade geometry was found to significantly modify the inclination angles where VIV was observed. In particular, the simulations of three of the tip designs showed a shifting of several degrees for the point where the maximum vibrations were recorded. For the specific tip geometry where only the sweep angle was taken into account, a total mitigation of the VIV was observed.


I. INTRODUCTION
Under certain conditions, in the presence of extreme winds or during maintenance inspections or during turbine erection, horizontal axis wind turbine rotors can adopt the so-called standstill configuration. In such a case, a pitch angle of 90 ○ is applied to the blades in order to aerodynamically break the rotor. For large modern horizontal axis wind turbine designs, the standstill configuration may pose a problem due to the high aeroelastic response of the blades. Indeed, several wind turbine manufacturers have unofficially indicated significant vibrations for their most recent designs. As these vibrations can lead to a potential structural failure, the research community is interested in giving a plausible explanation for the aeroelastic mechanisms that could be involved.
Even if the proneness of standstill wind turbine blades to experience flow induced vibrations was predicted in several works based on airfoils Pellegrino and Meskell, 2014;Zou et al., 2015;and Meskell and Pellegrino, 2019), the first research that studied a representative blade geometry was carried out by Heinz et al. (2016b). The authors presented a computational study of a 86 m long wind turbine blade, corresponding to the reference wind turbine (RWT) developed by the Technical University of Denmark (DTU) known as DTU 10MW RWT (Bak et al., 2013). It was shown that when the incoming flow led to an angle of attack in the vicinity of 90 ○ -due to a failure of the yaw system or its deactivation for maintenance-the structure could undergo extreme edgewise oscillations. The authors related the mechanisms driving these oscillations to the so-called vortex induced vibrations (VIVs) that rely on the synchronization of the flow shedding frequencies with the frequency of the structural response. Of particular importance for the triggering of VIV was the inclination of the freestream velocity. As for the present work, the inclination was defined as the relative angle between the incoming flow and the local blade axis at the root, being 0 ○ when these vectors were normal and positive when a velocity component from tip to root was introduced.
Indeed, the authors identified the inflow inclination as a correlating mechanism of the shedding frequencies along the span. In particular, considerable blade vibrations were measured for inclination angles higher than 20 ○ . Similar conclusions were made by the same authors in the subsequent work (Heinz et al., 2016a) based on a set of simulations performed on the 100 m long AVATAR blade (Lekou et al., 2015). Even if the aforementioned publications showed the susceptibility of modern wind turbine blades to undergo substantial flow induced vibrations at standstill, the mechanism driving the resulting aeroelastic response was not analyzed in depth. Indeed, the amplitude of vibrations included in both works presented a complex dependency with respect to the inclination angle, as well as to the considered wind speed and angle of attack. It is expected that a detailed analysis of the flow in such conditions will contribute to the progress in the understanding of the mechanisms influencing the blade response. This will also provide an answer to some of the design questions related to the existence of VIV, such as the impact of small geometrical changes on the resulting aeroelastic response of the blades. This concern is tightly related to the motivation of the current work. In fact, due to the required computational time, the number of simulations aiming at identifying the occurrence of VIV of new blade designs is extremely limited (if existent), that is, in the case of Smart Tip (2019), the wind energy project hosted the present research. In Smart Tip, the state of the art optimization techniques are being used in order to design a blade tip extension able to improve the power production capabilities of the rotor without increasing the blade loads. In the present work, the authors aimed at answering the following question: should the assessment of VIV be part of the optimization loop? In order to find a plausible answer, a phenomenological analysis was performed on the flow induced vibrations observed on the International Energy Association (IEA) 10MW RWT rotor blade (Bortolotti et al., 2019), which has a blade length of 96 m. This was achieved by means of several fluid structure interaction (FSI) simulations that couple a flow solver based on improved delayed detached eddy simulations (IDDES) with a multi-body structural solver. A range of high inclination angles was explored in order to identify the inflow conditions triggering important blade vibrations due to VIV. An additional set of simulations, based on the modified blade geometries accounting for tip extensions, was then carried out. A comparison of the flow patterns and aeroelastic responses of all the performed computations allowed us to elucidate some of the mechanisms involved in the VIV of wind turbine blades and to quantify the influence of the tip geometry.
The remainder of the present article is organized as follows: In Sec. II, a literature review is included, covering the most relevant publications concerning VIV and the aerodynamics of bluff bodies. Section III presents the numerical methods that were employed, followed by a description of the models that were generated, and the performed set of simulations. Section IV includes the results and discussion of the numerical study. Finally, the conclusions of the present work are given in Sec. V together with some potential directions for future work.

II. LITERATURE REVIEW
While the literature regarding VIV of wind turbines is scarce, several similarities can be established with the study of this phenomenon for other applications, as well as with traditional fundamental research lines. An overview of the most relevant VIV articles is given in this section. Special attention was paid in the identification of the different geometrical and inflow parameters that could influence the mechanisms involved in the VIV of wind turbine blades.
A. The canonical case: VIV of a circular cylinder The circular cylinder is, in words of Roshko, the quintessential bluff body (Roshko, 1993). It has served as a reference case for a considerable amount of experimental and computational research concerning flow induced vibrations. Particularly studied is the elastically mounted setup, where the cylinder is allowed to oscillate in the cross-flow direction. In such a configuration, large amplitudes of vibrations are measured when the shedding frequency of the non-moving cylinder is close to the frequency of the mounting. In this condition, the frequency of the flow loading can align with the frequency of motion once it is allowed to oscillate freely, a phenomenon that is often referred to as lock-in, vortex-resonance, or simply synchronization. Even if the transversal oscillations of a circular cylinder are probably the most elementary configuration where the VIV phenomenon is observed, it is far from being a simple study case. Among others, complex dependencies with respect to the Reynolds number, the mass ratio between cylinder and fluid, the damping ratio of the elastic mounting, or the surface roughness have been reported. To have an overview of the problem, the reader is referred to one of the regularly published review papers and monographs about the subject [for instance, Griffin et al. (1973), Morgenthal (2000), Williamson and Govardhan (2004), Païdoussis et al. (2011), Bearman (2011), and Gsell (2016)]. The complexity of the problem is further increased when considering additional aspects related to the studied geometry or the boundary conditions, as described in Subsections II B-II G.

B. Influence of the cross section
A wind turbine blade is a high-aspect ratio geometry, with airfoil-shaped cross-sectional shape, and is typically twisted and tapered. The mechanisms driving the VIV on a wind turbine blade may therefore substantially differ from those observed for circular cylinders and in fact also between different blade designs. This is because the geometric angle of attack will vary along the span according to the geometric twist, and the airfoil shaped cross-sectional shape will cause an asymmetry in the vortexshedding pattern due to the rounded leading edge and cusp trailing edge. Nakamura and Matsukawa (1987) investigated experimentally the VIV of several rectangular shape cylinders that were allowed to oscillate in the cross flow direction. The authors identified a synchronization region, showing a strong dependency on the aspect ratio of the cylinder. Ding et al. (2015) studied several simple crosssectional shapes undergoing VIV by means of numerical simulations. The study revealed a significant dependency of the observed VIV synchronization region with regard to the considered geometry. A similar conclusion was made from the experiments of Gomes and Lienhart (2013), where the effect of turbulence on the aeroelastic response was also addressed.

ARTICLE scitation.org/journal/phf
Apart from the influence that non-circular geometries could have on VIV, these cross sections are also prone to other types of flow induced vibrations. Several phenomena fall in this category, such as the low speed galloping (Nakamura and Hirata, 1989;Tamura and Itoh, 1999;and Tamura and Dias, 2003), the transversal galloping (Nakamura et al., 1991;Massai et al., 2018), or the torsional galloping (Robertson et al., 2003;Chen et al., 2012;and Fernandes and Armandei, 2014).

C. Influence of the inclination angle
As stated in the Introduction, the inclination angle has been identified as one of the most important parameters influencing the VIV experienced by wind turbines. The inclination angle is defined as the angle between the incoming flow and the structural axis, being 0 ○ when these vectors are normal and positive when a velocity component from tip to root is introduced. The effects of the inclination angle on the VIV of other structures have been extensively assessed in the literature, and some of the most relevant works are listed in this section. Bourguet and Triantafyllou (2014) compared the vibrations of an inclined cylinder with a non-inclined one. First, the authors assessed the wakes of a fixed configuration. While at 0 ○ of inclination, the shedding was parallel, and an oblique shedding was found at 80 ○ . Based on the analysis of the forces, the authors anticipated that the latter configuration was not susceptible to undergo VIV. However, the study of the same cases when the cylinder was allowed to oscillate led to a completely different scenario. A transition from oblique to a more parallel shedding was observed for the 80 ○ case. The frequency of this parallel shedding was found to be the structural frequency. This corresponded to a lock-in condition, similar to the classical VIV problem, but in this case, it was achieved through the reconfiguration of the wake due to the oscillation of the body. The aforementioned work was complemented by the study of Bourguet and Triantafyllou (2016). The authors aimed at identifying the mechanisms underlying the onset of the wake reconfiguration of the cylinder at 80 ○ . This was achieved by changing the magnitude of the incoming velocity while keeping the inclination angle as constant. The results showed that after a certain reduced velocity, the cylinder started to vibrate from rest. The authors characterized the reconfiguration of the wake by means of two frequency contributions. On one hand, the so-called Strouhal component that was related to the body wake was found even before large amplitudes of vibrations were observed. On the other hand, the appearance of the reconfiguration of the wake was accompanied by the raise of another frequency component, denoted as lock-in. This lock-in component was related to the body motion and corresponded to the structural frequency of the system.
It should be noted that inclined cylinders have been widely studied in connection to civil engineering applications, in particular to the study of wind-induced cable vibrations. These efforts tried to find a phenomenological answer for some of the important oscillations recorded for bridge cables, and it constitutes a research line by itself. The observed phenomena include rain-wind-induced vibrations (Matsumoto et al., 1988;Yeo and Jones, 2010), highspeed vortex excitation (Matsumoto, 1998;Matsumoto et al., 2001;Raeesi et al., 2012;and Jakobsen et al., 2012), and dry inclined cable galloping (Raeesi et al., 2012;Matsumoto et al., 2010).

D. Influence of the curvature
Modern wind turbine blades are characterized by a significantly curved axis due to the build-in out of plane prebend to provide clearance between the blades and the tower. Additionally, they may also account for a non-negligible in-plane sweep angle. These geometrical aspects are suspected to introduce important changes in the aeroelastic response of the structure, when compared to a straight design.
The characteristics of the flow passing a fixed curved cylinder substantially differ from the ones observed for the straight counterpart. Of particular relevance to the present research is the concave case that has been studied by several authors. For instance, the experiments of Seyed-Aghazadeh et al. (2015) revealed non-shedding regions behind a fixed concave cylinder-corresponding to the locations of high local curvature-together with the existence of a significant axial flow. The authors also pointed out that when the cylinder was allowed to vibrate, the shedding re-appeared all along the cylinder axis, being parallel to the local curvature. VIVs were finally recorded for this concave geometry, even if the maximum amplitude was attenuated with respect to the straight cylinder case. A similar conclusion regarding the mitigation of the maximum vibrations was made from the computations of Assi et al. (2014), and the suppression of the shedding at the high curvature regions was also observed in the experiments of Shang (2015) and the simulations of Jiang et al. (2018). The re-appearance of the vortex shedding was also found in the numerical study of De Vecchi et al. (2009), when a concave cylinder was forced to oscillate. Chaplin and King (2018) presented a set of experiments aiming at studying the VIV of a catenary riser with high concave curvature and very low bending stiffness. The authors observed a more complex structural response with respect to springmounted cylinders, accounting for broad-banded displacement spectra.
The simulations of Gallardo et al. (2014) on curved convex cylinders are also relevant to the present work, as they focused on turbulent wakes. The authors showed that the wake was mainly characterized by periodic shedding, with primary vortical structures oriented along the axis of the cylinder and secondary structures exhibiting different orientations. In the region of high local curvature, the coherence of the wake was broken due to the disappearance of the primary structures. The study of the time-average streamwise vorticity revealed the existence of two counter-rotating vortices, evolving along the axis of the cylinder. A broadband shedding spectrum was also observed.

E. Influence of the non-uniformity of the inflow
While the present work only considers a uniform incoming flow for simplicity, it should be noted that wind turbine blades are subjected to more complex inflow conditions. In particular, the incoming flow may account for a high turbulence intensity, together with a more or less pronounced shear profile due to the atmospheric boundary layer. Additionally, the studied blade may be in the wake of another blade or wind turbine. Bourguet et al. (2011) studied the VIVs of a long flexible cable in shear flow. The authors identified two distinct span regions with regard to their behavior in the frequency domain. The first one was the lock-in region, and it was characterized by a synchronization of the vibration and shedding frequencies.
The other region was referred to as the non-lock-in region and revealed a discontinuous cellular pattern, with frequencies differing from the vibration frequency. An insight into the influence of several shear profiles on the aeroelastic response of flexible cylinders is presented in the computations of Bourguet et al. (2013). The authors attributed to the considered profile the switching between narrowband and broadband lock-in. While in the narrowband case, the lock-in seemed to be localized and limited to a particular span region and the broadband lock-in spread to several span locations. The recent publication of Wang and Xiao (2016), which presented a validation of FSI simulations of a flexible riser in uniform and sheared inflow, revealed a considerable dependency of the predicted vibrations with respect to the velocity profile.
The mechanism of VIV is further complexified when the incoming flow is already perturbed by another bluff body, motivating a research line revolving around the so-called wake induced vibrations (WIVs). A well-studied academic setup for the study of WIV is two cylinders in tandem. When comparing the VIV of this configuration to the case of a single cylinder, substantial changes in the range of concerned velocities are expected (Bokaian and Geoola, 1984;Maeda et al., 1997;and Li et al., 2018). This includes the so-called infinite WIV, where the lock-in phenomenon is always observed beyond a certain wind speed.

F. Influence of the taper
As already mentioned, wind turbine blades are high aspect ratio tapered geometries. This implies a theoretical constant variation of the expected shedding frequencies in the spanwise direction, which should be taken into account while assessing VIV.
In Hover et al. (1998), experiments were performed in order to compare the vortex shedding of straight and tapered cylinders, being all stiff. Measurements showed that the tapered cylinder exhibited a wider region of wind speeds where loads were poorly correlated in the spanwise direction. Balasubramanian et al. (1998) presented a series of experiments for stiff tapered cylinders under uniform and shear inflow conditions. For both setups, cellular vortex shedding was observed, resulting in the clustering of shedding frequencies for particular regions of the span. Introducing a shear inflow, which offset the influence of the taper of the cylinder, resulted in a cancellation of the spanwise frequency variation along the span. When the influence of the shear inflow was added to the taper of the cylinder, a greater variation of the shedding frequencies along the span was observed. This work was followed by the experiments of Balasubramanian et al. (2001), where two cylinders accounting for positive and negative tapers were allowed to oscillate. Results revealed different lock-in regions and maximum recorded vibrations, depending on the taper sign and the consideration of uniform or shear inflow. The authors also pointed out that, in lock-in conditions, the wake was more organized due to the vibration of the cylinder. The recent computational studies of Kaja et al. (2016) concerning tapered cylinders showed that, for certain inflow conditions, the wake accounted for a combination of two distinct shedding patterns (i.e., 2S and 2P). Additionally, the studied tapered cylinder was found to have a wider range of velocities exhibiting lock-in, when compared to the uniform counterpart.

G. Influence of the tip
The publications presented in Secs. II A-II F revolved around the study of two-dimensional cases together with cylinders and cables. Due to the slenderness of these structures, the role played by the tip of the geometry is assumed to be relatively small. This explains why little or no effort was devoted to the study of the impact of the design of the tip on the resulting VIV. As concluded from the results of the present research, the case of a wind turbine blade is substantially different. Indeed, a strong influence of the tip geometry on the generated wake was found, especially for high inclination angles. Additionally, the tip region is the most aeroelastically active part of a wind turbine blade. Similar interactions were found in the literature for aircraft fuselage and missile aerodynamics, which has been studied for more than six decades. A brief discussion around the main findings concerning these topics is presented in this section.
Several early works characterized, by means of visualization techniques, the flow behind slender cylindrical bodies at different inclinations with respect to the incoming wind (Allen and Perkins, 1951;Thomson and Morrison, 1971). At low inclination angles, an unsteady flow regime was observed all along the body characterized by vortex streets. At high inclination angles (i.e., >25 ○ or >30 ○ ), no wake was identified in the vicinity of the tip of the body, and the appearance of two stationary symmetric vortices behind the body was reported for lower axial locations. These vortex lines were found to be straight, except for a curved region near the body. The aforementioned flow patterns for low and high inclination angles are often referred to as asymmetric and symmetric wakes, respectively. In this document, the terms unsteady and stationary will be adopted instead, as they can also be used in the context of non-cylindrical shapes. Depending on the inclination angle, both regimes can coexist in the flow. The review paper of Polhamus (1984) compiled the progress made by his contemporaries. The resulting side forces depended on the inclination angle, being almost negligible at 0 ○ and 90 ○ , and exhibiting a maximum at around 45 ○ . Additionally, the unsteadiness of the wake was found to be highly influenced by the geometry of the tip. This last statement was corroborated by the experiments of Zilliac et al. (1991), where the side forces were also found to be negligible for inclination angles higher than 60 ○ due to the stationary nature of the wake. In Levy and Hesselink (1995), the installation of small perturbations located at the tip of the body was found to be sufficient to trigger a change in the wake from stationary to unsteady. Ma and Liu (2014) presented the results of different numerical simulations at high inclination angles for an ogivecylinder body. The authors showed that for an inclination angle of 30 ○ , the wake close to the tip consisted of a pair of vortices, which could lead to large side-force fluctuations due to their oscillation. For an inclination angle of 40 ○ , this flow pattern was replaced by a set of vortical structures that were essentially steady. Similar shedding patterns were observed in the simulations of Ma and Yin (2018), which considered a hemisphere-cylinder body, and for delta wing geometries (Gursul, 2005;Luo and Ma, 2018), where the existence of leading edge vortices is well known. Jiang et al. (2015) studied the flow around an inclined prolate spheroid at 45 ○ , representative of a submarine geometry, by means of numerical simulations. The authors showed that the wake was dominated by a single coherent vortex tube, and the flow frequencies behind the body were found ARTICLE scitation.org/journal/phf to be significantly lower than those observed for the von Karman shedding. Modern research efforts concerning the study of the aerodynamics of this type of structure are also devoted to the implementation of control strategies to keep a stationary wake in order to avoid unfavorable unsteady loading. In this category fall the study of Liu et al. (2009), concerning the use of plasma actuators, and the work of Zhai et al. (2016), which explored the possibility of installing an oscillating flag at the tip. All the aforementioned publications were focused on the influence of the unsteady loads on the reliability of the structures and on the maneuverability of the vehicle, and no classical VIVs were reported for these applications. Nevertheless, it is worth noting that the wing rock, a phenomenon that can be categorized as flow induced vibration and that involves the presence of a stationary wake, has attracted the attention of several researchers. This phenomenon is characterized by a self-excited and large amplitude roll oscillation, and it has been observed for several fighter aircrafts at high angles of attack. The experiments of Ma et al. (2015) studied the wing rock by means of static, forced, and free oscillation tests of an aircraft model. By imposing the same roll motion that was recorded in a previous free oscillation test, the authors identified the generation of motion-induced unsteady vortices for span regions that recorded stationary wakes in static conditions. These unsteady vortices could lead, for some cases, to a limit-cycle oscillation of the roll angle. Additionally, the authors revealed the significant impact of the aircraft nose geometry on the recorded amplitude of oscillations. The implications of the model geometry (Ma et al., 2016) and Reynolds number (Ma et al., 2017) on the resulting oscillations were assessed by the same authors in subsequent studies. Degani et al. (2017) studied the wing rock phenomenon by means of computational methods, where both static and free oscillation conditions were assessed. Limit-cycle rotations of large amplitude were found for inclination angles of 35 ○ and 50 ○ . Additionally, an oscillation of the initially stationary vortical structures was observed when analyzing the flow of the free-to-roll simulations due to the motion of the structure.

III. METHODOLOGY
The numerical methods used for the present work consisted of a finite-element multi-body structural solver, a fluid solver, and a coupling framework that handled the communication between these through a staggered FSI approach. An overview of these tools is given in Sec. III A. This is followed by the description of the analysis that was carried out, stating also the scope and assumptions of the present study (Sec. III B). Finally, Sec. III C presents the numerical models that were used.

Structural solver: HAWC2
HAWC2 (Larsen and Hansen, 2015) is a finite-element multibody serial solver that allows us to perform multi-physics simulations of whole wind turbine assemblies both for onshore and offshore conditions. Among other components, the tool implements a hydrodynamic module as well as several standardized interfaces in order to account for the action of a wind turbine controller. In the present study, only the HAWC2 capabilities exclusively dealing with the problem of aeroelasticity were used (i.e., parts of its aerodynamics module and the aeroelastic solution). In this context, HAWC2 equips a built-in aerodynamics module based on an enhanced variant of the blade element momentum (BEM) method (Madsen et al., 2020). This relies on the definition of several sections along the blade, where the aerodynamic loads are calculated. The positions of these aerodynamic sections may differ from the discretization used for the structural definition of the problem, so a numerical integration is performed in HAWC2. Due to the limitations of BEM for standstill configurations (Pirrung et al., 2016), a higher fidelity aerodynamic solver based on computational fluid dynamics (CFD) was used in the present study (see Sec. III A 2). This was achieved by replacing the BEM-based loads by CFD-based loads during runtime, allowing us to re-use the aerodynamics section architecture of HAWC2, and thus minimizing the amount of modifications to perform in the code.
The framework that HAWC2 uses for the solution of multiphysics problems relies on an augmented form of the floating frame of reference (FFR) formulation [see, for instance, Shabana (2001)]. In FFR, the kinematics of the flexible body is described in a coordinate system, which translates and rotates relative to an inertial system. In this way, the blade of the present study was characterized by two complementary discretizations. On one hand, a set of linear isotropic Timoshenko beam elements was defined. On the other hand, an arbitrary number of bodies were also specified along the blade. These were connected through a set of constraints that allowed us to capture the geometrical non-linearities of the blade response under the presence of the CFD loading.
In the present study, HAWC2 solved the following equation of motion under the presence of constraints, where the set of unknowns included both the vector of generalized coordinates ⃗ u and the Lagrange multipliers ⃗ λ: (1) with M, C, and K being the inertia, damping, and stiffness matrices, respectively. The vector ⃗ f contains the external aerodynamic loads and ⃗ g the constraint equations. Finally, G is the Jacobian of the constraint equations with respect to the generalized coordinates.
In HAWC2, the solution of Eq. (1) is performed iteratively for every considered time instant based on a predictor/corrector approach. A Newmark-beta scheme is used for the time-integration (Newmark, 1959).

Flow solver: EllipSys3D
The three-dimensional CFD pressure-based solver EllipSys3D (Michelsen, 1992;and Sørensen, 1995) was used for the computation of the aerodynamic loads. EllipSys3D solves the Navier-Stokes equations under the hypothesis of an incompressible and isothermal flow. Assuming no external sources, these equations can be written in their differential form as Phys. Fluids 32, 065104 (2020); doi: 10.1063/5.0004005 32, 065104-5

ARTICLE scitation.org/journal/phf
where t refers to the time, ⃗ v is the velocity vector, ρ is the fluid density, and p is the pressure. τ refers to the shear stress tensor that under the hypothesis of a Newtonian fluid, it can be expressed as with μ being the dynamic molecular viscosity. EllipSys3D uses a finite volume discretization, and the solution of the Navier-Stokes equations is performed on a structured grid and in curvilinear coordinates. It accounts for moving mesh capabilities through the arbitrary Lagrangian-Eulerian formulation (Hirt et al., 1974), and several turbulence models have been implemented. For the simulations involved in this work, a zonal IDDES was employed (Gritskevich et al., 2012;Shur et al., 2008). This model uses a Reynolds-averaged Navier Stokes (RANS) model close to the blade surface and a large eddy simulation (LES) approach far away from the boundary layer. For the region involving the RANS method, the k-ω SST turbulence model of Menter (1994) was employed. EllipSys3D is second order in time and uses a third order accurate upwind scheme in the RANS region and a fourth order accurate scheme in the LES region. The solver implements a multiblock decomposition of the domain, allowing an efficient parallelization of the code through the use of the message passing interface (MPI) standard (Forum, 1994).

FSI coupling approach
The transfer of the aerodynamic loads from EllipSys3D into HAWC2 was performed once every time iteration. Consistently, the blade deflections predicted by HAWC2 were also transferred back to EllipSys3D so that they were accounted for while solving the Navier-Stokes equations. The coupling strategy between the aforementioned solvers was orchestrated by an in-house generic framework. This platform, referred to in the present document as DTU coupling, is a tool developed at the Wind Energy Department of DTU. It is based on the interpreted programming language Python (Van Rossum and Drake, 1995), and it can handle the execution and data transfer of an arbitrary number of numerical solvers. The current version of the DTU coupling supports the Linux operating system, and it was designed in order to work in a parallel environment through the use of MPI. As a prerequisite for the integration of each of the independent solvers in the DTU coupling, these need to be compiled as external libraries (i.e., Linux shared objects). Additionally, a set of Ccompliant interfaces should be implemented in each of the solvers. These interfaces are supported by most of the programming languages used for scientific codes (e.g., Fortran, C, and C++), and their objective is twofold. First, they allow us to independently trigger the different actions involved in the execution sequence (e.g., initial allocation of the variables, carrying out the solution for the current time iteration, and writing of output files at the end of the simulation). Second, the interfaces enable the data transfer with the DTU coupling during run-time by means of a collection of set and get methods. In practice, a Python module is also developed around each of the contributing solvers in order to build an object-oriented highlevel interface that facilitates the definition of the coupling strategies inside the DTU coupling. A schematic of the particular adaptation of the DTU coupling framework for the present study is illustrated in Fig. 1, where the Python modules developed around the HAWC2 and EllipSys3D shared objects are labeled pyhawc and pyellipsys, respectively.
The implemented FSI strategy was based on a loose coupling scheme that kept the built-in predictor-corrector strategy of HAWC2. A comprehensive description of the coupling sequence was already published by Heinz (2013), but a list summarizing the stages involved in every time iteration i is also included below.
1. HAWC2, predictor step: Prediction of the structural coordinates for the next time iteration⃗ u i+1 computed as with dt = ti +1 − ti being the time step. 2. DTU coupling: Transfer of⃗ u i+1 from HAWC2 to EllipSys3D. 3. EllipSys3D: Flow solution: (a) Application of⃗ u i+1 to the EllipSys3D model. This is achieved by interpolating the rotational and translational displacements of the beam nodes into the CFD surface grid using a spline based interpolation. (b) Re-adaptation of the inner CFD mesh with the use of a mesh deformation technique. For the present study, a generalized analytical approach was used. In particular, the deformation vector was smoothed based on a hyperbolic tangent function, taking as a variable the distance to the blade surface along the considered grid line.

Validation and previous studies
Both EllipSys3D and HAWC2 are commonly used tools in the wind energy research community, and they have been applied to a wide range of turbines and loading scenarios. The former solver was validated against the experimental results of two publicly available wind turbine models: the Mexico (Benchmann et al., 2011;Sørensen et al., 2016) and the NREL Phase VI (Sørensen et al., 2002). Validations of the employed turbulence model have also been published, including two-dimensional geometries  and a whole wind turbine rotor (Sørensen and Schreck, 2014). A validation of the standstill loads in pitch motion of the NREL Phase VI was also presented in Sørensen and Schreck (2011). The structural solver implemented in HAWC2 was verified and compared against higher fidelity models by Pavese et al. (2015). A recent publication also included a validation of the adopted FSI methodology in the context of the pull-release test of a wind turbine blade (Grinderslev et al., 2019). While no experimental data are available for the validation of the coupling framework when applied to large wind turbine rotors in standstill, several comparisons with lower fidelity models have been made in order to assess the appropriateness of the implementation. The user is referred to the monograph of Heinz (2013) for more details on the latter subject, where studies about the stability of the chosen FSI approach were also included.

B. Performed analysis
The 3-bladed wind turbine chosen for the present study, and referred to as IEA 10MW, corresponds to the offshore design presented in Bortolotti et al. (2019). With a length of 96 m, the IEA 10MW blade is assumed to be representative of state-of-the-art designs of horizontal axis wind turbines. This model is the result of a numerical optimization process, where the power production was maximized while keeping the blade loads at a certain level. The IEA 10MW is equipped with a series of airfoils of the FFA-W3 family (Bjorck, 1990). The thickness, twist, and chord length vary as a function of the radial location, the latter parameter having a maximum value of cmax = 6.02 m. The blade axis also has a considerable out-of-plane prebending, with the tip located more than 6 m of out-of-plane. When considering the whole wind turbine assembly, the blades also have a 3 ○ pre-cone angle and a 6 ○ tilt angle (see Fig. 2). However, these were disregarded in the present study for the sake of simplicity. The blade is made from composite materials, and its structural layout accounts for two main spars and three shear webs. In order to study the potential edgewise vibrations of the IEA 10 MW blade in standstill at an affordable computational time, a reduced set of inflow conditions was defined. Based on previous studies of the VIV of similar blades (Heinz et al., 2016b;2016a;, this was done by the following: • Fixing the magnitude of the freestream velocity U∞ to 18 m s −1 . This value represents a relatively high wind speed, but with a reasonable probability of occurrence.

FIG. 2.
Sketch with some of the geometrical parameters characterizing a wind turbine design.
• Fixing the orientation of the projection of U∞ into the local cross-sectional plane at the blade root. In particular, this value was set in such a way that the angle of attack at the tip of the undeflected blade geometry was 100 ○ . This is indeed one of the most favorable situations for the appearance of VIV based on the aforementioned studies. • Systematically varying the inclination angle I in the range (40 ○ , 70 ○ ). I is defined as the angle formed between the freestream velocity and the local blade axis at the root. It is null when both vectors are normal, and positive when a spanwise component from tip to root is introduced. Based on the aforementioned publications, important edgewise vibrations are expected for the values of I within the studied range. In order to improve the understanding of the flow, additional simulations were performed for I = 0 ○ and I = 10 ○ .
The set of inflow conditions described in the list above is illustrated in Fig. 3(a). In operational terms, they can be understood as a variation in the azimuthal position of the studied blade in the presence of a significant yaw misalignment, as depicted in Fig. 3(b). The chosen wind speed, 18 m s −1 , corresponded to chord Reynolds numbers in the range (1.0 ⋅ 10 6 , 5.6 ⋅ 10 6 ). For all the considered inclination angles, two simulations of the baseline geometry (BASELINE) were performed. The first one, labeled in this document as flexible, took into account the flexibility of the blade through the FSI coupling described in Sec. III A 3. The other computation assumed the blade as stiff, and corresponded to a pure CFD simulation carried out by EllipSys3D exclusively. Additionally, the influence of the blade tip on the resulting VIV was assessed by repeating this procedure with several blade geometry variations. In particular, four tip extensions were considered, with the same projected length into the rotor plane (i.e., 5% of the baseline blade projected spanwise length). The blade tip geometries were designed through a parameterization, which fitted an Akima spline through 3 spanwise points (Akima, 1970). The aforementioned points were located at 2.5% projected length inboard of the baseline tip, at the baseline tip, and at the extended tip. The chord at the location of the baseline tip was equal to the chord of the first inboard section. The chord at the tip of the extension was defined as 50% of the chord of the first inboard section, resulting in a more slender tip geometry compared to the baseline. The twist distribution was the same as 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 tip extensions considered represented the parametric space expected for the optimization and were characterized by smooth shapes representing specific geometrical design concepts:

Physics of Fluids
• EXT is a straight extension of the blade.
• SWEEP accounts for a backward sweep angle.
• DIHEDRAL (or simply DIHE) accounts for a dihedral angle, opposite to blade prebend. • COMBINED (or simply COMB) is a combination of SWEEP and DIHEDRAL.
A summary of the blades considered in the present work is included in Table I, and a detailed visualization of the tip geometries is depicted in Fig. 4.

CFD model
A CFD mesh for each of the blade variations considered in the present study was generated. This process was fully automatized in order to guarantee a similar resulting grid quality. The mesh generation accounted for the two consecutive steps described below, 1. Blade surface mesh generation: The openly available parametric geometry library (PGL) was used in a first step in order to create a structured grid of the blade geometry (Zahle, 2019). This mesh accounted for the local curvature of the blade axis, as well as for the tip geometry (if applicable). A total of 160 cells were used in the spanwise direction, and the   Fig. 5(a), and a detailed view around the blade tip is depicted in Fig. 5(b). 2. Volume mesh generation: A volume mesh was generated with the hyperbolic mesh generator Hypgrid (Sørensen, 1998), grown into a spherical domain, with outer boundaries located an approximate radius of 116cmax. A total of 256 cells were used in this normal direction. A boundary layer clustering was imposed with a first cell height of 10 −6 m, representing 1.7 ⋅ 10 −7 cmax and leading to y + values lower than unity. An additional refinement region was defined in order to increase the resolution close to the blade surface, extending up to 19cmax. Details of the cross-sectional mesh at half of the span are shown in Fig. 5(c) for the case of the BASE-LINE geometry. Figure 5(d) depicts the corresponding CFD domain boundaries together with the blade surface mesh for comparison.
The resulting mesh contained a total of 12.6 × 10 6 cells. The outer limit of the CFD domain was defined as an inlet (i.e., Dirichlet), except for a small opening where an outlet boundary condition was set (in the form of a zero normal gradient). The effect of the inclination angle was introduced by varying the wind speed direction imposed at the inlet, without modifying the position of the blade. A three level grid sequence was used during the simulation to speed up the development of the wake flow. Values of ρ = 1.225 kg m −3 and μ = 1.784 ⋅ 10 −5 kg m −1 s −1 were used. The flow was assumed to be fully turbulent, which is believed to be a good approximation at the high angles of attack considered in the present study. At the boundaries, the values of the specific turbulent dissipation rate ω = 10 6 s −1 and the turbulent kinetic energy k = 10 −2 m 2 s −2 were imposed.

Structural model
The blade was assumed to be clamped at the root. It was modeled by means of 19 Timoshenko beam elements, and a total of 10 bodies were used. For BASELINE, the mechanical properties of the beam elements of the flexible simulations corresponded to the original definition of the IEA 10MW blade. For the geometries accounting for a tip extension, these were estimated by stretching the baseline properties. This resulted in the same spanwise distribution compared to the tip region of BASELINE. The addition of the tip extensions did imply a slight change in the natural frequencies of the blade, as well as in the structural damping. The values of the first six modes for each of the blade geometry variants, obtained in a preliminary analysis performed with HAWC2, are listed in Table II.

Models integration
A total of 101 aerodynamic sections were defined along the blade, where the integrated CFD loads were injected into the HAWC2 solution. A time step of dt = 6 ⋅ 10 −3 s was chosen so that dtU∞/cmax = 1.8 ⋅ 10 −2 . Each simulation was run up to t = T = 45.3 s with TU∞/cmax = 135.
Preliminary simulations were performed, aiming at studying the sensitivity of some of the model parameters. A summary of the results for this study is presented in Appendix A. The computations presented in this document were run on the Jess high-performance computing the cluster, owned by DTU. Jess is equipped with a total 320 nodes, each of them with 20 cores running at 2.8 GHz. With the exception of the aforementioned sensitivity studies, all the simulations of the present work required a wall clock time of approximately 5 h when using 192 cores.

IV. RESULTS AND DISCUSSION
The results of the performed simulations are described as follows: Section IV A presents the CFD simulations of the different blade geometries in stiff conditions. In Sec. IV B, the results of the FSI simulations, which also accounted for the blade flexibility, are introduced.
The discussion revolves around the study of three main sets of results.
a. Wakes. The wake of the blades is visualized for the selected simulations. This is done through the use of the vorticity magnitude, the value of its y component, or the Q-criterion Qc. The latter quantity being defined as in Hunt et al. (1988), with Ω and S being the antisymmetric and symmetric parts of the velocity gradient tensor, respectively, Additional visualizations of the wakes around the tip were also performed in order to have a better insight into the influence of the inflow conditions and blade geometry on this region. In order to improve the readability of the manuscript, these are presented in the form of a complementary analysis in Appendix B. b. Blade loading. Details of the sectional blade loading in the edgewise direction are included for the selected simulations. These are presented in the form of spectra, computed by means of the power spectral density (PSD) based on the Welch method (Welch, 1967). Only the last 19 s of each simulation were taken into account for this computation. A linear detrend was applied to the signal, and a sampling period of 12 ms was used. The results of the PSDs were linearly scaled so that the peak height became an estimate for the root mean square (rms) amplitude. For completeness, most part of the time series involved in the computation of the PSDs are included in Appendix C.
c. Tip motion. For flexible simulations, details of the tip motion in the edgewise direction are given for the selected simulations. This is done by means of time series (in order to illustrate the amplification patterns) together with the corresponding peak-to-peak amplitude and the computed main frequency of the signal. Additionally, the mean values of the tip motion in other directions rather than edgewise are also presented.
A. Stiff simulations 1. Results for the BASELINE geometry Figure 6 shows the iso-surfaces of the Q-criterion for simulations of the BASELINE geometry at different inclination angles. The used instantaneous CFD flow field corresponds to the end of the simulation, where t = T = 45.3 s. For a better visualization, the vortical structures are colored based on the y component of the vorticity ωy. In particular, regions of positive ωy are displayed in red, while the negative vertical vorticity is shown in blue. At low inclination angles (i.e., I = 0 ○ and I = 10 ○ ), an unsteady turbulent wake dominated the flow. It was characterized by the shedding of disorganized vortical structures, as expected for the considered Reynolds number (Lienhard, 1966). When increasing the inclination angle, stationary vortices were observed for a certain extent of the outboard part of the blade (see the cases I = 40 ○ , I = 50 ○ , and I = 55 ○ ). In the present work, this region is referred to as the stationary wake. The leeward vortices of this stationary wake were convected away from the blade following the flow direction forming a set of straight vortical tubes. Consistently with previous experiences for other geometries, the vortex tubes exhibited a local curvature near the body surface (Allen and Perkins, 1951;Thomson and Morrison, 1971). At I = 70 ○ , the stationary wake was observed all along the blade span. This wake remained very close to the body surface, and the shedding only took place in the vicinity of the blade root. A close examination of the flow at high span allowed us to identify the presence of a pair of stationary vortex tubes emanating from the tip vicinity, irrespective of the considered inclination angle. The reader is referred to Appendix B 1 to get more details about this observation, where the wake topology is also further described.
A complementary visualization of the wakes for I = 0 ○ and I = 50 ○ is included in Figs. 7 and 8, respectively. Two-dimensional plots of ωy are depicted for different span locations and instants of time. In particular, four different time instants close to T were chosen. Figure 7 shows the distinct trailing and leading vorticities all along the blade span, corresponding to the shear layers. These shear layers eventually roll-up downstream, leading to detached wake vortices. The downstream distance where this vortex generation occurs varies with the spanwise distance y, being smaller at high span. From the comparison of the different time instants, it is clear that this process is inherently unsteady. The flow computed for the simulation at I = 50 ○ significantly differs from this pattern (Fig. 8). The steadiness of the leeward vortex tubes becomes evident by comparing the different time instants that are identifiable for all span locations. The shedding of vortical structures characterizing the unsteady wake only took place at low span (i.e., y = 17 m for this visualization), and it seemed to be constrained by the presence of the leeward vortex tubes. Figure 9 shows the PSDs of the sectional CFD loading in the edgewise direction, computed for different inclination angles for the BASELINE geometry. Several span locations are included as different curves. At 0 ○ of inclination, a broadband spectrum was computed for every span location. The lack of a dominating frequency could be related to the changes in the sectional geometry of the blade and, in particular, with the chord variation. This remark is in accordance with the previous observations made by Hover et al. (1998) for cylinders. The magnitude of the load fluctuations was found to be relatively low with respect to the other inclination angles studied, as also found for different geometries (Polhamus, 1984). The increase in I up to 10 ○ resulted in significant load fluctuations at outer most span locations. Dominating frequencies appeared for the PSDs in this region of the blade with similar values. This reveals an important load correlation in the spanwise direction, which has also been reported for other wind turbine blades (Heinz et al., 2016b;Horcas et al., 2019). The mechanism triggering this correlation is assumed to be the axial component of the freestream velocity. For the high inclination angles (i.e., 40 ○ and 50 ○ ), the existence of the stationary wake led to the suppression of edgewise unsteady loading for the high span region.

Physics of Fluids
To have a more quantitative assessment of the main frequencies of the CFD loading, a simple peak-identification algorithm was applied to the computed PSDs (Jones et al., 2001). The peak frequency fp showing the maximum value was then corrected by the neighboring points through the evaluation of a quadratic spline. The evolution of fp along the outer half of the blade for some of the simulations concerning the BASELINE geometry is depicted in Fig. 10.
To identify the presence of the stationary wake, the points were colored based on the PSD value of the selected peak, and a relatively low saturation value was used (i.e., 10 N m −1 ). Additionally, straight lines corresponding to the first edgewise frequency of the blade were included for reference. It should be noted that while this analysis presents some limitations due to the broadband nature of some of the PSDs and to their poor resolution in frequency, it was found to be a useful tool for the understanding of the flow. At 0 ○ inclination angle, a systematic increase in fp was found along the blade axis with significant fluctuations after approximately y = 75 m. As for the other simulations, a drop in the predicted peak frequency and PSD values was observed in the vicinity of the tip. This can be attributed to the presence of the pair of leeward vortex tubes emanated from the tip vicinity. For I = 10 ○ , distinct span regions of constant frequency were observed. This is known in the literature as a cellular pattern, and it is assumed to be triggered by the higher correlation of the flow in the spanwise direction. Additionally, considerable load fluctuations were recorded for the whole outer half of the blade. The fp predicted for both I = 0 ○ and I = 10 ○ simulations exhibited a considerable deviation with respect to the first edgewise frequency of the blade. Due to the existence of the stationary wake, the PSD values computed for I = 40 ○ , I = 50 ○ , and I = 55 ○ revealed negligible PSD levels. The value of the predicted frequency was then arbitrarily assigned by the peak finding algorithm without any physical meaning.
It is interesting to express the results of this peak identification methodology in terms of the reduced frequency, understood here as with U ref being a reference wind speed and l ref being a reference length. The U ref used in this study was the cross-sectional velocity (i.e., the projection of the freestream velocity vector into the local cross-sectional plane). This corresponds to the hypothesis known as the independence principle, stating that the flow around a cylindrical structure is solely characterized by the velocity component that is normal to its axis (Jones, 1947). According to this, the chosen l ref was the rejection of the chord vector into the cross-sectional velocity in such a way that the effect of the angle of attack was also taken into account. A graphical representation of both U ref and l ref parameters is included in Fig. 11. The fr distribution predicted for the BASE-LINE simulations in stiff conditions is depicted in Fig. 12. For 0 ○ of inclination, a similar reduced frequency was found for the whole outboard part of the blade. The value of this reduced frequency was very close to 0.16, which is the Strouhal number that has been previously reported in other wind turbine studies (e.g., Heinz et al., 2016b). At 10 ○ , the independence principle seemed to be partially violated because of the impact that the spanwise correlation had on fp. Since no clear vortical structures were observed for 40 ○ , 50 ○ , and 55 ○ , the levels of the PSD dropped and the predicted reduced frequencies had arbitrary values. Similar deviations from the independence principle were found in previous studies dealing with the effects of the flow inclination (Lucor and Karniadakis, 2003;Thakur et al., 2004;Hoang et al., 2015;Bourguet and Triantafyllou, 2016; and with the curvature of the geometry (Gallardo et al., 2014).

Influence of the tip geometry
At low inclination angles, the introduction of the tip extensions did not significantly alter the wake for most part of the blade span, and the resulting loads were found to be very similar. Figure 13 shows the reduced peak frequencies fr of the edgewise sectional loading for every geometry at I = 0 ○ . It can be observed that the independence principle also applies to the different extensions with a predicted Strouhal number centered around 0.16. The drastic decrease in the PSD levels observed for the BASELINE around the tip was also found for the EXT and SWEEP extensions. However, this phenomenon was not observed for the DIHE and COMB configurations. This can be attributed to the early shedding at the outermost span positions exhibited by those geometries, as described in Appendix B 2. A similar scenario was found when assessing high inclination angles, where the results of the BASELINE geometry exhibited a significantly developed stationary wake. Figure 14 shows the visualization of the Q-criterion at I = 50 ○ for the geometry variants with a tip extension. For the straight extension EXT as well as for the SWEEP geometry, the wake seemed to be qualitatively similar to the one computed for the BASELINE geometry [see  even if the generation of the leeward vortex tubes took place at different span locations. As reported for I = 0 ○ , the flow around the tips with a dihedral angle (i.e., DIHE and COMB) was found to significantly differ from the rest of the considered geometries. At high inclinations, the effect of the local curvature for those particular configurations seemed to lead to the early breakdown of the leeward vortex tubes at high span, resulting in a swirling flow that was convected away from the blade surface. A stationary wake was then found at lower span locations, analogously to the rest of the considered geometries. For the particular case of COMB, vortex shedding was also identified all along the blade span. This is further illustrated in Fig. 15, where contours of high vorticity at different sections and time instants are superposed to the solution computed at the end of the simulation. This significant sensitivity of the developed wake with respect to small geometrical variations of the tip is in line with the previous work, such as the ogive-cylinder body studies of Zilliac et al. (1991) and Levy and Hesselink (1995). While the blade loading computed for the stationary wake regions of the EXT and SWEEP geometries at I = 50 ○ was negligible, it was not the case for the configurations with a dihedral angle. Figure 16 shows the PSDs of the edgewise sectional loading at different sections for DIHE and COMB, expressed in terms of the reduced frequency. For the former configuration, the early breakdown of the leeward vortex tubes led to considerable sectional loading even at the region where the stationary wake was observed. The predicted loading was significantly higher for the COMB geometry due to the existence of blade shedding. For this geometry, a broad band spectrum accounting for a significant low frequency content was computed for the sections at y < 90 m. While clear main frequencies could not be established, an important contribution around 0.11 was observed. This value is slightly lower than the 0.16 predicted at I = 0 ○ . Figure 17 shows the visualization of the wakes of the different tip extensions when passing to I = 55 ○ by means of the iso-surfaces of the Q-criterion. The additional 5 ○ of inclination did not have a significant impact on the predicted flow for the considered geometries, with the exception of COMB. For this configuration, a clear stationary wake region was observed, and no shedding was reported. This led to an attenuation of the edgewise loading, in agreement with the observations made for the DIHE geometry at I = 50 ○ . Figure 18 depicts the time series of the blade tip position in the edgewise direction for several inclination angles. As for the z axis displayed in Fig. 4(c), an increase in the edgewise coordinate corresponded to a motion from the leading edge to the trailing edge. All the blade geometries are included, and only the last 30 s of every simulation are shown. For low inclination angles (i.e., 0 ○ , 10 ○ , and 40 ○ ), negligible edgewise vibrations were predicted for all the geometry variants. At I = 50 ○ and I = 55 ○ , a considerable increase in the vibration amplitude was recorded for several configurations, resulting in a limit cycle oscillation. Increasing the inclination angle up to 70 ○ led to the attenuation of the vibrations for all the performed simulations. The maximum peakto-peak amplitudes of the deflections in the edgewise direction are shown in Fig. 19(a), as a function of the inclination angle. These values were computed by means of the rms of every signal. For the BASELINE, a clear amplification region was observed between 47.5 ○ and 60 ○ with a maximum peak-to-peak amplitude of 2.3 m. This value is assumed to be high enough to entail a potential risk for the integrity of the considered blade. The EXT, DIHE, and COMB configurations showed a similar behavior, but the amplification region seemed to be shifted toward higher inclination angles. Additionally, the maximum peak-to-peak amplitude computed for these geometries was found to be lower than the one recorded for BASELINE. However, the predicted maximum deflections could slightly vary when considering longer simulation times, since the amplitudes for some of the simulations at 55 ○ were still increasing (see Fig. 18). No significant amplification of the vibrations was observed for the SWEEP configuration. Figure 19(b) shows the frequencies of the edgewise displacement in the range where high vibrations were identified (fx), computed based on the last periods of the signals. The results were normalized by means of the frequency of the first edgewise mode fn for every geometry, whose values were previously listed in Table II. All the aeroelastic frequencies were found to be close to the corresponding structural natural frequency (fx/fn ≃ 1), but no clear pattern identifying the onset of the amplification region was observed. It is worth mentioning that the amplification of the blade vibrations was restricted to the edgewise direction. For the particular case of the BASELINE geometry at I = 50 ○ , a maximum peak-to-peak value of 0.42 m was observed for the flapwise coordinate, and the torsion showed variations of less than 0.8 ○ .

Blade vibrations
For all the geometries, the decrease in the inclination led to a reduction in the blade bending at the tip. This is directly related to the magnitude of the cross-sectional velocity. Figure 20(a) quantifies this effect via the mean flapwise coordinate computed for every simulation. The decrease in the inclination angle also implied a reduction in the mean torsion at the tip, as depicted in Fig. 20(b).

Blade loading for the BASELINE geometry
At low inclination angles (e.g., 0 ○ , 10 ○ , and 40 ○ ), where the observed amplitudes of vibration were very small, the loading computed for the BASELINE flexible blade was found to be very similar to the one reported for the stiff counterpart. Figure 21 shows the PSDs of the sectional loading in the edgewise direction at I = 50 ○ , where a completely different scenario was found. Indeed, most of the sections exhibited a significantly high harmonic loading, indicating the existence of VIV. At this inclination angle, the PSDs of the different spanwise positions revealed the same main frequency with a non-negligible contribution of the corresponding subharmonics. These main loading frequencies corresponded indeed to the frequencies of motion fx, previously shown in Fig. 19(b). This fact points toward the existence of a lock-in between the blade deflection and the fluid excitation.

Physics of Fluids
As done for the stiff configuration, Fig. 22 depicts the evolution of the peak frequencies fp along the outer half of the blade for the FSI simulations. For inclination angles of 0 ○ and 10 ○ , the predicted peak frequencies were close to those already reported for the stiff counterparts. This implies that the shedding frequencies remained unaltered, being centered around reduced frequencies of approximately 0.16. For I = 40 ○ , the values of the PSD dropped, since the flow around the blade was found to be stationary. The aforementioned lock-in condition, where the peak frequency aligns with the frequency of motion, was found for all the inflow conditions that led to high amplitudes of vibrations (i.e., 47.5 ○ ≤ I ≤ 60 ○ ). This is exemplified in Fig. 22 by including the particular cases of 50 ○ and 55 ○ .

Blade loading of tip extensions
As for BASELINE, the loading computed at low I for the different tip extensions was not significantly affected by the consideration of the blade flexibility, so it is not documented here. Figure 23 shows the PSDs of the edgewise loading for two particular sections, where the results of all the geometry variants at I = 50 ○ were included. At the blade tip, only BASELINE (i.e., the geometry showing the highest amplitude of vibration at this angle) observed a clear lock-in. At y = 70 m, this condition was also found for the straight tip extension EXT, and broadband spectra were computed for the geometries accounting for a dihedral angle (i.e., COMB and DIHE). The latter statement is assumed to be related to the early shedding around the tip, already observed for the stiff configuration. No significant load fluctuations were observed for SWEEP at the outer half of the blade. Figure 24 includes the analogous plots for the case of I = 55 ○ . At the tip, all the simulations apart from SWEEP had a significant frequency content around the frequency of motion. For the case of COMB and DIHE, this was accompanied by additional contributions at higher frequencies. The lock-in of BASELINE, EXT, DIHE, and COMB was even more clear when studying the loads at y = 70 m. In summary, the study of the sectional loading for the different geometry variants establishes a direct connection between the lock-in condition and the high amplitudes of blade vibrations previously reported in Fig. 19(a). observed vortical structures were equivalent to those reported for the stiff configuration. The same remark applies for 40 ○ and 70 ○ that are two inclination angles accounting for a significant stationary wake region but where no lock-in was observed [Figs. 25(c) and 25(f)]. This is further illustrated in Fig. 26 for the case I = 40 ○ , where contours of high vorticity for both setups and at different sections are superposed. It can be seen that the spatial location of the flexible wake is slightly modified by the flapwise deflection, but the relative distance to the blade surface remained unaltered. For I = 50 ○ and I = 55 ○ , i.e., the inclination angles where lock-in was reported, a reconfiguration of the wake took place [Figs. 25(d) and 25(e)]. The stationary region computed for the stiff conditions was replaced by an unsteady system of vortical structures that were shed from the blade. This mechanism seemed to be synchronized all along the span, resulting in an almost parallel shedding with respect to the blade axis. A more detailed visualization of this process for the particular case of I = 50 ○ is depicted in Fig. 27 by means of the superposition of high vorticity contours for the stiff and flexible configurations. The shedding pattern corresponded to the 2S mode, with two single vortices being shed per oscillation cycle (Williamson and Roshko, 1988) (similarly to the traditional von Karman street).

Wakes of tip extensions
The similarity between the stiff and flexible wakes at low inclinations found for BASELINE was also reproduced when considering the rest of the geometry variants, so they are not documented here. Figure 28 depicts the wakes at I = 50 ○ for all the considered tip extensions by means of the iso-surfaces of the Q-criterion. The wake reconfiguration leading to parallel blade shedding, previously reported for BASELINE, was only observed for the EXT geome- try. This is in line with the high amplitudes of vibrations that were predicted for that configuration [see Fig. 19(a)] and with the corresponding harmonic loading (Fig. 23). The wakes of the other tip extensions corresponded indeed to the ones predicted for the stiff ARTICLE scitation.org/journal/phf geometry exhibited a broad band spectrum, without any identifiable main frequency. Figure 29 includes a visualization of the wakes when increasing the inclination angle up to 55 ○ by means of the iso-surfaces of the Q-criterion. In this case, all the geometry variants with the exception of SWEEP did exhibit a blade parallel shedding. This implies a wake reconfiguration with respect to the stiff conditions [see Fig. 17(a)], where an established stationary region was predicted for all the geometry variants. This parallel shedding is again in line with the high amplitude of vibrations predicted for all the blade shapes but for SWEEP and with the harmonic loading attached to them.

Concluding remarks about the coupling mechanism
The flexible results presented in this work reveal the existence of a lock-in phenomenon between the edgewise motion and the flow shedding. This lock-in was found to be a necessary condition in order to reach high amplitudes of edgewise vibration. This justifies the categorization of the studied aeroelastic problem as VIV. However, it should be clarified at this point that the lock-in observed in the present work significantly differs from the classical observations of this phenomenon in other applications, such as the two-dimensional cylinder in transverse oscillations. In this case, the lock-in appeared when the structural frequency was close to the natural shedding frequency. For the present study, no characteristic frequencies could be predicted by performing stiff simulations of the inflow conditions susceptible to VIV due to the existence of a clearly established stationary wake. This situation when a stationary wake is computed in stiff conditions but a parallel blade shedding is reported for the flexible case has been referred to as the wake reconfiguration in this work. The inflow conditions where the wake reconfiguration takes place were found to be significantly influenced by the tip geometry. Indeed, all the considered tip While the particularities of the mechanism triggering the lockin condition are not fully understood by the authors, the results of this work allow us to make some initial hypothesis that should be re-assessed in future studies. In this context, one of the plausible explanations for the appearance of the wake reconfiguration could revolve around the estimation of the natural shedding frequencies in the absence of a stationary wake. Indeed, if the appearance of this wake regime at high inclination angles is disregarded, one could still apply Eq. (7) with an imposed Strouhal number in order to assess how the shedding frequencies compare with the first edgewise frequency. The results of such academic exercise for the particular case of BASELINE are presented in Fig. 30. The imposed Strouhal number was 0.11 (corresponding to the reduced frequencies reported for the stiff DIHE at I = 50 ○ ), and the blade shape of the flexible configuration at the end of the simulation was used for the computation of U ref and l ref . It can be observed that for the inclination angles I = 50 ○ and I = 55 ○ , the predicted shedding frequencies for the high span region are very close to the first edgewise frequency (fp/fn ≃ 1), corresponding to the traditional lock-in condition.
Assuming that the tip region determines the overall blade response, as it can be inferred from the results of this work, the estimated frequencies in the tip vicinity could then be used as an indicator of the susceptibility of the system to undergo VIV. Figure 31 uses this assumption in order to quantify the effect of the different geometry variants by using the values of fp that were predicted 2 m inboard of the blade tip. Additionally, this figure presents the values of U ref and l ref involved in the computation of fp, as well as the evolution of the local inclination angle and the angle of attack (see Fig. 11). The comparison of BASELINE and EXT results points toward a delay, in terms of inclination angle, of the correspondence between the shedding frequency and edgewise frequency for the latter geometry. This can be explained by the chord distribution used for all the extensions (that assumed more slender tips), resulting in lower values of l ref . The effect of adding a dihedral angle to the straight extension, corresponding to the DIHE configuration, seemed to place the point where fp/fn ≃ 1 at even higher inclination angles. This could be explained by the higher U ref experienced by this Phys. Fluids 32, 065104 (2020); doi: 10.1063/5.0004005 32, 065104-21 ARTICLE scitation.org/journal/phf geometry, which was, in turn, related to the increase in the local inclination angle at high span locations. The consideration of a sweep angle in the SWEEP configuration resulted in an asymptotic behavior of the inclination dependency curve. This is related to the drastic reduction in l ref at high I values, which was related, in turn, to the decrease in the angle of attack seen by the high span locations. The SWEEP tip could then be seen as a shedding suppressing device for high inclination angles, since its sections tended to be aligned with the incoming flow. Finally, the results computed for the COMB configuration, which accounted for both sweep and dihedral angles, seemed to be dominated by the geometrical effects attached to the latter design variable. It is worth mentioning that the same conclusions could be made when assuming a Strouhal number of 0.16, which was the value predicted for low inclination angles, or when basing the analysis on the original geometries rather than on the deflected shapes. While the employment of such a simplified relation for the prediction of VIV regions is clearly limited, the exposed geometrical considerations could explain the shifting of the areas of lock-in identified in the present study for EXT, DIHE, and COMB geometries, as well as the absence of high amplitudes of vibrations for SWEEP. However, the fact that no natural shedding was observed for the stiff simulations is still questioning its applicability. Plausible explanations for this issue could revolve around the inherent limitations of the pure CFD simulations. In this context, it should be stated that the consideration of completely stiff blades prevented the mean deflection of the blades to come into play. This could have a substantial effect on the wakes predicted for the high inclination angles, which could eventually exhibit natural shedding, as it was predicted, for instance, for the COMB configuration at I = 50 ○ . To test this hypothesis, the stiff simulations should then be seen as an a posteriori assessment of the characteristics of the flow, rather than a preliminary analysis, as they were understood in the present study. Alternatively, the mechanism triggering the wake reconfiguration could be related to the appearance of motion-induced vortices that are obviously disregarded in the stiff simulations. While the role of these vortices in the VIV is still unclear, several authors have reported the impact that structural motion can have on the reconfiguration of the wake. Along these lines are the works of Tamura and Itoh (1999) and Tamura and Dias (2003) that explained the low speed galloping of rectangular cylinders by the motion-induced vortices generated by the structural displacements. The research of Bourguet and Triantafyllou (2014) and Bourguet and Triantafyllou (2016) on inclined cylinders also highlighted the limitations of pure CFD simulations for the prediction of VIV, which could only be understood by the reconfiguration of the wake when allowing the structure to vibrate. In the study of the concave cylinder included in Seyed-Aghazadeh et al. (2015), the authors reported the reappearance of the vortex shedding when the elasticity of the system was enabled. Finally, Ma et al. (2015) related the wing rock phenomenon to the reconfiguration of the stationary wake into an unsteady one. Assuming that the appearance of such motioninduced vortices is responsible for the lock-in condition identified in the present study, the necessary conditions for their generation remain the open question. Based on the results of this work, one of the potential factors leading to the generation of the motioninduced vortices could be related to the tip wake. In particular, all the stiff simulations accounting for inflow conditions located at the onset of the VIV region reported the presence of a pair of leeward vortex tubes located very close to the tip surface, as concluded from the detailed study found in Appendix B. With an increase in the inclination angle, these vortical structures remained attached to the blade surface, but their strength was significantly attenuated. Hence, the conditions for the appearance of the motioninduced vortices could be related to a trade-off between the proximity of the tip vortices and their strength, which could be already characterized with pure CFD simulations. The assessment of such hypothesis would require additional research. Since the aforementioned mechanism seems to be independent of the frequency of motion, this future work could include the evaluation of the sensitivity of the aeroelastic response with regard to the natural frequencies of the structure. Alternatively, the study of FSI simulations where a harmonic blade motion is imposed could also reveal valuable information.

V. GENERAL CONCLUSIONS
The present work addresses the problem of vortex induced vibrations (VIVs) of wind turbine blades in standstill. In order to progress in the understanding of this phenomenon, both a comprehensive literature review and a numerical study were carried out. The focus of this study was put in assessing the influence of geometrical modifications of the blade tip on the predicted aeroelastic response. Five different variants of a 96 m long wind turbine blade were used for this purpose. The first one corresponded to the original geometry of the IEA 10MW, and it was labeled BASELINE in the present work. The second one corresponded to a straight tip extension of 5% in the projected length, and it was referred to as EXT. The third variant considered in this study, SWEEP, was a tip extension with the same length, but it also accounted for a sweep angle of 20 ○ . The fourth extension, DIHE, consisted of a similar tip extension, but with a dihedral angle of 20 ○ . Finally, the COMBINED geometry variant accounted for a tip extension where both a dihedral angle and a sweep angle were taken into account.
The aeroelastic response of the different blade variants was assessed by means of fluid structure interaction (FSI) simulations. These FSI computations, labeled flexible, coupled a computational fluid dynamics (CFD) solver with a multi-body structural solver. A range of inclination angles I, understood here as the relative angle between the freestream velocity and the blade axis, were considered. The aim of such a parametric study was the identification of inflow conditions triggering the VIV phenomenon. In addition to the aforementioned FSI simulations, a set of pure CFD computations of the undeformed blade variants were also performed in order to have a better insight into the flow. These simulations were referred to as stiff in this document.
The results of the stiff simulations of the BASELINE geometry revealed a strong dependency of the wake with respect to the inclination angle. At low values of I, a clear shedding was identified for most part of the blade span, corresponding to the wake configuration labeled unsteady in this document. This led to a broadband edgewise loading spectrum with main frequencies centered around Strouhal numbers of approximately 0.16. For inclination angles higher than 40 ○ , the wake for most part of the blade span consisted of a system of stationary leeward vortex tubes. This wake configuration, referred to as stationary in this document, resulted in the complete mitigation of the dynamic blade loading. Similar conclusions could be made when studying the solutions of the stiff simulations for the geometries accounting for tip extensions. However, significant changes in the wakes at high inclination were reported, especially in the vicinity of the tip. Particularly relevant was the effect of the dihedral angle, which led to the early breakdown of the leeward vortex tubes emanated from the tip. Under certain inflow conditions, a generalized vortex shedding all along the span was also reported for the case of the COMBINED configuration, which resulted in a broadband edgewise loading with important contributions around a Strouhal number of 0.11. Additionally, the results of the SWEEP configuration at high inclination angles showed that one of the tip vortex tubes was rapidly convected away from the blade surface, and the detachment of the other one was considerably delayed.
The results of the flexible simulations of the BASELINE geometry showed significant edgewise vibrations for inclination angles between 47.5 ○ and 60 ○ . These were in the form of limit cycle oscillations, and a maximum peak-to-peak displacement of 2.3 m was predicted. Such an amplitude of vibration, if kept for a certain amount of time, is assumed to be high enough to entail a potential risk for the structural integrity of the blade. For the inflow conditions exhibiting high amplitudes of vibrations, an important harmonic loading was reported, where the main frequency corresponded to the frequency of motion. This lock-in condition between the fluid loading and the blade displacements was always accompanied by a reconfiguration of the stationary wake region into a parallel vortex shedding. This shedding was also observed at the tip vicinity, and it followed a 2S pattern. Since these observations revealed a clear relation between the vortex excitation and the structural displacements, the studied phenomenon was classified as a VIV phenomenon by the authors. Outside of the region where these VIVs were reported, the wakes predicted by the flexible and stiff simulations were found to be essentially equivalent. Similar lock-in regions were also found for EXT, DIHE, and COMBINED, but the range of concerned inclination angles seemed to be shifted toward higher values. For the particular case of SWEEP, no VIVs were reported for any of the inflow conditions considered for the present work. By assuming a common Strouhal number, the authors attributed these modifications of the VIV region to the geometrical particularities of each of the tips. In this way, the shifting toward higher inclination angles experienced by EXT was related to the lower chord of the airfoils located close to the tip. The dihedral angle introduced in DIHE was found to work in the same direction due to the increase in the cross-sectional velocities seen by the high span locations. The suppression of the VIV observed for SWEEP was assumed to be related to the angle of attack experienced by the airfoils of the tip vicinity, which tended to 0 for high inclination angles. For the COMBINED geometry, the effects introduced by the dihedral angle were found to be dominant with respect to the ones attached to the sweep angle. While these geometrical considerations could explain the modifications of the areas where VIV was observed, the mechanism underlying the wake reconfiguration was not completely understood by the authors. As a preliminary explanation, two of the limitations of the stiff simulations were pointed out. On one hand, the pure CFD computations did not take into account the mean deflection of the blade, which could have a non-negligible impact on the computed wake. On the other hand, a plausible cause for the wake reconfiguration could rely on the appearance of motion-induced vortices, which could be triggered depending on the proximity and strength of the leeward vortex tubes emanated from the tip. Additional research efforts should be conducted in order to check both hypotheses.
The contribution of the present work to the literature is twofold. On one hand, it constitutes the first study of wind turbine blades' VIV that includes a detailed analysis of the wakes. This allowed us to characterize the complexity of the problem from a flow perspective, putting into evidence the challenges of building a reduced order model targeting this application. On the other hand, the results of the present research showed that the introduction of different tip geometries has a significant influence on the resulting VIV. This observation implies that, for the optimization and design loops of large wind turbine blades, the susceptibility of every new candidate to undergo VIV should be re-assessed. In addition, it points toward the design of the tip being a potential solution for the alleviation of VIV on new or already existing blades.
The main limitations of the present study are related to the significant computational resources required for the type of analysis that was performed. This forced the amount of inflow conditions to be constrained to a particular subset. It was done by fixing the absolute value of the freestream velocity as well as the initial angle of attack of the different blade sections. While the conclusions made in the present work may hold for other values of these two parameters, supplementary simulations will be needed in order to make definite statements about the stability of the different tips with regard to VIV. In this context, it is also worth mentioning that the focus of this work was put in the vibrations at high inclination angles (i.e., >40 ○ ), while previous experiences have also reported considerable edgewise motion amplitudes for angles of around 20 ○ . Additionally, the simulations performed within the framework of the present study assumed a uniform inflow. It could then be interesting in the future to assess the impact of inflow turbulence in the predicted areas of VIV, as well as the influence of the atmospheric boundary layer. Furthermore, the scope of the present work was limited to the study of VIV of wind turbine blades, disregarding the other components of the machine. Recent publications have shown the plausibility of wind turbine towers to undergo this phenomenon (Livanos, 2018;Viré et al., 2019). The same numerical methods used in the present work could then be applied to the study of wind turbine towers' VIV as a future work. Finally, it should be emphasized that the present study was purely numerical due to the absence of relevant experimental data. While individual validations of the different codes involved in the FSI framework were previously performed, the generation of experiments targeting VIV of wind turbine blades would be important to corroborate the existence of the phenomena described in this document.

ACKNOWLEDGMENTS
This work was supported by the Smart Tip project, funded by Innovation Fund Denmark (J.nr. 7046-00023B).

APPENDIX A: SENSITIVITY STUDIES
Three independent studies were performed in order to assess the sensitivity of the numerical model described in Sec. III C when performing flexible simulations: • Time step sensitivity: Simulations accounting for different time steps dt were performed, and their results were compared with the original configuration (where dt was set to 6 ms). • Grid sensitivity: Additional meshes were created for this purpose based on the BASE grid resolution described in Sec. III C 1. The mesh labeled FINE was generated by doubling the resolution of the BASE grid in all directions. The COARSE1 mesh was generated by halving the resolution of the BASE grid in all directions, and the COARSE2 mesh was generated by repeating the coarsening procedure for the COARSE1 mesh. • Total time sensitivity: Longer simulations were run in order to study the evolution of the solution beyond t = 45.3 s.
The parameters describing each of the aforementioned simulations, together with the labels used in this section, are presented in Table III. The wall clock time W reported by the cluster Jess is also included. It should be noted that W was mainly dominated by the solution of the Navier-Stokes equations. A representative breakdown of the time spent by each of the components of the FSI coupling described in Sec. III A 3 is 99.05% for EllipSys3D, 0.9% for HAWC2, and 0.05% for the DTU coupling. Figure 32 depicts the time series of the tip edgewise position for each of the sensitivity studies. For brevity, only the results of the BASELINE geometry at I = 50 ○ are included. This corresponds to the case where the maximum VIVs were recorded, but similar observations were made for other configurations. The time step sensitivity study [ Fig. 32(a)] revealed similar amplitudes of vibrations when dt was in the range (3, 12) ms, with relative differences of less than 5%. A similar discrepancy was found when comparing the FINE and  Figure 32(c) shows that, with the total simulated time chosen for this particular case, the limit cycle oscillation was already reached.
In summary, the model parameters chosen for the present study offered a good trade-off between the accuracy of the results and the required computational resources. While small discrepancies in the predicted amplitudes of vibrations are expected, the chosen setup was found to be very reliable when distinguishing between VIV configurations from non-vibrating cases.

APPENDIX B: DETAILED STUDY OF TIP WAKES
This appendix presents a study of the wakes around the tip for the selected simulations. It provides complementary observations to those exposed in the main part of this manuscript, which does not include any detailed analysis of the flow around the tips. This study is ARTICLE scitation.org/journal/phf split into several sections, allowing the assessment of the influence of the tip extensions and the blade flexibility in an incremental manner. In this way, Appendix B 1 presents the tip wakes for the BASELINE geometry in the stiff configuration and describes how the topology is modified when varying the inclination angle. In Appendix B 2, the changes observed in the tip wakes due to the consideration of the different tip extensions is discussed. Appendix B 3 shows how the blade flexibility altered the tip wakes computed for the BASELINE geometry, relating the variations to the appearance of the lock-in condition. Appendix B 4 presents some of the observations made for the wakes of the tip extensions in the flexible configuration. Finally, some concluding remarks are included in Appendix B 5. For all the cases, the tip wakes are presented by means of the iso-surfaces of vorticity, which were clipped at a certain high span location (i.e., y value). Figure 33 depicts the tip wake of the BASELINE geometry in the stiff configuration for several inclination angles. At low inclination angles (i.e., I = 0 ○ and I = 10 ○ ), a part of the wake around the tip corresponded to the unsteady wake region. This area, labeled SH, is characterized by the periodic shedding of vortical structures. At the outer most span positions, the presence of a stationary pair of vortex tubes, which remained close to the blade surface, was also observed. These were formed from the roll-up of the vortex sheets around the leading edge and the tailing edge forming two distinct cores. These structures, labeled LE0 and TE0 in Fig. 33(a), had opposite signs of vertical vorticity due to the inversion of the roll-up direction. The presence of TE0 and LE0 even at I = 0 ○ could be attributed to the local curvature induced by the blade prebend, as illustrated in Fig. 4(b). Finally, a small vortex was also identified at the very tip. This tip vortex (TV) also remained stationary during the simulation. For high inclination angles (e.g., I = 40 ○ ), the onset of the unsteady wake region was shifted toward lower span locations, and all the observed structures around the tip were stationary. Together with LE0, TE0, and TV, two additional leeward vortex tubes generated from the trailing edge and leading edge of the blade were observed. These structures are labeled LE1 and TE1 in Fig. 33(c), respectively. The change in the span location where LE1 and TE1 were generated could be related to the considered inflow conditions, which assumed an angle of attack of 100 ○ at the tip. Additionally, geometrical aspects such as the blade sweep or the twist distribution could also play a role in this behavior. It should be remarked that the same process described for the generation of LE0, TE0, LE1, and TE1 was repeated all along the blade span, resulting in the region that was referred to as the stationary wake in the present work. The comparison of the results for the inclination angles I = 40 ○ , I = 50 ○ , I = 55 ○ , and I = 70 ○ allowed us to characterize the influence of the inclination angle on the aforementioned region. In particular, two different mechanisms could be attributed to an increase in I. On one hand, the leeward vortex tubes remained closer to the blade surface and for a longer span range. On the other hand, the strength of these vortices was gradually attenuated.

Effect of tip extensions in stiff conditions
The tip wakes behind the considered tip extensions were found to be substantially different. Figure 34  situation at low inclination angles (i.e., I = 0 ○ ). As for the BASE-LINE wake previously shown in Fig. 33(a), the leeward vortex tubes emanated from the tip of the EXT and SWEEP configurations prevented the outermost span locations to shed. This situation was changed when considering the dihedral angle (that was included in the DIHE and COMB cases). For these geometries, significant vortical structures were shed even at span locations very close to the tip, probably due to the increase in the local inclination angles for this region. The tip wakes of the different tip extensions for a high inclination angle (i.e., I = 55 ○ ) are depicted in Fig. 35. The wake of EXT was found to be equivalent to the one computed for the BASELINE geometry, with LE0 and TE0 remaining close to the tip geometry. For SWEEP, LE0 was convected away from the blade surface, and the detachment of TE0 was considerably delayed. The wake of DIHE accounted for the early breakdown of the vortex tubes around the tip probably due to the change in the local inclination angle induced by the dihedral angle. Additionally, both LE0 and TE0 seemed to be convected away from the blade surface faster than for the case of EXT. The same remarks apply for the COMB geometry variant. In addition, the early departure of LE0 was also observed for this case (since this geometry also accounted for a sweep angle).
The same vortical structures exhibited by the tip extensions at I = 55 ○ were observed at moderately lower inclination angles (e.g., 50 ○ ). Consistently to the BASELINE case, the distance between the leeward vortex tubes and the blade surface was slightly increased, and their strength was augmented.

BASELINE in flex conditions
As for the rest of the blade, the tip wakes for the BASELINE geometry in the flexible configuration were equivalent to the stiff counterparts for two different groups of simulations. On one hand, the tip wakes at low inclination angles exhibited the same topology for TE0, LE0, and TV, and the vortical structures observed in the shedding region SH were also found to be very similar. On the other hand, the topology of the tip wakes for high inclination angles where lock-in was not observed reproduced the patterns from the stiff configurations. It should also be remarked that for the latter case, the comparison of the results of different time instants revealed that the tip wakes were also completely stationary when considering the blade flexibility.
It is then interesting to have a better insight into the tip wakes for the simulations at high inclination angles where the lockin condition was indeed reported. This is illustrated in Fig. 36 for the particular case of I = 50 ○ , showing the wake at different time instants. The leeward vortex tubes LE0 and TE0, which remained close to the blade surface for the stiff case, were periodically released when enabling the blade flexibility. This release mechanism was synchronized with the structural displacement, following a 2S pattern. Hence, it is concluded that the vortical structures around the tip also contributed to the shedding characterizing the wake reconfiguration.

Effect of tip extensions in flexible conditions
The equivalence between the tip wakes of the flexible and stiff configurations, for the simulations that did not experience lock-in, was also observed for the tip extensions. It is then interesting to analyze the inclination angle 55 ○ , where all the considered geometries but SWEEP experienced lock-in. Figure 37 presents the wakes of the different tip extensions at selected time instants in order to highlight the effects of the wake reconfiguration. As already stated, the tip wake of the SWEEP geometry was found to be equivalent to the one computed in the stiff configuration, as previously shown in Fig. 35(b). For the other variants, the vortex tubes emanated from the tip vicinity were found to contribute to the generalized blade shedding. The release mechanism behind this unsteady process was found to be equivalent to the one reported for the BASELINE geometry, with two vortices being shed per period.

Summary of observations
The wake of the BASELINE stiff geometry always accounted for a stationary system of three vortices at the very tip: a leading edge vortex LE0, a trailing edge vortex TE0, and a small tip vortex TV. Depending on the considered inclination angle, the wake at lower span locations was characterized by either a series of stationary leeward vortex tubes emanated from the edges or the unsteady shedding of vortical structures. For both cases, the local modification of the blade geometry through the installation of a tip extension implied substantial changes in the wake. At high inclination angles, where some of the configurations were susceptible to undergo VIV, the introduction of a dihedral angle led to the early breakdown of LE0 and TE0. The effects of the introduction of a sweep angle were also noticeable in these conditions, implying the convection of LE0 away from the blade surface. The analysis of the tip wakes for the flexible simulations reporting the lock-in condition revealed that the tip also contributed to the generalized shedding, regardless of the considered geometry. Even if no clear relation between the tip wake topology and the appearance of the aforementioned phenomenon could be established, the stiff wakes of all the configurations that led to high amplitudes of vibrations did share a common feature. In particular, a pair of leeward vortex tubes, which remained very close to the blade tip surface, was observed. In this context, a possible mechanism for preventing SWEEP to undergo high amplitudes of vibrations could be related to the effects that the sweep angle had on LE0. For the rest of the geometries, the trade-off between the proximity and the strength of the LE0-TE0 system could be hypothesized as the steering mechanism behind the lock-in. However, no quantification of this premise was performed within the framework of the present study, requiring additional analysis.

APPENDIX C: TIME SERIES
This appendix includes most part of the time series that were used to generate the PSDs of the sectional loading included in the main part of this manuscript. A reference to each of these figures below was included in the legend of the corresponding spectra.

DATA AVAILABILITY
The data that support the findings of this study are available within this article.