Fluid – structure interaction simulations of a wind turbine rotor in complex flows, validated through field experiments

Aeroelastic simulations of a 2.3 MW wind turbine rotor operating in different complex atmospheric flows are conducted using high fidelity fluid – structure interaction (FSI) simulations. Simpler blade element momentum (BEM) theory based simulations are likewise conducted for comparison, and measurements from field experiments are used for validation of the simulations. Good agreement is seen between simulated and measured forces. It is found that for complex flows, BEM-based simulations predict similar forces as computational fluid dynamics (CFD)-based FSI, however with some distinct discrepancies. Firstly, stall is predicted for a large part of the blade using BEM-based aerodynamics, which are not seen in either FSI simulations or measurements in the case of a high shear. This leads to a more dynamic structural response for BEM-based simulations than for FSI. For a highly yawed and sheared flow case, the BEM-based simulations overpredict outboard forces for a significant part of the rotation. This emphasizes the need of validation of BEM-based simulations through higher fidelity methods, when considering complex flows. Including flexibility in simulations shows only little impact on the considered rotor for both FSI-and BEM-based simulations. In general, the loading of the blades increases slightly, and the rotor wake is almost identical for stiff and flexible FSI simulations. and yaw misalignments. This results in an unsteady loading of the rotor as the blades experience varying wind speeds during their rotation. In term, this leads to an unsteady structural response of the blades which changes the flow conditions once again through the interac-tions between fluid and structure. Aeroelastic phenomena, become increasingly relevant as wind turbines are expanding in size and blades become more flexible, due to mass restrictions. This increases the risk of instabilities from vortex induced vibrations (VIVs) and flutter, which in the worst case scenario can lead to structural failure. Additionally, designers need to consider the increased flexibility to avoid bending blades ever colliding with the tower, and to consider any increase in fatigue loads of the turbine. Less dramatic, but though very important, is also the effects of flexibility on power production, as flexible blades will bend and twist, changing the loads and thereby the driving forces of the turbine. This presents simulations of a 2.3 MW wind turbine rotor situated in simple to complex inflow conditions, using both BEM- and CFD-based aerodynamics. Methods are validated using experimental measurements from the DanAero field experiments. Additionally, the influence the rotor is investigated by to a multi-body finite element solver. cases complexities are studied Case 1: axi-symmetric with no tilt included, Case 2: with tilt Case and including tilt.

