the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Investigation of blade flexibility effects on the loads and wake of a 15 MW wind turbine using a flexible actuator line method

### Francois Trigaux

### Philippe Chatelain

### Grégoire Winckelmans

This paper investigates the impact of blade flexibility on the aerodynamics and wake of large offshore turbines using a flexible actuator line method (ALM) coupled to the structural solver BeamDyn in large-eddy simulations. The study considers the IEA 15 MW reference wind turbine in close-to-rated operating conditions. The flexible ALM is first compared to OpenFAST simulations and is shown to consistently predict the rotor aerodynamics and the blade structural dynamics. However, the effect of blade flexibility on the loads is more pronounced when predicted using the ALM compared with using the blade element momentum theory. The wind turbine is then simulated in a neutral turbulent atmospheric boundary layer with flexible and rigid blades. The significant flapwise and torsional mean displacements lead to an overall decrease of 14 % in thrust and 10 % in power compared to a rotor with no deformation. These changes influence the wake through a reduced time-averaged velocity deficit and turbulent kinetic energy. The unsteady loads induced by the rotation in the sheared wind and the turbulent velocity fluctuations are also substantially affected by the flexibility and exhibit a noticeably different spectrum. However, the influence of these load variations on the wake is limited, and the assumption of rigid blades in their deformed geometry is shown to be sufficient to capture the wake dynamics. The influence of the resolution of the flow solver is also evaluated, and the results are shown to remain consistent between different spatial resolutions. Overall, the structural deformations have a substantial impact on the turbine performance, loads, and wake, which emphasizes the importance of considering the flexibility of the blades in simulations of large offshore wind turbines.

- Article
(10963 KB) - Full-text XML
- BibTeX
- EndNote

Over the past few years, offshore wind energy has grown massively thanks to a drastic decrease in its levellized cost of energy. This significant reduction is primarily attributed to the increase in the rotor size, which tends to drive down both capital and operational expenditures. In fact, larger modern turbines typically exhibit a lower turbine cost per megawatt compared to their smaller counterparts. At the level of the wind farm, the decreased number of turbines also allows for reducing the cable requirements (Fingersh et al., 2006). The cost of maintenance per megawatt has also been reported to be a few percent lower per additional megawatt of rated power (Sørensen and Larsen, 2021). Furthermore, larger turbines typically have higher hub heights, which results in a more significant mean flow velocity across the rotor and hence an increased power extraction per unit area (Sørensen and Larsen, 2023). As a result, turbines with a rated power of 10–12 MW are now commonly used in offshore conditions, and the industry is rapidly developing 15–18 MW prototypes with blades much longer than 100 m.

However, as the rotor diameter increases, fundamental changes in design are necessary to constrain the total mass of the blades in order to reduce the centrifugal and gravitational loads acting upon them (Sieros et al., 2012). Therefore, recent blades are slender and made of lightweight composite materials, increasing their flexibility and potentially causing significant aeroelastic effects. The structural deformation of the blades modifies their aerodynamics by adding bending and twist, which can lead to substantial changes in the mean loads acting on the turbine. Additionally, the wind shear and turbulence create unsteady aeroelastic effects that impact the turbine lifetime. These variations can ultimately affect the flow past the turbine, modifying the wake evolution and recovery. A deeper understanding of the loads and wake of these large turbines is required to design reliable products and predict their interaction in a wind farm, which will allow for further reducing the cost of offshore wind energy.

The structural aspect is often neglected in simulations that use actuator line or disc methods to represent the flow past a turbine. While this approach was mostly valid for smaller turbines, the flow past large rotors is now noticeably affected by the changes in the aerodynamics. To evaluate these effects, some studies have considered the coupling of a flow solver to a structural solver using actuator methods as an interface. Storey et al. (2013) coupled OpenFAST to an actuator disc method to perform simulations of one or two turbines and to compare their results to measurements. Vitsas and Meyers (2016) used a modal structural model coupled to an actuator sector method to study two wind farm configurations. More recently, OpenFAST was also coupled to Nalu-Wind and PALM (Sprague et al., 2020; Krüger et al., 2022) to study the NREL 5 MW reference wind turbine. Hodgson et al. (2021) also performed validation of flexible actuator methods against high-fidelity blade-resolved simulations, highlighting the ability of the actuator line method (ALM) to accurately capture aeroelastic effects on rotors up to 5 MW (see also Sørensen et al., 2015; Hodgson et al., 2022). An elastic ALM based on the coupling between an ALM and an Euler beam model discretized using finite differences was also developed by Meng et al. (2018, 2020) and used for a pair of NREL 5 MW turbines and a wind farm. Spyropoulos et al. (2021) also developed a fluid–structure interaction framework based on the ALM in a flow solver for the compressible URANS equations. This framework uses a multi-body formulation with Timoshenko beam elements to solve the structural dynamics. It was compared in uniform flow against the blade element momentum (BEM) theory and a lifting line and was shown to accurately predict the blade deflections and the loads, also in yawed conditions. Della Posta et al. (2022, 2023) also considered a two-way coupling between a structural solver and an ALM, including the tower, for unsteady aerodynamics of the NREL 5 MW rotor. However, most of the aforementioned studies consider medium-scale rotors with simpler structural parameters than the new-generation turbines. Consequently, the effects of the nonlinearities induced by the large displacements are not considered. Additionally, the impact of the flexibility is rarely assessed by comparing the results to those obtained for rigid blades. Finally, most of the studies focus on the blade and rotor loads but do not assess the effect of the flexibility on the wake in depth.

This study investigates the effects of the deformation of the rotor on a very large reference wind turbine in realistic operating conditions. The open-source IEA 15 MW (Gaertner et al., 2020) turbine model is selected for its large rated power that is representative of modern large-scale offshore turbines and for its comprehensive documentation. Large-eddy simulations of the atmospheric boundary layer and the wake of the turbine are carried out using a flexible ALM coupled to the nonlinear structural solver BeamDyn from OpenFAST. Comparison between fully coupled aeroelastic simulations and simulations assuming rigid blades aims to provide a deeper understanding of the variation in the loads induced by the flexibility. Then, the effect of the blade displacement and load variation on the wake is assessed. The contribution of this study is therefore to provide further insights into the effect of the mean and unsteady deformations on the loads of the 15 MW turbine and to evaluate their impact on the wake.

This paper is structured as follows. The methodology is presented in Sect. 2, which describes the flow solver, the embedded flexible ALM, and the coupling with the structural solver. Then, Sect. 3 presents a comparison between results obtained with the flexible ALM and those obtained using OpenFAST simulations. Performed in both uniform and turbulent flows, this comparison also aims to verify the results obtained with the flexible ALM and assess the relevance of the various aerodynamic models for the simulation of large flexible wind turbines. The case of a wind turbine in realistic atmospheric conditions is then investigated in Sect. 4, where comparison between rigid and flexible rotors is carried out. Finally, the impact of the spatial resolution of the flow is also considered in Sect. 5.

## 2.1 Flow solver

The large-eddy simulation (LES) of turbulent flow is performed by solving the incompressible Navier–Stokes equations supplemented by a subgrid-scale (SGS) model, using an in-house massively parallel code (Duponcheel et al., 2014; Moens et al., 2018). Centred fourth-order finite-difference schemes are used for the spatial discretization. The temporal integration is carried out using a second-order Adams–Bashforth scheme (AB2). The Smagorinsky SGS model (Smagorinsky, 1963) is used, with the asymptotic value of *C*_{S}=0.027 corresponding to a LES where the grid size is much larger than the Kolmogorov scale (Meneveau and Lund, 1997). For a simulation including the ground effect, the wall model developed in Thiry (2017) is used to handle the connection of the LES to the ground.

The turbulent inflow consists of either synthetic turbulent fluctuations generated using the Mann algorithm (Mann, 1998) (used in this study for the comparison to the aerodynamic models of OpenFAST) or a neutral atmospheric boundary layer obtained using a co-simulation (Moens, 2018) (used for the actual study of the loads and wake). In the latter case, a simulation of a half-channel flow is carried out with periodic boundary conditions in the streamwise and lateral directions and slip at the top boundary. The flow is driven by a constant pressure gradient, chosen in combination with the ground roughness to obtain the prescribed hub velocity and turbulence intensity (TI). This simulation runs until statistical convergence is reached. Then, it runs concurrently with the main simulation that contains the wind turbines. A vertical plane of velocity is sampled from the channel simulation and is used as inflow for the main simulation, providing a converged flow equivalent to a neutral atmospheric boundary layer. The influence of the wind turbine on the flow is modelled using an actuator line method described hereunder.

## 2.2 Actuator line method (ALM)

The actuator line method (ALM) represents the effect of the blades on the flow by adding a forcing term to the LES equations (Sorensen and Shen, 2002). Here, the ALM is used with 2D force distribution and integral velocity sampling (Mikkelsen, 2004; Jha and Schmitz, 2018; Churchfield et al., 2017). The blades of the turbine are first discretized into a series of segments, with a control point placed at the centre of each segment, coinciding with the aerodynamic centre of the corresponding airfoil profile. The components of the aerodynamic force acting on each cross-section are computed using the effective flow velocity associated with the control point and the lift and drag coefficients of that airfoil. This effective velocity is obtained in two steps. First, the velocity at the flow solver grid points is linearly interpolated to the nodes of a 2D template centred on the control point. The template is depicted in Fig. 1 and lies in the airfoil plane, which is perpendicular to the actuator line, and is divided into cells of size *h*, equal to the flow grid spacing. Then, the effective velocity is taken as a weighted average of the velocities on the template, using Gaussian weights. Here, the Gaussian weights are obtained using a 2D Gaussian kernel of width *σ*. Its value is taken relatively to the flow grid spacing *h* and is set to $\mathit{\sigma}/h=\mathrm{2}$, similarly to Troldborg (2008). This procedure requires many interpolations, yet its computational cost is limited by grouping the interpolated velocities into a single buffer and performing the Message Passing Interface (MPI) communication once. This velocity sampling is also conservative, as the template used is Gaussian, and the sum of its weights is unity. The distribution of the force is also performed in two steps. First, the force evaluated at the blade control point is distributed on the cells of the 2D template using the same Gaussian weights as those used for the velocity. Then, the force on each cell of the template is redistributed on the closest flow grid points using a linear distribution kernel (as this keeps the distribution as local as possible). It should be noted that the use of integral sampling leads to some smoothing of the high-frequency velocity fluctuations. This effect is discussed in Sect. 3.2. The position of the control points is also updated to follow the time evolution of the deforming configuration.

