*3.1. Numerical CFD Simulations*

The computational domain adopted for the CFD simulations is depicted in Figure 6A. The domain was made in the conventional bullet shape, having a distance of 20 and 40 chords upstream and downstream of the airfoil, respectively. The dimensions of the bullet were in agreement with the most conservative suggestions that can be found in the literature. The choice of the bullet-shaped domain was due to the possibility of imposing only an inlet (at which values of the velocity vector, turbulence intensity, and turbulent viscosity ratio were assumed) and an outlet boundary (at which the value of the gauge pressure was assumed, and for all other quantities zero normal derivative was assumed), with benefits in terms of calculation stability. The same approach was followed successfully in [16]. The works by Balduzzi et al. [20] were taken as the main references in order to select the most suitable numerical settings for the solver. For the sake of completeness, an overview on the main settings of the simulation models is given in the following section.

**Figure 6.** (**A**) Computational domain; (**B**) detailed view of the mesh at the trailing edge of the smooth airfoil (**left**) and near the GF (**right**).

The commercial code ANSYS® FLUENT® was used in the two-dimensional form to solve the time-dependent unsteady Reynolds-averaged Navier–Stokes (URANS) equations in a pressure-based formulation. The fluid was air, modeled as an ideal compressible gas. Based on previous experience, the validation against experiments made use of the four equations Transition SST model. This choice was due to the very low Reynolds number achieved in experiments (max 180 k), which provoked a massive impact of transition phenomena. Conversely, in the sensitivity analysis based on the realistic conditions of Table 2, the authors decided to achieve the turbulence closure by means of Menter's SST *k-*ω model [21], which is a blend of *k*-ε and *k-*ω two-equation formulations. This was due to pretty higher Reynolds numbers, which made the transitional effects less relevant. Moreover, the same study cases were originally simulated with this turbulence model, which proved to be particularly effective and robust. The Coupled algorithm (non-segregated) was employed, where the Navier–Stokes equations set is directly solved through an implicit discretization of pressure in the momentum equations, with benefits in terms of robustness and convergence, as shown in [20]. The second order upwind scheme was used for the spatial discretization of the whole set of URANS and turbulence equations, as well as the bounded second order for time differencing to obtain a good resolution. To allow for the pitching movement of the airfoil, the domain was split into rotating and stationary parts. The interface was circular, having a diameter of 14C. To handle the coupling of the two domains, a general grid interface (GGI) was used. The computational mesh generated for the two domains was of the unstructured type, made with the native mesher of the ANSYS® package. Triangular elements were used throughout the domain, with a massive element refinement within the rotating region and an additional local refinement area around the GF, as shown in Figure 6B. In order to properly capture the flow behavior within the boundary layer, a 30-layer O-grid with prismatic elements was instead created around the airfoil and the GF. The first element height was always chosen so as to guarantee a dimensionless wall distance (*y*+) at the grid nodes of the first layer above the blade wall constantly lower than 1. According to the prescriptions of [22], the expansion ratio for the growth of elements

starting from the surface was kept below 1.05 to achieve good mesh quality. A mesh dependency study was carried out to ensure that the results were not affected by mesh density and quality. Regarding the grid created for the validation of experiments, please refer to [17] for any additional details. For the spatial discretization of the airfoils with GF three different meshes were created. The *y*<sup>+</sup> was lower than 1, thus the meshes differed by the number of elements along the airfoil surface and their size in the free stream area. Created meshes had around 160,000, 280,000, and 400,000 elements, respectively. The mesh sensitivity analysis was conducted for the NACA0021 airfoil with Gurney flap lengths equal to 1%, 2%, and 3% of the chord length, respectively (these values are indicated in Figure 7A–C as *GF1*, *GF2*, *GF3*). Different AoAs were tested, even though Figure 7 only reports the case at AoA = 0◦ at Re = 250 K for brevity.

**Figure 7.** Mesh quality influence on (**A**) lift, (**B**) drag, and (**C**) moment coefficients.

It can be observed that asymptotic convergence is reached as the mesh is refined. The mesh with a number of elements equal to 400,000 was considered as not affected by the mesh size and further computations were carried out using this mesh. The values of the coefficients presented in Figure 7 have been averaged over 500 time steps after getting a converged solution. As extensively discussed in [22], the influence of the time step size also has to be considered carefully in order to obtain the desired accuracy of the computational model. In order to perform reliable dependence studies of temporal discretization, both the one-side GF and the fish tail GF were tested, since they were thought to induce a quite different vortex shedding at the trailing edge for low AoAs (see [6]), thus leading to different characteristic Strouhal numbers. Tests were carried out using a time-step of 10−3, 10−4, and 10−<sup>5</sup> s, respectively. Upon examination of the results, it was apparent that a time-step of 10−<sup>4</sup> s was sufficient to achieve independent results (differences in the absolute value lower than 0.1%).