In the present study, fluid-structure interaction (FSI) simulations of a 2.3 MW NM80 wind turbine in complex flow are conducted through a coupling of the computational fluid dynamics (CFD) code EllipSys3D 1-3 and the structural solver from the aeroelastic code HAWC2. 4,5 These simulations are compared to pure HAWC2 simulations based on the lower fidelity Blade Element Momentum (BEM) aerodynamics for load calculations.
For atmospheric flow with complexities as yaw misalignment, shear and tilt, blade resolved FSI has not been studied extensively, as it is computationally heavy. Often, low fidelity BEM-based aeroelastic simulations will be conducted instead, as these are efficient and in many cases sufficiently accurate. It is however known, that the accuracy of BEM suffers from the need of airfoil input data and various corrections for 3D effects, stall, dynamic inflow etc. This creates the need of validation of complex flows from experiments and high fidelity methods like FSI.
Yu and Kwon 6 studied FSI of the NREL 5MW reference turbine under shear and yawed conditions separately. In their FSI approach, data between fluid solver and structural solver was only exchanged once per revolution, correcting the BEM loads with periodically converged CFD loads. The study found that especially the torsional degree of freedom is important to consider in FSI simulations, as the changing angle of attack from the torsional deformation alleviates the loading of the blades in the studied setup. This effect is further enhanced when including yaw and shear in the simulation which leads to more unsteady forcing on the rotor. In the phases of the revolution where the blade passes high velocity due to shear or enters the area where the blade rotates against the yawed inflow, the alleviating effect is seen to be highest due to the increased twist angle. Li et al. 7 likewise studied the NREL 5MW reference turbine in sheared and turbulent inflow including flexibility of the rotor. Here the FSI framework was based on a CFD solver coupled with a multi-body dynamics (MBD) structural solver.
Turbulence was imposed using the Mann turbulence box method. 8 The main conclusions of the study was that realistic atmospheric flow including shear and turbulence is important when designing large scale wind turbines. Additionally, the study concluded that inclusion of blade flexibility does not impact highly the wake behavior, whereas inflow turbulence have high impacts on wake diffusion. A comprehensive study was published by Santo et al., 9 where FSI is studied on a rotor in sheared inflow with the rotor structurally represented through shell elements. Various aspects are investigated such as the effect of tower shadow, tilt and yaw. Tower shadow is found to affect very locally around the azimuthal range when the blade passes the tower, and simulations with and without tower yielding very similar results in the remaining part of the revolution. Recently, Guma et al. 10 published a study looking into the aeroelastic response of the NM80 rotor, also studied in the present article, in turbulent inflow. Here, synthetic Mann box turbulence and the Delayed Detached Eddy Simulation (DDES) turbulence model were used to create and resolve turbulent structures in the wind flow. The fluctuating forces occurring on the blades, were used to calculate the fatigue damage on the blades by means of the so-called Damage Equivalent Loading (DEL). It was found that the DEL is mainly influenced by the turbulent inflow rather than the inclusion of flexibility. The study also investigated the impact of considering the entire wind turbine, including tower, compared to a rotor only, by modeling one-third rotor only (one blade with one-third hub), using periodic boundary conditions to resemble the full rotor effect. It was concluded that only small differences are found in average sectional load and deformations, whereas the DEL is more sensitive to the modeling method.
The objective of the present study is to investigate rotor aerodynamics in complex flows for different fidelity solvers along with influences of including flexibility. Cases with measurements available are chosen to validate results obtained from simulations, and highlight possible focus areas for future studies. The study is part of the IEA 29 Wind Task, 11 which has the purpose of comparing and developing simulation methods of various fidelities, and clarifying their limitations. Here, the studied rotor and inflow conditions are defined by the DanAero field experiments. 12 In the experiment the wind is measured at a nearby met mast. Loads based on pressure taps at four spanwise sections along the blade of the NM80 rotor are likewise available for comparison. The study presented in this article, serves as a continuation of the work presented in Grinderslev et al., 13 where the same flow cases were considered for pure CFD computations in another setup, comparing the CFD codes EllipSys3D and Nalu-Wind. 14 Here, it was found that the solvers are able to capture the complex flow to a satisfying degree. It was found, however, that care needs to be taken in terms of grid resolution from inlet to rotor to keep the described shear profile, as an expansion of grid cells will tend to smear out the profile. Additionally, it was found that omitting the constraint effect of the ground leads to a speed up of the flow, especially beneath the rotor, with the consequence of increasing forces on the blades passing the lower area. The experiences gained, are considered in the setups of the current study.
In the following section (Section 2) the methodology of the study is presented, describing the solvers and frameworks used along with the computational setups and flow cases. In the next section (Section 3), the results are presented for the three studied flow cases respectively focusing of resulting forces and flow response. Finally, the study is concluded in Section 4, recapping the results and discussing future applications.

| METHODOLOGY
In this section, the numerical methods and computational setups used in the study will be presented, along with the flow cases chosen for simulations.

| Flow solver
EllipSys3D is used to solve the incompressible Navier-Stokes equations in structured curvilinear coordinates using the finite volume method with a collocated grid arrangement. The code is parallel and highly scalable using Message Passing Interface (MPI) and multi-block decomposition, multi-grid method and grid sequencing. EllipSys3D has multiple convective schemes implemented, such as Central Difference (CDS), Second order Upwind (SUDS) and Quadratic Upstream Interpolation for Convective Kinematics (QUICK). For pressure correction, various algorithms are possible such as PISO, SIMPLE, SIMPLEC and variations hereof. Rhie-Chow interpolation is used to avoid odd/even pressure decoupling. Overset capabilities, including grid hole-cutting are implemented internally in the code. 15 Several turbulence models are implemented such as two equation Reynolds Averaged Navier-Stokes (RANS) models, k − ϵ and k − ω among others, hybrid models like Detached Eddy Simulations (DES), Delayed DES (DDES), Improved DDES (IDDES) and higher fidelity Large Eddy Simulations (LES) models.
Volume grid deformation is handled using a blend factor (0 ≤ F ≤ 1), handling how volume grid points are affected by surface grid deformation. Grid points close to the solid surface will replicate the full deformation of the surface by having a blend factor of 1 (F = 1), while cells further from the surface will move with a fraction of the deformation (F < 1). This is handled through the blend factor that depend on curve length along the grid line normal to the surface, either linearly or through a tanh function. The range of which the blend factor moves from 1 to 0 is tuned for the specific grid setup. When using the overset grid method, the deformation is only transferred to the grid group containing the rotor.
The code has been used thoroughly for many years for various test cases and was validated in, for example, the Mexico project 16,17 and for the Phase VI NREL rotor. 18,19 Recently, the code was validated for the specific case of the present NM80 rotor in atmospheric flow conditions by comparison with the CFD code Nalu-Wind 14 and measurements from the DanAero experiments. 13

