**1. Introduction**

Renewable energy sources have been gaining more importance in the last few decades as part of the strategies adopted by countries all over the world to limit the use of fossil fuel and fight pollution and global warming. In particular, wind energy has grown rapidly, resulting in the increasing size of modern wind turbines with the objective of reducing the cost of the produced energy [1]. Nevertheless, the adoption of more slender blades has also led to higher deflections during normal operation, and consequently, more interest toward the fluid–structure interaction (FSI) phenomenon in modern wind turbines. Recent works have computed that, while operating in design conditions, the flapwise deflection of a modern horizontal-axis wind turbine (HAWT) blade is in the order of 6–7% of its span [2–4]. Furthermore, their deflections have a sensible impact on the produced power, affecting its oscillation or introducing a reduction up to 6% [2–6].

However, due to the complexities involved, when the FSI of modern wind turbine blades is to be analyzed, simplified models are often adopted for this task. On the aerodynamic side of the problem, among others, blade element momentum theory (BEM) is widely used in the FSI modeling of wind turbines [7–9]. Despite its low computational cost, BEM theory is affected by many limitations, including the need to include tip-loss corrections to account for a blade of finite length [10]. Another class of widely used models are the actuator models, where the blades are represented by lines or surfaces exchanging momentum with the incoming wind flow [11–13]. This strategy has also been

used in FSI simulations of a multi-megawatt turbine, coupling it with a structural model based on non-linear beam theory [14].

In order to increase the level of detail in the extracted results, a higher computational cost is required. Computational fluid dynamics (CFD) simulations, performed by solving the Navier–Stokes equations on a computational grid rendering the geometry of the wind turbine, have recently been used in literature within FSI frameworks. In particular, Reynolds-averaged Navier–Stokes (RANS) models are often used to account for turbulence in the atmospheric wind flow [15–19]. This class of turbulence models are computationally cheap when compared to more complex turbulence models, which are also reported in recent literature within FSI simulations of wind turbines [2–4,20,21] but have an extremely high computational cost.

On the structural side of the FSI problem, the complex nature of modern composite blades [22] often leads to the adoption of simplified models. Among others, FSI works relying on multi-body dynamics [21] or using one-dimensional beam elements [14,15,19] can be found, whereas the implementation of detailed finite element models (FEMs), loyally reproducing the composite structure of modern blades, is still limited to only a few works [3,4,6,20] by two distinct research groups.

Despite its importance, the atmospheric boundary layer (ABL), namely the velocity gradient leading to higher wind speeds at higher heights, is very often neglected in FSI simulations of wind turbines. To the best of the authors' knowledge, only a very limited number of works account for it [6,21,23]. These works remark upon the importance of the ABL, highlighting how it induces oscillating loads, power, and deflections on modern blades.

Furthermore, given the aleatory nature of the wind, wind turbines always operate in unsteady and rapidly changing conditions. In particular, local oscillations in the wind speed, referred to as "gusts," lead to the varying performance and structural responses of the blades, which is particularly important in the fatigue life estimation of modern structures. Gusts can be of various shapes, sizes, lengths, and intensities [24], and can be induced by the terrain morphology, as well as by thermal or turbulent effects [25]. Despite their relevance in relation to large HAWTs, the majority of works currently available in the literature focus on vertical-axis wind turbines (VAWTs) or on small HAWTs. Wu et al. [26] carried out 2D RANS simulations to obtain the lift–drag polars of the sections of a VAWT in both steady wind and unsteady (i.e., affected by a gust) wind conditions and used them to obtain the performance of the turbine by means of a BEM-like strategy. Onol and Yesilyurt [27] analyzed a small VAWT by means of a 2D unsteady RANS (U-RANS) model, changing the approaching wind velocity according to the gust time variation prescribed by the IEC-64100 standards for wind turbines. Where 3D simulations are concerned, Bhargav et al. [28] performed CFD simulations of a VAWT while changing the inlet wind speed to follow a sinusoidal function but without considering the FSI.

Regarding to FSI simulations, the work of Timme et al. [29] analyzed the aeroelastic response of a straight wing impacted by a vertical gust. The shape of this gust corresponds to the "1 − cos" shape, also prescribed in the IEC-64100 standards. Also operating in an aeroelastic framework, Svacek and Horacek [30] analyzed the response of a flexibly supported airfoil to a vertical wind gust, mimicking the dynamics of a wing by means of 2D U-RANS simulations.

The available literature on the gust aeroelastic response to HAWT is extremely limited. First in time, the work of Younsi et al. [31] featured a BEM code in combination with a structural model reproducing the blade as being homogenously made of an elastic and isotropic material. A single blade was modelled, neglecting the ABL and simulating an extreme gust impacting on the whole blade. More recently, Castellani et al. [32] carried out both experimental and numerical works to investigate the response of a small HAWT (2-m diameter) to a periodical change in the incoming wind speed, aiming at the investigation of an optimal control strategy for such a machine. On the other hand, Ebrahimi and Sekandari [33] investigated the aeroelasticity of a multi-megawatt wind turbine to a sudden change in wind, coupling an unsteady vortex-lattice method (VLM) with a structural model relying on a modal approach. The entire turbine was subjected to a change of velocity magnitude and

direction, also accounting for the response of the control systems that react by changing the pitch of the blades and the yaw of the whole machine.

All the aforementioned works on HAWTs neglect the ABL and analyze a wind gust bigger than the analyzed turbine by means of simplified models. The investigated gust impacts on the whole rotor can be counteracted by the turbine control systems [33]. No work about the aeroelastic response of the blades of a large HAWT immersed in the ABL and locally impacted by a wind gust (i.e., soliciting only one blade) was found in the current literature. In this work, for the first time, two high-fidelity models, one for the CFD side and one for the FEM side of the FSI problem, were strongly coupled to dynamically analyze the aeroelasticity of the blade of a 100-m diameter HAWT immersed in the ABL and impacted by wind gusts of different intensities and morphologies. The proposed methodology is believed to be well-suited for advanced engineering applications, and in particular, to analyze the response of modern wind turbines to extreme load cases as part of the design process. This is an added value compared to the available literature, in which the aeroelastic analysis of a wind turbine blade attacked by a wind gust is scarcely reported.

The paper is structured as follows. The CFD model is presented in Section 2.1, the structural FEM model is described in Section 2.2, the coupling strategy is addressed in Section 2.3, followed by the gust model in Section 2.4. Section 3 contains the extracted results and their discussion, and the conclusions are drawn in Section 4.
