Influence of the installation of a trailing edge flap on the vortex induced vibrations of a wind turbine blade

This work explores the use of a trailing edge flap for the mitigation of the Vortex Induced Vibrations (VIV) of a wind turbine blade in standstill. Fluid structure interaction simulations are performed, coupling a computational fluid dynamics code with a multi-body structural solver. A publicly available 96 m long wind turbine blade is selected for the study, and labeled as IEA 10MW. The effects of a flap actuation in the range [ − 90 ◦ , 90 ◦ ] are analyzed, for a given set of inflow conditions, and compared to the aeroelastic response of the original IEA 10MW geometry. While the consideration of the flap does lead to reductions on the amplitude of the edgewise vibrations, these are only noticeable when operating outside of the design limits of the aerodynamic device. A simple non-dimensional analysis seems to indicate that the potential of the flap installation with regards to VIV mitigation is very limited, as the instability is only shifted towards different inflow conditions and the expected effects on the maximum amplitude of vibrations are minor.


Introduction
When horizontal axis wind turbines adopt the standstill configuration, due to for instance extreme wind or scheduled maintenance, the blades are rotated approximately 90°to aerodynamically brake the rotor. The standstill configuration entails a potential risk to the structural integrity of the blades, as current designs are presumed to experience large vibrations for particular configurations involving massive flow separation. This has motivated several research efforts, aiming at the characterization of the aeroelastic instabilities involved and the mitigation of their adverse effects. The main factor challenging these research activities is the complexity of this coupled problem. In particular, the application combines (i) high Reynolds numbers (ii) highly flexible structures and (iii) complex and variable inflow. Additionally and unlike other traditional applications involving aeroelastic instabilities in massive flow separation, it is still uncertain if a two-dimensional configuration could accurately reproduce the phenomenon. For all the aforementioned reasons, the design of ad hoc wind tunnel experiments requires large facilities and a careful aerodynamic and aeroelastic scaling. Moreover, the use of field data for the analysis of the phenomenon may be challenging, since it can be difficult to fully characterize the environmental conditions during an experimental campaign. While the acquisition of experimental data is essential for the proper characterization of the physics of standstill vibrations, several researchers have already approached the problem by means of numerical simulations, as a first exploration of the phenomenon. In particular, previous work * Corresponding author.
relied on the use of Fluid Structure Interaction (FSI) methods, which were previously validated for other relevant applications, and were able to predict significant vibrations under particular inflow conditions. The first of this category was the study of Heinz et al. (2016b), where the 86 m long wind turbine blade of the DTU 10 MW design (Bak et al., 2013) was studied. It was shown that when the incoming flow led to an angle of attack in the vicinity of 90°, due to e.g. a failure of the yaw system, the structure could undergo significant limit cycle oscillations in the edgewise direction. The authors related the appearance of these vibrations to the action of the vortex shedding, categorizing the phenomenon as Vortex Induced Vibrations (VIV). Similar conclusions were made in Heinz et al. (2016a), based on a set of simulations performed on the 100 m long AVATAR blade (Lekou et al., 2015). The work of Horcas et al. (2020) included a comprehensive study on the VIV of the IEA 10 MW rotor blade (Bortolotti et al., 2019), which has a blade length of 96 m, and characterized the problem from a flow and geometrical perspectives. In special, the aeroelastic instability seemed to be significantly affected by the relative angles of the incoming flow, by the blade curvature, and by the tip geometry. The authors pointed out that these particularities may make wind turbine VIV a relatively special application, when compared to other traditional fields such as bridge decks (Wang and Chen, 2022;Bai et al., 2021;Yang et al., 2021;Xue et al., 2021;Cao et al., 2020) or circular cylinders (Guo et al., 2021;Chen et al., 2020;Deng et al., 2011;Wang et al., 2019). The https://doi.org/10.1016/j.jweia.2022.105118 Received 30 March 2022; Received in revised form 27 July 2022; Accepted 27 July 2022 complexity of wind turbine blades vibrations in standstill was further explored in the recent simulations by Horcas et al. (2022), where a comprehensive analysis of the aeroelastic mechanism triggering the instability was carried out. The authors analyzed the relation between the transferred power and the wake topology for a range of inflow angles. It was concluded that the total transferred power, and therefore the amplitude of vibration, was the result of a combination of oblique shedding and parallel shedding. Additionally, the inflow angles leading to large amplitudes of vibrations were clustered around specific regions. This seemed to establish a relation between the inflow conditions leading to large amplitudes of motion and the natural tendency to shed of the different blade sections, even if the appearance of VIV was accompanied by a complex wake reconfiguration.
While a better understanding of the VIV mechanism could directly improve new machine designs with regards to the susceptibility of vibrations, it is also relevant to explore the use of mitigating solutions for existing blades. In particular, the present work addresses the potential benefits of a trailing edge flap installation. These aerodynamic devices are already investigated for wind turbine control applications, due to their load alleviation capabilities in operation, with early work on testing as rotor speed regulation devices (Miller et al., 1998), later on aeroservoelasticity (Montinari et al., 2018;Bernhammer et al., 2016;Bottasso et al., 2016;Barlas and van Kuik, 2010), higher fidelity numerical simulations (Heinz et al., 2011;Barlas et al., 2016;Jost et al., 2016) and full scale testing (Castaignet et al., 2014;Ai et al., 2019;Gomez Gonzalez et al., 2021). In the context of VIV mitigation, a preliminary numerical study based on FSI simulations of the AVATAR blade revealed a decrease on the amplitudes of vibration for particular flap angles and inflow conditions (Horcas et al., 2019). However, the underlying mechanisms behind this mitigation remained unknown, and the amount of simulated configurations was limited. A more comprehensive numerical campaign is included in the present work, and the implications of the flap actuation on the wake of the blade have been studied in detail. Additionally, the IEA 10 MW reference wind turbine was employed, in order to check if the observations made for the AVATAR blade could be replicated for other designs. In summary, the objective of the present work is to numerically explore if trailing edge flaps show a consistent impact on the aeroelastic response of wind turbine blades in standstill, and if their installation on current designs could be a suitable solution for VIV mitigation. That would represent the first assessment of the performance of an aerodynamic device, here understood as an add-on, for the control of such vibrations. As a first approach to the problem, the flaps would be assumed to be deployed during the whole simulation (i.e. considered to be static), and different relevant flapping angles are explored.
The remainder of the present article is organized as follows. Section 2 introduces the performed analysis, including details about the definition of the considered inflow conditions and the trailing edge flap geometry. Section 3 presents the numerical methods that were employed. This is followed by Section 4, where the numerical models that were built are described. Section 5 includes the results and discussion of the present numerical study. Finally, the conclusions are summarized in Section 6, together with some potential directions for future work.