| Aeroelastic solver
HAWC2 is an aeroelastic code used for calculating low fidelity BEM aerodynamics and structural responses of wind turbines. The structural part of the code is based on the multi-body formulation, accounting for non-linear effects of large deflections. Each structural component, that is, a blade or the tower, can be represented by a number of bodies assembled by linear Euler or Timoshenko beam elements. Subbodies are connected with constraint equations considering non-linearities. HAWC2 has a built-in aerodynamics module, that calculates aerodynamic forces using Blade Element Momentum (BEM) theory. To do this, airfoil polar input is needed to calculate forces on the blade, which is inserted through aerodynamic sections along the blade. The location of these sections is independent on the number of structural nodes, and will often be distributed with concentrations near root and tip to capture high gradients. From the aerodynamic sections, the loads are interpolated to the beam element nodes, to calculate the structural response. Multiple correction schemes are implemented to improve the BEM aerodynamics such as induction corrections from yaw, tip loss corrections, dynamic stall models, tower shadow effect, and much more; see Madsen et al. 5 HAWC2 is widely used by industry, and has been verified and validated in various studies, 5,20 considering the structural and aerodynamics aspects of the code respectively.

| FSI framework
The two codes are coupled, in a partitioned manner, through the Python framework, HAWC2CFD, originally created by Heinz et al. 21 and further developed by Horcas et al. 22 Through the loose coupling framework, the BEM aerodynamics module of HAWC2 is exchanged with the CFD loading instead.
Using predicted displacements and rotations from HAWC2, the CFD mesh is deformed, and a new flow field is found through EllipSys3D.
Loads based on pressure and viscous friction are integrated along intersections between the CFD surface mesh and planes perpendicular to the blade axis, defined by the aerodynamic sections of the HAWC2 model. Those loads are then communicated to the structural solver, so that the computed aeroelastic response is the result of the CFD loading. All communication between EllipSys3D and HAWC2 happens through the HAWC2CFD Python framework. In the work of Heinz et al., 21 a loosely coupled approach was found to be sufficient for wind energy related cases, due to the high mass ratio between turbine and air and is therefore used. This means, that for each main iteration each solver solves the time step only once, through internal sub-iterations, before advancing in time.
Studies involving the application of HAWC2CFD, for both operational and standstill configurations, include the studies of Heinz et al. 21,23 and Horcas et al. 22,24 The framework has further been validated with experiments through simulations of a pull-release test of a wind turbine blade in the large scale test facility of DTU. 25 The process of the framework between the main iterations can be described through the following steps: • The generalized displacements at the present time step are predicted by HAWC2 with second order accuracy, using kinematics from the previous time step.
• Displacements are send to EllipSys3D and the surface mesh is deformed while displacements are propagated into the volume mesh using a volume blend method.
• The Navier-Stokes equations are solved to calculate the flow field for the new time step through under-relaxed sub-iterations in EllipSys3D.
• Forces are computed and integrated on the CFD mesh surface at prescribed sections and communicated to HAWC2.
• Forces are applied at the aerodynamic sections of the HAWC2 model and actual deformations are calculated.
• Unless the solution has reached total simulation time, the simulation is advanced to the next time step and the procedure is repeated.

