A Parametric Study of Trailing Edge Flap Implementation on Three Di ﬀ erent Airfoils Through an Artiﬁcial Neuronal Network

: Trailing edge ﬂaps (TEFs) are high-lift devices that generate changes in the lift and drag coe ﬃ cients of an airfoil. A large number of 2D simulations are performed in this study, in order to measure these changes in aerodynamic coe ﬃ cients and to analyze them for a given Reynolds number. Three di ﬀ erent airfoils, namely NACA 0012, NACA 64(3)-618, and S810, are studied in relation to three combinations of the following parameters: angle of attack, ﬂap angle (deﬂection), and ﬂaplength. Results are in concordance with the aerodynamic results expected when studying a TEF on an airfoil, showing the e ﬀ ect exerted by the three parameters on both aerodynamic coe ﬃ cients lift and drag. Depending on whether the airfoil ﬂap is deployed on either the pressure zone or the suction zone, the lift-to-drag ratio, C L / C D , will increase or decrease, respectively. Besides, the use of a larger ﬂap length will increase the higher values and decrease the lower values of the C L / C D ratio. In addition, an artiﬁcial neural network (ANN) based prediction model for aerodynamic forces was built through the results obtained from the research.


Introduction
Wind turbine sizes have steadily been increasing over the past few years. Commercial offshore wind turbines with a maximum capacity of up to 6 MW are operational in the U.S., Europe, and China. Indeed, 10 MW turbines with diameters of 144 m and even 20 MW turbines with diameters of 240 m are currently under development, as Musial et al. [1] and Fichaux [2] et al. have reported. Barlas et al. [3] stated that one main focus is on developing new technologies capable of considerably reducing ultimate and fatigue loads in wind turbines. If an innovative blade design can result in a decrease in loading, Veers et al. [4] concluded in their study that the general cost will decrease, as rotor loads affect the loading of other components, such as the drive train and the tower.
According to the study of Yu et al. [5], two load control methods are widely used for commercial offshore wind turbines: collective pitch control and relatively advanced individual pitch control.
is deployed on the pressure side, the lift increases, and if the flap is deployed on the suction side, the lift decreases.
TEFs can be employed in two ways: discrete flaps require a moment on the hinge to achieve the required position, while flexible flaps present a smoother shape between the device and the airfoil, which increases their efficiency. The scheme of a trailing edge flap installed on an airfoil is shown in Figure 1.  Lackner M. et al. [25] presented the benefits of using flow control devices such as TEFs. Basualdo [26] showed that the use of variable geometry airfoils in wind turbine blades can lead to load alleviation. Andersen et al. [27] report the potential of flaps to alleviate 34% of the fatigue equivalent damage in flapwise loading, while Barlas et al. [28] find slightly lower values up to 27%. The success of the static angle testing of the modular blade flaps and the ability of the system to measure moment changes within the blade introduces numerous possibilities for further research [20].
A parametric study of trailing edge flap implementation on three different airfoils through an artificial neuronal network is presented in this article.
There are two main contributions made by this study: On the one hand the high number of simulations and different parameters combinations performed for three different airfoils, on the other hand, the use of an artificial neural network.
As shown in the methodology section, the large number of simulations performed is due to the high number of specifically designed combinations that have been carried out between different parameters (flap angle β and flaplength L F ). In addition, as it has also been specified in the same section, this research increases the analysis ratio of some parameters with respect to the state of the art.
For different angles of attack, the relationship between aerodynamic forces, flap position and flaplength has been presented for each airfoil using an original three-dimensional surface. This research, apart from adding new knowledge about the effects of TEFs applied to three airfoils sections, has compared the results corresponding to the different airfoils in each step.
An artificial neural network-based prediction model for aerodynamic forces built through the results obtained from the research is presented in this study. The relevance of this type of mathematical solutions for control tools has been increased for the last decade.
Finally, in order to consolidate the results of this study, these are complemented with data on the near wake region and pressure.