Performed analysis
The 3-bladed wind turbine chosen for the present study, and referred to as IEA 10 MW, corresponds to the offshore design presented in Bortolotti et al. (2019). With a length of 96 m, the IEA 10 MW blade is assumed to be representative of state-of-the-art designs of horizontal axis wind turbines. The IEA 10 MW is equipped with a series of airfoils of the FFA-W3 family (Bjorck, 1990). The airfoil thickness, chord length and twist vary as a function of the radial location. The blade axis also has a considerable out-of-plane prebending, with the tip located more than 6 m out-of-plane. The IEA 10 MW blade is made from composite  materials, and its structural layout accounts for two main spars and three shear webs.
To reduce the amount of computational resources needed, the present work studied the VIV phenomenon by simulating a single blade, assuming to be clamped at the root. The cone angle and the tilt angle, both related to the way the blade is assembled with respect to the tower axis, were also disregarded for simplicity.

Inflow conditions
The definition of the inflow configurations used in this study relied on the results of Horcas et al. (2020), where the present authors already assessed the stability of the IEA 10 MW blade in standstill. In particular, a clear range of inflow conditions leading to the amplification of the edgewise displacements of the baseline geometry, in the form of limit cycle oscillations, was identified (see Fig. 1). This was achieved by the variation of the so-called inclination angle in the range [40°,70°]. This variable is defined as the angle formed between the freestream velocity and the local blade axis at the root, being null when both vectors are normal and positive when a spanwise component from tip to root is introduced. In other words, the inclination angle corresponds to the vertical angle of the freestream velocity with respect to the blade axis (see Fig. 2). In wind turbine operational terms, the inclination angle is mainly related to the azimuthal position at which the considered blade is parked. Both the cited study and the present work considered the following inflow parameters as fixed: • Pitch angle = 100°: The pitch angle is understood as the orientation of the horizontal projection of the freestream velocity ∞, with respect to the rotor plane, as depicted in Fig. 2. can be seen as an estimation of the angle of attack seeing by the untwisted sections, assuming no blade deflection. In operational terms, the pitch angle is mainly related to the yaw misalignment refers to the spanwise projected length of the blade, and to the chord at a given section. experienced by the considered blade with respect to the incoming flow.
• Freestream velocity ∞ = 18 m s −1 : This value represents a relatively high wind speed, but with a reasonable probability of occurrence. For the IEA 10 MW blade, a value of 18 m s −1 corresponds to chord Reynolds numbers in the range [1.0⋅10 6 , 5.6⋅10 6 ].
It should be noted that the use of a uniform inflow represents a simplification of the targeted problem, as for instance (i) the atmospheric boundary layer would impose a spanwise wind speed dependency (ii) wind veer may be present or (iii) the inflow may be highly unsteady due to atmospheric turbulence. This hypothesis was adopted in order to allow for a simpler CFD setup, and to facilitate the subsequent analysis of the results (so that the aeroelastic instability was characterized by a single wind speed magnitude).