| Computational grids
In the present setup, only the wind turbine rotor is considered, and it is assumed that neglecting the tower will have a minor effect. The tower will mainly influence the inboard loading, while the outboard loading will be dominated by the velocity shear. This is supported by the findings of Guma et al. 10 of the same wind turbine.
To investigate rotors in atmospheric flow including shear and yaw, an overset grid approach offers good flexibility, as separate meshes can be constructed for each purpose and then overlapped, connected through interpolation algorithms. In EllipSys3D, flow information (velocity and pressure gradients) are communicated between the grids through donor and receiver cells using trilinear interpolation, and fluxes are corrected to ensure a divergence free flow field within the grid groups; see Zahle et al. 15 Three mesh groups are used and overlapped in this study; see Figure 1. Near the rotor, an O-O type mesh is grown from the blade surface mesh, extruding ≈15 m, discretized over 128 cells, normal to the surface using the inhouse grid generator HypGrid. 26 The first cell adjacent to the rotor surface has a height of 1 Á 10 −6 m to ensure a y + of less than 1. Each blade is represented through 128 grid points spanwise and 256 chordwise. The blade tip and grid around a blade section are presented in Figure 2.
Around the rotor mesh, a cylindrical disc mesh is constructed with pre-cut holes around the blades. This mesh rotates along with the rotor mesh, speeding up the hole-cutting algorithm, as the holes move along with the rotor. To allow for the rotor and disc to tilt, without cells moving outside the ground boundary, the last third of the disc mesh downstream is narrowed linearly.
To simplify the overset search of donor and receiver cells and speed up computations, all deformation from the rotor will be propagated to the rotor mesh in such a way, that only cells that are inside the hole region of the overlapping disc mesh deform. By this, there is no need for updating donor and receiver cells between the rotor and disc mesh, as the two rotate together. This choice necessitates the hole of the disc mesh to be far enough from the surface to leave room for the deformation of mesh cells without impairing the cell quality by introducing excessive skewness or stretching. In the presented grid, the holes are ≈17 m wide along the rotor axis, with the undeformed blades located in the middle of the holes. Deformation of the blades is then propagated to the rotor mesh ensuring cells further than ≈8.5 m from the initial surface are unaltered. To ensure this, volume mesh deformation is kept within 40% grid normal distance with a linearly decaying blend factor starting the decay at 15% grid normal distance. This, to ensure that the quality of the cells in the high gradient region (<15% grid normal distance) remain unstretched by moving the cells as rigid bodies, while the region of the rotor grid containing donor cells (>40% grid normal distance) does not move. Here, the meaning of grid normal distance is defined as the curve length along the grid line that protrudes normal to the surface. For the farfield domain a semi-cylindrical shape is chosen, as depicted in Figure 3, ensuring a good compatibility between overlapping grids. The farfield domain expands 2 km, ≈25D, in the flow direction, with the rotor placed 800 m, ≈10D, away from the inlet. Cells are stretched far from the rotor, F I G U R E 1 Grids used for simulations. (Left) side view, (right) front view. Red cells show receiver cells of overlapping grids. Blue, rotor grid; orange, disc grid; black, background grid. Entire background grid is not shown but compressed closer to resolve the wake. Across the flow direction, the domain radius is 500 m ≈6.25D in each direction, with the rotor placed in the center. Ground conditions are ensured 57.2 m ≈1.43D below the hub, using a symmetry boundary condition, which constrains the flow perpendicular to the boundary. Parallel to the flow, the symmetry condition works as a slip condition. This is chosen over the no-slip boundary condition or an ABL wall function, as no development of the prescribed flow is desired. By the symmetry condition, the input profile is maintained all the way from inlet to the rotor, without any significant change. The combined grids consist of a total of 39.3M grid cells, of which 14.2M cells are used for the rotor grid.

| Simulation setups
From the IEA29 Wind Task, 11 flow cases are defined based on field measurements from the DanAero experiments 12 measured on a met mast located ≈250 m from the considered wind turbine. A corresponding power law shear exponent α has been fitted to measurements to define the sheared inflow profiles. Yaw errors are based on the hub height measurements. Three flow cases with increasing complexity are studied, presented in Table 1 and depicted in Figure 4 along with DanAero measurements. The first case is an idealized axi-symmetric case, where no shear,

| EllipSys3D model
In the present study, unsteady RANS (URANS) simulations are conducted and compared to the DanAero averaged measurements. Sheared inflow is imposed through a power law, Equation (1), at the velocity inlet.
Turbulence is modeled through the k − ω SST model by Menter,27 and flow is assumed fully turbulent, as no transition model is included. The use of URANS limits the study to considering ensemble averaged effects on loads and flow only. If turbulent inflow or massive separation effects were to be studied the more computationally demanding DES models would be required. However, this is not the scope of the present work.
Air is described with density of 1.22 kg m −3 and a dynamic viscosity of 1.769 Á 10 −5 kg m −1 s −1 The convective terms are calculated through the QUICK scheme while an improved version of the the SIMPLEC algorithm 28 is used to couple the velocity and pressure. A time increment corresponding to 0.125 rotation per time step (ts) is used for all simulations corresponding to 1.7 Á 10 −3 s/ts for Case 1 and 1.29 Á 10 −3 s/ts for Cases 2 and 3.

