**3. Numerical Modelling**

A 2D numerical model of the DU-06-W-200 airfoil was implemented in Ansys-FLUENT® v2020 to obtain the aerodynamic coefficients numerically. The Reynolds-averaged Navier–Stokes equations were resolved in an incompressible approach using eddy–viscosity turbulence modelling. Both Spalart–Allmaras (S-A) and new generalized k-ω (GEKO) models were used. The one-equation S-A model [17] is widely used for external aerodynamic applications. Although it is known to provide reasonable solutions for flows with adverse pressure gradients and separation, its accuracy to predict separation is lower than optimal two-equation models such as k-ω omega SST and GEKO. In addition, all k-ω models in Ansys are implemented with a *y*+-insensitive wall treatment, avoiding the discussion concerning the optimal selection of wall formulations in k-ε models [18].

GEKO is a recent turbulence model framework (based on the ω-equation) which introduces free parameters into the equations. The main advantage is that relevant parameters can be decided and tuned by the user for given operative ranges, and without a negative impact on the basic model calibration. The main tuning parameter for the GEKO model is the coefficient CSEP, which controls the boundary layer separation, predicting a more aggressive detachment if its value is increased. In the case of airfoils, it is highly recommended to use a value between 2.0 and 2.5 [19]. Furthermore, the GEKO model has been executed also with the option for scale-adaptive simulation (SAS) activated, which deploys an improved URANS formulation for the resolution of the turbulent spectrum in unstable flow conditions. The SAS concept is based on the introduction of the von Kármán length scale into the turbulence scale equation, allowing the model to dynamically adjust to resolved structures in a URANS simulation, which results in an LES-like behavior in unsteady regions of the flow field (those with flow separation).

An extended domain, with a distance to the inlet of 12.5*c* and a distance to the outlet equal to 20*c* (domain size 32.5*c* × 25*c*), in line with typical values found in the literature, was considered accurate to avoid the effect of the boundaries on the development of the flow inside the domain region (see Figure 4). A C-mesh distribution has been employed around the airfoil, resulting in a [350 × 75] cell size for both pressure and suction sides of the airfoil. An averaged value of *y +=* 1.7 (at *Rec* = 200,000) has been achieved with the first mesh point located at roughly 0.05 mm from the wall. At the wake region, a structured mesh of [300 × 150] cells was also employed, resulting in 97,500 cells for the complete 2D model. Furthermore, an additional refined mesh with [525 × 150] nodes on the airfoil walls and 247,500 cells for the whole domain was also employed to check the solution sensitivity to the grid resolution.

The boundary conditions of the simulation domain are given in Figure 4 which includes details of the adopted mesh. A velocity inlet condition of 16.4 m/s was set at the domain inlet to match the Reynolds number (200,000) of the experimental measurements. Furthermore, up to 21 different angles of attack (*AoA*) were simulated to complete a detailed evolution of the aerodynamic coefficients, including negative and positive incidences: ±[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]. According to previous measurements, a turbulence intensity of 0.7% was fixed for a length scale one order of magnitude lower than the characteristic size of the test section in the wind tunnel. Both steady and unsteady simulations have been conducted, the latter necessary at high *AoA* for partially and fully detached flow conditions. A time-step size of 3 × <sup>10</sup>−<sup>4</sup> s was fixed in order to track the evolution of the vortex shedding with sufficient resolution. A time-averaged value of the airfoil coefficients was finally computed after periodically fluctuating regimes were achieved (typically, 50 times the airfoil chord flow-time).

The flow equations were discretized using the finite volume method with a secondorder scheme for momentum and turbulent variables. Second-order accuracy was also selected in the transport equation for the pressure correction. The discretization of the temporal terms (when necessary) corresponds to a bounded second-order implicit formulation. The SIMPLE algorithm was used for the pressure–velocity coupling for all studied cases. Spatial discretization regarding gradient terms was selected to be the least-squares

cell-based discretization. Finally, a convergence criterion of 10−<sup>6</sup> was fixed for the velocity components of the momentum equation, while a minimum threshold of 10−<sup>5</sup> was at least required for the rest of the equations. Simulations were performed using a four-node Intel Core i7-52820K at 3.3 GHz and 64 Gb RAM, with characteristic CPU times of 75 min for every execution (1 day of CPU time to complete the whole angular range) in the case of the refined mesh.

**Figure 4.** CFD simulation mesh detail views and boundary conditions.