Flap actuation
The trailing edge flap was dimensioned based on the results of Horcas et al. (2019), where different geometries were tested for the mitigation of the VIV of the AVATAR reference blade. The flap dimensions shown in Fig. 3 were found to be particularly effective, and they were adopted for this study. The flap is centered at 75% of the blade span, and has a spanwise length of 20% and a chordwise length of 20%. Different flap angles were explored in the present work, that were assumed to be fixed throughout the entire simulation. These angles were namely: 0°, ±10°, ±20°, ±30°, ±30°, ±70°, ±80°and ±90°. These values cover not only the expected operating limits of this technology, that usually lay in the range [−15°,+15°] Barlas et al., 2016;Gomez Gonzalez et al., 2021), but also larger angles. While this may lead to an unrealistic actuation, the larger angles were included in the present study in order to better characterize the effects of the flap on the aeroelastic response of the blade. Note that the baseline geometry of the IEA 10 MW corresponds to a flap actuation at 0°, and it is also referred to in this document as the CLEAN configuration.

Aeroelastic solver
HAWC2 (Larsen and Hansen, 2015) is a serial solver that can perform multi-physics simulations of whole wind turbine assemblies. Among other components, the tool implements an aerodynamics, structural multibody and hydrodynamics 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 aeroelastic capabilities were used. Its structural module uses a multibody formulation based on an augmented form of the floating frame of reference [see for instance Shabana (2001)], with a finite-element formulation to compute the elastic deflections of the bodies. For the present work, several sections were defined along the blade, where the aerodynamic loads are injected. 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.
Although HAWC2 uses linear Timoshenko beam elements for modeling the structural bodies, it can still capture large deflections (Gözcü and Verelst, 2020). This is achieved by representing the continuous structure by means of multiple bodies. These bodies are connected by a series of constraints, which satisfy the continuity of the structure. At every time step, HAWC2 solves the equation of motion under the presence of external loads and inertial loads, based on a predictor/corrector approach. For the present work, the external loads corresponded to the aerodynamic loads. The effects of gravity, which was aligned with the negative axis shown in Fig. 2, were also included as inertial loads. A Newmark-beta scheme was used for the time-integration (Newmark, 1959), where the coefficients and were set to 0.27 and 0.51, respectively.

