Fluid–Structure Interaction Simulations of a Wind Gust Impacting on the Blades of a Large Horizontal Axis Wind Turbine

The effect of a wind gust impacting on the blades of a large horizontal-axis wind turbine is analyzed by means of high-fidelity fluid–structure interaction (FSI) simulations. The employed FSI model consisted of a computational fluid dynamics (CFD) model reproducing the velocity stratification of the atmospheric boundary layer (ABL) and a computational structural mechanics (CSM) model loyally reproducing the composite materials of each blade. Two different gust shapes were simulated, and for each of them, two different amplitudes were analyzed. The gusts were chosen to impact the blade when it pointed upwards and was attacked by the highest wind velocity due to the presence of the ABL. The loads and the performance of the impacted blade were studied in detail, analyzing the effect of the different gust shapes and intensities. Also, the deflections of the blade were evaluated and followed during the blade’s rotation. The flow patterns over the blade were monitored in order to assess the occurrence and impact of flow separation over the monitored quantities.


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][3][4]. Furthermore, their deflections have a sensible impact on the produced power, affecting its oscillation or introducing a reduction up to 6% [2][3][4][5][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][8][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][12][13]. This strategy has also been 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.

Methodology
In this section, the adopted numerical model is described. It is built analogously to the model in Santo et al. [6]. Therefore, the methodology is shortly summarized here and the reader is referred to Santo et al. [6] for a more extensive description, validation, and sensitivity study of the chosen methods.

The CFD Model
On the CFD side of the employed FSI model, the domain illustrated in Figure 1 was used.
Energies 2020, 13,509 3 of 20 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 highfidelity 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.

Methodology
In this section, the adopted numerical model is described. It is built analogously to the model in Santo et al. [6]. Therefore, the methodology is shortly summarized here and the reader is referred to Santo et al. [6] for a more extensive description, validation, and sensitivity study of the chosen methods.