For aeroelastic simulations, the accuracy of the loads provided by the ALM is essential as they drive the structural dynamics. However, this accuracy is often questioned, especially at the blade tip, mostly due to the regularization of the forces that reduces the induced velocity. To mitigate this issue, it is possible to use either vortex-based smearing corrections (see Meyer Forsting et al., 2019; Dag and Sørensen, 2020; Martínez-Tossas and Meneveau, 2019; Kleine et al., 2023) or regularization kernels that lead to a more accurate prediction of the forces (Jha and Schmitz, 2018; Caprace et al., 2019). The second approach is employed here, using the 2D Gaussian regularization for both the evaluation of the effective velocity and the force distribution. This leads to a better estimation of the aerodynamic loads, hence decreasing the need for a smearing correction, especially when slender blades are considered. A more detailed description of the regularization used is provided in Trigaux et al. (2022), and its advantages over a 3D distribution are detailed in Caprace et al. (2019). Note that the shape of the tip vortices is affected by the type of mollification and that the use of a 2D mollification results in a sharp force gradient that can lead to numerical oscillations in some numerical methods. Here, these were verified to be of small amplitude and not to affect the flow dynamics away from these vortices and the loads on the blades. Here, the spatial resolution of the flow solver near the ALM is taken as 64 pts per turbine diameter D. Section 5 evaluates the influence of the resolution on the loads and the wake in the aeroelastic case.

## 2.3 Structural solver

The structural dynamics of the blades are simulated using the BeamDyn module from OpenFAST (Wang et al., 2017). BeamDyn implements nonlinear geometrically exact beam theory, which is particularly well suited to study slender blades made of composite materials as it accounts for the large displacements and allows for the coupling of degrees of freedom using full 6×6 cross-sectional mass and stiffness matrices. For the considered IEA 15 MW blade, this coupling mainly relates the flapwise and edgewise bending to the twist, which can substantially modify the angle of attack of the blade sections under classical loading conditions. The use of BeamDyn is justified by the importance of the nonlinear effects on large rotors (Manolas et al., 2015; Panteli et al., 2022). The structural data provided in Gaertner et al. (2020) are used. For this turbine, preliminary analysis indicates that the time step must be constrained by d*t*<0.005 s to correctly capture the dynamics of the first torsional mode. Each blade of the turbine is represented by one instance of the structural beam solver, for which the kinematics of the blade root is imposed. The simulations presented in this paper all consider a constant rotation speed, which decouples the blade root boundary condition from the other blades and drivetrain dynamics. The choice of a constant rotation speed is motivated by the need for exact comparison of the wake in different configurations, yet the absence of a model for the drivetrain dynamics impacts the rotor loads, as discussed in Sect. 4.2. Additionally, this approach allows us to use a larger time step than that typically used in OpenFAST simulations. BeamDyn is coupled to the ALM as described in the following section.

## 2.4 Coupling

The coupling scheme linking the flow and structural solvers is the second-order improved serial staggered (ISS) (Degroote, 2013), which is detailed in Fig. 2. This scheme evaluates the flow and the structural models at different time levels. The flow velocity is evaluated at time *t* and is used by the ALM to evaluate the aerodynamic forces based on the previous structural states $x(t-\mathrm{d}t/\mathrm{2})$ (the states include the blade deformation, position, and velocity). The structural response to these loads is then integrated in time, using sub-cycling (with constant forces) to meet the structural solver constraints on the time step. The forces are then re-evaluated by the ALM using the updated structural states and distributed on the flow mesh according to the deformed configuration.

The computational time required for solving the structural dynamics is much lower than that required by the flow solver. BeamDyn is very efficient thanks to the use of high-order spectral finite elements (Sprague and Geers, 2008). When the turbine operates at a constant rotation speed, the time step also remains relatively large, resulting in fewer necessary sub-steps. Additionally, here the structural dynamics of each blade is solved in parallel, leveraging the multiple CPUs dedicated to the flow solver that would otherwise remain in standby until the completion of the structural computation. As a result, the total overhead generated by using the structural solver is less than 5 % compared to simulations performed with rigid blades. More details on the implementation of the coupling between the ALM and the structural solver in parallel are provided in Appendix A.

The frames of reference in which the blade displacements, loads, and wake are presented are also defined according to Fig. 3. The global frame, in which the wake quantities are presented, is located at the tower bottom. The (${\widehat{\mathbf{e}}}_{x}$, ${\widehat{\mathbf{e}}}_{y}$, and ${\widehat{\mathbf{e}}}_{z}$) directions correspond to the streamwise, vertical, and lateral directions, respectively. The loads are presented in the blade root frame ${\widehat{\mathbf{e}}}_{\mathrm{b}}$. Its origin is located at the basis of the blade, and its orientation accounts for the shaft tilt angle, the blade rotation, and the cone angle. The displacements are also measured in the blade root frame and are considered to be the variation from the undeflected configuration. The torsion angle is defined along the ${\widehat{\mathbf{e}}}_{y,\mathrm{b}}$ axis such that a positive value of the angle corresponds to a nose-up rotation of the blade section.

The presented methodology allows us to simulate many rotations of wind turbines with flexible blades in LES at a reasonable computational cost (i.e. a few thousands of CPU hours for the simulations presented in the next section). In what follows, the methodology is compared to various aerodynamic models of OpenFAST before being applied to assess the effect of the flexibility on the IEA 15 MW turbine.

The methodology developed in the previous section is tested and compared against the results obtained using different aerodynamic models available in OpenFAST (v3.5.0). OpenFAST is a publicly available modular framework that uses physics-based models to perform aero-servo-elastic simulations of wind turbines under various conditions (Jonkman, 2013). The presented comparison aims to verify our methodology and point out the limitations of some aerodynamic models for large wind turbines.

Here, the IEA 15 MW reference wind turbine is considered with the following parameters. The blades have a length of 117 m, and the hub radius is 3 m, leading to a total diameter *D*=240 m. To avoid tower strike, the blades are pre-bent by 4 m, the shaft is tilted by an angle of −6°, and the rotor has a cone angle of 4°. The gravity loads are also included. All the parameters used for this model are provided in the turbine definition, except for the structural torsional damping factor, which was increased to avoid torsional instabilities when the turbine operates in near-rated highly loaded conditions.

The following investigation compares the results obtained using the flexible ALM to those obtained using two aerodynamic models from OpenFAST: the blade element momentum (BEM) theory model and a free vortex wake model (cOnvecting LAgrangian Filaments, OLAF) (Jonkman et al., 2015; Shaler et al., 2020). BeamDyn is used for the blade structural dynamics in all the considered cases.

## 3.1 Comparison in uniform flow

First, the case of a uniform inflow is considered. The inflow wind speed is set to *U*_{∞}=9 m s^{−1}, which corresponds to the region of the controller where the optimal tip-speed ratio (TSR) of 9 is maintained. Consequently, the constant turbine rotation speed Ω is set to 6.45 rpm. The ALM and OpenFAST simulations run for 200 s, which corresponds to ≃21.5 rotations, and the results are averaged over the last four rotations. For the ALM simulations, the computational domain size is $\mathrm{12}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{12}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{12}\phantom{\rule{0.125em}{0ex}}D$ (in the streamwise, vertical, and transverse directions, respectively). The turbine is located at the position $(\mathrm{4},\mathrm{6},\mathrm{6}\phantom{\rule{0.125em}{0ex}}D)$ to ensure a sufficient distance from the domain boundaries. This way, the blockage ratio is ≃0.5 %, which allows us to compare the ALM simulations to the OpenFAST aerodynamic models that do not account for the presence of lateral boundaries. A uniform resolution of 64 flow grid points per diameter is used, and the time step is set to ≃0.025 s to obtain 360 time steps per turbine rotation. For the OLAF model, 50 blade nodes are used in the analysis, and the time step is set to 0.155 s, which corresponds to an azimuthal increment Δ*θ*≃6*°* per time step, as suggested in the guidelines (Shaler et al., 2020). The circulation on the blade is obtained using the same airfoil polars as those used in the ALM. The near wake consists of 1024 vortex panels, which gives a wake length corresponding to approximately 17 rotations. There is no far-wake region. The last third of the wake is frozen to mitigate the effect of wake truncation. This wake length was verified to be sufficient in this case to provide a converged force distribution on the blades by comparing it to a simulation with a shorter near wake, using 600 panels extending for 10 rotations, which led to the same results. The vortex core is regularized using the Vatistas vortex model ^{1} with the optimized parameter option, which corresponds to using *r*_{c}=2Δ*s*, where Δ*s* is the mean spanwise spacing between the aerodynamic nodes. The size of the regularization parameter can noticeably affect the force distribution, especially near the tip. The regularization also evolves with time using a viscous diffusion model (here with the parameter *δ*=1000). For the BEM, 50 radial locations are used. Prandtl's tip-loss factor is used to model the tip and the hub losses (Moriarty and Hansen, 2005). The tangential induction is accounted for in the BEM equations, and the effect of the drag term is considered in both the axial and the tangential induction. The OpenFAST cases have a very small time step of 0.0005 s, which is necessary to ensure the convergence of BeamDyn. Interested readers can refer to the IEA 15 MW definition that contains example cases with OLAF and the BEM (Gaertner et al., 2020).

