*2.1. The CFD Model*

On the CFD side of the employed FSI model, the domain illustrated in Figure 1 was used.

**Figure 1.** Overview and sections of the CFD computational domain of the background for the ABL with indications of the used boundary conditions.

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

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*<sup>+</sup> 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.

**Figure 2.** (**left**) Overview of the body-fitted mesh around blades and supporting structures with the overset borders in red, and (**right**) sections of the mesh around each blade, taken at different spanwise positions.

The air flow was modelled as incompressible (air density constant equal to 1.225 kg/m3) 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):

$$
\mu(z) = \frac{v\_\*}{K} \ln \left( \frac{z + z\_0}{z\_0} \right) \tag{1}
$$

$$k = \frac{v\_\*^2}{\sqrt{\mathbb{C}\_{\mu}}},\tag{2}$$

$$\kappa(z) = \frac{v\_\ast^3}{K(z+z\_0)}.\tag{3}$$

In these equations, *K* and *C*<sup>μ</sup> are the von Karman constant (0.4187) and a constant of the *k*-ε model (0.09), respectively. Moreover, *v*<sup>∗</sup> is the friction velocity (a global index of the wind intensity), and *z*<sup>0</sup> 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*<sup>+</sup> *mod* and of the wall function constant *Emod*:

$$z\_{mod}^{+} = \frac{(z + z\_0)u\_f \rho}{\mu},\tag{4}$$

$$E\_{mod} = \frac{\mu}{\rho \mathbb{Z}\_0 \mu\_f}.\tag{5}$$

In these equations, <sup>ρ</sup> and <sup>μ</sup> are the air density (1.225 kg/m3) and its viscosity (1.7894 <sup>×</sup> <sup>10</sup>−<sup>5</sup> 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 m2/s2, 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.

## *2.2. 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.

**Figure 3.** Overview of the structural mesh in the outer shells and in the shear webs of each blade, with a detailed view of the composite layup in one of the 64,000 shell elements.
