Articles | Volume 6, issue 3
Wind Energ. Sci., 6, 627–643, 2021
Wind Energ. Sci., 6, 627–643, 2021

Research article 06 May 2021

Research article | 06 May 2021

Wind turbines in atmospheric flow: fluid–structure interaction simulations with hybrid turbulence modeling

Wind turbines in atmospheric flow: fluid–structure interaction simulations with hybrid turbulence modeling
Christian Grinderslev, Niels Nørmark Sørensen, Sergio González Horcas, Niels Troldborg, and Frederik Zahle Christian Grinderslev et al.
  • Department of Wind Energy, Technical University of Denmark, Risø Campus, 4000, Roskilde, Denmark

Correspondence: Christian Grinderslev (


In order to design future large wind turbines, knowledge is needed about the impact of aero-elasticity on the rotor loads and performance and about the physics of the atmospheric flow surrounding the turbines. The objective of the present work is to study both effects by means of high-fidelity rotor-resolved numerical simulations. In particular, unsteady computational fluid dynamics (CFD) simulations of a 2.3 MW wind turbine are conducted, this rotor being the largest design with relevant experimental data available to the authors. Turbulence is modeled with two different approaches. On one hand, a model using the well-established technique of improved delayed detached eddy simulation (IDDES) is employed. An additional set of simulations relies on a novel hybrid turbulence model, developed within the framework of the present work. It consists of a blend of a large-eddy simulation (LES) model by Deardorff for atmospheric flow and an IDDES model for the separated flow near the rotor geometry.

In the same way, the assessment of the influence of the blade flexibility is performed by comparing two different sets of computations. The first group accounts for a structural multi-body dynamics (MBD) model of the blades. The MBD solver was coupled to the CFD solver during run time with a staggered fluid–structure interaction (FSI) scheme. The second set of simulations uses the original rotor geometry, without accounting for any structural deflection. The results of the present work show no significant difference between the IDDES and the hybrid turbulence model. In a similar manner, and due to the fact that the considered rotor was relatively stiff, the loading variation introduced by the blade flexibility was found to be negligible when compared to the influence of inflow turbulence. The simulation method validated here is considered highly relevant for future turbine designs, where the impact of blade elasticity will be significant and the detailed structure of the atmospheric inflow will be important.

1 Introduction

As future wind turbines will have unprecedentedly long and flexible blades, the necessity of understanding the effects of aero-elasticity on the rotor performance and on its structural integrity increases. Along with this, large wind turbines interact with a larger part of the atmospheric boundary layer (ABL), often exceeding the height of the atmospheric surface layer. This also needs consideration in the design phase, as the rotor blades consequently experience a large variation of flow through each revolution, and flow cases which were not relevant to consider for past designs might occur. This needed knowledge can be obtained through high-fidelity methods, such as fluid–structure interaction (FSI) simulations, which model the coupled effects of both flow and structure. These simulations can further be used to develop and improve lower-fidelity engineering models used by wind turbine designers in industry.

FSI of wind turbines in atmospheric turbulent flow is not a widely studied topic, due to the computational costs of such simulations, especially when geometrically resolved wind turbines are modeled. Instead, a more efficient approach often chosen is the use of actuator lines/discs (Sørensen and Shen2002; Sørensen et al.2015). Here, the rotor is represented through body forces smeared in the computational fluid dynamics (CFD) grid, reducing the need of grid refinement significantly. An example of actuator-line-based FSI in turbulent flow is the work by Lee et al. (2013), where simulations of two aligned 5 MW wind turbines in a turbulent inflow modeled by large-eddy simulations (LES) were conducted. The structural response of the turbines was found through the FAST aero-elastic code (Jonkman and Buhl2005) to study the fatigue loading. It was found that especially the surface roughness and the rotor shadow effect had a large influence on the fatigue loading. As actuator lines merely represent the turbines through smeared forces, blade surface boundary layers and resultingly generated wake turbulence are not modeled. Likewise, the resulting shedding of vortices at the tips and roots is not highly resolved and improperly modeled. The far wake response is, however, sufficiently accurate when the inflow to the turbine has a high turbulence intensity (TI) (Troldborg et al.2015).

Looking at rotor-resolved CFD/FSI, using LES is still too computationally expensive for many practical applications. Instead, compromises are needed for the turbulence models. In the works by Santo et al., FSI for wind turbines, structurally represented through finite-element shells, were studied for steady ABL flows (Santo et al.2020a, b) using unsteady Reynolds-averaged Navier–Stokes (URANS), with the kϵ model. In Santo et al. (2020a), the effects of wind shear, yaw error, tilt, and tower shadow were all investigated, and it was found for instance that the introduction of yaw led to a decrease in blade deflection but a large increase in yaw moment on the hub. In Santo et al. (2020b), wind gusts were introduced by acceleration of the flow near the rotor top position. One conclusion found was that, for the setup and turbulence model used, a flow separation occurred when the velocity rapidly increased due to the gust, limiting the load increase avoiding any extreme deflections. To consider turbulent fluctuating flow, a popular alternative to LES is synthetic turbulence generators such as the method by Mann (1998). These methods efficiently create boxes of turbulent fluctuations, which can be used to create inflow for CFD simulations or inserted internally in the domain by additional body forces (Troldborg et al.2014). Along with this, a hybrid turbulence model like detached eddy simulation (DES) can be used to resolve the turbulence in the grid. This model combines the URANS approach for attached flow regions with LES in the separated regions. The use of synthetic turbulence is efficient as the modeling of turbulent fluctuations is fast, and the DES models need less grid resolution near the rotor than LES. The turbulence will, however, not be in balance with the CFD simulation shear as shown in Troldborg et al. (2014), and therefore the turbulence will change as it convects through the domain. Another drawback of this method is that the modeled turbulence is neutrally stratified and therefore cannot naturally handle atmospheric stability. Further, a potential problem of the synthetic turbulence methods is the assumption of homogeneous and Gaussian turbulence. Even though previous work (Berg et al.2016) has shown that the latter assumption does not significantly affect the loads on a wind turbine under normal conditions, one could easily come across cases where these assumptions do not hold. In Li et al. (2015) synthetic turbulence was used to study the geometrically resolved NREL 5 MW reference turbine in sheared and turbulent inflow including flexibility of the rotor. The FSI framework was based on a CFD solver coupled with a multi-body dynamics (MBD) structural solver, and turbulence was imposed at the inlet using the Mann turbulence box as input. The main conclusions of the study were that realistic atmospheric flow including shear and turbulence is important when designing large-scale wind turbines in terms of loading. Additionally, the study concluded that, for the specific turbine and flow cases, inclusion of blade flexibility does not highly impact the wake behavior, whereas inflow turbulence has high impacts on wake diffusion. Guma et al. (2021) recently published a study looking into the aero-elastic 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 is found that for low inflow velocities the DEL is mainly influenced by the turbulent inflow rather than the inclusion of flexibility, at least for the considered relatively stiff rotor.

An alternative more complex method for simulating geometrically resolved turbines in the ABL flow was proposed by Vijayakumar (2015). Here, a hybrid turbulence model was developed which combines spectral ABL LES simulations by Moeng (1984) with more feasible URANS-based kω SAS (Egorov et al.2010) simulations close to the rotor. By this combination, a large decrease in the necessary number of grid cells is achieved, as the URANS-based turbulence models the effect of all the turbulent scales. The model was applied to a single wind turbine blade in Vijayakumar et al. (2016), however, using pure CFD without a structural coupling.

In general, considering presently available high-performance computing capabilities, compromises are needed when doing high-fidelity aero-elastic modeling of wind turbines in atmospheric flow using FSI – either by reducing the rotor representation by actuator lines to allow LES simulations or instead by simplifying the turbulence modeling.

The objective of the present study is to move one step up the ladder of complexity by investigating rotor aerodynamics and aero-elasticity in turbulent LES inflow, using a novel turbulence model. The model is inspired by the one of Vijayakumar et al. (2016), combining the ABL turbulent flow modeling of the Deardorff LES model with the improved delayed detached eddy simulation (IDDES) engineering model near the rotor. This study considers blade resolved FSI simulations of the 2.3 MW NM80 wind turbine rotor including blade flexibility using a FSI coupling framework combining the CFD code EllipSys3D (Michelsen1992, 1994; Sørensen1995) and the structural solver from the aero-elastic code HAWC2 (Larsen and Hansen2007). For the specific rotor, measurements of inflow and blade loading are available for comparison with the computed results. The study is a continuation of Grinderslev et al. (2021), where FSI of the NM80 rotor was studied in various flow scenarios with sheared and yawed, albeit laminar, inflow using URANS turbulence modeling. The flow scenario studied here resembles one of these studied laminar scenarios, albeit with turbulent inflow through the novel hybrid turbulence model and an adjusted grid setup.

2 Methodology

In this section, the computational solvers are presented along with the simulation strategies such as FSI framework and precursor simulations. Further, the participating turbulence models will be introduced, to prepare for the discussion of the hybrid model. Finally, the computational grids used in the study are described along with the chosen simulation parameters.

2.1 Numerical methods

2.1.1 Flow solver

To solve the fluid flow, the DTU in-house CFD code EllipSys3D (Michelsen1992, 1994; Sørensen1995) is used. The code solves 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 the message passing interface (MPI) and multi-block decomposition, the 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 the solution of the pressure correction equation, various algorithms are implemented such as PISO, SIMPLE, SIMPLEC, and variations thereof. The Rhie–Chow interpolation is used to avoid odd–even pressure decoupling. Overset capabilities, including grid hole cutting, are implemented internally in the code (Zahle et al.2009).

Several turbulence models are implemented such as two-equation Reynolds-averaged Navier–Stokes (RANS) models, kϵ and kω among others; hybrid models like DES, DDES, and IDDES; and multiple LES models. In addition to these, a hybrid version of the LES and the IDDES model will be presented in this paper.

For FSI simulations, the deformation of grids is handled through a moving mesh method with a blend factor. The surface displacement is propagated along the grid lines normal to the surface, with a blend factor gradually diminishing to zero with increasing distance to the blade surface. The blending can be either linear in the distance to the blade surface or based on a tanh function. This ensures that mesh points in the vicinity of the blade surface are displaced as a solid body movement along with the blade, while points further away only move a fraction of the displacement. When using the overset grid method, the deformation is only transferred to the volume grid blocks containing the solid surface.

The code has been used extensively for a range of test cases and was validated in, e.g., the Mexico project (Bechmann et al.2011; Sørensen et al.2016) and for the Phase VI NREL rotor (Sørensen and Schreck2014; Sørensen et al.2002). Recently, the code was validated in Grinderslev et al. (2020b) for the specific case of the present NM80 rotor in atmospheric laminar flow conditions by comparison with the CFD code Nalu-Wind (Sprague et al.2019) and measurements from the DanAero experiments. Further, the FSI framework was used in Grinderslev et al. (2021) to simulate the coupled effects of the DanAero-inspired laminar wind flow and the structural response.

2.1.2 Aero-elastic solver

HAWC2 (Larsen and Hansen2007) is an aero-elastic code developed at DTU used for calculating blade element momentum (BEM) aerodynamics and structural responses of wind turbines. The structural part of the code is based on the multi-body dynamics formulation, accounting for non-linear effects of large deflections. Each structural component – i.e., a blade or the tower – can be represented by a number of bodies assembled by linear Euler or Timoshenko beam elements. Sub-bodies are connected with constraint equations considering non-linearities.

HAWC2 has a built-in aerodynamics module that calculates aerodynamic forces using BEM theory. As is common in BEM implementations, prediction of airfoil aerodynamic performance is based on pre-computed look-up tables of lift, drag, and moment, which are needed to calculate forces along the blade. Multiple correction schemes are implemented to improve the BEM aerodynamics, such as tip loss corrections, dynamic stall models, tower shadow effect, and many more; see Madsen et al. (2020).

HAWC2 is widely used by industry, and it has been verified and validated in Pavese et al. (2015) and Madsen et al. (2020), considering the structural and aerodynamics aspects of the code, respectively.

2.1.3 FSI framework

The codes EllipSys3D and HAWC2 are coupled, in a partitioned manner, through the Python framework, referred to as the DTU coupling, originally created by Heinz et al. (2016a) and further developed by Horcas et al. (2019) and García Ramos et al. (2021). Through the use of the coupling framework, the BEM aerodynamics module of HAWC2 is replaced by an interface to the EllipSys3D CFD code.

Using predicted displacements of nodes from HAWC2, the CFD mesh is deformed, and a new flow field is found through EllipSys3D. The loads predicted by the CFD solver are then applied to the HAWC2 structural model and a new deformation is found. All communication between EllipSys3D and HAWC2 happens through the DTU coupling framework. In Heinz et al. (2016a), a loosely coupled approach was found to be sufficient for wind-energy-related cases, due to the high mass ratio between the turbine structure and air, and is therefore used.

Studies involving the application of the FSI framework, for both operational and standstill configurations, include Heinz et al. (2016a, b) and Horcas et al. (2019, 2020). The framework has been validated with experiments through simulations of a pull–release test of a wind turbine blade in the large-scale test facility of DTU; see Grinderslev et al. (2020a). The process of the framework between the main iterations can be described through the following steps:

  • the displacements of the present time step are predicted by HAWC2 with second-order accuracy, using kinematics from the previous time step;

  • displacements are sent 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 and sent to HAWC2;

  • forces are interpolated to the aerodynamic sections of the HAWC2 model, and actual deformations are calculated;

  • unless the solution has reached the total simulation time, the simulation is advanced to the next time step and the procedure is repeated.

2.2 Turbulence modeling

A hybrid turbulence model has been developed to consider the dominant turbulence scales from the atmospheric boundary layer (ABL) down to the blade boundary layer (BBL), within the same simulation. To do this, the Deardorff one-equation LES turbulence model for ABL flows (Deardorff1972) is blended with the IDDES turbulence model (Shur et al.2008), which itself is a blend between URANS modeling in the BBL and LES in the separation region outside the BBL. The blending of the two models happens through the energy equation, which is solved for in both methods. In the Deardorff model, the transport equation of sub-grid-scale (SGS) energy e is solved, whereas the transport equation for total turbulent kinetic energy k is solved in the IDDES method. These energy expressions are blended through their respective terms of diffusion, convection, production, buoyancy, and dissipation using a smooth tanh blending function. By this, e of the Deardorff model coming towards the rotor is transformed into equivalent k of the IDDES, and vice versa in the wake region. In the following, the two models will be introduced, followed by a description of the blending for the hybrid model used in this study.

2.2.1 Deardorff large-eddy simulation model

In the Deardorff LES turbulence model (Deardorff1972, 1980), the turbulent eddy viscosity μt is calculated through the expression

(1) μ t = C k ρ l LES e .

Here, Ck is a constant of 0.1; ρ is the air density; and lLES is a mixing length scale, which for neutral stratification is set equal to ΔLES, which is the grid size, here defined as ΔLES=(dxdydz)1/3, with dx, dy, and dz being the grid spacing in the respective directions.

The SGS energy, e, is found by solving the following transport equation:

(2) D ρ e D t = - τ i j S i j + g θ 0 τ θ w , LES - C ϵ ρ e 3 / 2 l LES + x j μ + 2 μ t e x j ,

where g, t, and μ refer to the gravity, time, and molecular viscosity, respectively. Cϵ is equal to 0.93, and the buoyancy SGS fluxes τθi,LES=-μθθxi, with variable temperature θ and the eddy heat diffusivity being μθ=1+2lLESΔLESμt. θ0 is the surface reference temperature. The SGS stress tensor τij is defined as τij=-2μtSij using the strain rate Sij=1/2uixj+ujxi, with u being the velocity vector.

2.2.2 SST-based detached eddy simulation models

For the kω SST-based DES turbulence models, μt is found through the standard kω SST (Menter1993) approach, which is then altered in the dissipation term depending on the chosen DES model.

(3) μ t = ρ a 1 k max a 1 ω , Ω F 2 , with a 1 = 0.31

Here, k is a total turbulent kinetic energy, ω the specific dissipation rate, Ω the shear-strain rate, and F2 a limiting blending function. k and ω are found through the following two transport equations:

For k,

(4) D ρ k D t = - τ i j S i j + g θ 0 τ θ w , DES - ρ k 3 / 2 l ̃ + x j μ + μ t σ k k x j .

For ω,

(5) D ρ ω D t = γ ν t τ i j u i x j - ρ β ω 2 + x j μ + μ t σ ω ω x j + 2 1 - F 1 ρ σ ω 2 1 ω k x i ω x i .

Here F1 is a blending function, to shift between the standard kω model (F1=1) near the surface and the kϵ model (F1=0) within the boundary layer and further out, while σk, σω, β, and γ are parameters which themselves depend on F1. Finally β* and σω2 are constants. νt is the kinematic turbulent viscosity νt=μt/ρ. All constants and parameters as well as blending functions F1 and F2 can be found in the original work by Menter (Menter1993). In the original work by Menter, the buoyancy term is not considered. In the EllipSys3D version of the kω SST, however, it is considered as the second term of the right-hand side of Eq. (4), and the flux is found as τθi,DES=-μtθxi.

The length scale which appears in the k-equation serves to switch from URANS to LES mode and is defined as:

(6) l ̃ = min l k - ω , l DES ,

where lk-ω=k/(β*ω) and lDES is the length scale in the LES region. In the standard DES model (Spalart et al.1997; Travin et al.2004), lDES=CDESΔDES, where ΔDES=max(dx,dy,dz) and CDES is a F1-dependent parameter.

DES is known to be sensitive to sudden changes of grid refinements as grid-induced separation (GIS) can be introduced. Here, the modeled turbulent viscosity will drop instantly without the additional turbulence being resolved. It is also known to have a mismatch between the URANS and LES region, if used as a wall-modeled LES model. These issues are addressed in DDES (Menter et al.2003; Spalart et al.2006), IDDES (Shur et al.2008), and SIDDES (Gritskevich et al.2012) by using more advanced expressions for the length scale l̃.

2.2.3 Hybrid ABL–BBL model

In order to simulate the effect of turbulence on both ABL and BBL scales, a hybrid method is suggested where the Deardorff ABL LES model is blended together with the BBL DES models to avoid the need of excessive grid resolution in the BBL otherwise needed by LES. The blending is established through a blending function Fh, which is zero in the ABL region and one in the DES region, and then defining a hybrid turbulence kinetic energy k̃=Fhk+(1-Fh)e. Using these definitions, the energy equations Eqs. (2) and (4) are combined to give the following transport equation for k̃:

(7) D ρ k ̃ D t = - τ i j S i j + g θ 0 F h τ θ w , DES + 1 - F h τ θ w , LES - ρ k ̃ 3 / 2 F h l ̃ - C ϵ 1 - F h l LES + x j μ + μ t σ k F h + 2 1 - F h k ̃ x j .

The blending function Fh is defined as follows:

(8) F h = 0.5 - 0.5 tanh d w - R 2 δ blend .

Here, dw is the wall-normal distance; R is the wall distance to the location where Fh=0.5; and δblend is the transition distance between Fh=0.12 and Fh=0.88, where the blend is most rapid.

To allow the present method to work together with the kω model, an expression for ω is needed in the LES region. This expression is made through the standard kω turbulent viscosity expression, to ensure consistency through the blending regimes.

(9) μ t = C k ρ l LES k ̃ = ρ k ̃ ω LES ω LES = k ̃ C k l LES

Here, it is assumed that the blending from DES to LES happens in the region where F2=0, such that the viscosity limiter is inactive. A blended expression ω̃ is then found for the entire domain.

(10) ω ̃ = ω F h + ω LES 1 - F h

This allows the calculation of the turbulent viscosity similar to Eq. (3):

(11) μ t = ρ a 1 k ̃ max a 1 ω ̃ , F 2 Ω .

It is noted that in the Deardorff part of the model the turbulent viscosity, μt, is linearly proportional to the length scale, lLES, through ωLES; see Eq. (9). This needs to be considered if sudden changes are made to the grid resolution, as this will lead to a proportionally equal change to μt. This could for instance be the case with overset grids, as used in the present study, where sudden changes of grid resolutions are happening over the interface. This should in theory be fine, as the resolved turbulence adapts to the grid, but, as seen with the known GIS issues of the original DES model, the change in resolved turbulence, due to grid refinement, does not happen instantaneously. In the present study, this is handled by limiting the LES length scale ΔLES to the grid size of the background grid; see Sect. 2.3.1.

2.2.4 Turbulent inflow simulations

In this study, the turbulent flow of the atmospheric boundary layer is modeled through a LES precursor simulation using the Deardorff model. Here, a neutrally stratified wind profile is simulated and sampled for use as input in the successor simulation including the rotor.

In the successor simulation the hybrid LES–IDDES model is used. LES is used for turbulence modeling in the majority of the domain, except for the region close to the rotor. In this area, the IDDES model is utilized instead, which removes the LES grid requirement near the rotor surface.

The precursor conditions are approximating measurements from the DanAero field experiment (Bak et al.2010), where a met mast located ≈2.5 diameters from the considered rotor measured the wind field using cup anemometers at six points vertically, 17, 28.5, 41, 57 (hub height), 77, and 93 m. The data set from these cup anemometers is used to fit a corresponding neutral log-law wind profile to generate inputs for the Schumann–Grötzbach (SG) wall model (Schumann1975; Grötzbach1987) used in the simulation.

(12) U = u * / κ ln z / z 0 ,

where U is the wind speed, u* the friction velocity, κ the von Kármán constant (≈0.4), z the vertical coordinate, and finally z0 the roughness length. As a neutral stratification flow is modeled for simplicity, no temperature is modeled in the present study.

2.3 Simulation setups

2.3.1 EllipSys3D model

Air is described with a 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 a blend of the fourth-order central difference (CDS4) scheme in the LES area and the upwind QUICK scheme in the URANS part as described by Strelets (2001). An improved version of the SIMPLEC algorithm (Shen et al.2003) is used to couple the velocity and pressure. No transition model is applied, such that the blade boundary layer is assumed fully turbulent. A time increment corresponding to 0.125 rotation per time step is used for all simulations, corresponding to 1.29×10-3 s per time step. The rotation speed of the rotor is constant at 16.2 rpm, resulting in an effective Reynolds number of ≈6 million along the majority of the blade for the studied flow case.

Turbulence blending

To enable the hybrid turbulence modeling, a blending region must be defined. As mentioned in Sect. 2.2.3, a sudden grid refinement will create a sudden length scale change and thereby, if in the Deardorff LES region, a sudden change of turbulent viscosity. In the present setup with overset grids, it is therefore chosen to avoid the viscosity “jump” by keeping the LES length scale ΔLES at the background grid value. By this, the refinement does not change the dissipation length scale or the viscosity. Near the rotor, however, an IDDES region is prescribed depending on the wall distance. In this region, the refined mesh impacts the turbulent dissipation through ΔDES, as usual.

The length scale limit of the LES region caps the frequency range of the resolved turbulent kinetic energy to the background grid resolution. Close to the rotor, however, the small-scale detached flow is still captured through the IDDES model. For studies, with long distances from refinement to object or larger resolution differences between background and overlapping sub-grids, this strategy would likely not be optimal due to the capping of resolved frequencies being based on an unnecessarily large grid size.

In the present setup the blending between LES and IDDES happens with the midpoint (Fh=0.5) 8 m from the surface, with the majority of the blend happening over a distance of ±2 m from this point to ensure a smooth transfer from LES to IDDES; see Fig. 1.

Figure 1Blend factor, Fh. Red: IDDES region. Blue: Deardorff LES region. Isosurface of Fh=0.9 (≈6 m normal from surface).


Figure 2Concept of precursor to successor simulations along with domain sizes of conducted precursor and successor simulations. Colored contours showing flow velocity W.


2.3.2 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 (Bak et al.2013). This model was also used in a former study of FSI on the same rotor (Grinderslev et al.2021) comparing URANS FSI with BEM-based aero-elastic simulations for complex laminar flow scenarios. The blade has a prebend into the wind of ≈1.5 m at the tip. Each blade is structurally discretized into 22 bodies, each consisting of one Timoshenko beam element. As only the rotor is modeled in CFD, only blade flexibility is considered as well. This means that tower, nacelle, shaft, and hub are not active parts of the HAWC2 simulations. A total of 60 aerodynamic sections are distributed per blade, which are used for both BEM and CFD loads. For BEM aerodynamics, airfoil data are used, obtained during the original DanAero project through wind tunnel tests and corrected for 3D effects; see Bak et al. (2006, 2011). From Grinderslev et al. (2021), it is known that the airfoil data do not capture well the 3D effects and predict an earlier stall than seen in CFD or experiments. In this case, however, the BEM calculations are used for initialization only to get good estimate of the initial bending, and for that reason no further corrections to the airfoil data have been conducted. Dynamic stall corrections (Hansen et al.2004) and tip corrections (Glauert1935) are applied during the initializing BEM calculation. No controller is used, as a constant rotation speed of 16.2 rpm and pitch setting of −4.75 (increasing the angle of attack) are set. For simplicity, the yaw and tilt are omitted in the simulation setup. For the DanAero campaign used for comparison, a tilt of 5 and average yaw error of 6.01 were present, however.

2.4 Simulation method

2.4.1 Precursor simulation

For the precursor simulation, as the first step, the turbulent flow is developed by recycling the flow using periodic boundary conditions. This resembles the flow moving over a very long distance, building up the boundary layer, and producing the turbulence through shear production. In order to ensure a mean profile close to the desired measured wind velocity profile, the SG wall model is used. This forces the surface shear stress of the first adjacent cells to the ground to fit the log law. The flow is driven trough a constant pressure gradient calibrated to obtain the desired friction velocity and resulting velocity profile with a roughness length z0 of 0.73 m.

Initially, the grid sequencing scheme of EllipSys3D is utilized on three grid levels to speed up the simulation and reach a fully turbulent domain quickly. When the flow is fully turbulent and the mean flow profiles match the desired flow, planes consisting of velocity components U (horizontal perpendicular to mean flow direction), V (vertical and perpendicular to mean flow direction), and W (mean flow direction); pressure, P; and SGS kinetic energy, e, are sampled. The plane is centered in the cross-flow directions of 1000 m × 600 m with 4 m cell distances; see Fig. 2.

2.4.2 FSI simulation

The FSI successor simulation process is divided into phases depicted in Fig. 3.

Figure 3Process diagram of conducted simulations.


In the first phase, simulations without coupling to the structural solver are run to develop the flow and fill the domain with the sampled turbulent flow. In this phase, the grid sequencing scheme of EllipSys3D is used exploring coarser grids to minimize the simulation cost during spin-up.

Figure 4Grids used for simulations. (a) Side view, (b) 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.


When passing to the FSI framework, phase 2, HAWC2 is run for the same amount of revolutions using BEM aerodynamics corresponding to the mean flow profile to ensure compatibility in time between the solvers when coupling and obtaining a good guess of initial blade deformations. In phase 3 the coupling of EllipSys3D and HAWC2 is initiated with a smooth linear blending of forces over two revolutions to switch from BEM to CFD loading. This is done to avoid any large force jumps in the HAWC2 solver, to suppress undesired vibrations in the system. In the final phase, phase 4, pure CFD loads are used for calculation of structural response, and a full two-way coupling is simulated for the desired amount of revolutions.

2.5 Computational grids

2.5.1 Precursor simulation

The precursor domain is 4000×1300×8000 m (width × height × length) discretized 576×256×1920 cells divided into 8640 blocks of 323 cells. A total of ≈283 million cells are present in the precursor. The grid cells vary in size in the cross-flow directions to obtain higher resolution in the sampling area. In the sampling area, the cells are cubic with 4 m cell sides, while cells are slowly stretched towards the boundary sides and top. Periodic boundaries are prescribed on the vertical sides, while a symmetry condition is used on the top boundary, and the SG wall model is used for ensuring the Monin–Obukhov similarity law in the first cells adjacent to the wall by dictating the wall shear stress.

2.5.2 Successor simulation

For the rotor simulations, an overset grid method is utilized (Zahle et al.2009), as this allows for a stationary background grid including the ground, while a rotating grid can be used for the rotor. Flow information is then communicated by interpolation between the grids through donor and receiver cells within the overlapping region of the meshes.

In the present setup, only the rotor is considered, omitting the tower, hub, and nacelle, with a total of three overlapping mesh groups; see Fig. 4. The omission of these elements is expected to have a minor effect, supported by Guma et al. (2021), who studied the same rotor represented both as a rotor only and as a full turbine. Near the rotor, an O–O type mesh is grown from the blade surface, extruding ≈15 m, discretized with 128 cells from the surface using the grid generator HypGrid (Sørensen1998). 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 Fig. 5. The rotor diameter, D, is ≈80 m.

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. Thereby, the need of searching for hole, fringe, and donor cells between rotor and disc mesh for each time step is avoided, as the relations between the two meshes remain the same.

All deformation from the rotor is propagated to the rotor mesh in such a way that only cells that lie inside the hole region of the overlapping disc mesh deform. This is done to avoid deformation of the donor cells, keeping the interpolation coefficients between fringe and donor cells unaltered. Through this simplification, there is no need for updating communication tables for donor and receiver cells between the rotor and disc mesh as these also rotate together. This choice, however, necessitates that the hole of the disc mesh be far enough from the surface to leave room for the deformation of mesh cells without impairing the cell quality. In the present setup the holes are 17m wide in the rotor axis direction, where the main deformation is present, located with the undeformed rotor in the center. Displacements are propagated to the volume mesh, such that points within the inner 15 % of the grid curve length normal to the surface are moved as solid body motion to ensure no change of quality of the inner cells resolving the high gradient flow. Further out, from 15 to 40 % the volume blend factor is linearly decaying from 1 to 0, such that points from 40 % grid curve length and out are unchanged to avoid changes to donor cells. Note that communication tables and hole cutting still need an update between disc mesh and background mesh as the latter is static. The disc and rotor grids are similar to the setup used in Grinderslev et al. (2021); however in this study the background grid has changed to be suitable for LES simulations by using rectangular cells with low cell stretching in the area of focus.

Figure 5Near-rotor mesh at 25 m span and surface discretization at tip. Only every second line shown.


The background domain is a box of 1000×600×1350 m (12.5×7.5×16.9D) (width × height × length) using 352×256×640 cells, adding up to ≈58 million cells. A concentration of cells is present in the cross-flow directions around the rotor area down to 1 m side lengths; see Fig. 4 (right). Cells in the flow direction are kept constant at ≈1.4 m from the inlet to the rotor and 6 D behind it, before stretching towards the outlet. Boundary conditions are velocity inlet, outlet assuming fully developed flow, and symmetry conditions (slip walls) on sides and top boundaries. The ground has a no-slip wall condition, but with the SG wall model as in the precursor simulation. The rotor is placed ≈4.38D from the inlet, ≈6.25D from sides and top, and ≈12.5D from the outlet.

A total of 78 million cells are used for the combined setup.

3 Results

3.1 Precursor simulation

A total of 9750 s was simulated for the precursor simulation, of which the final 1000 s (equivalent to ≈270 rotor revolutions) was sampled for statistical postprocessing, in a period where the developed flow profile sufficiently matched the desired profile. The precursor was run on three grid levels with varying time steps: first the coarse period (δz=δy=δx16 m, Δt=1.0 s), then a medium period (δz=δy=δx8 m, Δt=0.5 s), and finally a fine period (δz=δy=δx4 m, Δt=0.25 s). The full precursor simulation was conducted on 1728 AMD EPYC 2.9 GHz processors on the computer cluster of DTU and lasted ≈45 wall clock hours. The sampling was conducted in the fine phase only as depicted in Fig. 6 along with the pre-multiplied spectra fS(f) at three different altitudes. As seen, the turbulence is well resolved with a decent inertial subrange following the Kolmogorov spectrum law with a decaying slope of -2/3.

Figure 6(a) Time series of wind speed, W, at three altitudes, approximately matching the rotor bottom, center, and top altitudes. (b) Spectra of wind speed time series (fine-resolution period only) using the Welch estimate.


From the sampling plane, depicted in Fig. 2, the wind speed profiles of W were extracted and horizontally and temporally averaged ±1D from the rotor position in the cross-plane direction as depicted in Fig. 7 (left). As seen, the relative difference of the averaged profile and the DanAero log-law fit match well with a maximum of 8 % at ≈14 m, which corresponds to only ≈0.5 m s−1 at the specific altitude. One difference to note, however, is the larger standard deviation, and thus turbulence intensity, of the sampled flow, with fluctuations that supersede the DanAero measurements. The complexity of fitting both the mean profile and turbulence intensity between measurements and LES simulations is high. In this specific case, the assumption of neutral stratification in the simulation, with no knowledge about stratification being available from the measurements, likely plays a role in the capabilities to match results. This was the best match obtained after multiple calibration attempts, considering both mean profile and turbulence intensity.

Figure 7(a) Horizontal and temporal average profile μW ±1σ (solid red and red patch), DanAero measurements and fitted log law (blue error bars and dashed). Horizontal averaging based on flow from ±1D from the rotor center on the sampled flow plane. (b) Relative error between log-law fit and μW profile.


Figure 8 depicts the resulting resolved and SGS flow shear stresses and resolved friction velocity, u*.

Figure 8Precursor results: (a) horizontal averages of resolved and SGS shear stresses at 9250 s and (b) resulting resolved friction velocity (u*)2=vw2+uv2.


3.2 Successor simulation

In the following, the results of the successor simulations are presented. First, the new turbulence model is compared to the same setup using only the IDDES turbulence model assuming a elastically stiff configuration. Further, results from simulations using the hybrid model with and without flexibility of the blades are presented to study the effect of the blade elasticity. For the initial phase 1 (see Fig. 3), simulations were conducted on 1189 AMD EPYC 2.9 GHz processors, while coupled simulations (phases 3 and 4 in Fig. 3) were conducted on 793 Intel Xeon 2.8 GHz processors. The initialization simulations took on the order of ≈26 wall clock hours, while the coupled simulations lasted for ≈180 wall clock hours per simulation. Figure 9, shows the Q criterion = 0.4 (Hunt et al.1988) of the flow, visualizing the turbulent structures up- and especially downwind of the rotor. As seen, the tip vortices in the wake are quickly broken up into smaller structures by the surrounding turbulent flow. This is also visible in Fig. 10, showing the flow velocity W at multiple downstream positions.

Figure 9Isobars of Q criterion = 0.4 colored with value of flow velocity W, ranging from −5.0 to 15.0 m s−1.


Figure 10Velocity W downstream of rotor.


Figure 11Instantaneous sections of flow velocity W for hybrid (a) and IDDES (b) simulations. Black dashed lines indicate the locations of the profiles presented in Fig. 12.


3.2.1 Impact of turbulence model

To study the impact of the presented turbulence model on the flow, simulations with the hybrid LES–IDDES blending enabled along with pure IDDES simulations are conducted. In the pure IDDES simulation, a slip wall condition is used on the terrain surface, contrary to the log law used for the LES–IDDES hybrid model. Simulations with and without the rotor present were simulated. In the empty setup, the hybrid model acts as a pure Deardorff LES model, as no blending region is defined. For all simulations, inflow is interpolated from the LES precursor planes to ensure identical inlet conditions. In the simulations comparing turbulence models, only the CFD code has been used, meaning that no flexibility of the blades is considered.

Firstly, the empty setups are presented in Fig. 11, showing the velocity component W at a vertical plane aligned with the flow direction intersecting the rotor center, for the simulations with Deardorff and IDDES turbulence modeling at the same time instance. From the planes, instantaneous velocity profiles and turbulence intensity profiles are extracted along the dashed lines, which are shown in Fig. 12. Both simulations show very comparable results. As seen, the velocity and TI profiles, extracted 96 m from the inlet, are practically identical, while a discrepancy is seen further downstream in the domain as a result of changing the turbulence and wall models. While both turbulence models perform close to identically for the present single-turbine-domain study, it could be speculated that this would not be the case when considering larger domains, for instance when studying multiple turbines at once. In that case, the differences in turbulence and wall modeling will likely result in different flow fields, due to the longer distances covered. Considering temperature effects on the flow might likewise reveal differences between the turbulence models. This could be in terms of temperature flux differences, along with the Deardorff model considering the stability in the mixing length scale, which is not the case in the IDDES model.

Figure 12Instantaneous sampled profiles at time = 129.6 s. (a, b) Flow velocity W, (c, d) turbulence intensity (TI). Samples extracted from planes close to the inlet (96 m from inlet) and far downstream (691 m from inlet); see dashed lines of Fig. 11.


Stiff simulations covering 35 rotor revolutions were also conducted with the two turbulence models, including the rotor in the simulations. Mean and standard deviations of azimuthal forces of the final 15 revolutions at two blade sections, near the mid-span and near the tip, are presented in Fig. 13. Only slight differences are seen in both mean and standard deviations between the two models, aligning well with what is seen in the empty domain simulations. As the incoming flow is not altered significantly by the choice of turbulence model, and the turbulence model near the blade is IDDES in both simulations, the forces are expected to be similar as well.

Figure 13Normal and tangential forces at 48 and 92 % blade length using hybrid or IDDES turbulence model. Temporal means and standard deviations based on the final 15 revolutions.


3.2.2 Impact of flexibility

To study the effect of the rotor flexibility, FSI simulations of flexible and stiff setups were performed. First, 35 revolutions were simulated through pure CFD, as presented before, followed by 25 revolutions with the FSI coupling enabled; see Fig. 3 for the FSI simulation process.

The following results are obtained using the hybrid turbulence model only, but similar results would be expected for pure IDDES simulations, based on the aforementioned findings. The effect of including the blade flexibility is assessed through the resulting blade displacements, torsion, and rotor loading.

Figure 14 depicts the tip displacement flap- and edgewise along with the resulting blade torsion at 60.1 and 95.3 % blade length. The tip displacement in the flapwise direction is ≈6 % of the blade length, with fluctuations up to ≈1 % due to the turbulent flow. Edgewise displacements are low and dominated by gravity, seen in the more regular pattern and low standard deviation.

Figure 14(a) Tip displacements flap- and edgewise. (b) Blade torsional deformation at 60.1 and 95.3 % blade length. Crosses represent instantaneous realizations.


Figure 15Integrated thrust and torque for stiff and flexible configurations.


Blade torsion is quite low as well, with less than 0.5 near the tip, increasing the angle of attack.

In Fig. 15, the integrated rotor thrust and torque are depicted, showing that in general only slight differences are seen by including flexibility. This is seen, with an increase in thrust of 1–5 %, while no significant change is seen in torque, other than a slight decrease in fluctuation amplitudes when including flexibility. The large- and small-scale fluctuations of the signals further indicate that the turbulent inflow has a higher impact on torque and thrust than the change seen from considering flexibility.

As mentioned, some differences are present in the simulation setup compared to the DanAero field experiment, being the omission of yaw, tilt, and tower along with the higher turbulence intensity of the generated flow.

Despite this disclaimer, the resulting forces at four sections of the blade are depicted in Fig. 16, showing the mean azimuthal pressure forces normal and tangential to chord for both flexible and stiff simulations along with the DanAero measurements. As seen, the forces agree well between the two simulations and the measurements, with the main differences being the lack of tower shadow at the inner sections, resulting in a drop of loading in the measurements not seen in the simulations (see Fig. 16, top row). The standard deviation of the forces is seen to be higher in simulations than for measurements, which is expected as the turbulence intensity of the sampled flow is higher than measurements as seen in Fig. 7.

Figure 16Normal and tangential forces for stiff and flexible simulations along with DanAero measurements. Sections at 33, 48, 76, and 92 % blade length.


The impact of including flexibility is quite small, and general observations are that normal and tangential forces respectively slightly increase and decrease when considering the flexibility of the rotor. This is expected for the NM80 rotor, which is quite stiff compared to modern wind turbines. The standard deviations of the forces due to turbulence are much higher than the difference between mean forces of stiff and flexible simulations. This shows that including turbulent inflow is more important than including flexibility, at least in the present rotor/flow case. In the simulations including the flexibility the standard deviations of the normal forces are up to 10 % of the mean near the tip and 15% near the root. For tangential forces this is even higher, with 24 % near the tip and 38 % near the root.

Spectral analysis of the resulting normal and tangential force signals at the 76 % blade length section are presented in Fig. 17 (left), showing the power spectral densities (PSDs) using the Welch estimate. As seen, both stiff and flexible simulations result in similar PSDs, with the main difference being the peak at the first edgewise frequency seen in the flexible signal. The majority of energy is found in the rotation frequency, 1P, and its harmonics. This is also the case when looking at the PSD of the tip displacement in the flap- and edgewise direction. Here, it is again the rotation frequency and its harmonics that dominate, along with a peak of the first edgewise mode.

Figure 17(a) PSD of normal and tangential forces at 78 % blade length for both stiff and flexible simulations. (b) PSD of tip displacement in flap- and edgewise direction.


4 Conclusions

This study investigates the phenomenon of aero-elasticity of wind turbines placed in atmospheric flow conditions, by means of high-fidelity numerical methods. Fluid–structure interaction (FSI) simulations of a 2.3 MW wind turbine rotor have been conducted using a novel turbulence model, blending the Deardorff large-eddy simulation (LES) model for atmospheric flows with the improved delayed detached eddy simulation (IDDES) model for the separated flow near the rotor boundary. Precursor simulations were conducted in a large domain in order to assure sampling of realistic turbulent atmospheric boundary layer (ABL) flow, matching well with the DanAero measurements, for the successor simulations.

As the first study, the hybrid model was compared to the pure IDDES turbulence model, by computational fluid dynamics (CFD) successor simulations of the turbulent ABL inflow with and without the rotor present. In empty simulations, this corresponded to a comparison between pure Deardorff LES and pure IDDES, while for rotor simulations the hybrid model used both Deardorff LES for the domain flow and IDDES for the near-rotor flow. It was found that there was no significant difference in either the flow or the rotor loading between the two methods, likely due to the short domain considered and assumptions omitting the Coriolis force and temperature effects.

Secondly, FSI simulations were conducted by coupling the CFD simulations to a structural solver. It was found that for the specific rotor, which is relatively stiff compared to modern turbines, only a small impact was found by considering the flexibility of the blades. A general increase of ≈1–5 % in total thrust was found, while the power producing torque was close to identical for stiff and flexible simulations.

Inflow turbulence on the other hand has a large influence on the rotor loading, with standard deviations as high as 15 % of the mean for normal forces and even higher tangentially. This emphasizes the importance of correct modeling of inflow turbulence.

5 Future studies

As has been shown in the present study, the developed hybrid turbulence model resulted in practically identical loading of the rotor to that of the IDDES model alone. Relevant future studies would be to investigate when this is not the case. This could for instance be simulations including stable/unstable stratification and/or the Coriolis force. Here, the IDDES model will probably be insufficient to capture the effects, as the model is calibrated for aerodynamics mainly and not ABL flows. The Deardorff LES model, however, is calibrated for such flows, and the mixing length scale depends on the stratification, as it is reduced for stable cases. The two models also model temperature effects differently as the flux (τθ,w) is based on different weights of the turbulent viscosity. Longer domains with multiple rotors could also be relevant, as there are time and distance for the two turbulence models to develop the flow differently, which is indicated by the results of the present study.

A relevant future study would likewise be to compare the method to more efficient BEM-based aerodynamics solvers with the precursor turbulence as input. Here, the CFD-based results could, if needed, be used to correct airfoil polars and calibrate the many correction models needed by BEM solvers to consider e.g. tip loss effects, dynamic inflow, and dynamic stall.

In terms of FSI, it would be natural to investigate more recent/future turbine designs, which are larger and much more flexible than the considered NM80 rotor. These rotors are at a higher risk of instability phenomena and operate in a larger part of the atmospheric boundary layer.

Code and data availability

The codes used to conduct the presented simulations are licensed and not publicly available. The structural model of the NM80 rotor is likewise not publicly available. However, result data can be shared upon request.

Author contributions

CG conducted all simulations, including grid generation, and implemented the presented turbulence model in the CFD code. Additionally CG did the analysis of the results and has been the main writer of the paper. NNS assisted with expertise in the code development and CFD setup, along with planning the study and analysis. SGH has supported in especially FSI setup along with planning and analysis. NT supported the development and implementation of the turbulence model. FZ assisted in the overset grid setup, and did development on the CFD code to accelerate simulations. Further, all authors contributed in writing and editing this paper.

Competing interests

The authors declare that they have no conflicts of interest.


The DanAero projects, from which experimental data were obtained, were funded partly by the Danish Energy Agency (EFP2007. Journal no. 33033-0074 and EUDP 2009-II. Journal no. 64009-0258) and partly by internal funding from the project partners.

Financial support

This study was conducted as part of the PhD project “Fluid-structure Interaction for Wind Turbines in Atmospheric Flow”, which was internally funded by the Technical University of Denmark (DTU) – Department of Wind Energy.

Review statement

This paper was edited by Johan Meyers and reviewed by two anonymous referees.


Bak, C., Johansen, J., and Andersen, P.: Three-dimensional corrections of airfoil characteristics based on pressure distributions (paper and poster), in: Proceedings (online), European Wind Energy Association (EWEA), 2006 European Wind Energy Conference and Exhibition, EWEC 2006, 27 February–2 March 2006, available at: (last access: 3 March 2021), 2006. a

Bak, C., Madsen, H., Paulsen, U. S., Gaunaa, M., Sørensen, N., Fuglsang, P., Romblad, J., Olsen, N., Enevoldsen, P., Laursen, J., and Jensen, L.: DAN-AERO MW: Detailed aerodynamic measurements on a full scale MW wind turbine, in: European Wind Energy Conference and Exhibition 2010, 2, Ewec 2010, 20–23 April 2010, Warsaw, Poland, 792–836, 2010. a

Bak, C., Troldborg, N., and Madsen, H.: DAN-AERO MW: Measured airfoil characteristics for a MW rotor in atmospheric conditions, in: Scientific Proceedings, European Wind Energy Association (EWEA), 4–17 March 2011, 171–175, available at: (last access: 3 May 2021), 2011. a

Bak, C., Madsen, H., Troldborg, N., and Wedel-Heinen, J.: DANAERO MW: Data for the NM80 turbine at Tjæreborg Enge for aeroelastic evaluation, Tech. rep., Technical University of Denmark, Denmark, 2013. a

Bechmann, A., Sørensen, N., and Zahle, F.: CFD simulations of the MEXICO rotor, Wind Energy, 14, 677–689, 2011. a

Berg, J., Natarajan, A., Mann, J., and Patton, E.: Gaussian vs non-Gaussian turbulence: impact on wind turbine loads, Wind Energy, 19, 1975–1989,, 2016. a

Deardorff, J. W.: Numerical Investigation of Neutral and Unstable Planetary Boundary Layers, J. Atmosl. Sci., 29, 91–115,<0091:NIONAU>2.0.CO;2, 1972. a, b

Deardorff, J. W.: Stratocumulus-capped mixed layers derived from a three-dimensional model, Bound.-Lay. Meteorol., 18, 495–527,, 1980. a

Egorov, Y., Menter, F. R., Lechner, R., and Cokljat, D.: The scale-adaptive simulation method for unsteady turbulent flow predictions. part 2: Application to complex flows, Flow Turbul. Combust., 85, 139–165,, 2010. a

García Ramos, N., Sessarego, M., and Horcas, S. G.: Aero–hydro–servo–elastic coupling of a multi-body finite-element solver and a multi-fidelity vortex method, Wind Energy, 24, 481–501,, 2021. a

Glauert, H.: Airplane Propellers, Division L, in: Aerodynamic Theory 4, Springer, Berlin, Heidelberg,, 1935. a

Grinderslev, C., Belloni, F., Horcas, S., and Sørensen, N.: Investigations of aerodynamic drag forces during structural blade testing using high-fidelity fluid–structure interaction, Wind Energ. Sci., 5, 543–560,, 2020a. a

Grinderslev, C., Vijayakumar, G., Ananthan, S., Sørensen, N., Zahle, F., and Sprague, M.: Validation of blade-resolved computational fluid dynamics for a MW-scale turbine rotor in atmospheric flow, J. Phys.: Conf. Ser., 1618, 052049,, 2020b. a

Grinderslev, C., Horcas, S., and Sørensen, N.: Fluid–structure interaction simulations of a wind turbine rotor in complex flows, validated through field experiments, Wind Energy,, in press, 2021. a, b, c, d, e

Gritskevich, M. S., Garbaruk, A. V., Schütze, J., and Menter, F. R.: Development of DDES and IDDES formulations for the kω shear stress transport model, Flow Turbul. Combust., 88, 431–449,, 2012. a

Grötzbach, G.: Direct numerical and large eddy simulation of turbulent channel flows, Encycloped. Fluid Mech., 6, 1337–1391, 1987. a

Guma, G., Bangga, G., Lutz, T., and Krämer, E.: Aeroelastic analysis of wind turbines under turbulent inflow conditions, Wind Energ. Sci., 6, 93–110,, 2021. a, b

Hansen, M., Gaunaa, M., and Madsen, H.: A Beddoes-Leishman type dynamic stall model in state-space and indicial formulations, Tech. rep., Risø National Laboratory, Risø, Denmark, 2004. a

Heinz, J., Sørensen, N., and Zahle, F.: Fluid–structure interaction computations for geometrically resolved rotor simulations using CFD, Wind Energy, 19, 2205–2221, 2016a. a, b, c

Heinz, J., Sørensen, N., Zahle, F., and Skrzypinski, W.: Vortex-induced vibrations on a modern wind turbine blade, Wind Energy, 19, 2041–2051,, 2016b. a

Horcas, S., Madsen, M., Sørensen, N., and Zahle, F.: Suppressing Vortex Induced Vibrations of Wind Turbine Blades with Flaps, in: Recent Advances in CFD for Wind and Tidal Offshore Turbines, edited by: Ferrer, E. M. A., Springer, Switzerland, 11–24, 2019. a, b

Horcas, S., Barlas, T., Zahle, F., and Sørensen, N.: Vortex induced vibrations of wind turbine blades: Influence of the tip geometry, Phys. Fluids, 32, 065104,, 2020. a

Hunt, J. C. R., Wray, A. A., and Moin, P.: Eddies, streams, and convergence zones in turbulent flows, Studying Turbulence Using Numerical Simulation Databases, NASA, USA, 193–208, 1988. a

Jonkman, J. M. and Buhl, M. L. J.: FAST User's Guide, Tech. rep., National Renewable Energy Laboratory, National Renewable Energy Laboratory Golden, Colorado, 2005. a

Larsen, T. and Hansen, A.: How 2 HAWC2, the user's manual, Tech. rep., Risø National Laboratory, Risø, Denmark, 2007. a, b

Lee, S., Churchfield, M. J., Moriarty, P. J., Jonkman, J., and Michalakes, J.: A numerical study of atmospheric and wake turbulence impacts on wind turbine fatigue loadings, T. ASME J. Sol. Energ. Eng., 135, 1–10,, 2013. a

Li, Y., Castro, A. M., Sinokrot, T., Prescott, W., and Carrica, P. M.: Coupled multi-body dynamics and CFD for wind turbine simulation including explicit wind turbulence, Renew. Energy, 76, 338–361,, 2015. a

Madsen, H. A., Larsen, T. J., Pirrung, G. R., Li, A., and Zahle, F.: Implementation of the blade element momentum model on a polar grid and its aeroelastic load impact, Wind Energ. Sci., 5, 1–27,, 2020. a, b

Mann, J.: Wind field simulation, Probabil. Eng. Mech., 13, 269–282,, 1998. a

Menter, F.: Zonal Two Equation Kappa–Omega Turbulence Models for Aerodynamic Flows, in: 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference,, 1993. a, b

Menter, F., Kuntz, M., and Langtry, R.: Ten years of industrial experience with the SST turbulence model, Heat Mass Transf., 4, 1–8, 2003. a

Michelsen, J.: Basis3D – A Platform for Development of Multiblock PDE Solvers, Tech. rep., Risø National Laboratory, Risø, Denmark, 1992. a, b

Michelsen, J.: Block Structured Multigrid Solution of 2D and 3D Elliptic PDE's, Tech. rep., Technical University of Denmark, Denmark, 1994. a, b

Moeng, C.: A Large-Eddy-Simulation Model For The Study OF Planetary Boundary-Layer Turbulence, J. Atmos. Sci., 41, 2052–2062, 1984. a

Pavese, C., Wang, Q., Kim, T., Jonkman, J., and Sprague, M.: HAWC2 and BeamDyn: Comparison Between Beam Structural Models for Aero-Servo-Elastic Frameworks, in: Proceedings of the EWEA Annual Event and Exhibition 2015, European Wind Energy Association (EWEA), EWEA Annual Conference and Exhibition 2015, 17–20 November 2015, Paris, France, 2015. a

Santo, G., Peeters, M., Van Paepegem, W., and Degroote, J.: Effect of rotor–tower interaction, tilt angle, and yaw misalignment on the aeroelasticity of a large horizontal axis wind turbine with composite blades, Wind Energy, 23, 1578–1595,, 2020a. a, b

Santo, G., Peeters, M., Van Paepegem, W., and Degroote, J.: Fluid–Structure Interaction Simulations of a Wind Gust Impacting on the Blades of a Large Horizontal Axis Wind Turbine, Energies, 13, 509,, 2020b. a, b

Schumann, U.: Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli, J. Comput. Phys., 18, 376–404, 1975. a

Shen, W., Michelsen, J., Sørensen, N., and Sørensen, J.: An improved SIMPLEC method on collocated grids for steady and unsteady flow computations, Numer. Heat Transf. Pt. B, 43, 221–239, 2003. a

Shur, M. L., Spalart, P. R., Strelets, M. K., and Travin, A. K.: A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities, Int. J. Heat Fluid Flow, 29, 1638–1649,, 2008. a, b

Sørensen, J. N. and Shen, W. Z.: Numerical Modeling of Wind Turbine Wakes, J. Fluids Eng., 124, 393–399,, 2002. a

Sørensen, J. N., Mikkelsen, R. F., Henningson, D. S., Ivanell, S., Sarmast, S., and Andersen, S. J.: Simulation of wind turbine wakes using the actuator line technique, Philos. T. Roy. Soc. A, 373, 20140071,, 2015. a

Sørensen, N.: General purpose flow solver applied to flow over hills, PhD thesis, Risø National Laboratory, Risø, Denmark, 1995. a, b

Sørensen, N.: HypGrid2D. A 2-d mesh generator, Tech. rep., Risø National Laboratory, Risø, Denmark, 1998. a

Sørensen, N. and Schreck, S.: Transitional DDES computations of the NREL Phase-VI rotor in axial flow conditions, J. Phys.: Conf. Ser., 555, 012096,, 2014. a

Sørensen, N., Zahle, F., Boorsma, K., and Schepers, G.: CFD computations of the second round of MEXICO rotor measurements, J. Phys.: Conf. Ser., 753, 022054,, 2016. a

Sørensen, N. N., Michelsen, J. A., and Schreck, S.: Navier–Stokes predictions of the NREL phase VI rotor in the NASA Ames 80 ft × 120 ft wind tunnel, Wind Energy, 5, 151–169,, 2002.  a

Spalart, P., Jou, W., Strelets, M., and Allmaras, S.: Comments on the Feasibility of LES for Wings, and on a Hybrid RANS/LES Approach, in: International conference, 1st Advances in DNS/LES: Direct numerical simulation and large eddy simulation, Greyden Press, USA, 1997. a

Spalart, P. R., Deck, S., Shur, M. L., Squires, K. D., Strelets, M. K., and Travin, A.: A new version of detached-eddy simulation, resistant to ambiguous grid densities, Theor. Comput. Fluid Dynam., 20, 181–195,, 2006. a

Sprague, M., Ananthan, S., Vijayakumar, G., and Robinson, M.: ExaWind: A multi-fidelity modeling and simulation environment for wind energy, in: Proceedings of NAWEA WindTech, 14–16 October 2019, Amherst, Massachusetts, USA, 2019. a

Strelets, M.: Detached eddy simulation of massively separated flows, in: 39th Aerospace Sciences Meeting and Exhibit, 8–11 January 2001, Reno, USA, 2001. a

Travin, A., Shur, M., Strelets, M., and Spalart, P. R.: Physical and numerical upgrades in the Detached-Eddy Simulation of complex turbulent flows, Fluid Mech. Appl., 65, 239–254, 2004. a

Troldborg, N., Sørensen, J., Mikkelsen, R., and Sørensen, N.: A simple atmospheric boundary layer model applied to large eddy simulations of wind turbine wakes, Wind Energy, 17, 657–669, 2014. a, b

Troldborg, N., Zahle, F., Réthoré, P., and Sørensen, N.: Comparison of wind turbine wake properties in non‐sheared inflow predicted by different computational fluid dynamics rotor models, Wind Energy, 18, 1239–1250,, 2015. a

Vijayakumar, G.: Non-steady dynamics of atmospheric turbulence interaction with wind turbine loadings through blade-boundary-layer-resolved CFD, PhD Thesis, The Pennsylvania State University, USA, 2015. a

Vijayakumar, G. J. B., Lavely, A., Jayaraman, B., and Craven, B.: Interaction of Atmospheric Turbulence with Blade Boundary Layer Dynamics on a 5 MW Wind Turbine using Blade-boundary-layer-resolved CFD with hybrid URANS-LES, in: 34th Wind Energy Symposium, 4–8 January 2016, San Diego, California, USA,, 2016. a, b

Zahle, F., Sørensen, N., and Johansen, J.: Wind Turbine Rotor-Tower Interaction Using an Incompressible Overset Grid Method, Wind Energy, 12, 594–619,, 2009. a, b

Short summary
This study investigates aero-elasticity of wind turbines present in the turbulent and chaotic wind flow of the lower atmosphere, using fluid–structure interaction simulations. This method combines structural response computations with high-fidelity modeling of the turbulent wind flow, using a novel turbulence model which combines the capabilities of large-eddy simulations for atmospheric flows with improved delayed detached eddy simulations for the separated flow near the rotor.