A comparison of the rotor thrust and power coefficients predicted by each method is provided in Table 1. For the rigid case, the thrust and power coefficients of the ALM and OLAF match well, whereas the BEM tends to predict a lower power coefficient. For the flexible case, the tip displacement is predicted similarly by all methods. These deflections are particularly significant in the flapwise direction, reaching a value of 12 m at the blade tip, which corresponds to around 10 % of the blade length. The torsion angle reaches −2.8° at the tip (corresponding to a nose-down rotation of the airfoil), which also significantly reduces the angle of attack of the blade section, resulting in a noticeable variation in the thrust and power coefficients when accounting for the blade flexibility. Specifically, all methods predict a 9 % decrease in the thrust coefficient. The variation in the power coefficient is closer to 5 % for the ALM and OLAF, whereas the BEM predicts a smaller reduction of 2 %. As a result, the difference in power coefficients between the methods is lower when flexibility is considered.

The aerodynamic forces acting on the rigid and flexible blade are displayed in Fig. 4. These forces are expressed in the blade root frame. For the rigid case, the forces obtained using OLAF and the ALM match very well. The only noticeable difference between the two methods is the small increase in the forces near the blade tip, which is more pronounced using OLAF. The latter is due to the regularization of the trailing vortex filaments, which decreases the induced velocity near the tip. To a smaller extent, the ALM also depicts a slight increase in the lateral forces (i.e. *F*_{z}, in the blade root frame) at the tip due to the mollification (i.e. the smearing of the shed vortices, Caprace et al., 2019). It is however less pronounced thanks to the use of the 2D mollification, without spanwise smearing. The force distribution predicted by the BEM is lower than for OLAF and the ALM, especially in the edgewise direction. This difference was also observed in other studies (Madsen et al., 2012; Perez-Becker et al., 2020) and can be explained by two factors. High thrust coefficients are somewhat challenging for the BEM as the limits of the validity of the theoretical model are reached, and the empirical Glauert correction is used. On the other hand, the mollification of the ALM and the regularization of the vortex core in OLAF also affect the induction, which can result in slightly higher forces (Caprace et al., 2019; Shaler et al., 2023).

When flexibility is considered, all methods predict a significant decrease in the loads on the outer part of the blade ($r/R>\mathrm{0.6}$). This decrease is more significant for OLAF and the ALM than for the BEM, especially in the edgewise direction. This difference arises from the fact that the BEM predicts less change in the axial induction due to the out-of-plane bending. Conversely, the effect of the torsion angle is captured similarly by all methods, which results in a similar load distribution in the flapwise direction. Additionally, on the inner part of the blade ($r/R<\mathrm{0.5}$), OLAF and the ALM show an increase in the loads. This phenomenon arises from the additional flapwise deformation of the blade, which modifies the location of the emission of the tip vortices, and was also observed to a smaller extent for the NREL 5 MW turbine (Dose et al., 2018; Trigaux et al., 2022). It is not captured by the BEM as the effect of the blade curvature on its aerodynamics is not modelled here (Fritz et al., 2022; Li et al., 2022). These findings are also consistent with the results predicted by the BEM and blade-resolved URANS simulation for the DTU 10 MW turbine presented in Sayed et al. (2019). This highlights the importance of using higher-fidelity methods when the aeroelasticity of large rotors is considered or at least considering the effect of out-of-plane bending on the blade aerodynamics (Fritz et al., 2022; Li et al., 2022), also for steady loads such as those presented in this section.

## 3.2 Comparison in turbulent flow

To verify the dynamic response of the structure under unsteady loads, a comparison between our framework and OpenFAST is also carried out with a turbulent inflow. The latter consists of the superposition of a uniform flow with *U*=9 m s^{−1} and a turbulent field generated using the Mann algorithm (Mann, 1998). The parameters of the turbulent field are the integral length scale, here set to *L*=88.5 m; the anisotropy factor Γ=3.9; and the turbulence intensity, set to TI =6 %. The turbulent fluctuations are generated in a 3D domain of size $\mathrm{8}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{2}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{2}\phantom{\rule{0.125em}{0ex}}D$, with a spatial discretization of 32 points per D. In OpenFAST, the turbulent fluctuations are directly interpolated from the 3D domain to the turbine blade nodes and added to the upstream velocity value. In the case of the ALM, the turbulent fluctuations are interpolated to the inflow plane and are then convected to the turbine location. The OpenFAST cases are kept identical to the cases with a uniform inflow. Only the time step at which OLAF emits a new panel is decreased to 0.1 s to better account for the turbulent fluctuations. The ALM case also remains identical, except that the distance between the turbine and the inflow plane is reduced to 2 *D* to reduce the variation in the turbulent fluctuations that occurs during their convection from the inflow to the turbine. The cases are run for 400 s to converge the statistics.

The dynamic structural deformation of the blade is considered for this comparison to assess whether the flexible ALM correctly captures the variations in the displacements at the relevant frequencies. Figure 5 depicts the power spectral density (PSD) of the tip displacement for the flapwise and edgewise translations and for the torsion angle. The PSDs are computed using the Welch algorithm with segments of duration equal to 50 s that overlap by half their length (Welch, 1967). The blade natural frequencies are also depicted. These were obtained using the mass and stiffness matrices provided by BeamDyn in the rotating undeformed configuration, without accounting for the coupling with the aerodynamic loads. It should be noted that, for a nonlinear beam solver such as BeamDyn, the blade natural frequencies are slightly affected by the deformation. Here, this difference was verified to be negligible such that the natural frequencies do not need to be updated using the blade deformed configuration. Very good agreement is observed between OpenFAST and the ALM. All the methods predict the same magnitude for the peak occurring at the 1P frequency for all three components. Additionally, the peaks of displacement occurring at the edgewise natural frequencies (E1 and E2) are present and reach the same value for all three solvers. The structural response close to the first torsional frequency (T1) is also similar for all methods. All three solvers predict a similar decrease in the frequency at which this peak arises compared to the T1 frequency. This shift is due to the coupling between the aerodynamic forces (that increase linearly with the angle of attack) and the torsion angle. The amplitude of the peak is however slightly smaller for the ALM. Indeed, the values of the PSD at all the higher frequencies (larger than 20P) obtained using the ALM are smaller than those predicted using OpenFAST. This is mostly related to the use of integral velocity sampling for the ALM, which leads to some smoothing of the fluctuations as the velocity is obtained as a weighted average. In contrast, OpenFAST interpolates the velocity fluctuations directly to the nodes of the aerodynamic solver (BEM or OLAF). However, the amplitude of these differences is limited, and the value of the PSD at the affected frequencies is much smaller than at the lower frequencies, which are predominant. Some additional differences in the lower-frequency part of the spectrum are related to changes in the turbulent field during its convection from the inflow to the turbine location.

In this section, rigid and flexible rotors are studied with sheared turbulent wind as inflow to evaluate the unsteady aeroelastic effects on the aerodynamic loads and the wake. To this end, the numerical setup consists of a domain size of $\mathrm{12}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{4}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{4}\phantom{\rule{0.125em}{0ex}}D$, with a spatial resolution of 64 grid points per *D*. The location of the turbine hub is $(\mathrm{2}D,{y}_{\mathrm{hub}},\mathrm{2}D)$, where *y*_{hub}=150 m. Although the distance between the turbine and the inflow is small, it has been verified to have a very small effect on the loads (approximately 1 % difference in thrust and 3 % difference in power when the distance to the inflow is doubled). Moreover, the same distance to the inflow is used for each considered case such that the comparison is not affected. To ensure a valid comparison between the rigid and flexible cases, the rotation speed at the blade root is set to Ω=0.75 rad s^{−1}, which corresponds to the design TSR of 9. The time increment of the flow simulation is set to obtain 360 time steps per turbine rotation. Additionally, there are four structural sub-steps per flow time step. A total of 200 rotations are simulated for each case, which corresponds to a physical time of 1675 s.

The inflow is a neutral ABL generated using the co-simulation technique (as described in Sect. 2.1). The numerical setup is described in detail in Appendix B. It was designed to obtain a time-averaged velocity of *U*_{hub}=10 m s^{−1} at hub height and with a TI close to 5 %. The obtained profiles of mean velocity and turbulence intensity are also provided in the Appendix, together with the methodology used to estimate the integral length scales, obtained as 106 m in the streamwise direction and 89.5 m in the lateral direction.

Three cases are compared in this section. The first case is the “rigid undeformed” rotor, which includes the upwind pre-bend of the blades. The second case is the flexible rotor, modelled using the flexible ALM that dynamically deforms based on the aerodynamic loads. The third case considers a rigid rotor that includes the mean deformation of the blades, noted as the “rigid deformed” case. The deformed geometry was obtained by time averaging the structural displacements obtained in the second case. This allows us to isolate the effect of the time-varying response of the blade on the rotor dynamics and wake.

This analysis is structured in three parts. First, the displacements of the blades are considered. Then, the distribution of the loads and their variation are presented. Understanding these aspects allows us to explain the differences observed in the wake, which are presented in the last part of the analysis.

## 4.1 Structural displacements