The CFD Model
On the CFD side of the employed FSI model, the domain illustrated in Figure 1 was used. An overset method (also called Chimera) approach was chosen to handle the rotation and deformation of the blades, resulting in a background grid and several component grids. Every grid was entirely made of hexahedral finite volumes and the state variables were stored in the cell centers. An overset method (also called Chimera) approach was chosen to handle the rotation and deformation of the blades, resulting in a background grid and several component grids. Every grid was entirely made of hexahedral finite volumes and the state variables were stored in the cell centers. Figure 1 shows the structured grid employed for the background. In its finest region, the cells were Energies 2020, 13, 509 4 of 20 cubic and had a 0.275-m edge. In this location, the analyzed HAWT was positioned, being five rotor diameters distant from inlet, symmetry sides, and top surfaces. The atmospheric pressure outlet was placed 15 diameters downstream of the rotor. These distances were chosen according to good practice guidelines for atmospheric flows in urban areas [34]. However, it is reported that remarkably smaller (approximately half) distances are found to be sufficient to minimize boundary influence [35].
A body-fitted hexahedral component mesh was constructed around each blade and around the tower and the nacelle, as depicted in Figure 2. The tower was considered rigid and the geometry of the tower was extracted from Harte and Van Zijl [36], being used for wind turbines of similar size. The hub height was 100 m and the nacelle had a length of approximately 12 m and a section of 5 m × 5 m. The right-hand side of Figure 2 shows different sections of the mesh around each blade, taken at different radial spanwise locations. The mesh on the blade walls was designed to have a y + in the log layer (between 30 and 300), while on the border of each body-fitted mesh, a mesh size comparable to the one used in the background was imposed.
Energies 2020, 13,509 4 of 20 Figure 1 shows the structured grid employed for the background. In its finest region, the cells were cubic and had a 0.275-m edge. In this location, the analyzed HAWT was positioned, being five rotor diameters distant from inlet, symmetry sides, and top surfaces. The atmospheric pressure outlet was placed 15 diameters downstream of the rotor. These distances were chosen according to good practice guidelines for atmospheric flows in urban areas [34]. However, it is reported that remarkably smaller (approximately half) distances are found to be sufficient to minimize boundary influence [35]. A body-fitted hexahedral component mesh was constructed around each blade and around the tower and the nacelle, as depicted in Figure 2. The tower was considered rigid and the geometry of the tower was extracted from Harte and Van Zijl [36], being used for wind turbines of similar size. The hub height was 100 m and the nacelle had a length of approximately 12 m and a section of 5 m × 5 m. The right-hand side of Figure 2 shows different sections of the mesh around each blade, taken at different radial spanwise locations. The mesh on the blade walls was designed to have a in the log layer (between 30 and 300), while on the border of each body-fitted mesh, a mesh size comparable to the one used in the background was imposed. The air flow was modelled as incompressible (air density constant equal to 1.225 kg/m 3 ) and the -(RANS) turbulence model was used. For this model, Richard and Hoxey [37] developed inlet profiles for the velocity (v), turbulent kinetic energy ( ) and its dissipation ratio ( ) as functions of the height from the ground ( . These profiles are obtained as analytical solutions of the transport equations of and and are reported in Equations (1)-(3): In these equations, and are the von Karman constant (0.4187) and a constant of themodel (0.09), respectively. Moreover, * is the friction velocity (a global index of the wind intensity), and is the aerodynamic roughness length, which provides a measure of how rough the ground wall is [38]. These last two parameters fully define the ABL profiles.
Nevertheless, in order to preserve these profiles while travelling through the computational domain, a modified formulation is necessary for the ground wall, as observed in different works [39,40]. Parente et al. [40,41] have therefore obtained modified wall functions, where the aerodynamic The air flow was modelled as incompressible (air density constant equal to 1.225 kg/m 3 ) and the k-ε (RANS) turbulence model was used. For this model, Richard and Hoxey [37] developed inlet profiles for the velocity (v), turbulent kinetic energy (k) and its dissipation ratio (ε) as functions of the height from the ground (z). These profiles are obtained as analytical solutions of the transport equations of k and ε and are reported in Equations (1)-(3): In these equations, K and C µ are the von Karman constant (0.4187) and a constant of the k-ε model (0.09), respectively. Moreover, v * is the friction velocity (a global index of the wind intensity), and z 0 is the aerodynamic roughness length, which provides a measure of how rough the ground wall is [38]. These last two parameters fully define the ABL profiles.
Nevertheless, in order to preserve these profiles while travelling through the computational domain, a modified formulation is necessary for the ground wall, as observed in different works [39,40].
Parente et al. [40,41] have therefore obtained modified wall functions, where the aerodynamic roughness length is directly included in the formulation of the non-dimensional wall distance z + mod and of the wall function constant E mod : In these equations, ρ and µ are the air density (1.225 kg/m 3 ) and its viscosity (1.7894 × 10 −5 kg/ms), respectively. The modified wall functions were therefore used on the ground wall, while on the walls of the wind turbine, standard wall functions were adopted.
All the investigated load scenarios were carried out at the turbine nominal operating point, where it is expected to operate for most of its life. Therefore, following input from the blade's manufacturer, the rotational speed of the machine was set to 1.445 rad/s and the wind velocity at the hub height was set to 8.5 m/s at the inlet, leading to a tip speed ratio of 8.5 using undisturbed conditions. The axis of rotation was perfectly aligned with the wind flow (i.e., no tilt angle, no yaw misalignment). Consistently, the friction velocity and the aerodynamic roughness length were set equal to 0.671 m/s and 0.5 m, respectively, to achieve the desired wind speed at the hub height. The inlet turbulent kinetic energy was 0.01512 m 2 /s 2 , leading to a turbulence intensity of 1.3% at the hub height.
The mass and momentum conservation equations were solved together using a pressure-based solver. For the momentum equation, a second-order upwind discretization scheme was adopted, while a first-order implicit scheme was selected for time discretization. All the simulations presented in this work featured the same numerical setup. The entire CFD model was implemented in Fluent 18.1 (Ansys Inc., Canonsburg, PA, USA).
According to the outcome of the mesh and time step sensitivity study carried out in Santo et al. [6], the selected mesh had a total of 55 million cells, 1.9 million belonging to each blade component mesh. The surface of each blade wall was divided into approximately 38,500 faces. The time step size was 0.01812 s, corresponding to 1.5 • of rotor rotation per time step. With these settings, the torque coefficient produced by the turbine, considering rigid blades and no gust, was computed to be 0.05221, which compares well with what the manufacturer provides for this operating point (0.0556).
It is remarked that the simulations presented in this work cannot be validated by experimental results, given the difficulty in reproducing the imposed gusts in a controlled way.

The Structural FEM Model
Each of the three blades was 50-m long and entirely made of composite material. In order to loyally mimic its structure, the blade was divided into approximately 64,000 shell elements with reduced integration. The elements were distributed on the outer mold layer of the blade and on its shear webs and shear caps ( Figure 3). In each shell element, a global layup orientation was defined, with respect to which, a ply orientation was defined in each ply of the composite stack in order to fully define the anisotropic properties of the used materials. The generation of the structural mesh is described in Peeters et al. [22] and the comparison of the computed eigenfrequencies with the values provided by the manufacturer is given in Santo et al. [6]. Gravity was also included in the model. Energies 2020, 13, 509 6 of 20

FSI Coupling
The CFD and the structural models outlined in the previous sections were coupled by means of Tango, an in-house code developed at Ghent University, resulting in a partitioned FSI simulation [42]. The Gauss-Seidel coupling algorithm was selected and three iterations were enough within each time step to reach a displacement absolute residual of about 5 mm on the fluid-structure interface. The fluid mesh was deformed according to what is prescribed by the structural solver on the fluid-structure interface. Each blade component mesh was deformed by means of a spring-based smoothing method, therefore adopting an arbitrary Lagrangian-Eulerian (ALE) formulation. Compared to other methods present in the literature [43], the adopted methodology can consistently preserve an appropriate value of the stretched cells in the region adjacent to each wall. This ensured a well-resolved near-wall flow at a reasonable computational cost.
Simulations were performed on 280 cores (10 nodes inter-connected by InfiniBand, each with 2 CPUs of the type 14-core Xeon E5-2680v4, 2.4 GHz) and approximately 10 days were necessary to perform a complete revolution. Starting the FSI simulation from the results of a CFD simulation, approximately 1.2 revolutions were necessary to reach the regime in time. Then, one full rotation was performed and analyzed. During this rotation, the loads, stress, and displacements of each blade could be univocally linked to its azimuth angle. The azimuth angle was set to 0 when the blade was horizontally positioned and in an upward motion. Therefore, a 90° azimuth angle corresponded to the blade pointing upward and 270° to the blade pointing downward.

Gust Model
Different gust shapes and intensities were imposed to impact one of the three blades. In order to investigate the worst-case scenario, the gusts were chosen to collide on the blade when the blade pointed upward (i.e., at a 90° azimuth angle); in this position, the blade was attacked by the highest wind speed, and therefore, the solicitations acting on it were also the highest. In this work, as in other works on wind gusts [29,44], the gusts were modelled as a local streamwise velocity increase to be superimposed on the local wind velocity field. Using to indicate the gust's additional velocity and the local wind velocity, the velocity in "gusted" conditions ′ can be expressed as: Notice that, at each time instant, all these velocities are functions of the position in the numerical domain. In order to minimize the computational time and numerical dissipation of the gusts [29,44], the gusts were added in the proximity of the turbine when the blade to be hit was positioned at a 60° azimuth angle, namely 30° and 20 time steps in advance with respect to the impact of the gust on the blade. They affected only a cylindrical region of 25 m in diameter (i.e., half of the blade span) and 12 m in length, whose axis was aligned with the wind direction. The frontal tip of this cylindrical region

FSI Coupling
The CFD and the structural models outlined in the previous sections were coupled by means of Tango, an in-house code developed at Ghent University, resulting in a partitioned FSI simulation [42]. The Gauss-Seidel coupling algorithm was selected and three iterations were enough within each time step to reach a displacement absolute residual of about 5 mm on the fluid-structure interface. The fluid mesh was deformed according to what is prescribed by the structural solver on the fluid-structure interface. Each blade component mesh was deformed by means of a spring-based smoothing method, therefore adopting an arbitrary Lagrangian-Eulerian (ALE) formulation. Compared to other methods present in the literature [43], the adopted methodology can consistently preserve an appropriate y + value of the stretched cells in the region adjacent to each wall. This ensured a well-resolved near-wall flow at a reasonable computational cost.
Simulations were performed on 280 cores (10 nodes inter-connected by InfiniBand, each with 2 CPUs of the type 14-core Xeon E5-2680v4, 2.4 GHz) and approximately 10 days were necessary to perform a complete revolution. Starting the FSI simulation from the results of a CFD simulation, approximately 1.2 revolutions were necessary to reach the regime in time. Then, one full rotation was performed and analyzed. During this rotation, the loads, stress, and displacements of each blade could be univocally linked to its azimuth angle. The azimuth angle was set to 0 when the blade was horizontally positioned and in an upward motion. Therefore, a 90 • azimuth angle corresponded to the blade pointing upward and 270 • to the blade pointing downward.

Gust Model
Different gust shapes and intensities were imposed to impact one of the three blades. In order to investigate the worst-case scenario, the gusts were chosen to collide on the blade when the blade pointed upward (i.e., at a 90 • azimuth angle); in this position, the blade was attacked by the highest wind speed, and therefore, the solicitations acting on it were also the highest. In this work, as in other works on wind gusts [29,44], the gusts were modelled as a local streamwise velocity increase to be superimposed on the local wind velocity field. Using v g to indicate the gust's additional velocity and v the local wind velocity, the velocity in "gusted" conditions v can be expressed as: Notice that, at each time instant, all these velocities are functions of the position in the numerical domain. In order to minimize the computational time and numerical dissipation of the gusts [29,44], the gusts were added in the proximity of the turbine when the blade to be hit was positioned at a 60 • Energies 2020, 13, 509 7 of 20 azimuth angle, namely 30 • and 20 time steps in advance with respect to the impact of the gust on the blade. They affected only a cylindrical region of 25 m in diameter (i.e., half of the blade span) and 12 m in length, whose axis was aligned with the wind direction. The frontal tip of this cylindrical region was positioned at an appropriate axial coordinate (approximately 1 m downstream of the blade tip's axial coordinate) in order to ensure the gust hit the targeted blade. This means that only in this cylindrical region, the gust velocity used in Equation (6) could be written as: In Equation (7), s and f are shape functions of the gust, depending respectively on the axial (a) and radial (r) coordinates inside the gusted cylindrical region. In particular, the function s sinusoidally ramped up from 0 to 1 in the first 3 m of the axial length of the cylinder, was kept constant and equal to 1 in the middle, and then sinusoidally ramped down from 1 to 0 in the last 3 m of the axial length. Figure 4 illustrates the initial position and shape of the imposed gusts, as well as the coordinates a and r.
Energies 2020, 13, 509 7 of 20 was positioned at an appropriate axial coordinate (approximately 1 m downstream of the blade tip's axial coordinate) in order to ensure the gust hit the targeted blade. This means that only in this cylindrical region, the gust velocity used in Equation (6) could be written as: In Equation (7), and are shape functions of the gust, depending respectively on the axial ( ) and radial ( ) coordinates inside the gusted cylindrical region. In particular, the function sinusoidally ramped up from 0 to 1 in the first 3 m of the axial length of the cylinder, was kept constant and equal to 1 in the middle, and then sinusoidally ramped down from 1 to 0 in the last 3 m of the axial length. Figure 4 illustrates the initial position and shape of the imposed gusts, as well as the coordinates and . For what concerns the gust amplitude , probabilistic analyses of wind gusts have shown that a gust amplitude of 5 m/s has more than an 80% probability of daily occurrence over central Europe when the wind conditions are similar to the ones chosen in this work [45]. Similarly, gusts exceeding this speed are commonly reported in Germany [46]. Also, more intense gusts (18 to 25 m/s) are fairly common in Europe, especially in coastal areas [25,47]. For these reasons, two gust amplitudes were used in this work, namely 5 and 10 m/s. Lastly, two shape functions were used in this work. In the first subsection, a novel gust shape function was proposed, imposing a local redistribution of the flow rate and no global change in it, providing a velocity deficit on its border and a velocity increase in its center. This gust shape was conceptually similar to the "extreme operating gust" from the 61400-IEC standards for wind turbines. Subsequently, in the second subsection, a consistently positive gust velocity was used, analogous to the "extreme coherent gust" from the 61400-IEC standards for wind turbines.

Results and Discussion
In this section, the results of the performed simulations are presented. In order to analyze the results, the blade was divided into 10 equally spaced sections, as depicted in Figure 5. These sections are addressed as "strips" and numbered according to the notation in the figure. For what concerns the gust amplitude A g , probabilistic analyses of wind gusts have shown that a gust amplitude of 5 m/s has more than an 80% probability of daily occurrence over central Europe when the wind conditions are similar to the ones chosen in this work [45]. Similarly, gusts exceeding this speed are commonly reported in Germany [46]. Also, more intense gusts (18 to 25 m/s) are fairly common in Europe, especially in coastal areas [25,47]. For these reasons, two gust amplitudes A g were used in this work, namely 5 and 10 m/s. Lastly, two shape functions f were used in this work. In the first subsection, a novel gust shape function was proposed, imposing a local redistribution of the flow rate and no global change in it, providing a velocity deficit on its border and a velocity increase in its center. This gust shape was conceptually similar to the "extreme operating gust" from the 61400-IEC standards for wind turbines. Subsequently, in the second subsection, a consistently positive gust velocity was used, analogous to the "extreme coherent gust" from the 61400-IEC standards for wind turbines.

Results and Discussion
In this section, the results of the performed simulations are presented. In order to analyze the results, the blade was divided into 10 equally spaced sections, as depicted in Figure 5. These sections are addressed as "strips" and numbered according to the notation in the figure. Furthermore, when presenting the results, the aerodynamic torque ( ) and axial force ( , aligned with the axis of rotation) acting on the blades was made non-dimensional by means of Equations (8) and (9), as is normally done in the wind energy community: where is the air density (1.225 kg/m 3 ), is the swept area of the rotor, and is its radius. The velocity is chosen to be the wind-free stream velocity at the hub height, namely 8.5 m/s.

Zero Net Flow Rate Gust
In this section, the effect of a gust corresponding to a local redistribution of the available flow rate is investigated. The parameter indicates the non-dimensional radial coordinate in the gusted cylinder, ranging from 0 (gust center) to 1 (gust border). Imposing C0 and C1 continuity at the border of the cylindrical region ( 1 = 0 and 1 = 0), maximum gust velocity at the center ( 0 = 1 and 0 = 0), and zero net flow rate ( = 0), the following polynomial expression is obtained for the radial shape function: The obtained function, compared to the "extreme operating gust" from the 61400-IEC standards (used as a function of space, as is done in De Nayer et al. [44]), produced a higher velocity increase, and thus more severe wind conditions, as depicted by Figure 6. Notice that since this gust shape does not modify the net mass flow rate in the affected area but only redistributes it, it is considered appropriate to be used on inlet boundaries in combination with incompressible solvers, without the necessity to correct the mass flow rate to preserve its steady value Furthermore, when presenting the results, the aerodynamic torque (T) and axial force (F, aligned with the axis of rotation) acting on the blades was made non-dimensional by means of Equations (8) and (9), as is normally done in the wind energy community: where ρ is the air density (1.225 kg/m 3 ), A is the swept area of the rotor, and R is its radius. The velocity v is chosen to be the wind-free stream velocity at the hub height, namely 8.5 m/s.

Zero Net Flow Rate Gust
In this section, the effect of a gust corresponding to a local redistribution of the available flow rate is investigated. The parameter r indicates the non-dimensional radial coordinate in the gusted cylinder, ranging from 0 (gust center) to 1 (gust border). Imposing C0 and C1 continuity at the border of the cylindrical region ( f (1) = 0 and f (1) = 0), maximum gust velocity at the center ( f (0) = 1 and f (0) = 0), and zero net flow rate ( 2π 0 1 0 f (r) r dr dϑ = 0), the following polynomial expression is obtained for the radial shape function: The obtained function, compared to the "extreme operating gust" from the 61400-IEC standards (used as a function of space, as is done in De Nayer et al. [44]), produced a higher velocity increase, and thus more severe wind conditions, as depicted by Figure 6.
Notice that since this gust shape does not modify the net mass flow rate in the affected area but only redistributes it, it is considered appropriate to be used on inlet boundaries in combination with incompressible solvers, without the necessity to correct the mass flow rate to preserve its steady value [44].
Depending on the radial position where the gust hits the blade, a different effect was found regarding its axial deformation. In order to assess where to hit the blade to obtain the highest effect, several simulations were carried out, positioning the center of the gust at a distance h (Figure 4) equal to 35 m, 40 m, and 45 m starting from the axis. Figure 7 summarizes the results of these simulations.
Energies 2020, 13, 509 9 of 20 is obtained for the radial shape function: The obtained function, compared to the "extreme operating gust" from the 61400-IEC standards (used as a function of space, as is done in De Nayer et al. [44]), produced a higher velocity increase, and thus more severe wind conditions, as depicted by Figure 6. Notice that since this gust shape does not modify the net mass flow rate in the affected area but only redistributes it, it is considered appropriate to be used on inlet boundaries in combination with incompressible solvers, without the necessity to correct the mass flow rate to preserve its steady value [44].  Depending on the radial position where the gust hits the blade, a different effect was found regarding its axial deformation. In order to assess where to hit the blade to obtain the highest effect, several simulations were carried out, positioning the center of the gust at a distance h (Figure 4) equal to 35 m, 40 m, and 45 m starting from the axis. Figure 7 summarizes the results of these simulations. Hitting the blade further from its axis of rotation (and thus from its constrained end) increased the lever of the increased axial force. At the same time, due to the tapering of the blade ( Figure 5), a smaller area was affected by the pressure increase and thus a smaller axial force increase was obtained. As shown in Figure 7, hitting the blade at 40 m (80%) of its span led to the highest axial deformation, irrespective of the chosen amplitude. For this reason, all the gusts analyzed in the remainder of this work were positioned at a 40 m height from the axis of rotation of the turbine. It was also noticed that, despite the blade always being hit by the gust at a 90° azimuth angle, a higher delay in its peak deflection was obtained when the gust was imposed further from its tip as a result of the higher portion of blade being displaced, and thus, of the higher inertia.
The axial force over the span of the blade was highly influenced by the wind gust. However, the differences on the lower 60% of the blade span (i.e., strips #1 to #6) were negligible, being smaller than 0.77% for = 5 m/s and smaller than 1.64% when = 10 m/s. Figure 8 shows the axial force evolution over each strip of the outboard 40% span of the blade, as well as the total axial force. Hitting the blade further from its axis of rotation (and thus from its constrained end) increased the lever of the increased axial force. At the same time, due to the tapering of the blade ( Figure 5), a smaller area was affected by the pressure increase and thus a smaller axial force increase was obtained. As shown in Figure 7, hitting the blade at 40 m (80%) of its span led to the highest axial deformation, irrespective of the chosen amplitude. For this reason, all the gusts analyzed in the remainder of this work were positioned at a 40 m height from the axis of rotation of the turbine. It was also noticed that, despite the blade always being hit by the gust at a 90 • azimuth angle, a higher delay in its peak deflection was obtained when the gust was imposed further from its tip as a result of the higher portion of blade being displaced, and thus, of the higher inertia.
The axial force over the span of the blade was highly influenced by the wind gust. However, the differences on the lower 60% of the blade span (i.e., strips #1 to #6) were negligible, being smaller than 0.77% for A g = 5 m/s and smaller than 1.64% when A g = 10 m/s. Figure 8 shows the axial force evolution over each strip of the outboard 40% span of the blade, as well as the total axial force.
of the higher portion of blade being displaced, and thus, of the higher inertia.
The axial force over the span of the blade was highly influenced by the wind gust. However, the differences on the lower 60% of the blade span (i.e., strips #1 to #6) were negligible, being smaller than 0.77% for = 5 m/s and smaller than 1.64% when = 10 m/s. Figure 8 shows the axial force evolution over each strip of the outboard 40% span of the blade, as well as the total axial force. The effect of the ABL is clearly seen in this figure, as well as in the ungusted condition, resulting in a higher loading when the blade pointed upward (90° azimuth angle) and lower loadings when it pointed downward (270° azimuth angle). The strips most sensibly influenced by the gust were #8 and #9 since the center of the gust was positioned exactly between these sections. These strips also provided an important contribution to the total bending moment, as they were located far from the axis of rotation. For both the amplitudes tested, the axial force over each strip resembled the distribution of velocity imposed for the gust, having a positive peak surrounded by two drops. However, it can be noted that the second drop in axial force was larger than the first one, especially for the case with = 10 m/s and on strips #8 and #9. This was due to the occurrence of flow separation in both cases, as illustrated in Figure 9. In this figure, the regions affected by separation were identified by marking the portions of the blade suction side where the tangential component of the wall shear stress was oriented according to the direction of rotation. The effect of the ABL is clearly seen in this figure, as well as in the ungusted condition, resulting in a higher loading when the blade pointed upward (90 • azimuth angle) and lower loadings when it pointed downward (270 • azimuth angle). The strips most sensibly influenced by the gust were #8 and #9 since the center of the gust was positioned exactly between these sections. These strips also provided an important contribution to the total bending moment, as they were located far from the axis of rotation. For both the amplitudes tested, the axial force over each strip resembled the distribution of velocity imposed for the gust, having a positive peak surrounded by two drops. However, it can be noted that the second drop in axial force was larger than the first one, especially for the case with A g = Energies 2020, 13, 509 11 of 20 10 m/s and on strips #8 and #9. This was due to the occurrence of flow separation in both cases, as illustrated in Figure 9. In this figure, the regions affected by separation were identified by marking the portions of the blade suction side where the tangential component of the wall shear stress was oriented according to the direction of rotation. On the root of the blade, a separation region was observed also when no gust was considered. In this region, the blade shape underwent a transition from a cylindrical root to an aerodynamically shaped body. For what concerns the outboard part of the blade, in the case with = 5 m/s, separation occurred only on a limited portion of the blade span and only on a restricted portion of the local chord length. On the other hand, when the gust amplitude was increased to 10 m/s, the separation area grew, expanding both in the spanwise and chordwise directions. This separation region was not reported during the first dip in the axial force. As soon as separation occurred, a sudden drop in the axial force was found. This was most intense on strip #9. As a result, the highest axial force was never reached at a 90° azimuth angle (Figure 8), but always a few degrees earlier. Furthermore, the second drop became longer in time and more intense in the axial force deficit with respect to the first one. This phenomenon strongly influenced the axial tip displacement, as shown in Figure 10. No sensible difference was observed in the tangential displacement because this was strictly related to gravity, as recognized in Santo et al. [6].  On the root of the blade, a separation region was observed also when no gust was considered. In this region, the blade shape underwent a transition from a cylindrical root to an aerodynamically shaped body. For what concerns the outboard part of the blade, in the case with A g = 5 m/s, separation occurred only on a limited portion of the blade span and only on a restricted portion of the local chord length. On the other hand, when the gust amplitude was increased to 10 m/s, the separation area grew, expanding both in the spanwise and chordwise directions. This separation region was not reported during the first dip in the axial force. As soon as separation occurred, a sudden drop in the axial force was found. This was most intense on strip #9. As a result, the highest axial force was never reached at a 90 • azimuth angle (Figure 8), but always a few degrees earlier. Furthermore, the second drop became longer in time and more intense in the axial force deficit with respect to the first one. This phenomenon strongly influenced the axial tip displacement, as shown in Figure 10. No sensible difference was observed in the tangential displacement because this was strictly related to gravity, as recognized in Santo et al. [6].
In this figure, the impact of the structural inertia was evident. The highest tip displacement was reached with a delay with respect to the axial force. Furthermore, the tip axial displacement immediately started to decrease when the outer border of the gust (where the velocity was decreased) impacted on the blade. Then, when the positive core of the gust hit the blade surface, the displacement started increasing again, and in the case with A g = 5 m/s, a higher tip displacement was achieved with respect to the ungusted configuration. In the case with A g = 10 m/s, the higher inertia of the blade due to its initial faster forward movement prevented the blade from reaching high peaks in its tip displacement, even if the axial force was increased, resulting in a maximum displacement lower than the case with A g = 5 m/s. In addition, when the flow separated (Figure 9), the tip displacement rapidly decreased, preventing the blade from reaching high displacements. This phenomenon was much more intense in the case with A g = 10 m/s since the area affected by the separation was larger (Figure 9), and consequently, the drop in the axial force was also more intense (Figure 8). When the blade moved out of the gust region, the tip showed an oscillatory motion, whose amplitude gradually decreased until the blade reached the (ungusted) regime condition again. occurred only on a limited portion of the blade span and only on a restricted portion of the local chord length. On the other hand, when the gust amplitude was increased to 10 m/s, the separation area grew, expanding both in the spanwise and chordwise directions. This separation region was not reported during the first dip in the axial force. As soon as separation occurred, a sudden drop in the axial force was found. This was most intense on strip #9. As a result, the highest axial force was never reached at a 90° azimuth angle (Figure 8), but always a few degrees earlier. Furthermore, the second drop became longer in time and more intense in the axial force deficit with respect to the first one. This phenomenon strongly influenced the axial tip displacement, as shown in Figure 10. No sensible difference was observed in the tangential displacement because this was strictly related to gravity, as recognized in Santo et al. [6].  A slightly different behavior was reported for the torque provided by the blade, as shown in Figure 11. Similarly to the axial force distribution, the torque differences over strips #1 to #6 did not exceed 1.43% for the A g = 5 m/s case and 3.22% for the A g = 10 m/s case.
Energies 2020, 13, 509 12 of 20 In this figure, the impact of the structural inertia was evident. The highest tip displacement was reached with a delay with respect to the axial force. Furthermore, the tip axial displacement immediately started to decrease when the outer border of the gust (where the velocity was decreased) impacted on the blade. Then, when the positive core of the gust hit the blade surface, the displacement started increasing again, and in the case with = 5 m/s, a higher tip displacement was achieved with respect to the ungusted configuration. In the case with = 10 m/s, the higher inertia of the blade due to its initial faster forward movement prevented the blade from reaching high peaks in its tip displacement, even if the axial force was increased, resulting in a maximum displacement lower than the case with = 5 m/s. In addition, when the flow separated (Figure 9), the tip displacement rapidly decreased, preventing the blade from reaching high displacements. This phenomenon was much more intense in the case with = 10 m/s since the area affected by the separation was larger (Figure 9), and consequently, the drop in the axial force was also more intense ( Figure 8). When the blade moved out of the gust region, the tip showed an oscillatory motion, whose amplitude gradually decreased until the blade reached the (ungusted) regime condition again.
A slightly different behavior was reported for the torque provided by the blade, as shown in Figure 11. Similarly to the axial force distribution, the torque differences over strips #1 to #6 did not exceed 1.43% for the = 5 m/s case and 3.22% for the = 10 m/s case. Similarly to what was observed for the axial force, the occurrence of separation was reflected in the fact that the highest peak in the torque was never reached at a 90° azimuth angle but always slightly in advance on strips #8 and #9. Differently, the difference in the two drops in torque was much smaller than what was observed for the axial force. This was due to the different effect of separation over the axial and tangential forces. When the flow angle increased, the tangential component of the lift and drag forces also increased, translating into an increase for the tangential force and a decrease for the axial force. At the same time, when separation occurred, the magnitude of lift and drag respectively decreased and increased, leading to a reduction for both the axial and tangential components (assuming that, as is typical, the lift-to-drag ratio was high enough). The two effects compensated for the tangential force and summed up for the axial force, leading to the observed difference in their dynamics. It is also remarked that the tip velocity induced by its transient deformation also sensibly affected the flow angle and the magnitude of the incoming relative velocity, as already observed in Santo et al. [6]. This was reflected in the regions around the 110° azimuth angle, where, on the most outboard strips, a higher torque contribution was found with respect to the ungusted condition as a consequence of the fast forward movement of the blade tip, which increased the incoming flow angle. Lastly, the total torque provided by the machine is provided in Figure 12. The peaks induced by the gust were comparable to the peaks induced by the tower-dam effect. It is also noted that a small effect was observed when the following blade went through the gusted Similarly to what was observed for the axial force, the occurrence of separation was reflected in the fact that the highest peak in the torque was never reached at a 90 • azimuth angle but always slightly in advance on strips #8 and #9. Differently, the difference in the two drops in torque was much smaller than what was observed for the axial force. This was due to the different effect of separation over the axial and tangential forces. When the flow angle increased, the tangential component of the lift and drag forces also increased, translating into an increase for the tangential force and a decrease for the axial force. At the same time, when separation occurred, the magnitude of lift and drag respectively decreased and increased, leading to a reduction for both the axial and tangential components (assuming that, as is typical, the lift-to-drag ratio was high enough). The two effects compensated for the tangential force and summed up for the axial force, leading to the observed difference in their dynamics. It is also remarked that the tip velocity induced by its transient deformation also sensibly affected the flow angle and the magnitude of the incoming relative velocity, as already observed in Santo et al. [6]. This was reflected in the regions around the 110 • azimuth angle, where, on the most outboard strips, a higher torque contribution was found with respect to the ungusted condition as a consequence of the fast forward movement of the blade tip, which increased the incoming flow angle. Lastly, the total torque provided by the machine is provided in Figure 12. Similarly to what was observed for the axial force, the occurrence of separation was reflected in the fact that the highest peak in the torque was never reached at a 90° azimuth angle but always slightly in advance on strips #8 and #9. Differently, the difference in the two drops in torque was much smaller than what was observed for the axial force. This was due to the different effect of separation over the axial and tangential forces. When the flow angle increased, the tangential component of the lift and drag forces also increased, translating into an increase for the tangential force and a decrease for the axial force. At the same time, when separation occurred, the magnitude of lift and drag respectively decreased and increased, leading to a reduction for both the axial and tangential components (assuming that, as is typical, the lift-to-drag ratio was high enough). The two effects compensated for the tangential force and summed up for the axial force, leading to the observed difference in their dynamics. It is also remarked that the tip velocity induced by its transient deformation also sensibly affected the flow angle and the magnitude of the incoming relative velocity, as already observed in Santo et al. [6]. This was reflected in the regions around the 110° azimuth angle, where, on the most outboard strips, a higher torque contribution was found with respect to the ungusted condition as a consequence of the fast forward movement of the blade tip, which increased the incoming flow angle. Lastly, the total torque provided by the machine is provided in Figure 12. The peaks induced by the gust were comparable to the peaks induced by the tower-dam effect. It is also noted that a small effect was observed when the following blade went through the gusted The peaks induced by the gust were comparable to the peaks induced by the tower-dam effect. It is also noted that a small effect was observed when the following blade went through the gusted Energies 2020, 13, 509 14 of 20 region, i.e., in the azimuthal range around 210 • . At this moment, the gust intensity in this region had lowered but had not disappeared.

The 1 + Cos Gust
In this section, a different gust shape function will be reported on, while the size and position of the gust is the same. The shape function reported on in this paragraph is given in Equation (11), where r indicates the non-dimensional radial coordinate in the gusted cylinder: This gust shape corresponds to the "extreme coherent gust" from the 61400-IEC, as used in similar works [26,29,44]. The present gust shape, when plugged into Equation (7), corresponds to a velocity increase in the whole region affected by the gust, contrary to the gust shape tested in the previous section. As already done for the zero net flow rate gust, two gust amplitudes will be used, namely 5 and 10 m/s. The total aerodynamic axial force and the contribution of the most outboard strips are illustrated in Figure 13.
Energies 2020, 13,509 14 of 20 region, i.e., in the azimuthal range around 210°. At this moment, the gust intensity in this region had lowered but had not disappeared.

The 1 + Cos Gust
In this section, a different gust shape function will be reported on, while the size and position of the gust is the same. The shape function reported on in this paragraph is given in Equation (11), where indicates the non-dimensional radial coordinate in the gusted cylinder: This gust shape corresponds to the "extreme coherent gust" from the 61400-IEC, as used in similar works [26,29,44]. The present gust shape, when plugged into Equation (7), corresponds to a velocity increase in the whole region affected by the gust, contrary to the gust shape tested in the previous section. As already done for the zero net flow rate gust, two gust amplitudes will be used, namely 5 and 10 m/s. The total aerodynamic axial force and the contribution of the most outboard strips are illustrated in Figure 13. Similarly to what was observed for the zero net flow rate gust, the maximum value in the axial force was always reached in advance of the 90° azimuth angle. Furthermore, a drop in the axial force was observed on each analyzed strip but not due to the gust shape function, indicating, also in this case, the occurrence of a flow separation, as depicted in Figure 14. Similarly to what was observed for the zero net flow rate gust, the maximum value in the axial force was always reached in advance of the 90 • azimuth angle. Furthermore, a drop in the axial force was observed on each analyzed strip but not due to the gust shape function, indicating, also in this case, the occurrence of a flow separation, as depicted in Figure 14. Similarly to what was observed for the zero net flow rate gust, the maximum value in the axial force was always reached in advance of the 90° azimuth angle. Furthermore, a drop in the axial force was observed on each analyzed strip but not due to the gust shape function, indicating, also in this case, the occurrence of a flow separation, as depicted in Figure 14. Strips #8 and #9 were the most affected by the imposed gusts. Flow separation induced the drops visible in Figure 13. In general, the area affected by the flow separation was broader compared to Figure 9, showing an expansion in both spanwise and chordwise directions. This was due to the more severe gust conditions imposed by the 1 + cos gust. As a result, the drop in axial force was more intense. In particular, when the gust amplitude was increased from 5 to 10 m/s, not only a larger separation region was obtained, but also a lower drop in the axial force. Importantly, the azimuthal range in which the axial force dropped below the ungusted case was wider (especially on strips #8, #9, and #10 in Figure 13) when the gust amplitude was doubled.  Figure 13. In general, the area affected by the flow separation was broader compared to Figure 9, showing an expansion in both spanwise and chordwise directions. This was due to the more severe gust conditions imposed by the 1 + cos gust. As a result, the drop in axial force was more intense. In particular, when the gust amplitude was increased from 5 to 10 m/s, not only a larger separation region was obtained, but also a lower drop in the axial force. Importantly, the azimuthal range in which the axial force dropped below the ungusted case was wider (especially on strips #8, #9, and #10 in Figure 13) when the gust amplitude was doubled.
The tip axial displacement was also heavily affected by the flow separation, as shown in Figure 15. It is remarked that, in the case with A g = 5 m/s, the highest tip axial displacement was obtained. Furthermore, in both cases, separation lowered the tip axial displacement in comparison with the ungusted configuration, acting as a protection mechanism against extreme deflection. This was clearly visible in the case with A g = 10 m/s, where separation prevented the blade from reaching high deflections and resulted in a maximum deflection lower than the case with A g = 5 m/s. In the latter case, the effect of the separation on the tip axial displacement was less evident and was restricted only to a marginal influence after a 120 • azimuth angle. When the blade surpassed the region affected by the gust, its tip continued oscillating around the ungusted deflection, until it went back into regime conditions by the end of the analyzed full revolution. The tip axial displacement was also heavily affected by the flow separation, as shown in Figure  15. It is remarked that, in the case with = 5 m/s, the highest tip axial displacement was obtained. Furthermore, in both cases, separation lowered the tip axial displacement in comparison with the ungusted configuration, acting as a protection mechanism against extreme deflection. This was clearly visible in the case with = 10 m/s, where separation prevented the blade from reaching high deflections and resulted in a maximum deflection lower than the case with = 5 m/s. In the latter case, the effect of the separation on the tip axial displacement was less evident and was restricted only to a marginal influence after a 120° azimuth angle. When the blade surpassed the region affected by the gust, its tip continued oscillating around the ungusted deflection, until it went back into regime conditions by the end of the analyzed full revolution.
Lastly, the torque provided by the blade hit by the gust is given in Figure 16. Lastly, the torque provided by the blade hit by the gust is given in Figure 16. It is remarked that, in the case with = 5 m/s, the highest tip axial displacement was obtained. Furthermore, in both cases, separation lowered the tip axial displacement in comparison with the ungusted configuration, acting as a protection mechanism against extreme deflection. This was clearly visible in the case with = 10 m/s, where separation prevented the blade from reaching high deflections and resulted in a maximum deflection lower than the case with = 5 m/s. In the latter case, the effect of the separation on the tip axial displacement was less evident and was restricted only to a marginal influence after a 120° azimuth angle. When the blade surpassed the region affected by the gust, its tip continued oscillating around the ungusted deflection, until it went back into regime conditions by the end of the analyzed full revolution.
Lastly, the torque provided by the blade hit by the gust is given in Figure 16. Similarly to what was explained in the previous section, the effect of separation on the torque was much less intense than on the axial force, resulting in very small drops below the ungusted condition. However, a sudden change in the trend was visible, especially on strips #8 and #9, around an 87° azimuth angle, due to the separation itself, which prevented the maximum torque being reached at a 90° azimuth angle. In the region downstream of the gust (from 110° azimuth angle onwards), the Similarly to what was explained in the previous section, the effect of separation on the torque was much less intense than on the axial force, resulting in very small drops below the ungusted condition. However, a sudden change in the trend was visible, especially on strips #8 and #9, around an 87 • azimuth angle, due to the separation itself, which prevented the maximum torque being reached at a 90 • azimuth angle. In the region downstream of the gust (from 110 • azimuth angle onwards), the differences observed in the torque, especially in the case with the highest gust amplitude, were attributed to the fast tip movement (Figure 15), which sensibly acted on the incoming flow angle and relative velocity magnitude.

Conclusions
In this work, the effect of a wind gust of various shapes and intensities on a modern horizontal axis wind turbine was investigated. A detailed and high-fidelity aeroelastic model was employed, implicitly coupling a computational fluid dynamics (CFD) solver based on an overset technique and a computational structural mechanics (CSM) solver, loyally reproducing the characteristics of the composite material of each blade. The gusts were superimposed on the atmospheric boundary layer.
First, the effect of a wind gust introducing a zero net flow rate was analyzed. This gust was conceptually similar to the "extreme operating gust" from the 61400-IEC standards but introduced a higher velocity increase. Results showed that an initial decrease in the blade loads and displacement was a consequence of the negative velocity increase imposed on the border of the gust. When the positive core of the gust impacted on the blade, the inertia of the structure caused a delay in the tip movement. Furthermore, despite the high peak reached by the aerodynamic axial force on the blade, flow separation over the span affected by the gust prevented the blade from reaching extreme deflections. Increasing the gust intensity, this protective effect was magnified by the broader area affected by the separation.
Subsequently, a gust analogous to the "extreme coherent gust" from the 61400-IEC standards was tested, introducing a consistently positive velocity increase. In this case, the peak tip deflection was shown to be higher than in the previous case as a consequence of the more severe wind conditions. Flow separation was also observed and affected a broader portion of the blade suction side, resulting in a fast reaction of the blade, whose tip underwent a fast axial movement. It was therefore concluded that, for each gust tested, flow separation acted as a protection mechanism and prevented the blade from reaching extreme deflections. For all the tested gusts, different dynamics were observed for the torque and axial force, especially when separation occurred.
Lastly, it can be concluded that the presented methodology allowed for the detailed investigation of the interaction between the blade structural response and the occurring wind gust: such data can be useful in the design stage. Among others, results can be used to better estimate the loading of the blade with respect to the meteorological data about the frequency, size, and intensity of the wind gusts in the site selected for the installation of the analyzed wind turbine.