the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Wind turbines in atmospheric flow: fluid–structure interaction simulations with hybrid turbulence modeling
Niels Nørmark Sørensen
Sergio González Horcas
Niels Troldborg
Frederik Zahle
In order to design future large wind turbines, knowledge is needed about the impact of aeroelasticity 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 highfidelity rotorresolved 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 wellestablished 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 largeeddy 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 multibody 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.
 Article
(9339 KB)  Fulltext XML
 BibTeX
 EndNote
As future wind turbines will have unprecedentedly long and flexible blades, the necessity of understanding the effects of aeroelasticity 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 highfidelity 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 lowerfidelity 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 Shen, 2002; 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 actuatorlinebased 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 largeeddy simulations (LES) were conducted. The structural response of the turbines was found through the FAST aeroelastic code (Jonkman and Buhl, 2005) 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 rotorresolved 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 finiteelement shells, were studied for steady ABL flows (Santo et al., 2020a, b) using unsteady Reynoldsaveraged 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 multibody 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 largescale 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 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 socalled 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 URANSbased 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 URANSbased 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 highperformance computing capabilities, compromises are needed when doing highfidelity aeroelastic 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 aeroelasticity 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 (Michelsen, 1992, 1994; Sørensen, 1995) and the structural solver from the aeroelastic code HAWC2 (Larsen and Hansen, 2007). 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.
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 inhouse CFD code EllipSys3D (Michelsen, 1992, 1994; Sørensen, 1995) is used. The code solves the incompressible Navier–Stokes equations in structured curvilinear coordinates using the finitevolume method with a collocated grid arrangement. The code is parallel and highly scalable using the message passing interface (MPI) and multiblock decomposition, the multigrid method and grid sequencing. EllipSys3D has multiple convective schemes implemented, such as central difference (CDS), secondorder 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 twoequation Reynoldsaveraged 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 Schreck, 2014; 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 NaluWind (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 DanAeroinspired laminar wind flow and the structural response.
2.1.2 Aeroelastic solver
HAWC2 (Larsen and Hansen, 2007) is an aeroelastic 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 multibody dynamics formulation, accounting for nonlinear 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. Subbodies are connected with constraint equations considering nonlinearities.
HAWC2 has a builtin aerodynamics module that calculates aerodynamic forces using BEM theory. As is common in BEM implementations, prediction of airfoil aerodynamic performance is based on precomputed lookup 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 windenergyrelated 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 largescale 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 secondorder 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 underrelaxed subiterations 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 oneequation LES turbulence model for ABL flows (Deardorff, 1972) 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 subgridscale (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 largeeddy simulation model
In the Deardorff LES turbulence model (Deardorff, 1972, 1980), the turbulent eddy viscosity μ_{t} is calculated through the expression
Here, C_{k} is a constant of 0.1; ρ is the air density; and l_{LES} is a mixing length scale, which for neutral stratification is set equal to Δ_{LES}, which is the grid size, here defined as ${\mathrm{\Delta}}_{\mathrm{LES}}=(\mathrm{d}x\cdot \mathrm{d}y\cdot \mathrm{d}z{)}^{\mathrm{1}/\mathrm{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:
where g, t, and μ refer to the gravity, time, and molecular viscosity, respectively. C_{ϵ} is equal to 0.93, and the buoyancy SGS fluxes ${\mathit{\tau}}_{\mathit{\theta}i,\mathrm{LES}}={\mathit{\mu}}_{\mathit{\theta}}\frac{\partial \mathit{\theta}}{\partial {x}_{i}}$, with variable temperature θ and the eddy heat diffusivity being ${\mathit{\mu}}_{\mathit{\theta}}=\left(\mathrm{1}+\frac{\mathrm{2}{l}_{\mathrm{LES}}}{{\mathrm{\Delta}}_{\mathrm{LES}}}\right){\mathit{\mu}}_{\mathrm{t}}$. θ_{0} is the surface reference temperature. The SGS stress tensor τ_{ij} is defined as ${\mathit{\tau}}_{ij}=\mathrm{2}{\mathit{\mu}}_{\mathrm{t}}{S}_{ij}$ using the strain rate ${S}_{ij}=\mathrm{1}/\mathrm{2}\left(\frac{\partial {\mathit{u}}_{i}}{\partial {x}_{j}}+\frac{\partial {\mathit{u}}_{j}}{\partial {x}_{i}}\right)$, with u being the velocity vector.
2.2.2 SSTbased detached eddy simulation models
For the k–ω SSTbased DES turbulence models, μ_{t} is found through the standard k–ω SST (Menter, 1993) approach, which is then altered in the dissipation term depending on the chosen DES model.
Here, k is a total turbulent kinetic energy, ω the specific dissipation rate, Ω the shearstrain rate, and F_{2} a limiting blending function. k and ω are found through the following two transport equations:
For k,
For ω,
Here F_{1} is a blending function, to shift between the standard k–ω model (F_{1}=1) near the surface and the k–ϵ model (F_{1}=0) within the boundary layer and further out, while σ_{k}, σ_{ω}, β, and γ are parameters which themselves depend on F_{1}. Finally β^{*} and σ_{ω2} are constants. ν_{t} is the kinematic turbulent viscosity ${\mathit{\nu}}_{\mathrm{t}}={\mathit{\mu}}_{\mathrm{t}}/\mathit{\rho}$. All constants and parameters as well as blending functions F_{1} and F_{2} can be found in the original work by Menter (Menter, 1993). 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 righthand side of Eq. (4), and the flux is found as ${\mathit{\tau}}_{\mathit{\theta}i,\mathrm{DES}}={\mathit{\mu}}_{\mathrm{t}}\frac{\partial \mathit{\theta}}{\partial {x}_{i}}$.
The length scale which appears in the kequation serves to switch from URANS to LES mode and is defined as:
where ${l}_{k\mathit{\omega}}=\sqrt{k}/\left({\mathit{\beta}}^{*}\mathit{\omega}\right)$ and l_{DES} is the length scale in the LES region. In the standard DES model (Spalart et al., 1997; Travin et al., 2004), l_{DES}=C_{DES}Δ_{DES}, where ${\mathrm{\Delta}}_{\mathrm{DES}}=\mathrm{max}(\mathrm{d}x,\mathrm{d}y,\mathrm{d}z)$ and C_{DES} is a F_{1}dependent parameter.
DES is known to be sensitive to sudden changes of grid refinements as gridinduced 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 wallmodeled 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 $\stackrel{\mathrm{\u0303}}{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 F_{h}, which is zero in the ABL region and one in the DES region, and then defining a hybrid turbulence kinetic energy $\stackrel{\mathrm{\u0303}}{k}={F}_{\mathrm{h}}k+(\mathrm{1}{F}_{\mathrm{h}})e$. Using these definitions, the energy equations Eqs. (2) and (4) are combined to give the following transport equation for $\stackrel{\mathrm{\u0303}}{k}$:
The blending function F_{h} is defined as follows:
Here, d_{w} is the wallnormal distance; R is the wall distance to the location where F_{h}=0.5; and δ_{blend} is the transition distance between F_{h}=0.12 and F_{h}=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.
Here, it is assumed that the blending from DES to LES happens in the region where F_{2}=0, such that the viscosity limiter is inactive. A blended expression $\stackrel{\mathrm{\u0303}}{\mathit{\omega}}$ is then found for the entire domain.
This allows the calculation of the turbulent viscosity similar to Eq. (3):
It is noted that in the Deardorff part of the model the turbulent viscosity, μ_{t}, is linearly proportional to the length scale, l_{LES}, 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 loglaw wind profile to generate inputs for the Schumann–Grötzbach (SG) wall model (Schumann, 1975; Grötzbach, 1987) used in the simulation.
where U is the wind speed, u_{*} the friction velocity, κ the von Kármán constant (≈0.4), z the vertical coordinate, and finally z_{0} 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 $\mathrm{1.769}\times {\mathrm{10}}^{\mathrm{5}}$ kg m^{−1} s^{−1}. The convective terms are calculated through a blend of the fourthorder 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 $\mathrm{1.29}\times {\mathrm{10}}^{\mathrm{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 smallscale 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 subgrids, 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 (F_{h}=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.
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 BEMbased aeroelastic 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 (Glauert, 1935) 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 z_{0} 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 crossflow 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.
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 spinup.
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 twoway coupling is simulated for the desired amount of revolutions.
2.5 Computational grids
2.5.1 Precursor simulation
The precursor domain is $\mathrm{4000}\times \mathrm{1300}\times \mathrm{8000}$ m (width × height × length) discretized $\mathrm{576}\times \mathrm{256}\times \mathrm{1920}$ cells divided into 8640 blocks of 32^{3} cells. A total of ≈283 million cells are present in the precursor. The grid cells vary in size in the crossflow 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ørensen, 1998). The first cell adjacent to the rotor surface has a height of $\mathrm{1}\times {\mathrm{10}}^{\mathrm{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 precut holes around the blades. This mesh rotates along with the rotor mesh, speeding up the holecutting 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.
The background domain is a box of $\mathrm{1000}\times \mathrm{600}\times \mathrm{1350}$ m ($\mathrm{12.5}\times \mathrm{7.5}\times \mathrm{16.9}$ D) (width × height × length) using $\mathrm{352}\times \mathrm{256}\times \mathrm{640}$ cells, adding up to ≈58 million cells. A concentration of cells is present in the crossflow 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 noslip wall condition, but with the SG wall model as in the precursor simulation. The rotor is placed ≈4.38 D from the inlet, ≈6.25 D from sides and top, and ≈12.5 D from the outlet.
A total of 78 million cells are used for the combined setup.
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 ($\mathit{\delta}z=\mathit{\delta}y=\mathit{\delta}x\approx \mathrm{16}$ m, Δt=1.0 s), then a medium period ($\mathit{\delta}z=\mathit{\delta}y=\mathit{\delta}x\approx \mathrm{8}$ m, Δt=0.5 s), and finally a fine period ($\mathit{\delta}z=\mathit{\delta}y=\mathit{\delta}x\approx \mathrm{4}$ 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 premultiplied spectra f⋅S(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 $\mathrm{2}/\mathrm{3}$.
From the sampling plane, depicted in Fig. 2, the wind speed profiles of W were extracted and horizontally and temporally averaged ±1 D from the rotor position in the crossplane direction as depicted in Fig. 7 (left). As seen, the relative difference of the averaged profile and the DanAero loglaw 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 8 depicts the resulting resolved and SGS flow shear stresses and resolved friction velocity, u^{*}.
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.
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 singleturbinedomain 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.
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 midspan 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.
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.
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 smallscale 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.
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.
This study investigates the phenomenon of aeroelasticity of wind turbines placed in atmospheric flow conditions, by means of highfidelity 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 largeeddy 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 nearrotor 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.
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 BEMbased aerodynamics solvers with the precursor turbulence as input. Here, the CFDbased 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.
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.
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.
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. 330330074 and EUDP 2009II. Journal no. 640090258) and partly by internal funding from the project partners.
This study was conducted as part of the PhD project “Fluidstructure Interaction for Wind Turbines in Atmospheric Flow”, which was internally funded by the Technical University of Denmark (DTU) – Department of Wind Energy.
This paper was edited by Johan Meyers and reviewed by two anonymous referees.
Bak, C., Johansen, J., and Andersen, P.: Threedimensional 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: https://orbit.dtu.dk/en/publications/threedimensionalcorrectionsofairfoilcharacteristicsbasedon (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.: DANAERO 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.: DANAERO 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: https://orbit.dtu.dk/en/publications/danaeromwmeasuredairfoilcharacteristicsforamwrotorinat (last access: 3 May 2021), 2011. a
Bak, C., Madsen, H., Troldborg, N., and WedelHeinen, 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 nonGaussian turbulence: impact on wind turbine loads, Wind Energy, 19, 1975–1989, https://doi.org/10.1002/we.1963, 2016. a
Deardorff, J. W.: Numerical Investigation of Neutral and Unstable Planetary Boundary Layers, J. Atmosl. Sci., 29, 91–115, https://doi.org/10.1175/15200469(1972)029<0091:NIONAU>2.0.CO;2, 1972. a, b
Deardorff, J. W.: Stratocumuluscapped mixed layers derived from a threedimensional model, Bound.Lay. Meteorol., 18, 495–527, https://doi.org/10.1007/BF00119502, 1980. a
Egorov, Y., Menter, F. R., Lechner, R., and Cokljat, D.: The scaleadaptive simulation method for unsteady turbulent flow predictions. part 2: Application to complex flows, Flow Turbul. Combust., 85, 139–165, https://doi.org/10.1007/s1049401092654, 2010. a
García Ramos, N., Sessarego, M., and Horcas, S. G.: Aero–hydro–servo–elastic coupling of a multibody finiteelement solver and a multifidelity vortex method, Wind Energy, 24, 481–501, https://doi.org/10.1002/we.2584, 2021. a
Glauert, H.: Airplane Propellers, Division L, in: Aerodynamic Theory 4, Springer, Berlin, Heidelberg, https://doi.org/10.1007/9783642914874_3, 1935. a
Grinderslev, C., Belloni, F., Horcas, S., and Sørensen, N.: Investigations of aerodynamic drag forces during structural blade testing using highfidelity fluid–structure interaction, Wind Energ. Sci., 5, 543–560, https://doi.org/10.5194/wes55432020, 2020a. a
Grinderslev, C., Vijayakumar, G., Ananthan, S., Sørensen, N., Zahle, F., and Sprague, M.: Validation of bladeresolved computational fluid dynamics for a MWscale turbine rotor in atmospheric flow, J. Phys.: Conf. Ser., 1618, 052049, https://doi.org/10.1088/17426596/1618/5/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, https://doi.org/10.1002/we.2639, 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, https://doi.org/10.1007/s1049401193784, 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, https://doi.org/10.5194/wes6932021, 2021. a, b
Hansen, M., Gaunaa, M., and Madsen, H.: A BeddoesLeishman type dynamic stall model in statespace 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.: Vortexinduced vibrations on a modern wind turbine blade, Wind Energy, 19, 2041–2051, https://doi.org/10.1002/we.1967, 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, https://doi.org/10.1063/5.0004005, 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, https://doi.org/10.1115/1.4023319, 2013. a
Li, Y., Castro, A. M., Sinokrot, T., Prescott, W., and Carrica, P. M.: Coupled multibody dynamics and CFD for wind turbine simulation including explicit wind turbulence, Renew. Energy, 76, 338–361, https://doi.org/10.1016/j.renene.2014.11.014, 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, https://doi.org/10.5194/wes512020, 2020. a, b
Mann, J.: Wind field simulation, Probabil. Eng. Mech., 13, 269–282, https://doi.org/10.1016/S02668920(97)000362, 1998. a
Menter, F.: Zonal Two Equation Kappa–Omega Turbulence Models for Aerodynamic Flows, in: 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference, https://doi.org/10.2514/6.19932906, 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 LargeEddySimulation Model For The Study OF Planetary BoundaryLayer 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 AeroServoElastic 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, https://doi.org/10.1002/we.2501, 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, https://doi.org/10.3390/en13030509, 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 RANSLES approach with delayedDES and wallmodelled LES capabilities, Int. J. Heat Fluid Flow, 29, 1638–1649, https://doi.org/10.1016/j.ijheatfluidflow.2008.07.001, 2008. a, b
Sørensen, J. N. and Shen, W. Z.: Numerical Modeling of Wind Turbine Wakes, J. Fluids Eng., 124, 393–399, https://doi.org/10.1115/1.1471361, 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, https://doi.org/10.1098/rsta.2014.0071, 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 2d mesh generator, Tech. rep., Risø National Laboratory, Risø, Denmark, 1998. a
Sørensen, N. and Schreck, S.: Transitional DDES computations of the NREL PhaseVI rotor in axial flow conditions, J. Phys.: Conf. Ser., 555, 012096, https://doi.org/10.1088/17426596/555/1/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, https://doi.org/10.1088/17426596/753/2/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, https://doi.org/10.1002/we.64, 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 detachededdy simulation, resistant to ambiguous grid densities, Theor. Comput. Fluid Dynam., 20, 181–195, https://doi.org/10.1007/s0016200600150, 2006. a
Sprague, M., Ananthan, S., Vijayakumar, G., and Robinson, M.: ExaWind: A multifidelity 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 DetachedEddy 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, https://doi.org/10.1002/we.1757, 2015. a
Vijayakumar, G.: Nonsteady dynamics of atmospheric turbulence interaction with wind turbine loadings through bladeboundarylayerresolved 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 Bladeboundarylayerresolved CFD with hybrid URANSLES, in: 34th Wind Energy Symposium, 4–8 January 2016, San Diego, California, USA, https://doi.org/10.2514/6.20160521, 2016. a, b
Zahle, F., Sørensen, N., and Johansen, J.: Wind Turbine RotorTower Interaction Using an Incompressible Overset Grid Method, Wind Energy, 12, 594–619, https://doi.org/10.1002/we.327, 2009. a, b