The PSD of the displacements is first considered in Fig. 6. In this section, the presented statistics of the blade displacements and loads are computed using the last 180 rotations of the simulations. The first 20 rotations are discarded as the wake induction has not yet converged. The power, thrust, and power spectra of the loads were verified to have reached statistical convergence after 20 rotations. The presented PSDs are computed using the Welch algorithm. The time series are divided into sub-segments of length equal to 22.5 rotations that overlap by half their length, resulting in a total of 15 segments. The PSDs are also computed for each blade and then averaged, resulting in smooth spectra. For all the components of the displacement, the 1P frequency is the most noticeable, as also noted for smaller turbines in Grinderslev et al. (2021). In the flapwise direction, the 1P frequency is excited by the rotation of the blade in the mean shear and in the spatially coherent turbulent structures. For the edgewise direction, the excitation of the 1P frequency is mostly attributed to gravity. In this direction, the spectrum reaches a much higher value at 1P than at the other frequencies, indicating the dominance of gravity over the aerodynamic loads. In the torsional direction, the peak at 1P is due to the coupling with the translational degrees of freedom and the effect of gravity on the deformed blade. The harmonics of the rotation frequency are also significant in all spectra, which is consistent with the results of Krüger et al. (2022). Concerning the blade structural frequencies, their excitation depends on the considered component. In the flapwise direction, the PSD exhibits no particular excitation at the flapwise (F1, F2, or F3) frequencies. Only a small peak is visible near the torsional frequency. The absence of a peak is explained by the high aerodynamic damping of the flapwise mode (Hansen, 2007). This damping is due to the fact that the blade flapwise displacement decreases the apparent streamwise velocity on the blade section, hence decreasing the angle of attack and the loads. In contrast, in the edgewise direction, distinguishable peaks of displacement occur at the edgewise natural frequencies (E1, E2) and close to the torsional frequency (T1). This is a result of the small aerodynamic damping in this direction, yet the amplitude of the peaks remains small compared to the one occurring at 1P. For the torsional direction, a peak occurs slightly below the T1 frequency, as observed in Sect. 3.2.

It is also interesting to consider the general shape of the displacement spectra. In the flapwise direction, the PSD follows a general downward slope that is constant throughout most of the spectrum. In the edgewise direction, a similar trend is noticeable before the 1P frequency and between the E1 and T1 frequencies. However, for the first harmonics of 1P, the general trend of the PSD remains constant. For the torsional direction, the general shape of the spectrum remains constant across the frequencies. Consequently, in this direction, some non-negligible high-frequency modes are observed in addition to the rotation frequency.

As the rotational frequency appears to be predominant, it is interesting to investigate the variation in the displacements along one rotation. The mean displacement and standard deviation are represented in Fig. 7 as a function of the azimuthal angle (with 0 representing the blade pointing upward). The values measured at each angle and over the simulated turbine rotations are also plotted; for each angle, they can be viewed as an approximation of the probability density function (PDF). In the flapwise direction, the minimal blade displacement is reached a few degrees after the blade's downward position. In this region, the PDF also has the smallest spread. Interestingly, the maximal mean flapwise displacement is not reached at the blade's upward position (i.e. 0°), where the flow velocity is maximal, but is rather constant from 0 to 90°. A similar behaviour was also observed for a 5 MW turbine by Santo et al. (2019). The PDF also exhibits a larger spread when the blade is in the upper-half region (i.e. between 270 and 90°). Overall, the flapwise deformation varies significantly along the rotation, oscillating between 12.4 and 13.6 m, primarily due to the mean shear of the flow. This displacement also represents more than 10 % of the rotor radius, which is quite significant, and illustrates the increased blade flexibility compared to smaller turbines. The PDF also reveals large variations around the mean displacement, attributable to the turbulent fluctuations. In contrast, the edgewise displacement presents a nearly perfect sinusoidal variation that is consistent with the gravity loads (Höning et al., 2024). The spread around the mean is also very limited, as the variation in the aerodynamic loads is small compared to that of the gravity. The torsional deformation follows the same phase as the edgewise displacement. This is due to the effect of the gravity loads and the flapwise deformation, which generates a torsion moment on the structural reference axis. Consequently, the torsion angle differs between the left and right parts of the rotor, which can lead to some additional load imbalance. The PDF also presents a larger width than for the edgewise displacement as a result of the higher sensitivity to the turbulent fluctuations.

The flapwise structural velocity is also considered as it substantially affects the angle of attack and hence the aerodynamic forces. In fact, the angle of attack *α* of a section of the blade can be approximated using the flow and structural velocities expressed in the blade root frame by neglecting the effect of the blade curvature on the angle of attack,

where *v*_{x} and *v*_{z} are the local flow velocities in the *x* (streamwise) and *z* (azimuthal) directions, *v*_{struct,x} and *v*_{struct,z} are the local structural velocities in these directions, *ϕ*_{struct,y} is the local torsion angle, and *β* is the local twist angle. Here, all these variables are expressed in the blade root frame. In general, the structural velocity in the *z* direction mostly consists of the rotational velocity Ω*r*, and the flow velocity is also negligible in comparison, leading to $\left({v}_{\mathrm{struct},z}-{v}_{z}\right)\simeq \mathrm{\Omega}r$. Additionally, since the normal velocity is also much smaller than the rotational velocity, $\mathrm{arctan}(\left({v}_{x}-{v}_{\mathrm{struct},x}\right)/(\mathrm{\Omega}r\left)\right)\simeq \left({v}_{x}-{v}_{\mathrm{struct},x}\right)/\left(\mathrm{\Omega}r\right)$. As a result, the expression for the angle of attack is simplified to

Consequently, the ratio ${v}_{\mathrm{struct},x}/\left(\mathrm{\Omega}r\right)$ has an effect on the angle of attack that is similar to that of the torsion angle. Figure 8 shows the mean flapwise velocity *v*_{struct,x} has a function of the blade azimuthal angle and the probability density function for each angle. The maximal positive velocity is 0.4 m s^{−1} and is reached in the left region of the rotor. The minimal velocity is −0.7 m s^{−1} and is reached when the blade is pointing downward. These values are more significant than those reported in Santo et al. (2019) for a 5 MW turbine in atmospheric conditions, which illustrates that the maximal flapwise velocity reached along the rotation likely increases with the rotor size. A large variation in the velocity occurs around the blade's downward position, which is attributed to the high mean vertical shear in the lower region of the rotor. The PDF of the flapwise velocity is also quite large around the mean, due to the turbulent fluctuations of the velocity. The figure also shows the contribution of the flapwise velocity ($-{v}_{\mathrm{struct},x}/\left(\mathrm{\Omega}r\right)$) to the variation in the angle of attack at the tip along a rotation. It is compared to the contribution of the torsional angle (${\mathit{\varphi}}_{\mathrm{struct},y}-{\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}}_{\mathrm{struct},y}$). The two effects are similar in magnitude and can result in a small variation in the angle of attack of up to 0.4°. Their phase is shifted by a quarter rotation, as the effect on the flapwise velocity mostly depends on the vertical shear, whereas the variation in the torsion follows the same phase as the edgewise displacements (and hence the gravity loads).

These displacements are not without consequences for the loads acting on the turbine, which are presented in the next section.

## 4.2 Comparison of the blade loads

This section presents a comparison of the loads obtained with the flexible and rigid rotors. Here, the focus is on the blade loads, presenting both the aerodynamic loads due to their impact on the wake and the structural loads due to their effect on the structural integrity. This choice also follows from the decision to maintain a constant rotation speed, primarily to facilitate the comparison between the wakes. In real-world conditions, the turbine controller adjusts the rotation speed and the blade pitch angle, and the three blades are dynamically coupled by their rigid connection with the turbine shaft, which modifies the resulting rotor loads, especially the rotor torque. To assess the impact of these changes, the simulations were also conducted with the controller. The loads and displacements of the individual blades were not noticeably affected and neither was the wake. However, the spectrum of the rotor torque was modified due to the coupling between the blades and the generator. Consequently, we chose to present the simulation at a constant rotation speed and to focus on the individual blade dynamics and the wake, which were confirmed to remain mostly unaffected by the controller.

First, the mean load distribution along the blade, depicted for each case in Fig. 9, is examined. The load distribution on the undeformed rotor clearly differs from that of the two cases that include the deformation. In particular, a significant reduction in the mean loads and their standard deviation is induced by the structural displacement on the outer part of the blade. This results in a substantial reduction in the power and thrust coefficients. The thrust coefficient is reduced by ≃ 14 %, going from *C*_{T}=0.80 for the undeformed rotor to *C*_{T}=0.69 for the deformed cases. Similarly, the power coefficient is reduced by 10 %, from *C*_{P}=0.50 to *C*_{P}=0.45. As anticipated, the rigid deformed and flexible cases exhibit a nearly identical mean force distribution. However, in the flexible case, the standard deviation presents a slight reduction, attributable to the additional compliance of the blade.

Although the radial distribution of the forces is the same for the rigid deformed and flexible cases, their distribution over the rotor plane varies. Figure 10 depicts the difference between the force distribution of the two cases. There is a notable difference at the bottom of the rotor in both the flapwise and the edgewise directions. In this region, the flexible rotor experiences more force due to the reduced flapwise deflection and the significant upstream blade velocity that increases the angle of attack, as observed in Fig. 8. In the normal direction, there is also a clear asymmetry between the left and right planes of the rotor. The flexible rotor presents higher forces compared to the rigid rotors on the left part of the plane due to the smaller torsion angle on this side. This asymmetry is not visible in the edgewise direction. The observed differences reach a maximum of around 10 % on the lower part of the rotor plane, which also results in a difference in the mean velocity deficit in the near wake, as discussed in the next section.

Figure 11 presents the PSD of the blade root moment in the two main directions. The spectra are given for the aerodynamic forces and the elastic structural response (which includes the inertia, centrifugal force, and gravitational force). Considering first the aerodynamic forces, one observes that the rotation frequency is still the most important component. The amplitude of the PSD at this frequency differs slightly between the cases, with the rigid undeformed case having the highest value and the flexible case the lowest. Large differences also appear close to the natural frequencies. The PSD is much lower close to the main bending frequencies (F1, E1) due to the dynamic response of the blade. In fact, at these frequencies, the blade responds significantly to the variations in the forces, and the resulting structural bending tends to smooth these changes. There is also a noticeable increase at the torsional frequency (close to T1) due to the subsequent modification of the angle of attack of the blade. The structural moment does not exhibit a decrease at the bending frequency, as the lower aerodynamic forces are compensated for by the inertial loads. In the flapwise direction, the spectrum of this moment is generally slightly smaller in the flexible case. There are also no specific peaks at the natural frequencies due to the high aerodynamic damping in this direction. In contrast, the edgewise root moment presents high peaks close to its natural frequencies. In this direction, the value of the PSD at the 1P frequency is also much higher than that of the aerodynamic loads, illustrating the importance of the gravity loads. The modification of the spectrum of the root moment can have a direct impact on the lifetime of the turbine and the damage equivalent load calculations, as also noted in Hodgson et al. (2021). It is therefore necessary to consider dynamic models when assessing the structural integrity of the rotor.