Aims and Methodology
Aiming to optimize airfoil aerodynamic performance, the purpose of this study is a numerical assessment of the aerodynamic characteristics of airfoil sections, specifically intended for wind turbine applications, equipped with TEFs. For this, a parametric study of three different airfoils, NACA 0012, NACA 64(3)-618, and S810, with TEFs installed has been carried out. Then, with the aim of facilitate future computational control tools for TEFs and aerodynamic predictions tools, all the results have been used to build a mathematical prediction model through an artificial neural network which is more detailed in Section 4.4.
The airfoils used, are part of two different airfoil families. On the one hand, those corresponding to the National Advisory Committee of Aeronautics (NACA) airfoils, NACA 0012 and NACA 64 (3)-618, and on the other hand a wind turbine airfoil from the National Renewable Energy Laboratory (NREL), S810. Figure 2 illustrates precisely the outline of the three aerodynamic profiles: The NACA 64(3)-618 however, belongs to the "6 series" of the NACA family. 6 denotes the series and indicates that this group is designed for greater laminar flow than the 4-digit series. The second digit, 4, is the location of the minimum pressure in tenths of chord (0.4c). The number in parentheses 3 indicates that low drag is maintained at lift coefficients 0.3 above and below the design lift coefficient (0.6) specified by the first digit after the dash in tenths. The final two digits specify the thickness in percentage of chord, 18% (see [31]).
In order to measure aerodynamic coefficients, and later, analyze them, a large number of 2D simulations are performed in this study.
TEFs actuate at different angles as flow control devices on airfoil sections and are capable of changing their aerodynamic characteristics significantly. As Menon et al. [33] explained, the increase or decrease in lift depends on the relative angle between airfoil and flap sections (actuation angle, β) and the angle of relative wind (angle of attack, α).
This study, apart from the effect of the angles α and β, will also analyze the effect of different flaplengths (LF) in the change of lift.
The NACA 64(3)-618 however, belongs to the "6 series" of the NACA family. 6 denotes the series and indicates that this group is designed for greater laminar flow than the 4-digit series. The second digit, 4, is the location of the minimum pressure in tenths of chord (0.4c). The number in parentheses 3 indicates that low drag is maintained at lift coefficients 0.3 above and below the design lift coefficient (0.6) specified by the first digit after the dash in tenths. The final two digits specify the thickness in percentage of chord, 18% (see [31]).
In order to measure aerodynamic coefficients, and later, analyze them, a large number of 2D simulations are performed in this study.
TEFs actuate at different angles as flow control devices on airfoil sections and are capable of changing their aerodynamic characteristics significantly. As Menon et al. [33] explained, the increase or decrease in lift depends on the relative angle between airfoil and flap sections (actuation angle, β) and the angle of relative wind (angle of attack, α).
This study, apart from the effect of the angles α and β, will also analyze the effect of different flaplengths (L F ) in the change of lift.
Procedure sketch showing the different combinations between flaplength, flap angle, and angle of attack, is shown in Figure 3.
Therefore, for each airfoil, a total number of 935 2D simulations were performed as a function of the angle of attack, the flap angle, and the flaplength. The three different airfoils NACA 0012, NACA 64(3)-618, and S810 have five flaplengths of 8%, 9%, 10%, 12%, and 14% of the airfoil chord c that are combined with eleven flap angles from β = −5° to β = 5°, in one degree steps, yielding a total of 55 meshes. The simulations were performed for seventeen angles of attack from α = −6° to α = 10°, and in one-degree steps.
Procedure sketch showing the different combinations between flaplength, flap angle, and angle of attack, is shown in Figure 3.  In the course of the European UPWIND project at the Institute of Aerodynamics and Gas Dynamics (IAG), Lutz et al. [34], analyzed an airfoil specifically designed for load alleviation purposes equipped with 10% of c trailing edge flap. Moreover, as proposed by Wei Jun Zhu [35], the trailing edge flaps accounts for approximately 10% of the airfoil chord length. With the purpose of adding new insight of the effect of different flaplengths in the change of the lift and drag, this study chose a range of different flap lengths close to the 10% of c.
The cases studied by Menon et al. [33] in a similar research, involved three different angles of actuation of the flap: −5 • , 0 • and 5 • . However, this study pretends to extend the previous study analyzing all the β angles between -5 • and 5 • , in one degree steps, so that the mapping charts shown later in the results section are supported by a greater number of points and are, therefore, of a higher quality.
The range of α angles (−6 • till to +10 • in one degree steps) used for the simulations of this study has been the necessary and sufficient to obtain the minimum and maximum values of the lift-to-drag ratio C L /C D , experiencing minimum and maximum figures respectively for these limits, as will be seen later in the results section.
According to Ju et al. [36], a high lift-to-drag ratio and a high lift coefficient are one of the principal aerodynamic requirements that a wind turbine airfoil should satisfy. So, this study has focused on the maximum C L /C D values, while stall regions will be studied in a following paper.
Computational simulations of the NACA 0012, the NACA 64(3)-618, and the S810 airfoils without a TEF were performed and validated against the experimental data provided by Douvi and Margaris [37], Abbot [19], and Reuss Ramsay et al. [38], respectively. These experimental data were obtained under similar boundary conditions to those of this study.