| HAWC2 model
The structural model used for the NM80 turbine was created and validated internally at DTU Wind Energy as part of the original DanAero project. 29 The blade has a prebend into the wind of ≈1.5 m at the tip, and a tilt of 5 is applied for Cases 2 and 3. Each blade is structurally discretized into  32 and tip corrections are applied through the commonly used Prandtl tip correction as described by Glauert. 33 No controller is used, as a constant rotation speed and pitch setting is set for each case, as presented in Table 1. The wind flows for BEM calculations are defined similarly to the CFD input, with the prescribed power law profiles and yaw settings.

| Simulation method
For all three cases, simulations with and without flexibility are conducted to assess its influence on loads and wind flow. The simulation process is divided in phases which are depicted in Figure 5.
F I G U R E 4 Conducted cases. Solid lines represent simulation input, and markers with errorbars represent measurement means ± one standard deviation As a first step (phase 1), pure CFD simulations without coupling to the structural solver, are run to develop the flow until convergence in total thrust and power (>40 revolutions). When passing to the HAWC2CFD framework, HAWC2 is initially run for the same amount of revolutions to ensure compatibility in time between the solvers (phase 2). In the following phase (phase 3), the CFD initialization is restarted to run along with the HAWC2 simulation. In this phase, the calculated BEM loads are linearly blended with the CFD loads over two revolutions in the calculations of deformations. This is done to avoid any large force jumps in the HAWC2 solver, such that no undesired vibrations are introduced to the system.
After the blending, pure CFD loads are used for calculation of structural response and a full two-way coupling is simulated (phase 4). This will run until the total simulated time is reached (5-10 extra revolutions). In terms of computational effort, phase 1 is the most costly part of the simulations due to the many revolutions simulated using CFD. Here, 1 to 2 weeks have been needed on 300 AMD 2.9 GHz processors. The second phase is quickly finished in a manner of minutes on one processor, as only HAWC2 is run. Phases 3 and 4 have similar cost per revolution simulated as phase 1, as the coupling to HAWC2 only yields a marginal extra cost due to the loose coupling scheme, however less revolutions are needed.

| RESULTS
In this section, forces, structural responses and flow behavior for all flow cases with both stiff and flexible blades are presented. Both CFD-and BEM-based aerodynamics are considered and compared. The cases are presented individually, with a step-wise increase of flow complexity between the cases. Azimuth angle is defined positive when moving clockwise looking at the rotor towards the downstream direction, starting at 0 when the studied blade points upwards. The torsional deformation is based around the mid-chord axis of the blade. Sectional pressure forces will be presented by means of their normal and tangential components with respect to the chord line of the blade (F n,c and F t,c ). Note that only the pressure contribution to the sectional CFD loads are shown, to better compare to the pressure based measurements. Additionally, the rotor integral loads, that is, torque and thrust (including the viscous contribution from the CFD), will be provided.

| Case 1-Axi-symmetric flow
For the first case with idealized flow, the forces, structural response and flow response will be axi-symmetric, except for the influence of gravity.
In the following, a brief description of the resulting forces and structural response will be presented to give an impression of the simulation capabilities before moving to the more complex flows.

F I G U R E 5 Process diagram of conducted simulations
F I G U R E 6 Normal and tangential pressure forces along blade, Case 1 (axi-symmetric)

Forces:
Normal and tangential forces along the blade for both CFD-and BEM-based simulations with and without flexibility are presented in Figure 6. As seen, all simulations agree well, with only small discrepancies, also with DanAero measurements. Inboard, stall is occurring in the CFD simulations, which leads to the higher discrepancies here, where 3D effects are also significant. In the figures, only the azimuthal position of 0 is shown, which due to the axi-symmetry is close to identical with the other azimuthal positions with only minor differences due to gravity. Including flexibility in CFD-based simulations, lead to marginal differences in forces, whereas the BEM-based flexible simulation result in a slight increase of normal force outboard. This results in an increase of the total thrust for both CFD (+≈1.0%) and BEM (+≈1.7%) based aerodynamics, as seen in Table 2. For torque however, CFD predicts an increase of ≈0.59%, whereas BEM predicts a decrease of ≈0.33% when including flexibility of the blades. The BEM-based simulations predict ≈5% higher total thrust and ≈7% higher total torque than CFD.