Similarly to the displacement, the variation in the blade root moment of the aerodynamic forces along the rotation is further investigated, as the 1P frequency was shown to be the most significant. The latter is represented in Fig. 12. The aerodynamic forces, rather than the structural response, are considered to remove the influence of gravity, which would be dominant on the edgewise moment.

The mean flapwise moment acting on the rigid turbine is 21 % higher than that obtained in the cases that include the deformation, which have a comparable mean value. The amplitude of the deviation along the rotation, defined as $({M}_{\mathrm{max}}-{M}_{\mathrm{min}})/{M}_{\mathrm{mean}}$, is however similar in both rigid cases (14.0 % for the rigid undeformed case and 13.2 % for the rigid deformed case). This deviation is reduced to 9.3 % in the flexible case, mostly due to a higher minimal value when the blade is in its downward position. In the edgewise direction, the mean value of the root moment is 12 % higher in the rigid undeformed case. The deviation is 31.5 % in this case, compared to 23.7 % in the rigid deformed case and 17.5 % in the flexible case.

Interestingly, there is also a noticeable phase lag between the rigid and the flexible cases. While the rigid case predicts a minimal moment at around *θ*=195*°*, the flexible case predicts a minimum at *θ*=215*°*. This phase lag also applies to the rotor torque and ultimately to the power. It could also have an effect on load alleviation control strategies that use the root bending moment as input, such as individual pitch control.

## 4.3 Impact on the turbine wake

In this section, the impact of the flexibility on the wake is considered. The wake can be affected in three ways. Firstly, the changes in the thrust coefficient *C*_{T} affect the total velocity deficit behind the turbine. Secondly, the variation in the load distribution along the blade span, described in Sect. 4.2, changes the shape of the velocity deficit. If the loads decrease more smoothly near the tip, the strength of the tip vortices also decreases. Finally, the displacements of the blade, described in Sect. 4.1, affect the location of the emission of the tip vortices. Specifically, the tip vortices are generated further downstream due to the flapwise bending, and the flexibility affects the distance between individual tip vortices, potentially modifying their interaction. This section aims to understand whether these effects have a significant impact on the wake, its stability, and the resulting velocity deficit.

First, the wake statistics are considered. The statistics are computed over the last 150 rotations to ensure that the wake is developed up to the end of the computational domain before averaging. The discarded first 50 rotations correspond to approximately 1.5 flow-through times, which was verified to be sufficient to reach a statistically converged flow. The averaging window of 150 rotations then corresponds to approximately 4.4 flow-through times. Figure 13 illustrates the time-averaged velocity deficit and turbulent kinetic energy (TKE) profiles. The differences between the profiles of the mean velocity deficit are comparable to the changes in the mean load distribution, which are the main source of differences between the wakes. Behind the tip of the turbine, the deficit is more pronounced in the rigid undeformed case compared to the deformed and flexible cases. Additionally, there is a steeper gradient of mean velocity behind the tip in this case. The TKE is also noticeably higher in the near-tip region, which results in a faster turbulent diffusion of the wake. For the rigid undeformed case, the velocity deficit observed 3 and 4 *D* behind the turbine is smaller but more diffuse than in the deformed cases. In the far wake (6 *D* behind the rotor), the velocity deficit and TKE profiles are almost equivalent. The time-averaged velocity deficit and TKE are very similar between the rigid deformed and flexible cases, which indicates that the unsteady variations around the mean deformation have a limited effect on the time-averaged wake statistics.

Some differences between the near wake of the rigid deformed and flexible cases can nevertheless still be observed. Figure 14 presents the difference in the mean velocity and TKE over vertical planes behind the rotor. The difference in the mean velocity behind the rotor ($x/D=\mathrm{0.5}$) is very close to that observed for the blade normal forces on the rotor area (see Fig. 10). The small changes in the local thrust result in differences in the mean velocity deficit that reach up to 2.5 % of the upstream velocity. The most substantial difference occurs behind the bottom of the rotor, where the velocity is significantly higher behind the rigid deformed turbine than behind the flexible turbine. Noticeable differences also exist behind the left and right parts of the rotor. At a distance $x/D=\mathrm{2}$, these differences start to vanish due to the wake destabilization that increases the turbulent mixing. The TKE exhibits a different behaviour, as it is globally lower in the rigid deformed case compared to the flexible case. The largest differences are observed for *θ*=90*°* and 270*°*, which also correspond to the angle with the highest variance of flapwise and torsional displacements. The increase in the TKE in the flexible case is attributed to the unsteady variations in the loads and position that further contribute to the fluctuations in the flow. However, these fluctuations remain moderate compared to the incoming turbulence of the wind and therefore have a limited impact on the wake behaviour.

A snapshot of an instantaneous out-of-plane vorticity field behind the rotor is depicted in Fig. 15 on a vertical plane. We stress that the turbulent inflow (provided by the co-simulation of the ABL) is at exactly the same time *t* for the three cases (i.e. it is synchronized), which allows us to compare fine details in the wake. The magnitude of the tip vortices behind the rigid rotors is significantly higher for the rigid undeformed rotor. This results from the higher gradient of the forces observed in the near-tip region for the undeformed blades. Additionally, these tip vortices are emitted further upstream compared to those of the deformed cases. As a result, the destabilization of the wake occurs at a smaller distance from the rotor. The rigid deformed and the flexible cases exhibit very similar vortical structures, even in the finer details. The wakes are almost identical, which further confirms the minimal influence of the unsteady fluctuations in the blade deformations on their mean deformation.

The streamwise velocity is measured at various locations in the wake for each case to compare the flow statistics. Four locations are selected, as depicted in Fig. 15. The power spectral density of the streamwise velocity (obtained using the Welch algorithm) is then depicted for each case in Fig. 16. The first location is 1 *D* upstream of the hub in the induction zone (Probe 1), and the three PSDs are virtually identical, which confirms that the differences in the downstream flow can be fully attributed to the changes in the turbine aerodynamics. The flow velocity is also sampled behind the blade tip at two positions. In the very near wake (0.5 *D* behind the rotor, Probe 2), the blade passing frequency is clearly noticeable for the deformed and flexible cases, whereas it is not visible in the rigid undeformed case. This can be related to the fact that the tip vortices are well separated at this location in the deformed and flexible cases, while they already merged into a vortex sheet in the rigid undeformed case. Additional differences can be observed further in the wake (1 *D* behind the rotor, Probe 3), where the rigid undeformed case presents a much higher PSD at the frequencies lower than 1P. This increase is due to the destabilization of the vortex sheet, which results in lower-frequency velocity variations. In contrast, the rigid deformed and flexible cases do not exhibit such an increase, indicating that the wake remains stable at this location. Finally, the velocity is measured further downstream, 3 *D* behind the hub (Probe 4), where the influence of the hub vortices is minimal. At this location, the three PSDs are similar. For each location, the differences between the rigid deformed and the flexible case are small, confirming the limited impact of the unsteady variations on the wake.

The wake recovery rate is measured for each case by computing the power available in the flow *P*_{w} on the rotor-swept area *A*_{rot}. The latter is obtained as a function of the streamwise coordinate *x* as

The variation in the power (normalized by the power *P*_{0} at the inflow) is depicted in Fig. 17. The available power in the wake of the rigid undeformed case is lower in the near wake due to the higher thrust coefficient. It also presents a faster recovery rate since it reaches the same level as the two other cases for $\mathrm{2}<x/D<\mathrm{6}$. However, in the far wake, it is also slightly lower than in the rigid deformed and flexible cases. These two cases again present almost identical statistics.

This analysis shows the necessity of including the mean deformation to accurately capture the near wake but also predicts a minimal effect of the unsteady variations of the displacements and loads around their mean. The spatial resolution of the presented simulations is 64 points per D, corresponding to 3.75 m, which is larger than the variation in the displacements along the mean deformed configuration. To further assess the impact of the small variations, it is important to verify that the resolution of the flow simulations is sufficient. The next section therefore assesses the effect of the spatial resolution on the loads and near-wake dynamics.

In this section, the results obtained with a spatial resolution of 64 grid points per diameter are compared to those of a higher-resolution simulation to assess whether the small deviations in loads and displacement around the mean are sufficiently resolved. The case of a flexible rotor is therefore reconsidered with a grid resolution twice as fine, corresponding to 128 grid points per *D*. The time step is also divided by 2 to maintain the Courant–Friedrichs–Lewy (CFL) number constant. The mollification of the ALM forces is kept constant relative to the grid size (i.e. $\mathit{\sigma}/h=\mathrm{2}$) and is thus 2 times smaller for the highest resolution. The resolution of the precursor simulation is similarly increased, and the smaller scales of the turbulence are obtained by restarting the ABL simulation at the new resolution and converging it during ≃25 flow-through times. The simulation with the resolution of 64 points per D is also performed using the same precursor simulation, which is interpolated to the coarser resolution at the inflow plane. This ensures a valid comparison between the two simulations.

The mean force distribution on the blades is depicted for both resolutions in Fig. 18. The two distributions and their variances are very close. The thrust coefficient *C*_{T} is 0.687 for the simulation with 64 points per D, and it increases slightly to 0.692 for 128 points per D (+0.7 %). Similarly, the power coefficient increases from 0.439 to 0.445 (+1.3 %). The general shape of the distribution remains consistent between the two simulations, also near the tip, which suggests a sufficient discretization of the blade and the wake.