Computational Configuration and Procedure
CFD techniques have been employed, in order to obtain the aerodynamic features of the trailing edge flap (TEF) on an airfoil. Nowadays, non-commercial and proprietary CFD codes are used to reproduce most physical problems relatively well. In this work, the open source code OpenFOAM [39] was used to simulate the effects of different airfoil chords, depending on the TEF length with different flap angle positions installed on three different airfoils: NACA 0012, NACA 64(3)-618, and S810.
The potentialFoam solver which solves potential flow problems was used to generate starting fields, in order to speed up the convergence process. This solver is a suitable tool to generate the Symmetry 2020, 12, 828 6 of 30 initial conditions for more advanced solvers, such as simpleFoam, the one used in the present work. The simpleFoam solver has been applied for steady-state, incompressible and turbulent flows using the RANS (Reynolds averaged Navier-Stokes) equations. Throughout the calculations, this solver uses the k-ω SST (shear stress transport) turbulence model developed by Menter [40], due to its robust performance in detached flows, as reported by Kral [41] and Gatski [42]. The turbulence model is a combination of two models: Wilcox's k-ω model for near wall regions and the k-ε model for the outer region and in free shear flows. Second-order discretization schemes were employed in the computational simulations.
Python code was used for the generation of an optimized O-grid 2D mesh with a finite sharp trailing edge around the different types of airfoil under analysis. The coordinate points defining the airfoils, some from the airfoil coordinates database [43] and from Timmer et al. [44], were interpolated, in order to improve the surrounding mesh quality. The computational domain generated, shown in Figure 4, was a function of the O-grid radius, the flap angle and length, the airfoil geometry and the corresponding parameters such as chord and maximum thickness. The cell number in the chord-wise direction was around 382 and in the normal direction 400. The outer boundary, according to the criteria followed by some authors [45], can be 40 chords away from the airfoil.
For this study, the patches of particular interest are inlet (inflow), outlet (zero gradient), and wall (airfoil). For the current computations, the turbulence intensity has been set to 5%, and the mixing length has been estimated as a quarter of the airfoil chord length.
After having performed a sensitivity analysis, the number of cells of the O-grid 2D mesh was around 120,000. The python code output was fully compatible with the OpenFOAM blockMesh mesh utility, and the structured meshes were generated with proper stretching functions, yielding low expansion rates and dimensionless wall distance (y+ ≈ 1) for the Reynolds number under consideration. Figure 5 shows the mesh around the trailing edge flap for its two extreme deployment positions, installed on one of the airfoils in this study. Each position of the TEF requires the generation of a new mesh for each airfoil. As proposed in Wei Jun Zhu [35], the trailing flap is designed as a flexible part which bends its sides when it flaps. Fincham et al. [46] used a similar concept in their study. However, in our particular model shown in Figure 1, for simplification, the flap shape is not modified, it is considered as a rigid body. The sub-chord length trailing flap rotates an angle β as a rigid body for each case with similar merging conditions for all the flap installation angles β, as a positive direction of the flap motion is considered counter clockwise rotation. The cell number in the chord-wise direction was around 382 and in the normal direction 400. The outer boundary, according to the criteria followed by some authors [45], can be 40 chords away from the airfoil.
For this study, the patches of particular interest are inlet (inflow), outlet (zero gradient), and wall (airfoil).
The wall condition is a non-slip condition, constant zero velocity and no gradient pressure. The variables of particular interest which require boundary conditions (BC) are the fluid quantities velocity and pressure. The side patches of the domain have been defined as slip conditions.
For the current computations, the turbulence intensity has been set to 5%, and the mixing length has been estimated as a quarter of the airfoil chord length.
After having performed a sensitivity analysis, the number of cells of the O-grid 2D mesh was around 120,000. The python code output was fully compatible with the OpenFOAM blockMesh mesh utility, and the structured meshes were generated with proper stretching functions, yielding low expansion rates and dimensionless wall distance (y+ ≈ 1) for the Reynolds number under consideration. Figure 5 shows the mesh around the trailing edge flap for its two extreme deployment positions, installed on one of the airfoils in this study. Each position of the TEF requires the generation of a new mesh for each airfoil. As proposed in Wei Jun Zhu [35], the trailing flap is designed as a flexible part which bends its sides when it flaps. Fincham et al. [46] used a similar concept in their study. However, in our particular model shown in Figure 1, for simplification, the flap shape is not modified, it is considered as a rigid body. The sub-chord length trailing flap rotates an angle β as a rigid body for each case with similar merging conditions for all the flap installation angles β, as a positive direction of the flap motion is considered counter clockwise rotation. expansion rates and dimensionless wall distance (y+ ≈ 1) for the Reynolds number under consideration. Figure 5 shows the mesh around the trailing edge flap for its two extreme deployment positions, installed on one of the airfoils in this study. Each position of the TEF requires the generation of a new mesh for each airfoil. As proposed in Wei Jun Zhu [35], the trailing flap is designed as a flexible part which bends its sides when it flaps. Fincham et al. [46] used a similar concept in their study. However, in our particular model shown in Figure 1, for simplification, the flap shape is not modified, it is considered as a rigid body. The sub-chord length trailing flap rotates an angle β as a rigid body for each case with similar merging conditions for all the flap installation angles β, as a positive direction of the flap motion is considered counter clockwise rotation. The orthogonality and skewness of the mesh of the three airfoils have been studied for the three airfoils and for three flap angles β. We chose the extreme values β = −5 and +5 degrees and the one with β = 0 degrees. The Table 1 summarizes the results of both mesh quality parameters. A general rule is that the maximum skewness for a mesh in most flows should be kept below 0.95, which is the case in the present study as the maximum skewness is 0.87, for the NACA 64(3)-618 case (see Table 1): As Belamadi et al. [47] used in a similar and recent research, the simulations considered a Reynolds number of Re = 10 6 . The resulting flow speed is 15.11 m s −1 . The Reynolds number is defined in Equation (1): where ρ is the density, µ the absolute viscosity, c the reference length, which in the case of an airfoil is the chord length, and U ∞ the free stream velocity. Their values are shown in Table 2. The Reynolds number is known, so the high (y) that the first row of grid cells should have is y = 3.65 × 10 −5 m. Figure 6 shows the wall y+ distribution over the NACA0012 airfoil for two flap angle configurations β = 5 • and −5 • , angle of attack of α = 0 • and flaplength of 8% of c.