Flow solver
In the present study, the three-dimensional Computational Fluid Dynamics (CFD) solver EllipSys3D (Michelsen, 1992(Michelsen, , 1994Sørensen, 1995) was used to solve the incompressible Navier-Stokes equations, and to compute the sectional aerodynamic loading of the blade. This was subsequently introduced in the aeroelastic solution of HAWC2.
EllipSys3D is based on a finite volume discretization, and the equations are solved on a structured grid in curvilinear coordinates. The flow variables are stored in a collocated grid arrangement, using Rhie-Chow interpolation (Rhie and Chow, 1983) to avoid odd/even pressure decoupling. The solution follows a pressure-based strategy, where the Poisson equation is handled by means of an iterative predictor/corrector method, relying for the present study on the implicit SIMPLE approach (Patankar and Spalding, 1972). EllipSys3D is secondorder accurate in time, by using a three-point backward scheme to advance the solution. The solution of the momentum equations in every direction is decoupled through the use of a red/black Gauss-Seidel point solver. EllipSys3D implements moving mesh capabilities through the arbitrary Lagrangian-Eulerian formulation (Hirt et al., 1974), and several turbulence models have been implemented. For the simulations of the present work, a hybrid model combining the Reynolds-Averaged Navier Stokes (RANS) model and the Large Eddy Simulation (LES) approach was used, in order to find the balance between resolved turbulent scales and computational time. In particular, the Improved Delayed Detached Eddy Simulation (IDDES) model was employed (Shur et al., 2008;Gritskevich et al., 2012). This turbulence model can be seen as an enhancement of the traditional Detached Eddy Simulations (DES) method (Spalart et al., 1997), aiming at preventing several adverse phenomena related to hybrid turbulence models such as grid induced separation and the logarithmic layer mismatch. For the region involving the RANS method, the k-SST turbulence model of Menter (1994) was employed. In EllipSys3D, the diffusive terms are discretized using a central differencing scheme. In this study, the third order accurate upwind scheme QUICK was used for the discretization for the convective terms in the RANS region (Leonard, 1979), while a fourth order central differencing scheme was employed in the LES region. EllipSys3D implements a multiblock decomposition of the domain, allowing for an efficient parallelization of the code through the use of the Message Passing Interface (MPI) standard (Forum, 1994). A multigrid solution strategy based on a five-level additive Schwarz is employed in order to accelerate the convergence of the pressure-correction equation.

FSI coupling approach
The transfer of the aerodynamic loads from EllipSys3D into HAWC2 was performed once every time step. Consistently, the blade deflections predicted by HAWC2 were also transferred back to EllipSys3D. These were interpolated into the CFD surface mesh by means of a spline based interpolation, and subsequently propagated into the volume grid through a mesh deformation algorithm. In particular, the surface 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. The implemented loose two-way FSI coupling kept the same predictor-corrector strategy of the original HAWC2 aeroelastic solver. A comprehensive description of the coupling sequence can be found in Heinz (2013), including details about the transferred variables and their integration in each of the solvers.

Previous validation efforts
Several validation campaigns were carried out around the Ellip-Sys3D flow solver, targeting applications such as flow over complex terrain (Koblitz et al., 2014;Cavar et al., 2016) or airfoil simulations, where a two-dimensional variant of the code was normally used. To the latter group belong the airfoil catalogue presented by Bertagnolio et al. (2001), the studies of passive vortex generations conducted by Baldacchino et al. (2016), or the deformable trailing edge validation campaign presented by Gomez Gonzalez et al. (2018). Concerning blade-resolved simulations EllipSys3D was validated against the experimental results of three wind turbine models in operation, namely the Mexico (Benchmann et al., 2011;Sørensen et al., 2016), the NREL Phase VI (Sørensen and Schreck, 2014;Sørensen et al., 2002) and the NM80 (Grinderslev et al., 2020c). Standstill simulations using EllipSys3D in imposed motion were validated in Sørensen and Schreck (2011), where comparisons of the loads for a full blade under massive flow separation were carried out. The structural solver of HAWC2 was verified and compared against higher fidelity models in Pavese et al. (2015). Validations of the aeroelastic solution of HAWC2 have been also carried out, including the study of operational configurations and flutter-like instabilities Volk et al., 2020). Finally, it is worth mentioning that the described EllipSyS3D-HAWC2 coupling was validated for the pull-release test of a wind turbine blade (Grinderslev et al., 2020a) and for a full flexible rotor in operation (Grinderslev et al., 2020b).

CFD model
The meshing strategy employed for the present work consisted in two subsequent steps. First, a surface mesh was created with the openly available Parametric Geometry Library (PGL) (Zahle, 2020). This grid had a total of 256 cells in the spanwise direction, and included a tip consisting of 4 blocks of 32 × 32. The chordwise direction was discretized with 256 cells, with 16 cells laying on the trailing edge. A refinement of the spanwise resolution around the limits of the flap geometry was also considered, to smoothly accommodate the deflection of the blade surface (see Fig. 4(a)). Secondly, a spherical volume grid was generated with the hyperbolic mesh generator Hypgrid (Sørensen, 1998), with outer boundaries located 7.3 blade lengths away from the surface, in order to generate a O-O topology. A total of 256 cells were used when marching in the normal direction. A boundary layer clustering was imposed, with a first cell height of 1⋅10 −6 m, to target + values lower than unity. An additional refinement region was defined in order to increase the resolution close to blade surface, extending up to 1.2 blade lengths. Fig. 4(b) depicts half of the CFD outer domain, together with the blade surface mesh for reference. The resulting volume grid contained 576 blocks of 32 × 32 × 32 cells, totaling 18.9 million cells and 20.7 million nodes. The boundary conditions at the far field were defined as an inlet (i.e. Dirichlet), except for a small opening where an outlet boundary condition was set (i.e. zero normal gradient).  The effect of the inclination angle was introduced by changing the wind direction at the inlet, keeping the blade in the upwards position. During the simulation, a three level grid sequence was used to speed up the development of the wake flow. The air density was fixed to 1.225 kg m −3 , and the dynamic viscosity was set to 1.784⋅10 −5 kg m −1 s −1 . 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, a value of 1 ⋅ 10 6 s −1 was imposed for the specific dissipation and 1 ⋅ 10 −2 m 2 s −2 for the turbulent kinetic energy.

