Transition turbulence model calibration for wind turbine airfoil characterization through the use of a Micro-Genetic Algorithm

The aerodynamic characterization of airfoils is of crucial importance for the design and optimization of wind turbines. The present paper tries to provide an engineering methodology for the improvement of the accuracy and reliability of 2D airfoil computational fluid dynamics models, by coupling the ANSYS Fluent solver and a Micro-Genetic Algorithm. The modeling strategy provided includes meshing optimization, solver settings, comparison between different turbulence models and, mainly, the calibration of the local correlation parameters of the transition turbulence model by Menter, which was found to be the most accurate model for the simulation of transitional flows. Specifically, the Micro-Genetic Algorithm works by generating populations of the missing local correlation parameters. In doing so, it is possible to search for the minimization of the error in lift calculations. For each specific Reynolds number, the calibration was carried out only at the Angle of Attack where the lift drop occurs and the airfoil completely stalls. This new idea allowed for a relatively rapid and good calibration as demonstrated by the experimental–numerical comparisons presented in this paper. Only the experimental stall angle and the relative lift coefficient were, therefore, necessary for obtaining a good calibration. The calibration was made using the widely known S809 profile data. The correlation parameters, obtained as so, were subsequently used for testing on the NACA 0018 airfoil with satisfactory results. Therefore, the calibration obtained using the S809 airfoil data appeared to be reliable and may be used for the simulation of other airfoils. This can be done without the need for further wind tunnel experimental data or recalibrations. The proposed methodology will, therefore, be of essential help in obtaining accurate aerodynamic coefficients data. This will drastically improve the capabilities of the 1D design codes at low Reynolds numbers thanks to the possibility of generating accurate databases of 2D airfoil aerodynamic coefficients. The advantages of the proposed calibration will be helpful in the generation of more accurate 3D wind turbine models as well. The final objective of the paper was thus to obtain a fine and reliable calibration of the transition turbulence model by Menter. This was specifically made for an accurate prediction of the aerodynamic coefficients of any airfoil at low Reynolds numbers and for the improvements of 3D rotor models.