Reynolds Number [-] ρ [kg m −3 ] μ [kg m −1 s −1 ] c [m]
[m s −1 ] 10 6 1.225 1.8375 × 10 −5 1 15.11 The Reynolds number is known, so the high (y) that the first row of grid cells should have is y = 3.65 × 10 −5 m. Figure 6 shows the wall y+ distribution over the NACA0012 airfoil for two flap angle configurations β = 5° and −5°, angle of attack of α = 0° and flaplength of 8% of c. In addition, the wall y+ parameter of the computations of the three airfoils have also been studied for two flap angle configurations β = 5° and β=−5°. Table 3 summarizes the results obtained in all cases for the non-dimensional wall y+ parameter. Angle of attack with α = 0° and flaplength of 8% of c. In addition, the wall y+ parameter of the computations of the three airfoils have also been studied for two flap angle configurations β = 5 • and β = −5 • . Table 3 summarizes the results obtained in all cases for the non-dimensional wall y+ parameter. Angle of attack with α = 0 • and flaplength of 8% of c. Simulations were performed on a personal server-clustered parallel machine with Intel ® Core i7-6700 CPU 3.40 GHz × 8 cores and 32 GB RAM on 64 bits Linux. Each domain corresponding to an O-grid mesh was automatically divided into eight subdomains to be solved in parallel, thereby reducing the simulation time.
The C D on the airfoil has been calculated by the application of Equation (2) in the far field planes normal to the flow direction and applied to the streamwise velocity in different wake rakes (WR), according to Beans et al. [48] and Young [49].
where A is the area determined by the chord length and the span, A = b × c. In our cases, 2D simulations have been carried out, and the drag coefficient C D has been calculated per meter of airfoil span.
The lift coefficient C L was determined by the OpenFOAM inbuilt code, which determines the value by means of the Equation (3): Symmetry 2020, 12, 828 9 of 30