Structural model
The blade was modeled by means of 19 Timoshenko beam elements and 10 bodies. The mechanical properties of the beam elements were obtained from the definition of the IEA 10 MW blade (Bortolotti et al., 2019). The natural frequency and the logarithmic decrement of the first four structural modes of the blade in standstill configuration, and in the presence of gravity, are listed in Table 1. These were obtained in a preliminary analysis, also performed with HAWC2.

Trailing edge flaps
The trailing edge flaps were modeled as a local modification of the CLEAN mesh. This was achieved by introducing the flap geometry in the blade surface mesh using the in-house deformation library, FFDlib [Madsen, 2020, Chapter 5] which allows for deformation of CFD meshes using the Free-From Deformation (FFD) approach (Sederberg and Parry, 1986). The basic FFD principle is to embed a geometry in a rubber-like box (gray, transparent box in Fig. 5) defined by a fixed number of control points (CPs). One can then carry out well-defined operations on the CPs such as a twist of certain box sections as seen in Fig. 5(d). The deformation will then travel inwards throughout the box and the embedded geometry will experience a similar deformation. The internal box space in FFDlib is mathematically described by a trivariate Bernstein polynomial [Madsen, 2020, Eq. 5.1]. The FFD approach has numerous salient features such as (i) exact mesh representation, (ii) analytical gradients, and (iii) -continuity control. Most of the flap angles were included during the simulation, so that the inner grid was re-adapted based on mesh deformation, thus allowing e.g the use of dynamic flapping in future studies. However, for the particular flap angles of ±70°, ±80°and ±90°, the input volume mesh of the CFD solver was re-generated with the considered flap geometry, in order to prevent negative volume cells from occurring when dealing with large local rotations. Fig. 5 shows the FFD box fitted to the CLEAN blade and how a trailing edge flap for a −90°flap activation can be generated with the FFD box. As is seen in Fig. 5(d), the −90°flap angle is gradually imposed on FFD box sections in the chordwise direction towards the trailing edge, where the −90°are reached. This ensures a −90°flap angle is imposed at the trailing edge, while the resulting geometry mimics a deformable flap.
The trailing edge geometries at the flap center (i.e. 75% of the projected length of the blade in spanwise) are depicted in Fig. 6(a). These geometries are shown based on a coordinate system which origin corresponds to the leading edge, and are normalized by the chord of the CLEAN blade at that location . The corresponding cross-sectional meshes for the particular cases of the CLEAN configuration, and the flap actuations at ±90°, are included in Figs. 6(b), 6(c), and 6(d). For clarity purposes, only 1 out of 4 grid lines are shown for all meshes. It can be noted that actual deformation is only visible for the aft 15% of the chord in Fig. 6(a) although the FFD box (gray, transparent in Fig. 5) indeed covers the entire indicated flap area seen in Fig. 3 embedding all CFD mesh points in the last 20% of the chord. The explanation for this apparent discrepancy is the 2 -continuity condition which was imposed in FFDlib in the chordwise direction to ensure a smooth transition to the non-deforming CFD mesh region. As seen in Fig. 6(a) the resulting deformation is indeed exhibiting a very smooth transition to the non-deforming region ensuring no sharp kinks occur in the surface.

Models integration
A total of 101 aerodynamic sections were defined along the blade, where the integrated CFD loads were injected into the HAWC2 solution. These sections were set to be normal to the local orientation of the blade axis, and were equally spaced in terms of its length (Fig. 7). A time step length of 6 ⋅ 10 −3 s was chosen, based on the conclusions made in Horcas et al. (2020) for a similar mesh under the same inflow conditions. Each simulation started from rest, without any structural excitation rather than the application of the gravity loads at the first time step, and the freestream velocity was assumed in the farfield as an initial condition. The computations were run up to 45.3 s, a total simulated time that was defined based on the available computational resources for the present work. Preliminary simulations were performed, aiming at assessing the sensitivity of the model with respect to this total simulated time and with regards to the mesh resolution.
A summary of this study, that included simulations of the CLEAN geometry and several flap activation configurations, is included in Appendix A.

Results and discussion
The presentation and discussion of the results are split into three different parts. In Section 5.1, an overview of the predicted amplitudes of motion is given. Section 5.2 then presents the wake topology, for selected vibrating and non-vibrating simulations. Finally, Section 5.3 includes additional comments about the characterization of the coupling mechanism.