Structural response:
Only a small torsional deformation of the blades is seen radially, as depicted in Figure 7.

| Case 2-Highly sheared flow
For the more complex Case 2, with high shear, tilt and a low yaw error, the blades pass varying wind speeds during their revolutions. Dominantly in this case, is the sheared wind flow, creating large differences in the free stream wind velocity. Figure 8 shows the Q-criterion 34 of the flow illustrating the resolved tip vortices, which move faster at high altitudes due to the sheared inflow.
On top of this, the yaw and tilt result in the relative wind speed increasing around 180 and 90 and and decreasing around 0 and 270 , compared to a non-yawed and non-tilted rotor.

Forces:
For Case 2, larger differences are seen between CFD-and BEM-based simulations, and influences of including flexibility are seen as well. Figure 9 shows the azimuthal development of normal and tangential forces of the 78% blade length section, along with measurements from the DanAero campaign. In general a good agreement is seen between simulations and measurements. However, the BEM-based simulations drop in especially normal forces when the blades pass the upper positions, which is not observed in neither CFD nor measurements. Tangentially, all simulations overpredict the forces compared to measurements. This overprediction is higher using CFD than BEM as the tangential forces also drop at the upper position. One reason for this discrepancy, likely has to do with the integrated pressure forces from measurements being very sensitive to the distribution of the pressure taps, which was quite coarse at the trailing edge.
The discrepancy between forces calculated through BEM-and CFD at 0 and 90 is present along almost the entire blade length as depicted in Figure 10 for 0 azimuth. For the BEM-based aerodynamics, the angle of attack reaches the stall region, which lowers the lift on the blade. This is not the case with CFD nor in measurements. The likely reason for this discrepancy is, that the input airfoil data defines the stall domain too early, despite being corrected for 3D effects. The 3D corrections described by Bak et al., 31 show that the corrections done mainly influence the sections of the inboard part of the blade. The stall seen here, however, occurs on a large part of the blade from ≈30% blade length and outwards.
As CFD is not based on predefined airfoil data, the stall is not predicted for those simulations.
Integration of forces show that including flexibility increases the total thrust and reduces the power producing torque. This is the case for both BEM-and CFD-based simulations, however with CFD predicting higher thrust and torque due to the stall discrepancy. As seen in Table 3, by including flexibility an increase of 2.0% and 0.8% in mean thrust and a decrease of 0.8% and 4.4% in mean torque for CFD-and BEM-based simulations, respectively.
Structural response: For this highly loaded case, a large relative difference is seen in the torsional deformation presented in Figure 11 (left). As seen, the BEMbased simulations twist more than the CFD-based simulations, however still with small angles below 1 . The BEM-based simulations predict approximately twice the torsional deformation at the tip, which further enhances the angle of attack and the tendency to stall.
The discrepancies between the BEM-and CFD-based simulations are also seen in the blade deformation, depicted using the tip displacements in Figure 11 (right). Both simulations consistently predict an overall deformation of the mean tip displacement in the flapwise direction of F I G U R E 8 Instantaneous iso-surfaces of Q-criterion = 0.1 and velocity magnitude contours of Case 2 F I G U R E 9 (left) Normal and tangential pressure forces at ≈78% during one revolution. (right) Rotor with and without displacement. Case 2 with high shear ≈2.4 m. This implies that the prebend is neutralized and the blade is even bend 1 m further along the wind direction, as shown in Figure 9. For BEM-based simulations, where stall is occurring in some part of the azimuth, the tip displacement is seen to fluctuate several times per revolution, whereas the CFD-based simulation results in a 1P motion of the blade. Edgewise deformations are governed by gravity, and match between CFD-and BEM-based simulations. The edgewise deformation is slightly smaller than for Case 1, due to the applied pitch of 4.75 and tilt of 5 changing the stiffness in the gravity direction.