Validation
The validations are represented in two different types of plots for each airfoil: (a) coefficient of lift against angle of attack C L -α; and, (b) drag coefficient to lift coefficient C D -C L in Figures 7-9. Each pair of plots for each airfoil is represented by a different figure. Symmetry 2020, 12, x FOR PEER REVIEW 9 of 30 between experimental data and CFD does exist. In this last validation of airfoil S810, as in the first validation of NACA 0012, due to the high number of simulations performed and the limited volume of experimental data, the volume of CFD data is larger than the experimental data.
Symmetry 2020, 12, x FOR PEER REVIEW 9 of 30 between experimental data and CFD does exist. In this last validation of airfoil S810, as in the first validation of NACA 0012, due to the high number of simulations performed and the limited volume of experimental data, the volume of CFD data is larger than the experimental data.  Figure 7a corresponds to the validation performed for the symmetric airfoil NACA 0012. The C L -α curves show that the lift CFD data at a Reynolds number of Re = 10 6 are coincident with the experimental ones [37]. Figure 7b, corresponding to the representation of the C D -C L curves, shows an over estimation of the drag CFD data with respect to the experimental ones. In order to highlight the trend of CFD data more clearly, the volume of CFD data used for the validation of the NACA 0012 airfoil is higher than the volume of the experimental points.
In the case of the NACA 64(3)-618 airfoil, the CFD data were obtained at a Reynolds number of Re = 10 6 , but the experimental ones [19] were obtained at a Reynolds number of Re = 3 × 10 6 . This difference might explain an additional deviation in C D -C L curves. Therefore, Figure 8b shows that the CFD data for the drag are overestimated with respect to the experimental data. Figure 8a shows almost coincident C L -α curves for the CFD and the experimental data.  The following validation corresponds to airfoil S810, the CFD and the experimental data [38] of which are compared in Figure 9. In this case, CFD and experimental data were obtained at a Reynolds number of Re = 10 6 . Figure 9a shows coincident C L -α curves for the CFD and experimental data along the entire curve except for the last two points of the experimental curve. C D -C L curves of Figure 9b shows an over estimation of the CFD data with respect the CFD data. Therefore, a concurrency between experimental data and CFD does exist. In this last validation of airfoil S810, as in the first validation of NACA 0012, due to the high number of simulations performed and the limited volume of experimental data, the volume of CFD data is larger than the experimental data.
Taking into account the concurrency between CFD and experimental C L -α curves, and in spite of the minimum differences between CFD and experimental C D -C L curves, we can safely confirm the validation of the CFD model used in this study.

4.1.
Lift-to-Drag Ratio C L /C D as a Function of the Angle of Attack, α, for Different Flap Angles β Lift-to-drag ratio, C L /C D , as a function of the angle of attack, α, is used to show the results obtained from computational simulations. As explained in the previous section, simulations were performed for seventeen angles of attack, by combining eleven flap angles with five flaplengths. A flaplength of 8% of chord length c is chosen in this section. The curves corresponding to the rest of the flaplengths 9%, 10%, 12%, and 14% of c are included in Appendix A.
As explained in the introduction, the TEF concept is based on lift increase, when the flap is deployed on the pressure side (negative β angles), and lift decrease when the flap is deployed on the suction side (positive β angles) of the airfoil. It can therefore be deduced from the following Figures 10-12, corresponding to the profiles for NACA 0012, NACA 64(3)-618, and S810, respectively, that results are in concordance with the aerodynamic results expected when studying a TEF inserted on an airfoil. As the flap is deployed towards the pressure area of the airfoil, approaching higher values of negative β angles, the lift increases. As the flap is deployed towards the suction zone of the airfoil, approaching higher values of positive β angles, the lift decreases.
As shown in the following Figures 10-12 of the C L /C D curves, the C L /C D ratio of each airfoil for the same angle of attack, α, increases and it decreases as a function of the β angle in the same proportion with respect to the curve, β = 0 • .
For a flap angle of β = −5 • , the range of α angles of attack chosen for this study, allow obtaining a maximum lift-to-drag ratio C L /C D for the three airfoils studied. As shown in Figures 10-12        As the positive and negative angle of attack, α, increases, the drag increases in a larger proportion than the lift, in such a way that when reaching the highest positive angle of attack, α = 10 • , and at the highest negative angle of attack α = −6 • , the C L /C D ratio values for different β angles are almost overlapping. These near overlaps are shown in Figures 10 and 12, corresponding to airfoils NACA 0012 and S810, respectively, where the drag increases in the same proportion for α positive angles and for negative ones in the case of NACA 0012, and almost in the same proportion in the case of S810. However, as shown in Figure 11, drag increases substantially as the positive value of α increases for the NACA 64(3)-618 airfoil and less so when the negative value of α increases.
As shown in Figure 10    Each plot corresponds to one of the three mentioned angles β. For each β angle, the C L /C D ratio of four intermediate angles of attack, α = −3 • , α = 0 • , α = 3 • and α = 6 • , will be compared for different flaplengths of 8%, 9%, 10%, 12%, and 14% of c. Figure 13, corresponding to the symmetric airfoil NACA 0012, shows the variations of C L /C D of four angles of attack as a function of the different flap lengths for the angles β = 0 • , β = −5 • , and β = 5 • , respectively.
In Figure 13a, corresponding to the flap position β = 0 • , no variations of the C L /C D ratio are observed as a function of flaplength. This is because for β = 0 • , the variation of flaplength does not vary the geometry of the airfoil as the flap is an integrated part of it. However, Figure 13b corresponding to the position of flap β = −5°, clearly shows an upward trend of the CL/CD ratio for all of the four intermediate angles of attack as the flap length increases. This trend indicates that when the flap is deflected towards the lower surface in its extreme position, β = −5°, a position in which the highest lift force is obtained, the maximum CL/CD ratio is obtained with a 14% flap length. Figure 13c, corresponding to the position of flap β = 5°, shows a downward trend of the CL/CD ratio as the flap length increases. This downward trend indicates that when the flap is deflected towards the upper surface at its extreme position, β = 5°, a position in which the least lift force is obtained, the minimum CL/CD ratio is obtained with a flaplength of 14% of c. Therefore, a greater flaplength maximizes the lift in the positions in which the highest lift force is obtained, and minimizes the lift in the positions in which the lowest lift force is obtained.   Figure A4b, corresponding to the NACA 64(3)-618 airfoil, with a value of α = 6 • , very close to stall values, shows that when the flap is deflected towards the lower surface in its extreme position, β = −5 • , as the flap length increases, the C L /C D ratio decreases slightly, instead of increasing. Figure 14, corresponding to the symmetric airfoil NACA 0012, at four intermediate angles of attack, α = −3 • , α = 0 • , α = 3 • and α = 6 • , shows the C L /C D ratio variation as a function of the flap angle β and the different flaplengths in the form of a three-dimensional surface. As the flap is deployed towards the pressure area of the airfoil and negative β increases, the lift increases, and the C L /C D ratio also increases. If the flap is deployed towards the suction zone and positive β increases, the lift decreases, and the C L /C D ratio also decreases. In addition, a greater flaplength maximizes the C L /C D ratio in the positions in which the highest lift force is obtained, and minimizes the C L /C D ratio in the positions in which the lowest lift force is obtained. β and the different flaplengths in the form of a three-dimensional surface. As the flap is deployed towards the pressure area of the airfoil and negative β increases, the lift increases, and the CL/CD ratio also increases. If the flap is deployed towards the suction zone and positive β increases, the lift decreases, and the CL/CD ratio also decreases. In addition, a greater flaplength maximizes the CL/CD ratio in the positions in which the highest lift force is obtained, and minimizes the CL/CD ratio in the positions in which the lowest lift force is obtained. The same three-dimensional type of surfaces for NACA 64(3)-618 and S810 non-symmetrical airfoils are shown in Figures A6 and A7 of Appendix C. The drag increases at a value of the angle of attack, α, quite close to stall values and with the flap deflected towards the pressure zone, in such a way that, as the flaplength increases at these extreme α angles, the CL/CD ratio increases less, shows no increase at all, and can even decrease. The same three-dimensional type of surfaces for NACA 64(3)-618 and S810 non-symmetrical airfoils are shown in Figures A6 and A7 of Appendix C. The drag increases at a value of the angle of attack, α, quite close to stall values and with the flap deflected towards the pressure zone, in such a way that, as the flaplength increases at these extreme α angles, the C L /C D ratio increases less, shows no increase at all, and can even decrease.