Amplitudes of motion
Some of the performed simulations led to significant vibrations in the edgewise direction, in the form of limit cycle oscillations. The frequency of these oscillations corresponded to the frequency of the first edgewise mode (i.e. 0.67 Hz). Fig. 8 depicts the peak-to-peak amplitude of motion at the tip, for all the performed computations. For the CLEAN geometry (black curve), a clear VIV amplification region was found between = 45°and = 65°. This region is consistent with the previous work of Horcas et al. (2020) on the baseline geometry of the IEA 10 MW, and previously shown in Fig. 1. Similar amplitudes of vibration were found when employing relatively low positive flap angles ( Fig. 8(a)). On the other hand, the flap angles of 70°, 80°a nd 90°showed a noticeable mitigation of the maximum amplitude. In addition, the onset of the vibration seemed to be shifted towards higher inclination angles. This effect was also observed when actuating the flap with negative angles (Fig. 8(b)), in a relatively symmetrical manner. Consistently with the results for the AVATAR 10 MW wind turbine shown in Horcas et al. (2018), the negative flap actuation was found to be more effective with regards to the decrease of the amplitude of vibrations. This is further illustrated in Fig. 9 by means of the maximum peak-to-peak amplitude of motion (i.e. among all the simulated inclination angles), as a function of the considered flap angle.    Fig. 10 illustrates the different wake topologies found in the present work, through selected configurations at = 50°. The wake topology for the flap angle of −90°, shown in Fig. 10(a), is representative of all the simulations leading to low amplitudes of vibration (≲1 m peak-topeak). It corresponded to the stationary wake described in Horcas et al. (2020), and was dominated by series of stationary vortex tubes, shed from the leading edge and trailing edge in an alternate way. Fig. 10(b) depicts the wake for the CLEAN geometry, that led to a high amplitude of vibration for this inflow (i.e. 2.3 m peak-to-peak). It can be observed that the wake topology was substantially different. The stationary wake was re-configured into the so-called unsteady wake (Horcas et al., 2020), with a generalized vortex shedding for most part of the blade, including the spanwise region where the flap was installed. This shedding led to a harmonic fluid loading, which fundamental frequency corresponded to the frequency of motion. A considerable slant angle was observed in the release of the vortical structures for the inner half of the blade, so that they were released in an oblique manner with regards to the surface. For a better insight on the differences between the stationary wake and the unsteady wake the reader is referred to Appendix B, where animations of the two discussed simulations are included. The wake for the flap actuating at 20°, where an amplitude of motion of 2.4 m was computed, was found to be equivalent to the CLEAN wake (Fig. 10(c)).

Wake visualization
For vibrating configurations with lower amplitudes of motion the wake topology was relatively similar, but the detachment of the vortical structures from the blade surface seemed to start at lower spanwise locations. This is illustrated in Fig. 10(d) for the particular case of a flap actuating at −30°, that predicted a peak-to-peak amplitude of vibration of 1.5 m.

About the coupling mechanism
The study of the wake, shown in Section 5.2, highlighted the complexity of the flow for the targeted application. Indeed, the simulations showing large amplitudes of vibration manifested a predominantly oblique shedding, which took place for most part of the blade. While the re-configuration of the steady wake into this complex shedding pattern should be the focus of future studies, one of the observations concerning the flow around the flap region already allowed to progress in the understanding of the VIV phenomenon. In particular, the analysis of the wakes revealed that for the simulations where the flap actuation had a mitigation effect on the amplitudes of vibration, the detachment of the vortical structures from the blade surface was shifted towards lower span locations. This was not only observed for the particular inclination angle shown in Fig. 10, but also for other values of . To quantity the effects of this spanwise shifting on the coupling mechanism, the integrated average power in the edgewise direction was used, here defined as: where 1 and 0 define the time window for the computation of , and were chosen to correspond to the last five periods of the frequency of motion. The variable denotes each of the 101 aerodynamic sections that were defined along the span of the blade. In this way, corresponds to the sectional loading in the edgewise direction, anḋto the corresponding velocity of motion of the blade axis.
The distribution of can be understood as a measure to identify the regions contributing to the oscillations, as positive gradients highlight the sections that transfer energy from the flow to the structure, and vice versa. Fig. 11 depicts the integrated average power as a function of the vertical coordinate , for selected simulations at = 54°. Note that the value of at the blade tip indicates the total average power transferred, so that it is directly related to the amplitude of the limit cycle oscillation. In this way, the curves shown in blue correspond to the flap angles leading to the largest attenuation of the vibrations, when compared to the CLEAN geometry. All the included configurations showed a region of monotonic increase, starting around mid span. The flap angles that significantly mitigated the motion seemed to shift this region towards lower . However, at the spanwise location where the flap geometry started (i.e. = 61.8 m), all the simulations predicted a similar value of . Indeed, the region steering the differences on the total average power corresponded to the spanwise range where the flap was installed. So that, as the flap actuation became closer to −80°, the distribution of along the flap length became more and more flat. This shows that the effect of the flap installation on the coupling mechanism remained relatively local in spanwise. Additionally, this observation seems to be consistent with the spanwise shifting on the release of the vortical structures attributed to the flap actuation, and previously shown in Fig. 10. In view of that, it seems reasonable to think that a necessary condition for the VIV to arise is the triggering of vortex shedding around the flap region. To further explore the plausibility of this hypothesis, a non-dimensional analysis was carried out. This relied on the establishment of a Strouhal law around the section of the blade located at the flap center. Inspired by the previous work of Horcas et al. (2020), the notion of reduced frequency was used, here defined as:

=
(2) with being the frequency of motion (i.e. 0.67 Hz). is an effective wind speed. Following the independence principle (Jones, 1947), was chosen to be the horizontal projection of the freestream S.G. Horcas et al. velocity ∞ cos( ). stands for an effective chord. was assumed to be the projected length with regards to for a planar cut located at the center of the flap, consistently to previous work on 2-D shapes [e.g. Lind and Jones (2015)]. For reference, an illustration of and is included in Fig. 12(a), and the distribution of the latter variable with regards to the flap angle is depicted in Fig. 12(b). Fig. 13(a) shows the peak-to-peak amplitudes of motion for all the performed simulations, as a function of . When compared to Fig. 8, where the inclination angle was used instead, the results seem to better collapse. Indeed, the maximum amplitudes of vibration are in the range of of [0.12,0.14], irrespective of the considered flap angle. This range is consistent with the Strouhal numbers found for wind turbine airfoils (Skrzypiński et al., 2014;Pellegrino and Meskell, 2013), indicating that the studied instability could keep some similarities with the classical lock-in phenomenon observed for 2-D shapes in transverse oscillation. More in particular, the appearance of the instability could be related to the natural tendency of the sections around the flap region to shed at the structural frequency. Additionally, the fact that a common range was found could indicate that the shifting on the values leading to vibrations is directly related to the modification of . That would imply, in turn, that a trailing edge flap is unable to completely suppress VIV, so that its effects are limited to the modification of the amplitudes of vibration. However, the presented results only show a significant mitigation of the CLEAN vibrations when operating at extremely high flap angles (either positive or negative), which are indeed unrealistic for such aerodynamic device. Additionally, the maximum amplitudes of vibration attributed to the flap activity and the shifting of the VIV amplification region in terms of seemed to be interrelated. In fact, a shifting towards higher seemed to decrease the maximum amplitude of motion, and vice versa. This effect is suspected to be linked to the change on the values of , as shown in Fig. 13(b). In other words, the fact that a flap can have a Fig. 13. Non-dimensional study for all the simulations of the present work, with the colors representing the flap angles: (a) peak-to-peak amplitude of the tip motion in the edgewise direction, as a function of the reduced frequency , and (b) maximum peakto-peak amplitude per flap angle, as a function of the mean effective speed weighted by the amplitudes of vibration̄. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) limited action on the value of could also restrain its impact on the amplitudes of vibration.

Conclusions
The objective of the present work was to assess the potential of the installation of a trailing edge flap for the mitigation of Vortex Induced Vibrations (VIV) of a horizontal axis wind turbine blade in standstill. Fluid Structure Interaction (FSI) simulations of a representative modern blade design were performed, coupling a Computational Fluid Dynamics (CFD) solver with a multi-body finite element code. For a set of inflow conditions, the aeroelastic response of the baseline geometry (referred to as CLEAN) was compared to the results found when installing a flap actuating in the range [−90°,90°]. The flap geometry was introduced in the simulations through the manipulation of the CFD surface mesh, using a free-form deformation approach.
Overall, and consistently to previous research based on a different blade design, the consideration of the flap resulted in a relative decrease of the amplitudes of vibration in the VIV region, when compared to the CLEAN configuration. However, the mitigation was only significant for extremely large flap angles (either positive or negative), falling outside of the expected operating range for this aerodynamic device. The results of the present research shed light on the basic mechanism underlying the vibrations mitigation. In particular, it was found that the installation of a trailing edge flap seemed to shift the release of vortical structures towards lower spanwise locations. This wake modification, in turn, seemed to considerably affect the amount of power that was transferred all along the flap length.
A simple non-dimensional analysis showed that the regions of instability did not disappear, but were simply shifted. This shifting was found to be consistent with the impact that the flap geometry has on the effective chord length, indicating that it has an upper limit. Additionally, there seemed to be a clear link between the shifting of the unstable conditions and the amplitude decrease when installing a flap. In other words, it indicated that the potential of a trailing edge flap with regards to the mitigation of vibrations is limited, at least when operating at a constant angle. Future work could be devoted to explore if the effects of a dynamic flap actuation on the wake could lead to a different scenario with regards to the suppression of the instability. Additionally, it could be relevant to explore other solutions, maybe involving the operational conditions of the machine, such as allowing the blade to freely rotate during the simulation to model the idling configuration. Finally, it should be emphasized that the present study was purely numerical. While several validations of the presented CFD and FSI methods were previously performed, a targeted standstill blade vibrations experiment would be necessary to support the observations made.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The data that support the findings of this study are available upon reasonable request to the corresponding author.