Flow response:
As the flow in Case 2 is mainly dominated by shear, the tip vortex rings are convected faster at higher altitudes than near the ground, which is visible in the vertical section of the flow field in Figure 12. Horizontally, the yaw error of 6.02 pushes the wake deficit to the side, where passing blades will point downstream, which takes approximately eight radii distance to even out downstream as depicted in Figure 12. No significant impact is seen on the flow pattern when including flexibility of the rotor, as presented in Figure 12, showing the velocity deficits at various upstream and downstream locations averaged over six revolutions.

| Case 3-Highly yawed and sheared flow
For the most complex Case 3, with high shear, tilt and high yaw error it is not obvious in advance what will dominate the loading of the rotor. While the blade passes the low altitudes, the high yaw error increases the relative wind speed, compensating for the lower free stream velocity from shear. Figure 13 shows how the vortices are convected away from the rotor following the yawed flow direction.

Forces:
Due to the many contributions to the flow asymmetry by shear, tilt and yaw error, the forces seen by the blade through one revolution are asymmetric as well. This is exemplified in Figure 14 showing the normal-and tangential-to-chord forces along the blade, along with Figure 15 showing the azimuthal variation of a section located at ≈78% blade length.  Both flexible BEM-and CFD-based simulations agree well with measurements in normal forces, whereas larger discrepancies are seen in tangential forces, however, not in the same azimuthal ranges. Inclusion of flexibility leads to slightly larger mean values and amplitudes of forces in both BEM-and CFD-based simulations. Taking the section at 78% blade length as an example, mean and standard deviation of forces over one revolution are found and presented in Table 4. Here, DanAero statistics are based on the mean values presented with standard deviation error bars in Figure 15. As seen, including flexibility results in increasing means and standard deviations for normal forces, while tangential loads increase in standard deviation only.
Statistics over one revolution only tells the story of general levels, but nothing about the actual variation during the revolution. Looking closer at the variation in Figure 15, it is found, that normal forces peak in the range of 240 -280 azimuth depending on simulation method, with CFDbased simulations peaking earlier than BEM-based simulations.
Looking at the normal force distributions along the blade, depicted in Figure 14, at 90 azimuth the largest difference between BEM-and CFD-based FSI are observed. Here, the BEM-based simulations predict ≈25% higher normal forces than the CFD-based simulations and measurements at the 92% blade section. Tangential forces along the blade differ more between the simulations and measurements; see Figure 14. As with normal forces, the BEM-based simulations overpredict the forces at 90 azimuth, and also at 0 azimuth. However, it seems that the CFD-based simulations overpredict the forces at 180 and 270 azimuth, whereas the BEM-based simulations agree more with measurements. Integration of forces show that including flexibility increases the total thrust as was also seen for Case 2; however, torque is slightly increased for CFD while reduced for BEM; see Table 5. In this case, BEM-based simulations predict higher thrust and markedly higher torque due to the aforementioned overpredicting of forces.
Inclusion of flexibility is seen to result in an increase of 1.9% and 2.2% in mean thrust for CFD-and BEM-based simulations respectively, and an increase of 0.34% but a decrease of 4.8% in mean torque. This shows that the present BEM simulations have higher sensitivity to inclusion of flexibility.

Structural response:
In the present highly yawed and sheared Case 3, the forces and motions are not varying smoothly in a 1P fashion as it was in Case 2, which was dominated by shear only. Instead, forces and thereby also structural response seem to have a 2P variation. This is due to the relative  Figure 16. After passing 180 , the relative wind speed increases as the sheared free stream wind speed increases and the motion is against the yawed wind direction. This results in a peak of tip displacement at 240 -280 depending on aerodynamics method, as was also seen with the forces in Figure 15. Through the first phase of the revolution, BEM-based simulations predict an ≈0.25 m larger flapwise tip displacement, caused by the overshooting of forces seen in this azimuthal range depicted in Figure 14. The change in twist is again seen to be small, and predicted higher by BEM than CFD aerodynamics.