The variations in the aerodynamic forces and structural root moment are also depicted in Fig. 19. The presented PSDs match perfectly until the E1 frequency. For higher frequencies, the loads predicted using 128 points per D are higher. This is due to the additional unsteadiness induced by the smaller turbulent scales that are resolved up to a higher frequency when at high resolution and the smaller mollification size that reduces the smoothing induced by the integral velocity sampling. The magnitude of the PSD at these high frequencies is however much lower than that of the small frequencies; hence it is unlikely that this additional part of the spectrum significantly affects the fatigue of the blades and the wake. The peaks of the PSD located at the E2 and T1 frequencies are also obtained in both cases, indicating that the main components of the structural response are not significantly impacted by the resolution.

Figure 20 presents the instantaneous vorticity field behind the rotor with the two resolutions. The same velocity fields obtained at a resolution of 128 points per D and at the same time *t* are used at the inflow of both simulations. The only difference is that the inflow is interpolated from 128 to 64 points per D at the inflow plane of the coarser simulation. At higher resolution, the tip vortices have a smaller core size and hence a higher vorticity at their centre (since they essentially have the same circulation in both cases). The hub vortices are also better resolved. The shed vortex sheet behind the outer part of the blade is also distinguishable. However, the larger structures of the wake are similar in both cases. This results in virtually identical velocity deficits as depicted in Fig. 21. The TKE has a higher value behind the hub and in the near-tip regions, which results from the sharper mollification and the smaller scales captured by the simulation. However, this does not affect the mean velocity deficit.

The PSD of the streamwise velocity at various locations in the domain is also compared for both resolutions. The locations where the flow velocity is sampled are depicted in Fig. 20, and the obtained spectra are shown in Fig. 22. The first spectra are obtained using the streamwise velocity measured 1*D* upstream of the hub (Probe 1). At this location, the PSD obtained for the simulation with 64 points per D decreases much faster at higher frequencies, as expected from the smaller LES cut-off frequency. There are also some slight differences below the cut-off frequency due to the interpolation at the inflow from the ABL with 128 points per D resolution to the LES with 64 points per D resolution and the variation occurring during the convection of the turbulence between the inflow and the sampling location. In the near wake, the PSD of the velocity measured 0.5 *D* behind the blade tip (Probe 2) is very similar between the two resolutions (except at the highest frequencies). When measuring 1 *D* downstream of the blade tip (Probe 3), the peak of the PSD obtained for the simulation with 128 points per D is more pronounced near the 3P frequency due to the more distinct tip vortices, whereas the tip vortices have merged into a vortex sheet in the simulation with 64 points per D. Further in the wake, 3 *D* behind the hub (Probe 4), the spectra are in good agreement at the lower frequencies and only show deviations at the higher frequencies.

This work presented an investigation of the effect of the flexibility of the blades on the loads and the wake for the large IEA 15 MW reference wind turbine. A flexible actuator line method (ALM), coupled to the nonlinear structural solver BeamDyn, is used in large-eddy simulations (LES) to compute the unsteady aerodynamic forces exerted by the turbulent inflow on each blade and their influence on the flow. This methodology leverages the LES capabilities to realistically simulate the unsteadiness of the inflow and the wake produced by the turbine in order to quantify the effect of the structural blade deformations on the loads and the wake.

The results obtained using the flexible actuator line method are first compared to those obtained using OpenFAST. In uniform flow, the free vortex wake model of OpenFAST and the present flexible ALM predict very similar distributions for the forces, which are slightly higher than those predicted using a BEM in the tangential direction. The effect of the deformation of the blades on the time-averaged loads is also predicted similarly by OLAF and the ALM. The magnitude of this effect is larger than that predicted by the BEM, which indicates that higher-fidelity aerodynamic models are necessary to account for the effect of large displacements on the blade aerodynamics. Comparisons carried out using a turbulent inflow without mean shear (generated using the Mann algorithm) show that the power spectral densities (PSDs) of the blade displacements are very similar for all three methods, which confirms that the flexible actuator line model is able to correctly capture the structural dynamics.

The displacement and loads of the blades are then investigated in a turbulent atmospheric boundary layer (generated using a co-simulation). Notable flapwise and torsional displacements are observed, which significantly modify the aerodynamics of the blades. The resulting thrust coefficient decreases by 14 % and the power coefficient by 10 % compared to the rigid undeformed case. Including the mean deformation while keeping the rotor rigid allows us to correctly recover the mean load distribution along the blade, which leads to thrust and power coefficients equivalent to those obtained using the flexible actuator line method. However, the variation in the loads during the blade rotation is different, especially close to the blade's downward position, due to the high mean shear of the flow. Moreover, the aerodynamic loads and the structural response at the blade natural frequencies are noticeably affected by the flexibility. We also show that the variation in the angle of attack due to the flapwise blade velocity is significant and of the same order of magnitude as that due to the blade torsion.

The wake behind the turbine is then considered. The results show that it is essential to include the mean deformation of the blade to obtain the correct mean velocity deficit and profile of turbulent kinetic energy. Once the mean deformation is added, the differences between the rigid deformed and the flexible cases are limited. In the near wake, the time-averaged velocity obtained in the flexible case is up to 4 % slower at the bottom of the rotor-swept area than in the rigid deformed case. The turbulent kinetic energy is slightly higher over the rotor area, which results from the unsteady variations in the loads and displacements. However, these changes do not significantly affect the global wake behaviour.

Finally, the effect of the LES flow solver resolution on the wake is investigated by increasing the spatial discretization to 128 grid points per diameter. The mean load distribution is comparable to that obtained using the resolution with 64 grid points per diameter. The PSD of the root moment is also identical up to the first edgewise frequency. For the higher frequencies, a difference in magnitude is observed due to the higher-frequency fluctuations captured by the LES (due to the twice higher cut-off wavenumber of the grid). A snapshot of the wake vorticity field shows tip vortices with a smaller core size, and hence a higher vorticity at their centre, and smaller turbulent scales in the flow. However, the time-averaged velocity deficits are largely similar between the two simulations. The TKE is higher behind the blade root and tip at high resolution, yet this does not impact the turbulent wake diffusion rate. As a result, it is unlikely that a higher resolution would significantly increase the impact of the unsteady displacements and load variations on the wake. It would however be necessary if the spectrum of the loads must be obtained at high frequencies.

The impact of the structural deformations on the loads and wake of a large 15 MW turbine is thus significant. These effects consist of a reduction in the thrust and power coefficients, a variation in the load distribution over the blade, and a modification of the spectrum of the blade root moments. These changes, in turn, substantially affect the near wake by altering the profiles of the velocity deficit and turbulent kinetic energy. In the relatively steady operating conditions considered in this paper (i.e. constant time-averaged wind velocity, turbulence intensity, and rotation speed), the proper mean load distribution and wake behaviour can be obtained using rigid blades but only when adjusting their geometry to incorporate their mean structural deformation. Here, the latter was obtained by time averaging the structural deformation produced by our fully coupled unsteady simulations, but it could also be obtained in a preprocessing step, using alternative aeroelastic simulation tools. This then removes the need for the unsteady coupling of the structural solver and the actuator line of the LES solver, thereby reducing the computational cost slightly. However, the unsteady variations in the mean would still differ compared to those of the fully coupled simulation, resulting in a different spectrum of the loads. It is hence recommended to include the structural dynamics of the blades if the unsteady variation in the loads is of interest, particularly for a blade structural integrity assessment or to assess the benefits of load alleviation control methods.

The following Appendix describes in more detail the strategy employed for parallelizing the computations related to the ALM and the structural solver. A separate subcommunicator is defined for each wind turbine, which is split from the main communicator based on the turbine location. Any process containing a grid point closer to the hub than a distance *D* is included in the turbine communicator. As the subcommunicators of each turbine generally comprise a distinct set of processors, the tasks related to the ALM and the structural solver can be performed concurrently for each turbine. The different steps performed on the local set of processes are shown in Fig. A1. The subcommunicator is first used to sample the velocity within a restricted region of the physical domain, thereby decreasing the communication overhead. The aerodynamic forces are computed on the root processor using the sampled velocities. Subsequently, the forces are broadcasted to the structural solver, which executes one instance of BeamDyn per blade on different processors. The time integration of the blade dynamics is therefore performed in parallel, and the resulting root moments of each blade are sent to the root process. The latter computes the drivetrain dynamics and the controller commands at each structural sub-step. The beam solver only takes the hub azimuth and the pitch angle (and their time derivatives) as input and only needs to provide the forces and moments on the hub as output. Consequently, the overhead of the added communications required to solve the blade dynamics in parallel is largely overcome by the benefits of the parallelization. After completion of the sub-steps, the blade displacements are broadcasted to all the processors of the subcommunicator to be used for the re-evaluation of the forces and their distribution. The following parallelization allows us to keep the computational cost of the ALM and its coupling to BeamDyn reasonably low and provide an efficient scaling in the case of multiple turbines.

This section presents the statistics of the neutral ABL that is used as inflow for the wind turbines. It is obtained by simulating a turbulent half-channel flow of height *H*=4 *D*; hence *H*=960 m. The flow is driven using a pressure gradient $(-\mathrm{d}p/\mathrm{d}x)/\mathit{\rho}=\mathrm{1.148}\times {\mathrm{10}}^{-\mathrm{4}}$ m s^{−2} (with *ρ*=1.225 kg m^{−3} for air). The lateral boundary conditions are periodic. The upper boundary condition is a no-slip condition. At the ground, a wall model is used. The latter imposes the shear stress at the wall using the law of the wall for a rough wall,

