*3.2. Pitching Movements and Conventions*

The cycloidal motion of each blade on board a Darrieus VAWT generates a continuous variation of the angle of attack on the blade itself as a function of the relative positioning of the chord and the oncoming wind. This, in turn, leads to a variable intensity of the relative speed and to discontinuous forces exerted by the airfoils.

For the conventions used in the study, please refer to Figure 8. The overall torque *T* produced by the blade is given by Equation (2), where *L* and *D* are the lift and drag forces, respectively, and *M* it the moment around the blade-spoke connection point, which does not always correspond to the aerodynamic centre [23].

$$T = (L\sin\alpha - D\cos\alpha)R + M\tag{2}$$

**Figure 8.** Forces acting on the turbine blade.

The dependency of the AoA on the relative positioning between the blade and the wind, as well as on the force really exerted by the blade, which induces a reduction of the oncoming wind, lead to the well-known variation trends of the AoA that are non-symmetrical between the upwind and the downwind halves of the revolution and characterized by very steep variation rates. As a result, recent research works (e.g., [17]) showed that simulating a blade in pure pitching motion is not sufficient to give reliable estimation of the functioning of cycloidal motion.

To this end, in the present study, the AoA trends depicted in Figure 4 were applied directly to the blades. An average value of the Reynolds number (calculated based on the actual ones taken from CFD) was imposed at the inlet boundary, while the actual forces (lift, drag, and pitching motion) in motion were reconstructed point-by-point by means of the relative speed value also presented in the same figure, in order to have a more precise estimation. By doing so, the variation of the airfoil performance with the Reynolds number is unfortunately neglected, but this cannot be avoided in the case of a pitching approach like the one presented here. However, since each simulation is carried out with each specific equivalent TSR, this variation is thought to be sufficiently small to not compromise the accuracy of the results. In order to compare airfoils and turbines working in different conditions, the introduction of dimensionless coefficients is needed. The main coefficients used in the following of the study are the torque coefficient (Equation (3)) and the power coefficient (Equation (4)):

$$C\_T = \frac{T}{\frac{1}{2}\rho c^2 V\_\infty^2} \tag{3}$$

$$\mathcal{L}\_p = \frac{P}{\frac{1}{2}\rho A V\_{\infty}^3} \tag{4}$$

#### *3.3. Radial Basis Functions (RBFs) Interpolation for Data Reduction*

The parametric analysis carried out on the airfoil performance with different GFs provided only a finite set of points. In order to obtain a continuous response surface, the multivariate radial basis functions (RBFs) interpolation was applied. This method is very efficient for interpolation of large scattered data sets. It also has some drawbacks connected with unstable solutions and fast growth of the computational cost for large data series and also non-negligible error connected with Runge's phenomenon at the boundary of the domain. Drawbacks notwithstanding, it seems to be a suitable choice for the considered case of data reduction [24].

A multivariate function Φ is called radial (RBF) if its value at each point depends only on its distance from the base point, what is written in mathematical notation as:

$$\Phi(r) = \Phi(\|r - r\_0\|)\tag{5}$$

where · is the Euclidian norm in the *R<sup>n</sup>* space and *r*<sup>0</sup> is the vector describing the position of the base point. The radial function based interpolation takes the form of a linear combination of base functions attached to all *N* collocation nodes giving following equation:

$$\mathbf{u}(r) = \sum\_{i=1}^{N} \alpha\_i \Phi(\|r - r\_0\|) \tag{6}$$

where α<sup>i</sup> is an unknown interpolation coefficient. The values of interpolation coefficients can be found by collocating the interpolation function of Equation (6) and then solving the resultant linear set of equations which can be briefly written as:

$$\alpha\_{\mathbf{i}} = \Phi^{-1} \cdot \mathbf{u}(r) \tag{7}$$

where the interpolation (or collocation) matrix is computed as:

$$\Phi = \begin{Bmatrix} \Phi\_{\bar{i}\bar{j}} \end{Bmatrix} = \begin{Bmatrix} \Phi(\|r\_{\bar{i}} - r\_{\bar{j}}\|) \end{Bmatrix} \text{where } \bar{i}, \bar{j} = 1, \dots, N \tag{8}$$