C L /C D Ratio as a Function of the Flap Angle β and the Different Flap Lengths L F
Having studied the influence of the flap angle, β, and the flaplength, L F , on the C L /C D ratio for certain angles of attack, α, the next step will be to build a prediction model of aerodynamic forces taking into account all the parameters studied throughout this research. This model is based on an artificial neural network (ANN).

Modeling of the CFD Results with an Artificial Neural Network
Taking advantage of the results previously generated, the purpose of this section is to make a prediction model that could facilitate future computational control tools for TEFs and aerodynamic predictions tools. Bernhammer et al. [50] demonstrated the load reduction potential of smart rotors by designing an individual flap controller. In addition, the prediction model presented below could be an interesting contribution to the blade element momentum (BEM) theory. BEM theory with airfoil data is a widely used technique for prediction of wind turbine aerodynamic performance, being the reliability of the airfoil data an important factor for the prediction accuracy of aerodynamic forces and power, see Yang et al. [51].
An artificial neural network (ANN) has been designed to model the C L /C D ratio of the three different airfoils as a function of the flap angle β and the different flap lengths L F . Moreover, the ANN contains the model of the C L /C D ratio for different values of the angle of attack, α of the incoming airflow. This model was built using Matlab [52] mathematical software.
The selected topology for the ANN is the MultiLayer Perceptron with BackPropagation (MLP-BP), which is known for its good characteristics to model every surface. In the work of Aramendia et al. [53] a MLP-BP neural network is designed to model the aerodynamic behavior of a DU91W(2)250 airfoil with different length gurney flaps (GFs). A MLP-BP is presented by Saenz-Aguirre et al. [54] to store the data corresponding to the matrix Q(s,a) of a reinforcement learning algorithm for the yaw angle control of a wind turbine.
The MLP-BP designed in this paper presents one input layer with 3 neurons, corresponding to each one of the inputs (α, L F and β), one hidden layer with 25 neurons and one output layer with 1 neuron, corresponding to the output of the ANN (the C L /C D ratio). The training process of the ANN has been carried out with a data set of 220 tuples and the distribution of the data has been set to 90% for the training, 5% for the validation, and 5% for the test.
After the training, the values of the regression coefficient (R) in the test and the values for the mean squared error (MSE) have been obtained. Both are indicators of a correct training process and they are shown below in Table 4. The ANN of this study and the CFD results for the airfoil NACA 0012 are shown in Figure 15. Following the procedure of a similar study carried out by Urkiola et al. [55], the model has been represented as a surface. Black dots represent the CFD results. The same comparative representations for NACA 64(3)-618 and S810 non-symmetrical airfoils are shown in Figures A8 and A9 of Appendix D. Finally, with the purpose of adding more information to the results previously exposed, the results of the streamwise velocity u (ux due of α = 0°), the wake turbulence kinetic energy TKE and the pressure coefficient Cp have been studied for a select number of cases of the three airfoils NACA 0012, NACA 64(3)-618 and S810.