Flow response:
As seen on Figure 17, the tip vortex rings convected downstream of the rotor are close together in the area where the blade points downstream compared to the area where the blade points upstream. This is due to the asymmetric velocity deficit behind the rotor resulting in unbalanced convection of the vortex structures. As depicted in the velocity profiles in Figure 17, no significant effect is seen in resulting flow by including the flexibility of the rotor. This was also seen in the study by Li et al. 7 on the NREL 5MW turbine.
T A B L E 4 Mean (μ) and standard deviation (σ) of forces over one revolution for section at 78% blade length, presented in Figure 15 μðF  Three cases with increasing complexities are studied being firstly, Case 1: axi-symmetric flow with no tilt included, secondly, Case 2: highly sheared flow with tilt and a small yaw error and finally Case 3: highly yawed and sheared flow including tilt.
For Case 1, both simulation methods compare well with each other and measurements. A ≈5% difference is found in total thrust and ≈7% difference is found in rotor torque. Discrepancies are mainly seen inboards and near the tip of the blades. Structurally, BEM predicts ≈5% larger tip displacement than CFD-based simulations.
For the complex Cases 2 and 3 some distinct discrepancies are however found. For highly sheared inflow in Case 2, with the largest free Despite BEM-based aerodynamics being in good agreement with CFD and measurements in the simple case, the study emphasizes the need of validation for BEM-based aerodynamics, when simulating complex flow cases. In the highly sheared case, a better prediction would be likely found if airfoil data was corrected better for 3D effects on the stall onset, possibly by use of CFD simulations. It is more difficult to assess the exact source of discrepancy of the highly yawed Case 3.
Inclusion of flexibility is found to slightly increase the forces seen by the blade, both for CFD-and BEM-based simulations. This is likely due to the prebend being neutralized, increasing the swept area of the rotor. The torsional stiffness of the rotor is high and deformation is mainly seen in translatoric displacements and not so much in torsion. However, the slight torsional deformation seen increases the angle of attack which also explains some of the load increase seen when including flexibility.
Investigating the flow resulting from the CFD-based simulations likewise show very little influence by the inclusion of flexibility in the simulations as no significant differences are seen in the wake deficits. One limitation of the present study, is the quite high stiffness of the considered rotor, reducing its quality as representative for new and future wind turbines with blades spans of 2-3 times the considered rotor. This leads to a need of conducting similar simulations of larger turbines, to do general conclusions on the influence of flexibility. Additionally, including turbulent inflow would be beneficial, to assess the impact on fluctuating flow which on large blades might be close to the frequencies of the first blade modes.

ACKNOWLEDGEMENTS
This study was funded by the Technical University of Denmark -Department of Wind Energy, through the internally funded PhD project "Fluid-structure Interaction for Wind Turbines in Atmospheric Flow." The DanAero projects, from which experimental data were obtained, were funded partly by the Danish Energy Authorities (EFP2007. Journal nr.: 33033-0074 and EUDP 2009-II. Journal nr. 64009-0258) and partly by eigenfunding from the project partners. Finally, the authors would like to thank the EUDP for funding the project "Participation in IEA Task

PEER REVIEW
The peer review history for this article is available at https://publons.com/publon/10.1002/we.2639.

DATA AVAILABILITY STATEMENT
The structural model of the considered wind turbine is not publicly available, and the structural and fluid solvers are licensed. The flow data that support the findings of this study are available on request from the corresponding author. In order to investigate the sensitivity of the simulations to grid resolution, the grid sequencing of EllipSys3D was utilized. By this it is possible to compare two grid levels, one with blocks of 16 × 16 × 16 cells (coarse grid) and one with blocks of 32 × 32 × 32 cells (fine grid) on the same setup and geometry. Here, the fine setup is the same as described in Section 2.2, which was used for studies presented. The coarse grid setup consists of 4.9M cells, while the fine setup consists of 39.3M cells.
In the present sensitivity study, mean torque and thrust is compared between grid levels for all three flow cases and presented in Table A1.
Relative difference is calculated with fine grid as benchmark, as this is the grid designed for overlap matching.
As seen, for cases 1 and 2, both thrust and torque agree well, in both grid levels, with a maximum error of 2.08% in torque for Case 1. For Case 3, a larger discrepancy is seen with -10.8% error for rotor torque. Along with this, large fluctuations were seen in the torque and thrust time signals for the coarse grid setup, which were not present for the fine setup, neither in the former cases regardless of grid level. This is likely due to the overset connectivity, which is not as good in this coarse case, where the large yaw is present between rotor/disc grids and background domain grid. Considering computational resources, the much improved overset connectivity of the fine grid, and the thrust and torque signals becoming stable without big changes to mean values, the fine grid is considered sufficient for use in all cases. Note: Mean (μ) and standard deviation (σ) over five revolutions. Relative difference based on means with fine grid as baseline