with *κ*=0.385 and *y*_{0}=0.0016 m. The computational domain is $\mathrm{8}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{4}\phantom{\rule{0.125em}{0ex}}D\times \mathrm{4}\phantom{\rule{0.125em}{0ex}}D$ and is discretized using a uniform structured mesh with a resolution of 64 grid points per *D*. The time step is set to d*t*=0.05 s, and the simulation is run for 584 flow-through times ${t}^{*}=t{U}_{\text{hub}}/\left(\mathrm{8}D\right)$ before it reaches convergence. The following statistics are then computed over ≃100 flow-through times.

The obtained profiles of the mean streamwise velocity and the turbulence intensity are shown in Fig. B1. The mean velocity profile indeed corresponds to a logarithmic profile with *u*_{τ}=0.332 m s^{−1}, and the turbulence intensity decreases linearly with height. The mean velocity of 10 m s^{−1} is obtained at the hub, as well as a TI of ≃5 % when computed using the three velocity components. There are some small oscillations in the TI near the ground resulting from the large gradients at the wall. However, these oscillations are located at the bottom of the domain and do not affect the rotor region. The auto-correlation functions of the streamwise velocity, ${\mathit{\rho}}_{u,u}^{x}$ and ${\mathit{\rho}}_{u,u}^{z}$, in the streamwise and lateral directions, respectively, are also depicted in Fig. B2. These are used to compute the integral length scale in the streamwise and lateral directions, ${L}_{u}^{x}$ and ${L}_{u}^{z}$, obtained by integration (Stanislawski et al., 2023):

where $\widehat{\mathit{\delta}}$ is the distance where the correlation function is smaller than 0.05. This leads to ${L}_{u}^{x}$ = 106 m = 0.44 *D* and ${L}_{u}^{z}$ = 89.5 m = 0.37*D*.

The LES flow solver BigFlow is not publicly available. The presented simulation results are available upon request.

FT, PC, and GW conceptualized the research and established the methodology. FT developed the software for the coupling and ran and post-processed the simulations with assistance from PC and GW during the investigation process. FT prepared the original draft with contributions from all authors, which was then reviewed by PC and GW. GW provided funding.

The contact author has declared that none of the authors has any competing interests.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

The authors would like to thank Matthieu Duponcheel and Maud Moens for their advice and assistance during multiple stages of this research.

This research has been supported by the Belgian Federal Government under the Energy Transition Fund (PhairywinD project). Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under grant no. 2.5020.11, and by the Walloon Region. The present research also benefited from computational resources made available on Lucia, the tier-1 supercomputer of the Walloon Region, with infrastructure funded by the Walloon Region under grant no. 1910247.

This paper was edited by Jens Nørkær Sørensen and reviewed by two anonymous referees.

Caprace, D.-G., Chatelain, P., and Winckelmans, G.: Lifting Line with Various Mollifications: Theory and Application to an Elliptical Wing, AIAA Journal, 57, 17–28, https://doi.org/10.2514/1.J057487, 2019. a, b, c, d

Churchfield, M. J., Schreck, S. J., Martinez, L. A., Meneveau, C., and Spalart, P. R.: An Advanced Actuator Line Method for Wind Energy Applications and Beyond, 35th Wind Energy Symposium, 9–13 January 2017, Grapevine, Texas, 1998, https://doi.org/10.2514/6.2017-1998, 2017. a

Dag, K. O. and Sørensen, J. N.: A New Tip Correction for Actuator Line Computations, Wind Energy, 23, 148–160, https://doi.org/10.1002/we.2419, 2020. a

Degroote, J.: Partitioned Simulation of Fluid-Structure Interaction, Archives of Computational Methods in Engineering, 20, 185–238, https://doi.org/10.1007/s11831-013-9085-5, 2013. a

Della Posta, G., Leonardi, S., and Bernardini, M.: A Two-Way Coupling Method for the Study of Aeroelastic Effects in Large Wind Turbines, Renew. Energ., 190, 971–992, https://doi.org/10.1016/j.renene.2022.03.158, 2022. a

Della Posta, G., Leonardi, S., and Bernardini, M.: Large Eddy Simulations of a Utility-Scale Horizontal Axis Wind Turbine Including Unsteady Aerodynamics and Fluid-Structure Interaction Modelling, Wind Energy, 26, 98–125, https://doi.org/10.1002/we.2789, 2023. a

Dose, B., Rahimi, H., Herráez, I., Stoevesandt, B., and Peinke, J.: Fluid-Structure Coupled Computations of the NREL5MW Wind Turbine by Means of CFD, Renew. Energ., 129, 591–605, https://doi.org/10.1016/j.renene.2018.05.064, 2018. a

Duponcheel, M., Bricteux, L., Manconi, M., Winckelmans, G., and Bartosiewicz, Y.: Assessment of RANS and Improved Near-Wall Modeling for Forced Convection at Low Prandtl Numbers Based on LES up to Re_{τ} = 2000, Int. J. Heat Mass Tran., 75, 470–482, https://doi.org/10.1016/j.ijheatmasstransfer.2014.03.080, 2014. a

Fingersh, L., Hand, M., and Laxson, A.: Wind Turbine Design Cost and Scaling Model, Tech. rep., National Renewable Energy Lab. (NREL), Golden, CO, USA, https://doi.org//10.2172/897434, 2006. a

Fritz, E. K., Ferreira, C., and Boorsma, K.: An Efficient Blade Sweep Correction Model for Blade Element Momentum Theory, Wind Energy, 25, 1977–1994, https://doi.org/10.1002/we.2778, 2022. a, b

Gaertner, E., Rinker, J., Sethuraman, L., Zahle, F., Anderson, B., Barter, G. E., Abbas, N. J., Meng, F., Bortolotti, P., Skrzypinski, W., Scott, G., Feil, R., Bredmose, H., Dykes, K., Shields, M., Allen, C., and Viselli, A.: IEA Wind TCP Task 37: Definition of the IEA 15-Megawatt Offshore Reference Wind Turbine, Tech. rep., National Renewable Energy Lab (NREL), Golden, CO, USA, https://doi.org/10.2172/1603478, 2020. a, b, c

Grinderslev, C., Sørensen, N. N., Horcas, S. G., Troldborg, N., and Zahle, F.: Wind turbines in atmospheric flow: fluid–structure interaction simulations with hybrid turbulence modeling, Wind Energ. Sci., 6, 627–643, https://doi.org/10.5194/wes-6-627-2021, 2021. a

Hansen, M. H.: Aeroelastic Instability Problems for Wind Turbines, Wind Energy, 10, 551–577, https://doi.org/10.1002/we.242, 2007. a

Hodgson, E. L., Andersen, S. J., Troldborg, N., Forsting, A. M., Mikkelsen, R. F., and Sørensen, J. N.: A Quantitative Comparison of Aeroelastic Computations Using Flex5 and Actuator Methods in Les, J. Phys. Conf. Ser., 1934, 012014, https://doi.org/10.1088/1742-6596/1934/1/012014, 2021. a, b

Hodgson, E. L., Grinderslev, C., Forsting, A. R. M., Troldborg, N., Sørensen, N. N., Sørensen, J. N., and Andersen, S. J.: Validation of Aeroelastic Actuator Line for Wind Turbine Modelling in Complex Flows, Frontiers in Energy Research, 10, 864645, https://doi.org/10.3389/fenrg.2022.864645, 2022. a

Höning, L., Lukassen, L. J., Stoevesandt, B., and Herráez, I.: Influence of rotor blade flexibility on the near-wake behavior of the NREL 5 MW wind turbine, Wind Energ. Sci., 9, 203–218, https://doi.org/10.5194/wes-9-203-2024, 2024. a

Jha, P. K. and Schmitz, S.: Actuator Curve Embedding – an Advanced Actuator Line Model, J. Fluid Mech., 834, R2, https://doi.org/10.1017/jfm.2017.793, 2018. a, b

Jonkman, J.: The New Modularization Framework for the FAST Wind Turbine CAE Tool, in: 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 202, https://doi.org/10.2514/6.2013-202, 2013. a

Jonkman, J. M., Hayman, G. J., Jonkman, B. J., Damiani, R. R., and Murray, R. E.: AeroDyn V15 User's Guide and Theory Manual, NREL Draft Report, 46 pp., https://www.nrel.gov/wind/nwtc/assets/pdfs/aerodyn-manual.pdf (last access: 16 August 2024), 2015. a

Kleine, V. G., Hanifi, A., and Henningson, D. S.: Non-Iterative Vortex-Based Smearing Correction for the Actuator Line Method, J. Fluid Mech., 961, A29, https://doi.org/10.1017/jfm.2023.237, 2023. a

Krüger, S., Steinfeld, G., Kraft, M., and Lukassen, L. J.: Validation of a coupled atmospheric–aeroelastic model system for wind turbine power and load calculations, Wind Energ. Sci., 7, 323–344, https://doi.org/10.5194/wes-7-323-2022, 2022. a, b

Li, A., Gaunaa, M., Pirrung, G. R., and Horcas, S. G.: A computationally efficient engineering aerodynamic model for non-planar wind turbine rotors, Wind Energ. Sci., 7, 75–104, https://doi.org/10.5194/wes-7-75-2022, 2022. a, b

Madsen, H. A., Riziotis, V., Zahle, F., Hansen, M. O. L., Snel, H., Grasso, F., Larsen, T. J., Politis, E., and Rasmussen, F.: Blade Element Momentum Modeling of Inflow with Shear in Comparison with Advanced Model Results, Wind Energy, 15, 63–81, https://doi.org/10.1002/we.493, 2012. a

Mann, J.: Wind Field Simulation, Probabilistic Eng. Mech., 13, 269–282, https://doi.org/10.1016/S0266-8920(97)00036-2, 1998. a, b