The radial function form has to be chosen adequately with respect to the considered problem. Thus, to find the most suitable function, a dedicated analysis was carried out for different common types of radial functions. The functions are shown in Table 3.

**Table 3.** Utilized radial basis functions.


The shape parameter *c*, which appears in the definition of different radial basis functions, is a parameter that controls the shape of the basis function and hence the size of the region of influence of a given basis function around the collocation point. The higher the value of shape parameter, the bigger is the region of influence of the basis function; unfortunately, this also causes deterioration of conditioning of the approximation problem.

It needs to be pointed out that the interpolation matrix Φ is symmetric and positively defined, hence it is always invertible. However, by incrementing the number of nodes and consequently the number of base functions, the matrix conditioning becomes worse. The condition number of a matrix measures error is introduced by the finite arithmetic of computations on computers [25]. The matrix condition number for inversion is given by following equation:

$$(\Phi) = \|\Phi\| \|\Phi^{-1}\|\tag{9}$$

where is the conditioning number. In case of an interpolation based on the RBFs, the conditioning number value of the interpolation matrix is strictly connected with the shape of applied interpolation function. In the case of the functions presented in Table 3, this is controlled by the shape parameter [26]. For every interpolation problem it is possible to find optimal value of the shape parameter value; see [26]. To this end, two additional simulations, regarding the geometrical parameters of the airfoil thickness and the Gurney flap length, which had not been covered during the case studies, were done for each interpolation data set. Further, based on the additional simulations, the root main square (RMS) method was used to assess the interpolation accuracy, which is given by the following equation:

$$RMS = \sqrt{\frac{\sum\_{j=1}^{N} \left(\overline{f}(p\_j) - f(p\_j)\right)}{N}} \tag{10}$$

where *p* is the vector of points, which in the interpolation corresponds with the collocation points, *f* indicates the interpolation function, and *f* is the numerical values of the original function. Figure 9 shows the values for the matrix of conditioning number and resulting RMS error as a function of the shape parameter.

**Figure 9.** Shape parameter influence on: (**A**) matrix conditioning number, (**B**) root mean square (RMS) error, (**C**) maximum error.

#### **4. Results and Discussions**

## *4.1. Experimental Assessment of the Numerical Model*

As discussed, an extensive validation of the numerical model was carried out preliminary upon comparison with experimental data from the Hermann Föttinger Institute (HFI) of the TU Berlin [17].

The test case was a NACA0021 airfoil under a Reynolds number of 180 k. The numerical model was tested both in terms of static polars and in dynamic pitching motion. For the sake of brevity, only the results with a one-sided GF (2.5%c) are reported here: Figure 10 displays the comparison of static polars and Figure 11 the performance in sinusoidal movement within an AoA range between 10 and 30 deg and a reduced frequency *k* of 0.05.

$$k = \frac{\omega}{\mu\_0} \cdot \frac{c}{2} = \frac{\pi f c}{\mu\_0} \tag{11}$$

**Figure 10.** Comparison experimental results and CFD computations for static flow around airfoil equipped with GF at different values of angle of attack.

**Figure 11.** Comparison experimental results and CFD computations for dynamically changed angle of attack.

It has to be remarked that the numerical data in Figure 10 were run in unsteady RANS mode for several incidence angles (see Figure 12). Figure 12, in particular, shows the contours of normalized velocity, i.e., the local velocity divided by the undisturbed one. In particular, as common practice, the use of unsteady simulations was mandatory for high AoAs after a stall, where the analysis of the residuals showed issues in convergence history, characterized by the intense fluctuations of the calculated aerodynamic forces already shown by [6].

**Figure 12.** Normalized velocity contours for 0◦, 2◦, 4◦ and 12◦, 14◦, 16◦ angles of attack.

In thte case of GFs, the use of unsteady simulations was also necessary for very low AoAs (between 0◦ and 4◦), where the GF itself produces some vortices that detach alternatively from the corners and then are convected downstream. In those cases, the timestep for the simulation was based on the previous experience of [6]. Upon examination of both comparisons, sound agreement can be noted overall between experiments and simulations, proving the effectiveness of the method.