Appendix A. Sensitivity studies
Before performing the numerical campaign presented in this work, a series of preliminary analysis were performed, aiming at assessing the sensitivities of several CFD parameters with regards to the aeroelastic response predicted by the FSI simulations. A summary of the findings for the most relevant parameters, namely the mesh discretization and the total simulated time, are included in this appendix. It should be noted that, while only selected results are shown (both in terms of flap actuation and inflow), similar conclusions could be drawn for other configurations.

A.1. Mesh discretization
Additional grids were created to assess the sensitivity of the results with regards to the mesh discretization. The topologies of these additional meshes were based on the grid resolution described in Section 4.1, here labeled as BASE. In this context, the FINE mesh was generated by doubling the resolution of the BASE grid in all direction. The COARSE1 mesh was generated by halving the resolution of the BASE grid in all direction, and the COARSE2 mesh was generated by repeating the coarsening procedure for the COARSE1 mesh. A summary of the final number of mesh elements for each of the grid variants is presented in Table A.2. For reference, a representative wall clock time reported by the cluster Jess, that was used for performing all the simulations of the present study, is also included.
Jess is a high-performance computing cluster owned by the Technical University of Denmark (DTU). It is equipped with a total of 320 nodes, each of them with 20 cores running at 2.8 GHz. Fig. A.14 depicts the time series of the blade tip position in the edgewise direction, for a selected configuration of the mesh sensitivity study, that accounted for a flap actuation. In particular, the case where the flap operates at −20°is shown, assuming an inclination angle of 54°. The results show a very small decrease in the amplitude when passing from FINE to BASE, with a relative difference of −1.7%.
For the type of analysis included in the present work, this difference was found to be acceptable. Note that passing from BASE to COARSE1 also showed a relatively low amplitude difference (i.e. 3.6%). However, as the mesh sensitivity study could not be performed for every single configuration included in the present work, it was decided to keep the finer BASE resolution. The effect of the mesh coarsening was noticeable when switching to COARSE2. The relative difference for that particular mesh, when compared to COARSE1, was of 8.5%.

A.2. Influence of total simulated time
While a fixed total simulated time of 45.3 s was assumed for all the computations of the present work, selected configurations were re-run for a longer time.    those found at the end of the simulation. However, the peak-to-peak amplitudes of motion for both CLEAN and the flap at −30°seemed to be under estimated when being evaluated at 45.3 s, relatively to the final stabilized values (Table A.3). In other words, the buildup time for these vibrating configurations was found to be larger than the total simulated time chosen. This does introduce a certain degree of uncertainty to the peak-to-peak amplitudes shown in the present work. Nevertheless, it was assumed that the use of a total simulated time of 45.3 s did not compromise the identification of the most unstable situations. Or, equivalently, that there was a correlation between the buildup time and the maximum amplitude at the limit cycle oscillation. For the CLEAN configuration, the animation of the Z component of the flow velocity reveals a generalized release of vortical structures, starting at the outboard part of the blade. Besides the near-root region, the release pattern was found to be synchronized with the frequency of motion. This becomes even more clear when looking at the Y vorticity for a planar cut crossing the center of the flap, that was included for completeness. Indeed, the animation of the vorticity shows a 2 shedding pattern, with a pair of vortices being shed per period. Note that, probably due to the high inclination of the inflow, the shedding was not completely parallel to the surface, but happened in an oblique manner for the inner half of the blade. Near the root, the shedding took place in a de-correlated way with regards to the motion, and a broadband spectrum was computed for the loading at that region.

Appendix B. Wake animations
For the simulation accounting for a flap actuation of −90°, the wake was mainly characterized by series of stationary vortices. These did not lead to a clear alternate release of vorticity, as it can be observed in the sectional cut of vertical vorticity. However, the small motion attributed to this simulation seemed to imply the oscillation of these structures. It should be noted that, in the near-root region, a broadband shedding was also observed, consistently to the results found for the CLEAN configuration.

Appendix C. Supplementary data
Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jweia.2022.105118.