Manolas, D. I., Riziotis, V. A., and Voutsinas, S. G.: Assessing the Importance of Geometric Nonlinear Effects in the Prediction of Wind Turbine Blade Loads, J. Comput. Nonlin. Dyn., 10, 041008, https://doi.org/10.1115/1.4027684, 2015. a

Martínez-Tossas, L. A. and Meneveau, C.: Filtered Lifting Line Theory and Application to the Actuator Line Model, J. Fluid Mech., 863, 269–292, https://doi.org/10.1017/jfm.2018.994, 2019. a

Meneveau, C. and Lund, T. S.: The Dynamic Smagorinsky Model and Scale-Dependent Coefficients in the Viscous Range of Turbulence, Phys. Fluids, 9, 3932–3934, https://doi.org/10.1063/1.869493, 1997. a

Meng, H., Lien, F.-S., and Li, L.: Elastic Actuator Line Modelling for Wake-Induced Fatigue Analysis of Horizontal Axis Wind Turbine Blade, Renew. Energ., 116, 423–437, https://doi.org/10.1016/j.renene.2017.08.074, 2018. a

Meng, H., Li, L., and Zhang, J.: A Preliminary Numerical Study of the Wake Effects on the Fatigue Load for Wind Farm Based on Elastic Actuator Line Model, Renew. Energ., 162, 788–801, https://doi.org/10.1016/j.renene.2020.07.140, 2020. a

Meyer Forsting, A. R., Pirrung, G. R., and Ramos-García, N.: A vortex-based tip/smearing correction for the actuator line, Wind Energ. Sci., 4, 369–383, https://doi.org/10.5194/wes-4-369-2019, 2019. a

Mikkelsen, R. F.: Actuator disk methods applied to wind turbines, PhD thesis, Technical University of Denmark, PhD thesis at Technical University of Denmark, https://orbit.dtu.dk/en/publications/actuator-disc-methods-applied-to-wind-turbines (last access: 16 August 2024), 2004. a

Moens, M.: Large Eddy Simulation of wind farm flows: improved actuator disk model and investigation of wake phenomena, PhD thesis at Université Catholique de Louvain, http://hdl.handle.net/2078.1/207779 (last access: 16 August 2024), 2018. a

Moens, M., Duponcheel, M., Winckelmans, G., and Chatelain, P.: An Actuator Disk Method with Tip-Loss Correction Based on Local Effective Upstream Velocities, Wind Energy, 21, 766–782, https://doi.org/10.1002/we.2192, 2018. a

Moriarty, P. J. and Hansen, A. C.: AeroDyn Theory Manual, Tech. rep., National Renewable Energy Lab, Golden, CO, USA, https://doi.org/10.2172/15014831, 2005. a

Panteli, A. N., Manolas, D. I., Riziotis, V. A., and Spiliopoulos, K. V.: Comparative Study of Two Geometrically Non-Linear Beam Approaches for the Coupled Wind Turbine System, J. Wind Eng. Ind. Aerod., 231, 105231, https://doi.org/10.1016/j.jweia.2022.105231, 2022. a

Perez-Becker, S., Papi, F., Saverin, J., Marten, D., Bianchini, A., and Paschereit, C. O.: Is the Blade Element Momentum theory overestimating wind turbine loads? – An aeroelastic comparison between OpenFAST's AeroDyn and QBlade's Lifting-Line Free Vortex Wake method, Wind Energ. Sci., 5, 721–743, https://doi.org/10.5194/wes-5-721-2020, 2020. a

Santo, G., Peeters, M., Van Paepegem, W., and Degroote, J.: Dynamic Load and Stress Analysis of a Large Horizontal Axis Wind Turbine Using Full Scale Fluid-Structure Interaction Simulation, Renew. Energ., 140, 212–226, https://doi.org/10.1016/j.renene.2019.03.053, 2019. a, b

Sayed, M., Lutz, Th., Krämer, E., Shayegan, Sh., and Wüchner, R.: Aeroelastic Analysis of 10 MW Wind Turbine Using CFD–CSD Explicit FSI-coupling Approach, J. Fluid. Struct., 87, 354–377, https://doi.org/10.1016/j.jfluidstructs.2019.03.023, 2019. a

Shaler, K., Branlard, E., and Platt, A.: OLAF User's Guide and Theory Manual, Tech. rep., National Renewable Energy Lab (NREL), Golden, CO, USA, https://doi.org/10.2172/15014831, 2020. a, b

Shaler, K., Anderson, B., Martínez-Tossas, L. A., Branlard, E., and Johnson, N.: Comparison of free vortex wake and blade element momentum results against large-eddy simulation results for highly flexible turbines under challenging inflow conditions, Wind Energ. Sci., 8, 383–399, https://doi.org/10.5194/wes-8-383-2023, 2023. a

Sieros, G., Chaviaropoulos, P., Sørensen, J. D., Bulder, B. H., and Jamieson, P.: Upscaling Wind Turbines: Theoretical and Practical Aspects and Their Impact on the Cost of Energy, Wind Energy, 15, 3–17, https://doi.org/10.1002/we.527, 2012. a

Smagorinsky, J.: General Circulation Experiments with the Primitive Equations. I. The Basic Experiment, Mon. Weather Rev., 91, 99–164, 1963. a

Sørensen, J. N. and Larsen, G. C.: A Minimalistic Prediction Model to Determine Energy Production and Costs of Offshore Wind Farms, Energies, 14, 448, https://doi.org/10.3390/en14020448, 2021. a

Sørensen, J. N. and Larsen, G. C.: Assessment of Energy Production and Costs Associated with a Massive Exploitation of Wind Power in the North Sea, J. Phys. Conf. Ser., 2505, 012013, https://doi.org/10.1088/1742-6596/2505/1/012013, 2023. a

Sorensen, J. N. and Shen, W. Z.: Numerical Modeling of Wind Turbine Wakes, J. Fluid. 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. R. Soc. A, 373, 20140071, https://doi.org/10.1098/rsta.2014.0071, 2015. a

Sprague, M. A., Ananthan, S., Vijayakumar, G., and Robinson, M.: ExaWind: A Multifidelity Modeling and Simulation Environment for Wind Energy, J. Phys. Conf. Ser., 1452, 012071, https://doi.org/10.1088/1742-6596/1452/1/012071, 2020. a

Sprague, M. A. and Geers, T. L.: Legendre Spectral Finite Elements for Structural Dynamics Analysis, Commun. Numer. Meth. En., 24, 1953–1965, https://doi.org/10.1002/cnm.1086, 2008. a

Spyropoulos, N., Papadakis, G., Prospathopoulos, J. M., and Riziotis, V. A.: Investigating the Level of Fidelity of an Actuator Line Model in Predicting Loads and Deflections of Rotating Blades under Uniform Free-Stream Flow, Appl. Sci.-Basel, 11, 12097, https://doi.org/10.3390/app112412097, 2021. a

Stanislawski, B. J., Thedin, R., Sharma, A., Branlard, E., Vijayakumar, G., and Sprague, M. A.: Effect of the Integral Length Scales of Turbulent Inflows on Wind Turbine Loads, Renew. Energ., 217, 119218, https://doi.org/10.1016/j.renene.2023.119218, 2023. a

Storey, R. C., Norris, S. E., Stol, K. A., and Cater, J. E.: Large Eddy Simulation of Dynamically Controlled Wind Turbines in an Offshore Environment, Wind Energy, 16, 845–864, https://doi.org/10.1002/we.1525, 2013. a

Thiry, O.: Investigation of wall shear stress models for large eddy simulation, PhD thesis at Université Catholique de Louvain, http://hdl.handle.net/2078.1/184590 (last access: 16 August 2024), 2017. a

Trigaux, F., Chatelain, P., and Winckelmans, G.: A Flexible Actuator Curve Model for Aeroelastic Simulations of Wind Turbines in Atmospheric Boundary Layers, J. Phys. Conf. Ser., 2265, 022050, https://doi.org/10.1088/1742-6596/2265/2/022050, 2022. a, b

Troldborg, N.: Actuator line modelling of wind turbine wakes, PhD thesis at Technical University of Denmark, https://orbit.dtu.dk/en/publications/actuator-line-modeling-of-wind-turbine-wakes (last access: 16 August 2024), 2008. a

Vitsas, A. and Meyers, J.: Multiscale Aeroelastic Simulations of Large Wind Farms in the Atmospheric Boundary Layer, J. Phys. Conf. Ser., 753, 082020, https://doi.org/10.1088/1742-6596/753/8/082020, 2016. a

Wang, Q., Sprague, M. A., Jonkman, J., Johnson, N., and Jonkman, B.: BeamDyn: A High-Fidelity Wind Turbine Blade Solver in the FAST Modular Framework, Wind Energy, 20, 1439–1462, https://doi.org/10.1002/we.2101, 2017. a

Welch, P.: The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging over Short, Modified Periodograms, IEEE T. Acoust. Speech, 15, 70–73, https://doi.org/10.1109/TAU.1967.1161901, 1967. a

^{1}

This model goes back to the regularization by Rosenhead (1935) of point vortices, where $\mathrm{1}/r=r/{r}^{\mathrm{2}}$ is replaced by $r/({r}^{\mathrm{2}}+{r}_{\mathrm{c}}^{\mathrm{2}})$.

- Abstract
- Introduction
- Methodology
- Comparison to the BEM theory and a free vortex wake model
- Investigation of rigid and flexible rotors in atmospheric conditions
- Investigation of the impact of the spatial resolution
- Conclusions
- Appendix A: Parallel implementation of the ALM and structural solver
- Appendix B: Statistics of the simulated atmospheric boundary layer
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Methodology
- Comparison to the BEM theory and a free vortex wake model
- Investigation of rigid and flexible rotors in atmospheric conditions
- Investigation of the impact of the spatial resolution
- Conclusions
- Appendix A: Parallel implementation of the ALM and structural solver
- Appendix B: Statistics of the simulated atmospheric boundary layer
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References