Introduction and background
The design and optimization of wind turbine rotors are usually made using mono-dimensional codes, which allow for very fast responses. However, the widely used 1D blade element momentum models (BEM), as well as the other simplified codes, need lift and drag coefficients for a wide range of Reynolds numbers to provide accurate and reliable results. The experimental airfoil data can be obtained only in wind tunnels and they are thus highly expensive to obtain, both in terms of time and money. For this reason, many CFD models have been developed by the researchers to allow for a quick prediction of the aerodynamic coefficients [1][2][3][4]. Furthermore, the increasing interest towards the low Reynolds numbers, in ''mini'' and ''micro'' wind turbine applications, has led to the modeling of an aerodynamic behavior. This behavior is characterized by the further complication of the laminar to turbulent boundary layer transition. Indeed, the correct modeling of the transitional flows and, specifically, the effect of the laminar instability on the flow separation, is a key step for obtaining an accurate prediction of the airfoil aerodynamic coefficients at low Reynolds numbers. It is widely known, in fact, that the laminar boundary layer is more sensitive to separation than the turbulent boundary layer. Moreover, the onset of laminar bubbles is a phenomenon that must be taken into account in airfoil modeling [5][6][7][8][9]. Therefore, the use of a fully turbulent model as well as a non-optimized transition turbulence model lead, in general, to a delayed separation, which occurs at AoAs higher than the real one. As a consequence, the lift coefficient is usually overestimated while the drag coefficient is underestimated. This may not be acceptable for reliable 1D codes, in which further simplifications are carried out to take into account the 3D effects. The 2D aerodynamic coefficients must, therefore, be as accurate as possible [10].
The ideas presented in this paper come out of the necessity to generate databases of 2D airfoil aerodynamic coefficients for the use in the in-house 1D BEM model and to obtain a calibrated transition turbulence model, suitable for full CFD wind turbine rotor simulations. Previous works by the authors [11][12][13] have demonstrated an excellent accuracy of the four equations transition SST model, developed by Menter et al. [14][15][16]. Furthermore, scientific literature underlines and demonstrates that, the use of transition turbulence models, at low Reynolds numbers [10,[17][18][19][20][21][22][23], is essential in obtaining an adequate CFD modeling of airfoil and rotor aerodynamic behaviors.
The original model by Menter [14][15][16] was based on three local correlation parameters, which control the transition onset inside the boundary layer. These parameters were empirically calibrated using flat plate and gas turbine blade experimental data. The functions related to this optimization are still proprietary and thus not given in references [14][15][16]. Many authors have proposed different ways to find the missing functions [24][25][26][27][28]. Other authors successfully used the transition model by Menter for the solution of various low Reynolds aerodynamic problems [18,21,[29][30][31][32][33][34][35]. During the last few years, the transition model by Menter has been extensively studied and has demonstrated good accuracy in the cases in which the laminar to turbulent boundary layer transition played an important role. Most of the research has focused on the deduction of the missing functions, basing their calibration using flat plate test cases. The promising results have encouraged the investigation of the performance of the transitional models for 2D airfoil and 3D wind turbine fluid dynamic modeling. The c-Re h model appears to be suitable for this research field as well. However, the model has proven to be highly sensitive to Reynolds number, geometrical shapes and turbulence parameters, therefore, necessitating for a fine calibration case by case. The developers themselves recognized that, for airfoil and wind turbines, a specific calibration would be necessary [36].
The objective of this work was, therefore, to provide a rapid and reliable strategy for the development of CFD models, which were able to capture the average steady behavior of the flow-field around airfoils at low Reynolds numbers. This was the aim in obtaining accurate lift and drag coefficients for a wide range of low Reynolds numbers. The results might subsequently be extended to 3D rotors as well, with good accuracy and reliability [11][12][13]. The calibration of the missing parameters was made by coupling the ANSYS Fluent solver with a Micro-Genetic Algorithm (lGA), developed by the authors in MATLAB. The idea, leading to this study, was to investigate the possibility of fine tuning the empirical correlations, by directly specifying the values inside the solver for only a specific AoA and to verify the effects of this tuning over the entire range of AoAs. The calibration process was, therefore, carried out for the typical wind turbine airfoils S809 for which a large set of experimental data was available. The optimal parameters, which were found for the S809, were then tested by simulating the NACA 0018 airfoil, finding very satisfactory agreement with experimental data. Moreover, the empirical correlations were optimized for Reynolds number between 300,000 and 1,000,000. Polynomial functions were, therefore, found to interpolate the parameter trends as a function of the Reynolds number. This allows for the calculation of aerodynamic coefficients at any desired Reynolds number and for any airfoil within the Reynolds number range of optimization. No parameter recalibrations or further experimental data are necessary.

Transition and turbulence modeling
The laminar to turbulent boundary layer transition is a highly complex phenomenon; particularly important at low Reynolds numbers (Re \ 2 million). The causes leading to the onset of the transition have not yet been fully clarified. However, four types of transition mechanisms have been recognized as the main factors responsible for the first order effects [37]. From an engineering point of view these explanations are more than sufficient and allow for a good comprehension and modeling of the phenomenon. The types of mechanisms are the natural transition, the bypass transition, the separation-induced transition and the wake induced transition [37,38]. It is obvious that an accurate transition model must take into account all the aforementioned mechanisms.
The formulation proposed by Menter et al. [14][15][16]36] is based on the widely known SST k-x turbulence model with two additional transport equations. The first is an intermittency equation used to trigger the transition process. The second equation is formulated in terms of the transition onset Reynolds number Re ht . The transport variables are forced to follow the values provided by the experiments. The proposed transport equations do not directly model the physics of the transition as it is contained in the experimental correlations. Therefore, the fact that the model is based on empirical correlations makes it suitable for any kind of transition mechanism but only if appropriate values can be provided [14]. This stresses the need for an accurate calibration of the local parameters for airfoil and wind turbine simulations.
The c-Re h model formulation is widely presented in the aforesaid literature, therefore, only a brief review, on what is of particular interest to this paper, is reported.
The transport equation for the intermittency c is: The production and destruction transition sources are provided by the terms P and E. Specifically: In (2) and (3) the term F length controls the transition length while Re hc is the critical Reynolds number where the intermittency first starts to increase in the boundary layer. Therefore, Re hc is the point where the model is activated to match both Re ht and F length . Clearly this must occur upstream of the transition Reynolds number, and therefore, Re hc \ Re ht .
The transport equation for the transition momentum thickness Reynolds number is: f Re ht is the transported value of the boundary layer momentum thickness Reynolds number. Therefore, the model contains three empirical correlations (6), which must be accurately provided to close the model: These local correlation variables were calibrated in the original model by implementing functions obtained through the use of experimental data, carried out for flat plates and gas turbine blades [14-16, 30, 36]. As previously seen, many authors proposed different ways to obtain the functions with satisfactory results. However, the Fluent solver allows for a direct specification of the three missing correlation using a user-defined function (UDF) written and compiled in C language. The advantage of a direct specification of the values is related to the aforementioned sensitivity of the model to different Reynolds numbers. In this way a precise calibration, made case by case, may be obtained and will be useful for thorough research on wind turbine airfoil transition effects. The value of Re ht depends on the turbulent intensity TI and on the Thwaites' pressure gradient k. The correlating function for Re ht is given in the original paper by Menter [15,16] and reads: The momentum thickness Reynolds number was, therefore, estimated in this paper, using the inlet boundary conditions which were TI = 0.1% and k = 0 [30].
Once the value (8) is known, only F length and Re hc must be calibrated through the proposed lGA.

CFD airfoil models
The first step for an accurate and reliable CFD model is the definition of the computational domain. In this work, a C-type domain was developed. The C-type domain allows for a rapid change of the inlet velocity components without changing the mesh. In this way, both different Reynolds numbers and AoAs were easily set.
The computational domain was the same for both S809 and NACA 0018 airfoils and is shown in Fig. 1 along with the dimensions in terms of chord length (c).
To obtain the best balance between computation time and accuracy, a grid independence study was carried out for both S809 and NACA 0018 profiles. The mesh was developed using ANSYS Icem. Full structured grids were implemented, gradually increasing the number of nodes on the profile as well as on the direction orthogonal to the airfoil and on the wake edge. Specifically, three grid refinements were tested so as to have a y ? \ 1 for all the grids. In Table 1 the specific grid features are presented. In Fig. 1 the boundaries are set as velocity inlet and pressure outlet. The velocity inlet boundary type allows for a fast implementation of the velocity magnitude and components so that the AoA (from 0°to 89°) and the Reynolds number of interest were easily implemented.
As aforesaid, the calibration process was carried out for the AoA in which the lift coefficient drops in stall. For this reason the grid independence study was done for both the airfoils and for all the Reynolds numbers of interest at this AoA and using the default transition turbulence model. For each case, several other AoAs were simulated to check the validity of the grid study for a wide range of conditions. The best compromise was always found with Grid 2 which ensured a good accuracy and a noticeable restraint of the calculation time. In Table 2 the lift and drag coefficient comparison, related to the grid independence study, are presented. In Fig. 2 details of the structured C-type mesh for the S809 airfoil are shown (Grid 2). The accurate resolution of the boundary layer and the uniformity of the grid are clearly evident.
For all the simulations, the Fluent solver was set as a steady state, pressure-based coupled solver with an absolute velocity formulation. The fluid was air with a density equal to 1.225 [kg/m 3 ] and a dynamic viscosity equal to 1.7894 10 -5 [Pa s]. Through the velocity inlet boundary condition, the velocity components were imposed along with turbulent intensity (TI = 0.1%) and turbulent viscosity ratio (TVR = 10). This was done by following the literature suggestions given [30,36]. This lower TI corresponds to that of the typical low turbulence level wind tunnels and it implies a natural or a separation-induced transition mechanism.
A least squares cell-based scheme was used for the gradient spatial discretization while second-order upwind schemes were used for pressure, momentum and turbulence equations. It was very important to set a correct value for the Courant number, specifically in full separated conditions at AoAs [ 28°. In the Fluent steady coupled solver, the Courant number allows for the control of the sub-iteration time step. In this way, even if the steady solver is used, a transient sub-iteration is performed and the time scale of the simulation can be adapted to the flow time scale. In attached conditions and in incipient and moderate stall conditions, the Courant number was set equal to 25 thus allowing for a rapid convergence. At AoAs [ 28°the unsteadiness, generated by the stall, become noticeable. In doing so, the Courant number must be reduced to adapt the temporal scale in such a way that an accurate average of the lift and drag coefficient can be obtained. Three iteration monitors were activated. The criterion to check the convergence was to achieve constant residuals below 10 -3 for all the equations. Constant trends for the lift and drag coefficient monitors were checked as well. A full multi grid (FMG) initialization was implemented to provide an optimal initialization, which resulted in a faster convergence. The RANS equations were solved using the aforesaid four equations transition turbulence To provide the local correlation parameters F length and Re hc a User Defined Function (UDF) was written and compiled in C language. The UDF was used by the lGA for the calibration process as well. Substantially, for each lGA simulation, MATLAB re-wrote the UDF text file with the new generated parameters and the new UDF file was subsequently read by Fluent which performed the new simulation, launched directly inside MATLAB. In this way, the calibration process was fully automated for each Reynolds number.

Micro-Genetic Algorithm details and optimization procedure
The Simple Genetic Algorithms (SGAs) have been widely used in the last decades in engineering and industrial optimization problems. They have been applied in many different fields such as automatic programming, machine learning, economics and ecology [39,40].
Unfortunately, when the GAs are coupled with any CAE tool (e.g. Fluid Dynamics codes), the computational time, needed for the evaluation of the objective function, greatly increases. Subsequently, the ''good'' solution may not be found in a reasonable time.
The Micro-Genetic Algorithms (lGAs) have been implemented to deal with this issue. It has been shown that they have been able to find the near-optimal solution with much less function evaluations, compared to a SGA [41][42][43].  The step-by-step procedure for the lGA implementation proposed by Krishnakumar [41] is presented below: 1b.
A population of size 5 is randomly generated; 2b.
The fitness is evaluated and the best string is determined. It has been labeled as string 5 and carried to the next generation (elitist strategy); 3b.
The remaining four strings for reproduction, based on a deterministic tournament selection strategy, are chosen; 4b.
Crossover with P c = 1 is applied; 5b.
The mutation rate equal to 0 is set; 6b.
Nominal convergence is checked. If converged, proceed to step 1b; 7b.
Go to step 2b.
Krishnakumar suggested the use of a deterministic tournament selection strategy due to the low number of individuals. In fact, it would be foolish to consider the law of averages and, consequently, to use the usual stochastic selection procedures [44].
Concerning step 6b, after the population converges to a prescribed measure, it is important to introduce new strings based on either genotype convergence of phenotype convergence, or based on the occurrence of a fixed number of generations. Thanks to this re-initialization of the population, new genetic material is introduced after every convergence. Obviously, the re-initialization phase does not deal with the best string, which must always be carried to the next generation.
Following the procedure proposed by Krishnakumar a lGA was developed basing on a SGA which was implemented in MATLAB [45]. A deterministic tournament selection routine was written to choose the four individuals to fill the mating pool for each generation. The individuals were mated with P c = 1 and their offspring were not subjected to any mutation. A control loop was used to reinitialize the population every five generations. A subroutine was implemented to further enhance the quality of the genetic available material. The subroutine task was to store all the individuals created thus for and to generate new ones by randomly choosing from within the parameter space, thus avoiding useless copies.
Since the optimization problem of the present paper deals with a uni-modal function, the proposed lGA was tested to solve the following function [40]: with À5:12 x; y; z 5:12 and dx ¼ dy The result of the best-so-far fitness for the in-house GAs is reported in Fig. 3 (left) while Fig. 3 (right) shows the result obtained by Krishnakumar [41] for the same test function. Specifically, in Fig. 3, the test function (9) was evaluated using the in-house GAs (lGA and SGAs) and those developed by Krishnakumar. The comparisons show a very similar trend and a faster convergence for both the lGAs. Moreover, Fig. 3 demonstrates the validity of the in-house lGA compared to that of Krishnakumar. The main difference between the in-house GAs and those developed by Krishnakumar is the type of encoding. Indeed, while Krishnakumar used a binary string of length 30 to map the variables, in the in-house algorithm a realcoded GA was used. The advantages of this type of representation are the absence of difference between the genotype and phenotype, a faster convergence and a reduction of the deception probability. The lines in Fig. 3 appear to be jagged due to the fact that the GA optimization was conducted for 25 different random starts and an ensemble average of the best string so far was calculated.
Once tested, the lGA was coupled to ANSYS Fluent to perform the optimization of the missing parameters. Specifically, a routine, which read the parameters Re hc and F length , randomly chosen by the lGA, was implemented for writing the new UDF that must be read by the CFD code. The same routine allowed for starting the 2D simulation in batch mode and for extracting and storing the results in terms of C l and C d . Furthermore, the lift and drag coefficients, calculated by each Fluent simulation, were compared by the routine with the specific experimental values. In this way, the genetic material, provided to the lGA was improved and the error was minimized. The algorithm was stopped after 80 generations, which corresponded to 400 function evaluations. This represented a very good tradeoff between accuracy and low computation time compared to the SGAs.
However, making 400 evaluations means that Fluent must be run 400 times. Therefore, having previously established the number of iteration to reach the convergence, each lGA run lasted at least 7-8 days.
The calibration procedure may be detailed as follows. First a Fluent case is implemented with the specific velocity and AoA settings, following what is reported in the previous section. A test run is carried out so that the journal file can be written and the number of iterations to reach the convergence can be fixed. At this point, the case file is associated to the lGA in MATLAB and the calibration process is launched. The lGA generates a population of five individuals for both Re hc and F length , which are automatically written in the UDF. The Fluent solver is thus started in batch by reading the journal file and the simulation proceeds until the prescribed number of iterations is reached. Now, the converged lift and drag coefficients are read by the lGA from the Fluent text files and are compared to the experimental values. Once the five simulations, related to the five individuals, are completed, the lGA generates a new population taking into account the error between the numerical and experimental data. The best couple of individuals that is the one that minimize the errors on lift and drag coefficients, is carried into the new population. The lGA continues with this logic until the 80 populations are generated.
The innovative idea of this work is to choose the AoA where the lift curve completely falls in stall and to calibrate the missing local correlation parameters only at this AoA. The rationale behind this idea is that the stall point can be accurately captured only if the boundary layer is correctly modeled. In transitional flow conditions, this occurs only if the laminar to turbulent transition is accurately predicted. This means that, for airfoils, the natural or the separationinduced transition must be reproduced in a precise way. Therefore, the calibration of the missing parameters at the stall AoA would lead to a good optimization of the transitional model without the need to use the expensive and difficult experimental analysis of the boundary layer. Clearly, a small amount of experimental data, about the stall angles and the relative lift coefficients, were needed to appropriately calculate the error in the lGA. However, once having calibrated the model for some Reynolds numbers, the optimal parameters may be extrapolated for any other Reynolds number within the range of optimization, by interpolating the correlation variables with polynomial functions.

Results and discussions
To check the reliability of the proposed optimization, the calibration process was first carried out for the S809 airfoil, for which a large set of widely validated experimental data was available. The correlation parameters, obtained as so, were subsequently interpolated depending on the Reynolds number, which was between 300,000 and 1 million. The optimal parameters were then tested by simulating the NACA 0018 airfoil. A good agreement between numerical and experimental data was found. The comparison with the default SST transition and the fully turbulent SST k-x model results, further demonstrate the good reliability of the proposed methodology.
In Fig. 4 the calibrated missing parameters for the S809 airfoil are shown as a function of the Reynolds number. Furthermore, the interpolating functions, which are able to fit the optimal values of F length and Re hc , at the desired Reynolds number, are highlighted. These are third-order polynomial functions and are defined as follows: Re hc ¼ 3:9592e À16 Re 3 À 9:598e À9 Re 2 þ 6:884e À4 Re þ 984:0408 The experimental data from Colorado State University (CSU), Ohio State University (OSU) and Delft University of Technology (DTU) [46] were available for Reynolds numbers 300,000, 500,000, 650,000 and 1 million. The calibration was, therefore, carried out at AoA = 17°for Re equal to 300,000, 500,000 and 650,000, and at AoA = 20°f or Re equal to 1 million. These were the stall angles as can be seen in Figs. 5 and 6.
The comparisons between experimental and numerical results are shown in Figs. 5 and 6 for the S809 airfoil at the specific Reynolds numbers.
Analyzing the comparisons in Figs. 5 and 6, one can clearly see an excellent improvement of the predictive capability of the optimized transition turbulence model. In the attached flow region, where the lift coefficient trend is linear, all the turbulence models demonstrate high accuracy, as expected. In the incipient stall region, between 7°-8°and 16°-18°of AoA, the calibrated transition turbulence model proves to be very accurate, with a maximum relative error at.
Re 500,000 and Re 650,000 less than 8% at about 8°of AoA. The stall angle and the relative C l are well predicted for all the Reynolds numbers. The default transition model, instead, is not as accurate, particularly in the flattened lift region before the deep stall, where a noticeable overestimation is detectable. The stall angle is delayed and the relative C l is nearly twice the experimental value. The fully turbulent SST k-x model is even less accurate with an overestimated prediction, higher than the default transition model. The stall angle is 4°/6°higher than the experimental value except for Reynolds number equal to one million. In this case, as the flow becomes more turbulent, the discrepancies between the models are reduced. However, the experimental data for Re 1 million are affected by a slight  uncertainty in the stall angle as evidenced in Fig. 6 (right). The drag data are not presented here due to the fact that only the experimental pressure drag was available; therefore there were no comparable data at low AoAs. However, a good agreement was also found at medium and high AoAs. The reason for the improvements on the prediction capabilities of the calibrated transition model lies in the improved modeling of the laminar behavior of the boundary layer (laminar bubbles) and, specifically, on a much accurate evaluation of the laminar to turbulent boundary layer transition due to laminar bubbles. At lower Re, when the boundary layer should be entirely laminar, the effects of the thick shape of the profile lead to an earlier separation-induced transition which is typical of the S809 airfoil. The proposed calibration of the local variables allows for an adequate prediction of this phenomenon thus resulting in a good agreement between numerical and experimental data. This is evidenced in Fig. 7 where the comparison between the calculated skin friction coefficients for the suction side of the S809 airfoil at AoA 0°and Re 500,00 is shown. Up to x/c = 0.4 the three models predict an identical trend of the C f . This is the laminar region and here within, the models have the same SST k-x low-Re formulation. After this point, the particular shape of the airfoil leads to the generation of a laminar separation bubble which forces the laminar to turbulent boundary layer transition. The discrepancies in this area, between x/ c = 0.4 and x/c 0.6, are evident. The SST k-x model overestimates the skin friction due to excessive turbulent kinetic energy production, thus resulting in a sudden transition. Both the transition models, instead, correctly predict the drop of the skin friction due to the laminar bubble. However, the optimized model shows a slightly larger extension of the bubble which is crucial for the correct prediction of stall angles and lifts and drag coefficients. Indeed, as the AoA increases, the laminar bubble moves upstream, thus a larger extension will correctly result in an earlier separation. A much larger extension is avoided thanks to the minimization of the error in the lGA. This is due to the fact that a larger extension should result in lower lift coefficients and thus in higher errors.
The above is supported by the post-processing images presented hereinafter. In Fig. 8 the contours of intermittency are shown. According to the skin friction distribution in Fig. 7, the intermittency inside the boundary layer is similar for both the transition models. However, the calibrated model (left) presents a slightly delayed and much gradual intermittency transition from 0 to 1, compared to the default model (right). This denotes that the laminar bubble is a slightly extended and the turbulent boundary layer develops downstream.
The images in Fig. 9 present the contours of turbulent kinetic energy for the calibrated (a), the default (b) transition model and for the fully turbulent SST k-x model (c). Only slight differences can be detected between default and calibrated model results. Mostly, only a small increase of the turbulent kinetic energy production in the calibrated model (a) is present, which is related to the slightly larger extension of the laminar bubble. According to the skin friction chart in Fig. 7, the fully turbulent SST k-x model (c) shows a much higher production of turbulent kinetic energy and this happens upstream compared to (a) and (b) cases. The magnitude of turbulent kinetic energy in (c) is lower than in (a) and (b) due to the lack of the bubble prediction, but the extension and the location are noticeably different.
The improvements, thanks to the implemented calibration, are evident in Fig. 10 where the velocity magnitude contours for the three models are plotted at the experimental stall angle (AoA = 17°) for the S809 airfoil at Re 500,000. Although the differences in the skin friction, intermittency and turbulent kinetic energy prediction appear to be very limited between (a) and (b), these differences are more than sufficient in obtaining a correct separation prediction. Indeed, in Fig. 10a a full separation When simulating a full 3D wind turbine or when using CFD models for the calculation of the aerodynamic coefficients for use in 1D BEM models, such an error is not acceptable. Therefore, the proposed calibration of the transition turbulence model will drastically improve the reliability of both BEM and CFD 3D models at low Reynolds numbers.
In Fig. 11 details of the turbulent kinetic energy contours at the leading edge (AoA 17°and Re 500,000) of the S809 airfoil are shown. In (a) the production of turbulence is much larger due to the full flow separation but the magnitude is lower. In (b) the higher magnitude of turbulence production causes the boundary layer to remain attached for a long extension. In (c) the turbulence magnitude is even higher and the flow results further attached. This is due to the widely known stability of the turbulent boundary layer against the tendency to separation, compared to the laminar boundary layer. A correct laminar bubble modeling is, therefore, a key step for an adequate prediction of the flow-field around airfoils and wind turbines.
The reliability of the proposed calibration was finally tested on the NACA 0018 airfoil. The local correlation parameters were calculated using Eqs. (11) and (12), at the Reynolds numbers for which the experimental data were available. The CFD 2D model was developed through the set-up presented in ''CFD airfoil models''. The results are shown in Fig. 12, where the experimental lift coefficient trends for the NACA 0018 are compared to the CFD simulation results. Also in this case, the calibrated, the default and the fully turbulent model results are presented. The improvements for the calibrated model are quite remarkable.
Specifically, at Re 300,000 (Fig. 12a) the behavior of the turbulence models appears to be similar to that of the S809 (Fig. 6). The calibrated transition model slightly underestimates the C l values in the incipient stall region, from 8°to 16°of AoA, but predicts in a satisfactory manner the stall angle and the relative lift coefficient. The default transition model instead slightly overestimates the lift values in the same region and shows a delayed stall angle by 2°of AoA. The SST k-x model in this case is not as accurate, even in the attached flow condition where the transition models are instead very accurate. Furthermore, the fully turbulent model presents the same stall angle delay as the default transition model. However, the advantages of the calibration are strikingly evident. The experimental data used for the calibration and the comparison were those of Timmer [47].
At higher Reynolds numbers, Re 1 million (Fig. 12b), the turbulence model behavior is different from that of the S809. The good agreement with the experimental data of Timmer [47] for the calibrated transition model further confirms the reliability of the proposed optimization. Only a slight overestimation is detectable in the incipient stall region. The default model and the SST k-x model in this case present a premature drop of the lift coefficient and the usual overestimation in the incipient stall region. The reason for this is not clear at the moment and will be further investigated. A similar trend is evident in Fig. 12c at Re 700,000. Here the numerical data are compared to three different sets of experimental data by Timmer [46], by Sheldal [48] and by Jacobs [49]. The striking differences between the experimental data in incipient and deep stall conditions make the level of uncertainty of this kind of experimental measures clearly understandable. However, the calibration was made using the data of Timmer and again shows a good improvement of the transition model predictive capabilities.

Conclusions
The present paper describes the methodology, put forth by the authors, for the optimization of the local correlation parameters which control the boundary layer transition inside the transition turbulence model developed by Menter. Unlike the numerous calibrations, proposed in the scientific literature, which have tried to reproduce the default correlating functions which had been missing, the present methodology was created specifically for the development of 2D CFD airfoil models. The idea was to directly optimize the two missing parameters, F length and Re hc , by coupling the ANSYS Fluent solver and a micro-genetic algorithm, specifically developed by the authors for this application. The use of the experimental lift coefficient at the stall AoA for the S809 airfoil in the lGA is already sufficient in obtaining a good calibration for the whole range of AoAs at the specific Reynolds number. The rationale behind this idea is thoroughly presented in the paper and demonstrates highly accurate results. In this way, the calibrated local correlation parameters, presented in this paper, may be used for the simulation of other airfoils within the range of Reynolds numbers, which were between 300,000 and 1 million. This will, therefore, make the wind tunnel experiments superfluous.
Several tests were carried out on other airfoils. For the NACA 0018, the results appeared to be accurate enough; therefore the proposed calibration may be considered a powerful tool for reliable airfoil simulations in transitional flow conditions. Furthermore, this will allow for the calculation of 2D aerodynamic coefficients for any airfoil and will help in the development of improved 1D BEM wind turbine models. The research on the use of different airfoils was in fact often limited due to the lack of experimental data at the desired Reynolds numbers. The use of the polynomial interpolating functions, obtained in this paper, allows for the generation of databases of accurate aerodynamic coefficients at the desired transitional Reynolds number. Certainly, this will be of great importance for future research on wind turbines. Moreover, the calibrated correlation parameters may be used to improve the predictive capabilities of CFD 3D wind turbine models. Specifically in the cases in which the boundary layer transition plays an important role and the rotor performance are mostly influenced by transitional effects.

Compliance with ethical standards
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.