Streamwise Velocity
For the same angle of attack α = 0°, six different cases have been chosen for each airfoil. With the intention of covering the largest possible area of data, the minimum and maximum flaplengths of 8% of c and 14% of c with the three extreme flap positions, β = −5°, β = 0° and β = 5° have been combined.
The results corresponding to the streamwise velocity and wake turbulence kinetic energy have been checked in the following wake locations: x/c = 1.05, x/c = 1.25 and x/c = 1.5. Taking into account the concurrency between CFD results and the ANN, we can confirm that the ANN model could be a good option for future computational control tools for TEFs and aerodynamic predictions tools.
Finally, with the purpose of adding more information to the results previously exposed, the results of the streamwise velocity u (u x due of α = 0 • ), the wake turbulence kinetic energy TKE and the pressure coefficient Cp have been studied for a select number of cases of the three airfoils NACA 0012, NACA 64(3)-618 and S810.

Streamwise Velocity
For the same angle of attack α = 0 • , six different cases have been chosen for each airfoil. With the intention of covering the largest possible area of data, the minimum and maximum flaplengths of 8% of c and 14% of c with the three extreme flap positions, β = −5 • , β = 0 • and β = 5 • have been combined.
The results corresponding to the streamwise velocity and wake turbulence kinetic energy have been checked in the following wake locations: x/c = 1.05, x/c = 1.25 and x/c = 1.5.
The data corresponding to the NACA 0012 airfoil are shown below in this results section. The data corresponding to the airfoils NACA 64(3)-618 and S810 are shown in Appendices E and F.
The results of the streamwise velocity and the wake turbulence kinetic energy, shown in Figures 16 and 17

Turbulence Kinetic Energy
Then, as shown in Figure 17, the results of the normalized wake turbulence kinetic energy are displayed in the same format in which the streamwise velocity has been previously shown.
The normalized turbulence kinetic energy distributed in different wake positions for three extreme flap positions and two flaplengths is shown in Figure 17a. If the curves corresponding to

Pressure Coefficient
As shown in Figure 18, the curves of the pressure coefficient corresponding to the flap position β = 0° are coincident for the flaplengths of 8% of c and 14% of c. The curves corresponding to the flaplength of 14% of c are the ones with the greatest area of pressure.
The curves corresponding to the positions of β = 5° and β = −5° are symmetric with respect to the pressure curve β = 0°. Besides, for the same flaplength, the curves of β = 5° and β = −5° are coincident. As mentioned on previous occasions, this is a phenomenon that happens due to the symmetry of the airfoil NACA0012. As shown in Figures A14 and A15 of Appendix F, this does not happen for airfoils NACA 64(3)-618 and S810.
The slight disturbance appeared between x/c = 0.25 and x/c = 0.5 reflects the laminar-turbulent transition zone.

Conclusions
A computational study of three airfoils, NACA 0012, NACA 64(3)-618, and S810 installed with trailing edge flaps has been performed by means of 2D computational fluid dynamic simulations at a Reynolds number of Re = 10⁶. In this work, the open source code OpenFOAM [39] was used for the simulations. The simpleFoam solver, combined with the k-ω SST turbulence model, was applied for As shown in Figure 16a, due to the reduction of the influence of the airfoil, the curves of normalized streamwise velocity furthest from the trailing edge, x/c = 1, are those that present a lesser alteration. This happens indifferently from the flaplength. The measurements of these alterations can be obtained in Figure 16b,c. The same effect for the airfoils NACA 64(3)-618 and S810 is shown in Figures A10 and A12 of Appendix E.
As shown throughout the results section, also for the streamwise velocity, a greater effect is generated with a flap of 14% of c. The effect of a flaplength of 14% of c is also greater for the airfoils NACA 64(3)-618 and S810.
Because the NACA 0012 is a symmetrical airfoil, the effect generated by the two extreme flap positions of β = 5 • and β = −5 • is absolutely symmetric with respect to the horizontal axis. This does not happen for the NACA 64(3)-618 and S810 asymmetrical airfoils.

Turbulence Kinetic Energy
Then, as shown in Figure 17, the results of the normalized wake turbulence kinetic energy are displayed in the same format in which the streamwise velocity has been previously shown.
The normalized turbulence kinetic energy distributed in different wake positions for three extreme flap positions and two flaplengths is shown in Figure 17a. If the curves corresponding to different wake positions are superimposed for each flaplength of 8% of c and 14% of c, enlargements of the disturbed areas are shown in Figure 17b,c, respectively.
As for the streamwise velocity, the curves furthest from the trailing edge are those that present a lesser alteration. Also, for turbulence kinetic energy, a greater effect is generated with a flaplength of 14% of c. The effect generated by the flap positions β= 5 • and β = −5 • is also symmetric with respect to the horizontal axis for airfoil NACA 0012.
The results corresponding to the airfoils NACA 64(3)-618 and S810 are shown in Figures A11 and A13 of Appendix E.

Pressure Coefficient
As shown in Figure 18, the curves of the pressure coefficient corresponding to the flap position β = 0 • are coincident for the flaplengths of 8% of c and 14% of c. The curves corresponding to the flaplength of 14% of c are the ones with the greatest area of pressure.
The curves corresponding to the positions of β = 5 • and β = −5 • are symmetric with respect to the pressure curve β = 0 • . Besides, for the same flaplength, the curves of β = 5 • and β = −5 • are coincident. As mentioned on previous occasions, this is a phenomenon that happens due to the symmetry of the airfoil NACA0012. As shown in Figures A14 and A15 of Appendix F, this does not happen for airfoils NACA 64(3)-618 and S810.
The slight disturbance appeared between x/c = 0.25 and x/c = 0.5 reflects the laminar-turbulent transition zone.

Conclusions
A computational study of three airfoils, NACA 0012, NACA 64(3)-618, and S810 installed with trailing edge flaps has been performed by means of 2D computational fluid dynamic simulations at a Reynolds number of Re = 10 6 . In this work, the open source code OpenFOAM [39] was used for the simulations. The simpleFoam solver, combined with the k-ω SST turbulence model, was applied for steady-state, incompressible, and turbulent flow using RANS equations.
Having previously validated the CFD model, the procedure followed in this work has been schematically shown in The results add new knowledge about the effects of TEFs applied to three airfoil sections specifically intended for wind turbine application. This research can provide useful data for the research community in developing new blade designs and aerodynamic forces control strategies for wind turbine rotors. The obtained results are in accordance with the aerodynamic results expected when studying a TEF inserted on an airfoil.
As the flap is deployed towards the pressure area of the airfoil, reaching higher values of negative β angles, the lift increases. On the contrary, as the flap is deployed towards the suction zone of the airfoil, reaching higher values of positive β angles, the lift decreases. In addition, the results for all the three airfoils have shown that the greatest flaplength studied here maximizes the C L /C D ratio in the positions in which the highest lift force is obtained, and minimizes the C L /C D ratio in the positions in which the lowest lift force is obtained.
To show the results on which this last conclusion is based, athree-dimensional surface representing C L /C D ratio variation as a function of the flap angle β and the different flaplengths has been addressed.
The On the one hand, taking these results into account, aerodynamic forces reduce the blade of a wind turbine by deflection of the flap towards the suction zone, such that when positive β increases, this would be an acceptable option. On the other hand, the deflection of the flap towards the pressure zone, when positive β decreases, would be an acceptable means of increasing the output power of the turbine. The recommended flaplength in both cases would be the largest of them all, in this study 14% of c.
The prediction model built to obtain the C L /C D ratio variation of the airfoils as a function of flap angle, β, and flap length, L F , based on an ANN has been shown in Section 4.4. As the numerical and graphical results show, this approach might represent a good option to facilitate future computational automatic-control tools for TEFs installed on airfoils.
Complementary data about streamwise velocity and turbulence kinetic energy in different points of the wake region as well as pressure coefficient data added to this study can give more information and consolidate the results previously exposed.

Appendix B
(c) (d) Figure A3. Lift-to-drag ratio CL/CD curves as a function of angle of attack, α, for different angles, β.