Design and Optimization of a Multi-Element Hydrofoil for a Horizontal-Axis Hydrokinetic Turbine

Hydrokinetic turbines are devices that harness the power from moving water of rivers, canals, and artificial currents without the construction of a dam. The design optimization of the rotor is the most important stage to maximize the power production. The rotor is designed to convert the kinetic energy of the water current into mechanical rotation energy, which is subsequently converted into electrical energy by an electric generator. The rotor blades are critical components that have a large impact on the performance of the turbine. These elements are designed from traditional hydrodynamic profiles (hydrofoils), to directly interact with the water current. Operational effectiveness of the hydrokinetic turbines depends on their performance, which is measured by using the ratio between the lift coefficient (CL) and the drag coefficient (CD) of the selected hydrofoil. High lift forces at low flow rates are required in the design of the blades; therefore, the use of multi-element hydrofoils is commonly regarded as an adequate solution to achieve this goal. In this study, 2D CFD simulations and multi-objective optimization methodology based on surrogate modelling were conducted to design an appropriate multi-element hydrofoil to be used in a horizontal-axis hydrokinetic turbine. The Eppler 420 hydrofoil was utilized for the design of the multi-element hydrofoil composed of a main element and a flap. The multi-element design selected as the optimal one had a gap of 2.825% of the chord length (C1), an overlap of 8.52 %C1, a flap deflection angle (δ) of 19.765◦, a flap chord length (C2) of 42.471 %C1, and an angle of attack (α) of –4◦.


Introduction
For applications requiring high lift forces at large angles of attack and low flow velocities, the traditional hydrofoils tend to produce separation of the flow near the trailing edge, causing a decrease in the hydrofoil performance. Multi-element hydrofoils constitute a proper alternative since the lift force is increased, leading to a rise in the camber when operating at a high angle of attack and a delay of the flow separation near the trailing edge. This delay in the flow separation on the deflected flap element is achieved by introducing a slot ahead of the flap for boundary layer control [1].
Multi-element hydrofoils are applied across a wide range of engineering disciplines. Most of the modern aircrafts use multiple flaps and slots to increase both the surface area and the camber during the critical takeoff and landing stages of the flight [2]. In turn, in the automotive industry, cars competing in Formula One races use multi-element foils to increase the down-force produced by the rear wing [3].
Horizontal-axis turbines are generally utilized for wind turbines, focusing on the use of multi-element airfoils to eliminate the losses of performance due to the requirement of thick airfoils in the inboard section of a wind turbine blade to resist structural, aerodynamic, and gravity loads resulting from the operation of the turbine [4]. Indeed, results reported by Ishaan Sood [5] demonstrated the capability of a turbine with multi-element blades to produce higher maximum power coefficient (C Pmax ) than a wind turbine with a traditional airfoil under equal operational conditions; for the same nominal power, the multi-element blade was proved to have higher efficiency and, thus, higher performance than that achieved with a traditional hydrofoil. On the other hand, the designed multi-element blade achieved the nominal power with a smaller blade radius, leading to material and manufacturing cost savings. The multi-element airfoil designed by Ragheb and Selig [4] offers an increase in the relationship between the lift and drag coefficients (C L /C D ) when compared to those included in the Delft University family of wind turbine airfoils [6,7] for C L values ≤ 1.7.
Hydrokinetic systems (HKs) are a class of zero-head hydropower systems by which the energy is extracted from the kinetic energy of river, tidal, and ocean currents in contrast to the potential energy associated with falling water [8][9][10][11][12]. HKs produce electricity for 24 h a day, as long as the running water is available, with minimal costs associated with the produced energy and minimal or mitigatable environmental impact in comparison with the large hydroelectric power plants [13,14]. The methodologies used for the design and analysis of HKs are similar to those ones utilized for wind turbines, which take advantage of the wind kinetic energy [8][9][10][11][12]. The turbine employed can be characterized by its rotational-axis orientation (i.e., horizontal-and vertical-axis turbines) with respect to the water flow direction. It is important to note that the axial turbine usually has high efficiency, self-starting capability, lack of torque fluctuation, and high-speed operation in comparison with a vertical-axis one.
Currently, there are several hydrokinetic technologies under development [8][9][10][11][12]. A number of researchers have carried out several studies focused on various aspects of the technoeconomic feasibility, design optimization, reliability, and location of the turbine [13][14][15][16][17][18][19][20]. Additionally, the augmentation application, anchoring, and environmental monitoring have been studied [10][11][12]. With respect to the design optimization, efforts have generally focused on the maximization of the power coefficient; i.e., on the maximization of the fraction of power in water that can be extracted by the hydrokinetic turbine. The optimization strategy has a direct impact on the blade shape due to the hydrofoil shape contributes to the generation of lift by creating suction on the hydrofoil upper surface [14][15][16][17][18][19]. Numerical analysis has exploited different aspects related to the blade geometric characteristics in order to increase the turbine efficiency. Furthermore, the use of the diffuser-enhanced concept has been explored in many experimental and commercial turbine models. The studies reported a power increase with the use of a diffuser. From the authors' knowledge, only one work has been reported to consider the use of a double-blade hydrofoil for generating the maximum lift [21]. It is important to note that there are not applications of multi-element hydrofoils for hydrokinetic appliances. However, several numerical studies on multi-element airfoil configurations for wind turbines have been reported [22][23][24].
The design of multi-element hydrofoils and airfoils is a complex optimization process due to the geometry of several elements that provides a significant number of independent parameters to be optimized. In fact, during the optimization process, both the position and the shape of the different hydrofoil parameters must be considered. Physical parameters, such as the angle of attack and the fluid velocity, are also crucial factors to be taken into account. Therefore, the number of multi-element airfoil or hydrofoil design variables could rapidly increase by a considerable number of factors to be optimized [25].
For the design of multi-element airfoils and hydrofoils, surrogate models are the most used tools [26]. Several researchers have successfully utilized these models to optimize the performance of wind turbines [27]. Surrogate models are models that can assess the new design objective(s) based on existing training samples. Additionally, there are statistics or approximation mathematical methods [28]. Surrogate models can be used to predict the design objective(s) and optimize the operating conditions throughout the use of an optimizer to find the optimum design and, finally, to utilize the high fidelity approach to verify the proposed objective(s) [29][30][31][32].
Different mathematical models have been usually used to construct surrogate models [26]. Kriging or Gaussian process regressions are usually applied for the multi-objective design of multi-element airfoils [33][34][35]. The Kriging model enables the search process to be efficient, resulting in a drastic reduction of the computational time required [34]. Moreover, it is a non-bias method to construct the surrogate model [35]. Therefore, for creating the surrogate model, the Kriging mathematical model was decided to be used in the current study, which is based on the work developed by Forrester and coworkers [28].
Concerning the optimization schemes, they have been highlighted to range from Tabu search [29,31,32] to genetic algorithm (GA) [34,36], which are the most often used. In addition, Benini and coworkers [36] obtained satisfactory results when using Matlab ® R2019a software (R2019a, MathWorks Inc, Natick, MA, USA) [37], which offers a general vision of a multi-objective genetic algorithm (MOGA) that is available in the global optimization toolbox in Matlab (gamultiobj function of Genetic Algorithm and Direct Search ToolboxTM), which in turn uses a controlled, elitist genetic algorithm that is a variant of NSGA-II [38,39] to create a set of points on the Pareto front.
Under this scenario, the objective of this study is the design and optimization of a multi-element profile for a horizontal-axis hydrokinetic turbine based on a traditional hydrofoil. For this purpose, the optimization problem is described in detail, as well as the surrogate model proposed to solve the optimization problem. Additionally, the results obtained from the proposed surrogate model and their convergences are presented and discussed. Finally, the multi-element hydrofoil selected as optimum is compared to the traditional one that served as the basis.

Description of the Optimization Problem
This study is focused on optimizing a multi-element hydrofoil based on a blade geometry of a hydrokinetic turbine of 1 kW designed by using the Blade Element Momentum Theory (BEM), with a water velocity (V 1 ) of 1.5 m/s, a tip speed ratio (TSR) of 6.325, a power coefficient (C P ) of 0.4382, a transmission efficiency (η) of 70%, a turbine rotor diameter of 1.58 m, and three blades of 0.79 m radius (R) [8,39,40]. For the turbine design, several airfoils analyzed like hydrofoil were used, including S805, S822, Eppler 420, Eppler 421, Eppler 422, Eppler 423, Eppler 857, Wortmann FX 74-CL5-140, Wortmann FX 74-CL5-140 MOD, Douglas/Liebeck LA203A, Selig S1210, Selig S1223, and UI-1720 airfoils. In Table 1, the main characteristics of each of these aerodynamic airfoils to be evaluated are listed.
For the hydrodynamic analysis of the airfoil, considered as hydrofoil, JavaFoil software (2.20, Developer: Martin Hepperle, DLR, German Aerospace Center, 38108 Braunschweig, Germany) was used. For the hydrokinetic turbine blade profiles, small angles of attack (α) are usually selected where C L and C D are high and low, respectively. These coefficients depend on V 1 and, hence, on the Reynolds number, since when the viscosity forces are greater than the inertial ones, the friction effects are increased, affecting the velocities, the pressure gradient, and the lift generated by the hydrodynamic profile [40,41]. Therefore, for the design of the hydrokinetic turbine blade, profiles with a high C L /C D ratio were chosen to be studied. In addition, the selection of a profile with a large section and a considerable thickness to withstand the mechanical forces induced during the operation of the blade is required, avoiding incurring in a drag increase [39][40][41].
The profiles were analyzed for a Reynolds number characteristic of hydrokinetic turbines, which was equal to 750,000. The chord length (C) remained unitary in traditional and multi-element profiles. C L and C D were analyzed every 1 • , within a wide α range. This procedure was carried out for the whole set of the analyzed profiles. The comparison between the diagram of C L and C D versus α for the airfoils used as hydrofoils obtained by numerical methods is illustrated in Figure 1. For the hydrodynamic analysis of the airfoil, considered as hydrofoil, JavaFoil software (2.20, Developer: Martin Hepperle, DLR, German Aerospace Center, 38108 Braunschweig, Germany) was used. For the hydrokinetic turbine blade profiles, small angles of attack (α) are usually selected where CL and CD are high and low, respectively. These coefficients depend on V1 and, hence, on the Reynolds number, since when the viscosity forces are greater than the inertial ones, the friction effects are increased, affecting the velocities, the pressure gradient, and the lift generated by the hydrodynamic profile [40,41]. Therefore, for the design of the hydrokinetic turbine blade, profiles with a high CL/CD ratio were chosen to be studied. In addition, the selection of a profile with a large section and a considerable thickness to withstand the mechanical forces induced during the operation of the blade is required, avoiding incurring in a drag increase [39][40][41].
The profiles were analyzed for a Reynolds number characteristic of hydrokinetic turbines, which was equal to 750,000. The chord length (C) remained unitary in traditional and multi-element profiles.
CL and CD were analyzed every 1°, within a wide α range. This procedure was carried out for the whole set of the analyzed profiles. The comparison between the diagram of CL and CD versus α for the airfoils used as hydrofoils obtained by numerical methods is illustrated in Figure 1. From the obtained results, Eppler 420 and Selig S1223 hydrofoils with α equal to 16° and 12°, respectively, were determined to be the best profiles. The Selig S1223 hydrofoil resulted to have a CLmax (2.7950) higher than that associated with the Eppler 420 hydrofoil (2.572); nevertheless, Eppler profile exhibited higher CLmax/CD (47.77) ratio compared to Selig profile (39.59). In addition, Eppler From the obtained results, Eppler 420 and Selig S1223 hydrofoils with α equal to 16 • and 12 • , respectively, were determined to be the best profiles. The Selig S1223 hydrofoil resulted to have a C Lmax (2.7950) higher than that associated with the Eppler 420 hydrofoil (2.572); nevertheless, Eppler profile exhibited higher C Lmax /C D (47.77) ratio compared to Selig profile (39.59). In addition, Eppler profile was thicker and capable of supporting higher hydrodynamic loads during its operation; therefore, this profile was chosen for the blade design of a hydrokinetic turbine of 1 kW. For this case, the average Reynolds number (789749.677) along the blade and the average relative velocity (V rel = 5.517 m/s) were considered to determine C, which was equal to 0.1773 m. V rel and C parameters were used in the CFD simulations of the assessed hydrofoil. From a numerical analysis in CFD, a maximum C L /C D ratio of 39.05 was obtained (being C L and C D equal to 1.42 and 0.036, respectively) under a α value of 3 • for the Eppler 420 hydrofoil. This maximum C L /C D ratio was taken as an optimization parameter during the multi-element hydrofoil analysis. It was hypothesized that there were multi-element hydrofoil configurations able to provide better C L /C D ratios in comparison with those ones produced by the traditional hydrodynamic profiles.
In order to define the optimal geometric configuration of a multi-element hydrofoil from the selected Eppler 420 hydrofoil, an optimization procedure was conducted in this study. The optimization of the multi-element hydrofoil required that the relative position of the different elements was varied. For this purpose, the "gap-overlap definition" was decided to be used. The "gap-overlap definition" utilizes three variables to define the flap positions: "gap", "overlap", and "deflection angle" [34]. The "gap" is defined as the vertical distance between the trailing edge of the main element and the flap; and it is always, by definition, a positive value. In turn, the "overlap" refers to a measurement of the elements overlapping, which is determined along the stowed configuration chord line. The overlap has a positive value when the elements do overlap, whereas a negative value of the overlap indicates an increase in the element separation. The measurements of the gap and overlap were provided in terms of a percentage of the chord length of the main element (gap = % C 1 , overlap = % C 1 ). The third parameter of the definition of the gap-overlap is the "deflection angle" (δ), which is defined as the angle between the main element and the flap chord; δ a positive value when the rotation flap is mobilized in a clockwise direction [29][30][31][32][33][34]42]. In addition to the mentioned variables, the flap chord length (C 2 ) was considered as a variable, which was given as a percentage of the main element chord (C 2 = % C 1 ), like the gap and the overlap. It is noteworthy that α and fluid velocity are also variables to be taken into account. All the parameters considered for optimization purposes are shown in Figure 2. therefore, this profile was chosen for the blade design of a hydrokinetic turbine of 1 kW. For this case, the average Reynolds number (789749.677) along the blade and the average relative velocity (Vrel = 5.517 m/s) were considered to determine C, which was equal to 0.1773 m. Vrel and C parameters were used in the CFD simulations of the assessed hydrofoil. From a numerical analysis in CFD, a maximum CL/CD ratio of 39.05 was obtained (being CL and CD equal to 1.42 and 0.036, respectively) under a α value of 3° for the Eppler 420 hydrofoil. This maximum CL/CD ratio was taken as an optimization parameter during the multi-element hydrofoil analysis. It was hypothesized that there were multielement hydrofoil configurations able to provide better CL/CD ratios in comparison with those ones produced by the traditional hydrodynamic profiles.
In order to define the optimal geometric configuration of a multi-element hydrofoil from the selected Eppler 420 hydrofoil, an optimization procedure was conducted in this study. The optimization of the multi-element hydrofoil required that the relative position of the different elements was varied. For this purpose, the "gap-overlap definition" was decided to be used. The "gapoverlap definition" utilizes three variables to define the flap positions: "gap", "overlap", and "deflection angle" [34]. The "gap" is defined as the vertical distance between the trailing edge of the main element and the flap; and it is always, by definition, a positive value. In turn, the "overlap" refers to a measurement of the elements overlapping, which is determined along the stowed configuration chord line. The overlap has a positive value when the elements do overlap, whereas a negative value of the overlap indicates an increase in the element separation. The measurements of the gap and overlap were provided in terms of a percentage of the chord length of the main element (gap = % C1, overlap = % C1). The third parameter of the definition of the gap-overlap is the "deflection angle" (δ), which is defined as the angle between the main element and the flap chord; δ a positive value when the rotation flap is mobilized in a clockwise direction [29][30][31][32][33][34]42]. In addition to the mentioned variables, the flap chord length (C2) was considered as a variable, which was given as a percentage of the main element chord (C2 = % C1), like the gap and the overlap. It is noteworthy that α and fluid velocity are also variables to be taken into account. All the parameters considered for optimization purposes are shown in Figure 2. In this work, a C2 in the range from 30% to 75% of C1, a α between -5° and 20°, a δ between 10° and 30°, a gap between 1% and 5% of C1 and an overlap (ovl) in the range from -5% to 20% of C1 were proposed to be used. The limits of the values were obtained from previous experiences and recommendations found in the literature [22]. In this work, a C 2 in the range from 30% to 75% of C 1 , a α between -5 • and 20 • , a δ between 10 • and 30 • , a gap between 1% and 5% of C 1 and an overlap (ovl) in the range from -5% to 20% of C 1 were proposed to be used. The limits of the values were obtained from previous experiences and recommendations found in the literature [22].
The optimization objectives of this study were to increase the lift and reduce the drag simultaneously. It is noteworthy that the lift was measured from C L , which is defined by Equation (1): where ρ ∞ and U ∞ are the fluid density and velocity, respectively; S refers to the hydrofoil surface; F ⊥U ∞ is a component of the fluid dynamic force orthogonal to the upward flow direction. In turn, the drag was measured by C D , which is defined as expressed in Equation (2): where F is the component of the fluid dynamic force in the upward flow direction. Cavitation inception was assumed to occur on the hydrofoil when the local pressure on the blade section was below the vapor pressure of the fluid, and it was predicted from the pressure distribution [43]. The cavitation number (σ) was defined as described in Equation (3): where P 0 is the absolute pressure (P 0 = P A + ρgh), P A is the atmospheric pressure and ρgh is the gauge pressure, which is calculated as the product among the water density (ρ), the distance between the free water surface and the hydrokinetic rotor center (h) and the local gravity (g); P V is the vapor pressure at the flow temperature; and V refers to the fluid velocity. The pressure coefficient (C pre ) was defined as expressed by Equation (4): It is important to note that cavitation inception can be predicted from the pressure distribution, since cavitation occurs when the local pressure (P L ) is equal to P V or when the minimum negative pressure coefficient (|minC pre |) is equal to σ [44]. By reviewing atmospheric pressures and temperatures of the regions of interest [45], in addition to previous works [46], a σ value equal to 4 was established. The value of σ was included within the surrogate model as a nonlinear restriction.
The multi-objective optimization problem can be mathematically defined as represented by Equation (5), which was subjected to the set of restrictions defined according to Equation (6). maxC L . or min − C L . minC D

Multi-Element Hydrofoil Optimization Framework
For the optimization of the multi-element profile, V rel and C were taken as the parameters for executing 2D CFD simulations using Ansys Fluent software [47] with the k-ω SST turbulence model. This turbulence model is commonly used for hydrokinetic turbine modelling [48][49][50][51][52][53][54] because the referred model has demonstrated higher performance for complex flows, including adverse pressure gradients and flow separations, as occurs in horizontal-axis hydrokinetic turbines. The k-ω SST turbulence model offers an improved prediction of adverse pressure gradients in the near wall regions when compared to the standard k-ω and k-ε models [48]. The flowchart concerning the optimization methodology used for the design of the Eppler 420 hydrofoil, considered as a multi-element hydrofoil, is shown in Figure 3. The optimization methodology is based on studies reported in the literature [27,30,44,45] and consists of the steps described below.

1.
Initial sampling plan. The substitute model must be first trained using a series of initial simulations, whose evaluation is expensive. These start points are defined by a design of experiment (DoE) technique and should be kept to a minimum. For the current study, a Latin Hypercube Sampling (LHS) was used with 100 points optimized according to the Morris-Mitchell criterion to ensure a uniform distribution of the sample points in the design space [28]. LHS can cover the whole design space to randomly sample and effectively simulate the sample output [55]. The design space is the set of all possible combinations of the design variables that are involved in the multi-element hydrofoil design. methodology is based on studies reported in the literature [27,30,44,45] and consists of the steps described below. 1. Initial sampling plan. The substitute model must be first trained using a series of initial simulations, whose evaluation is expensive. These start points are defined by a design of experiment (DoE) technique and should be kept to a minimum. For the current study, a Latin Hypercube Sampling (LHS) was used with 100 points optimized according to the Morris-Mitchell criterion to ensure a uniform distribution of the sample points in the design space [28]. LHS can cover the whole design space to randomly sample and effectively simulate the sample output [55]. The design space is the set of all possible combinations of the design variables that are involved in the multi-element hydrofoil design. 2. CFD simulations. For the CFD simulations of the multi-element hydrofoil, was set at 5.517 m/s and α was varied within the previously defined range for the 2D CFD simulation, using  In order to build a surrogate model, a limited number of CFD simulations was used. The surrogate model can be explored through evolutive algorithms to find the optimal solutions [56]. 3. Mathematical model. In surrogate-based optimization, the surrogate replaces one or more of the objective functions, and the search for the optimum is, therefore, carried out throughout the surrogate model. It must be noted that the surrogate model has to be previously constructed based on a limited, but carefully chosen, number of runs of the original sampling plan function [56]. In the current study, it was decided to use Gaussian processes or Kriging models [28], which take the data of the parameters and the results of the CFD simulations to create a surrogate model. In the design space, the set of non-dominated solutions of the surrogate model lies on a surface, which is commonly known as the Pareto front. Non-dominated solutions are those ones in which superior solutions do not exist within the design space. There are two popular ways of constructing Pareto sets. The first approach combines the optimization criteria into a single objective function; for this purpose, thresholds and penalty functions are often used, as well as weights for linear combinations of the design parameters. The second way for constructing Pareto sets is by using population-based search schemes by means of utilizing algorithms developed for this purpose. In such schemes, a set of designs is worked on concurrently, which evolved toward the final Pareto set in one process. For this, designs are compared to each other and progressed whether they are of high quality and are widely spaced apart from other competing designs. Moreover, an explicit weighting function is not usually required by the referred schemes to combine the objective functions of interest [28,56]. Once the surrogate-based optimization has provided the set of Pareto-optimal solutions, the most intuitive and simplest way of testing the validity of the surrogate is by running additional analyses of the high-fidelity model of the expensive function on representative points of the Pareto set and comparing the outcomes to the predictions. When the comparison is not a satisfactory one, an update of the surrogate is required. The update of the surrogate can be conducted by a simply re-calibration of the surrogate to the set of the original sampling points, plus the additional sampling points. Steps 4, 5, and 6 explain this process [56]. 4. Search. From the surrogate model, new design points are created by using a genetic algorithm (GA). For this propose, the multi-objective GA of the gamultiobj function of the Matlab ® software [37] is used. The goal of the algorithm is to find a set of optimal solutions along the Pareto front for a combination of criteria. The initial population size was equal to 20. This number was chosen In order to build a surrogate model, a limited number of CFD simulations was used. The surrogate model can be explored through evolutive algorithms to find the optimal solutions [56].

3.
Mathematical model. In surrogate-based optimization, the surrogate replaces one or more of the objective functions, and the search for the optimum is, therefore, carried out throughout the surrogate model. It must be noted that the surrogate model has to be previously constructed based on a limited, but carefully chosen, number of runs of the original sampling plan function [56].
In the current study, it was decided to use Gaussian processes or Kriging models [28], which take the data of the parameters and the results of the CFD simulations to create a surrogate model. In the design space, the set of non-dominated solutions of the surrogate model lies on a surface, which is commonly known as the Pareto front. Non-dominated solutions are those ones in which superior solutions do not exist within the design space. There are two popular ways of constructing Pareto sets. The first approach combines the optimization criteria into a single objective function; for this purpose, thresholds and penalty functions are often used, as well as weights for linear combinations of the design parameters. The second way for constructing Pareto sets is by using population-based search schemes by means of utilizing algorithms developed for this purpose. In such schemes, a set of designs is worked on concurrently, which evolved toward the final Pareto set in one process. For this, designs are compared to each other and progressed whether they are of high quality and are widely spaced apart from other competing designs. Moreover, an explicit weighting function is not usually required by the referred schemes to combine the objective functions of interest [28,56].
Once the surrogate-based optimization has provided the set of Pareto-optimal solutions, the most intuitive and simplest way of testing the validity of the surrogate is by running additional analyses of the high-fidelity model of the expensive function on representative points of the Pareto set and comparing the outcomes to the predictions. When the comparison is not a satisfactory one, an update of the surrogate is required. The update of the surrogate can be conducted by a simply re-calibration of the surrogate to the set of the original sampling points, plus the additional sampling points. Steps 4, 5, and 6 explain this process [56].

4.
Search. From the surrogate model, new design points are created by using a genetic algorithm (GA). For this propose, the multi-objective GA of the gamultiobj function of the Matlab ® software [37] is used. The goal of the algorithm is to find a set of optimal solutions along the Pareto front for a combination of criteria. The initial population size was equal to 20. This number was chosen by multiplying the number of free variables (5 parameters in the current study) by a factor of 4 [36,38]. The total number of generations defined in GA was equal to 100 [34].

5.
Evaluation of new designs. When obtaining the optimal design points of the GA [30,57,58], the three design points of the Pareto front with the highest C L , C D , and C L /C D ratio were evaluated in CFD. This process tends to improve the quality of the surrogate model, and it is useful for reducing a set of candidates prior to further CFD analysis [30,57,58]. For the design point with the best C L /C D ratio, additional CFD studies were carried out by varying α in the integer values close to the α given by GA for the referred design point up to a maximum C L /C D ratio of the studied geometry configuration is achieved. 6.
Addition of new design points. Once the results of the CFD simulations of the new design points are obtained, the data are added to the initial sampling to create a new surrogate model and an optimization cycle until the stop criterion is met [28,57]. The purpose of this step is to add points for creating a new surrogate model providing a more optimal objective function. 7.
Stop criterion. During this stage, the same number of new designs points than those ones considered in the initial sampling plan were assessed. Therefore, a total of 200 CFD simulations were considered in order to find the optimal design point that defines the geometric configuration of the multi-element hydrofoil. A proper hydrofoil for the hydrokinetic turbine application must have a high C L /C D ratio for improving the performance, and a high C pre (lower suction) on the suction side to prevent cavitation. After 200 iterations, the optimized multi-element hydrofoil was defined by the best design point of the last Pareto front (Pareto optimal front) that achieved the optimization requirements (maximum C L and minimum C D ), which were subjected to the considered constraints.

Results and Discussion
Through the established surrogate model, a Pareto front was constructed by using GA, as illustrated in Figure 5. The two axes represent the two objective functions that must be simultaneously minimized; therefore, the value of C L was reported as negative. In the figure, the results concerning the initial sampling, the design suggested by the surrogate model, the starting design (Eppler 420 hydrofoil), and the selected multi-element design based on the C L /C D ratio are depicted.
From Figure 5, it can be observed that few of the initial designs contribute to the Pareto front, and some of them grant a better C L /C D ratio than those ones corresponding to the starting Eppler 420 hydrofoil. Additionally, the designs supplied by the surrogate model contributed to the Pareto front with new designs that fill the gaps in the Pareto front of the initial sampling plan and move it forward.
In general terms, the surrogate model produced solutions that did not infringe the cavitation restriction (|min C pre | ≤ 4), being only five designs proposed by the surrogate model that did not meet this restriction. Figure 6 shows the pressure coefficient distributions (C pre ) for the traditional Eppler 420 hydrofoil and for the original and optimized multi-element hydrofoil. C pre values were observed to be higher in the multi-element hydrofoil than those ones in the traditional hydrofoil. In fact, C pre is highly decreased on the upper surface near the leading edge of the traditional Eppler 420 hydrofoil in comparison with the values obtained for the multi-element hydrofoil; then, the flow might not provide enough kinetic energy to withstand the adverse pressure gradient downstream and it would separate. All the tested hydrofoils satisfied the cavitation constraint [11]. From Figure 5, it can be observed that few of the initial designs contribute to the Pareto front, and some of them grant a better CL/CD ratio than those ones corresponding to the starting Eppler 420 hydrofoil. Additionally, the designs supplied by the surrogate model contributed to the Pareto front with new designs that fill the gaps in the Pareto front of the initial sampling plan and move it forward.
In general terms, the surrogate model produced solutions that did not infringe the cavitation restriction (|min Cpre|≤4), being only five designs proposed by the surrogate model that did not meet this restriction. Figure 6 shows the pressure coefficient distributions (Cpre) for the traditional Eppler 420 hydrofoil and for the original and optimized multi-element hydrofoil. Cpre values were observed to be higher in the multi-element hydrofoil than those ones in the traditional hydrofoil. In fact, Cpre is highly decreased on the upper surface near the leading edge of the traditional Eppler 420 hydrofoil in comparison with the values obtained for the multi-element hydrofoil; then, the flow might not provide enough kinetic energy to withstand the adverse pressure gradient downstream and it would separate. All the tested hydrofoils satisfied the cavitation constraint [11]. The traditional Eppler 420 hydrofoil and the initial and optimized multi-element hydrofoils are shown in Figure 7. In Table 2, the starting Eppler 420 hydrofoil is compared to the multi-element design selected as optimal, which exhibited a gap of 2.825% C1, an overlap of 8.52% C1, a δ of 19.765˚, a C2 of 42.471% C1, and a α of -4˚. Based on the CL/CD ratio, the selected multi-element design had a  The traditional Eppler 420 hydrofoil and the initial and optimized multi-element hydrofoils are shown in Figure 7. In Table 2, the starting Eppler 420 hydrofoil is compared to the multi-element design selected as optimal, which exhibited a gap of 2.825% C 1 , an overlap of 8.52% C 1 , a δ of 19.765 • , a C 2 of 42.471% C 1 , and a α of -4 • . Based on the C L /C D ratio, the selected multi-element design had a better performance than that achieved by the Eppler 420 (i.e., the traditional one). The traditional Eppler 420 hydrofoil and the initial and optimized multi-element hydrofoils are shown in Figure 7. In Table 2, the starting Eppler 420 hydrofoil is compared to the multi-element design selected as optimal, which exhibited a gap of 2.825% C1, an overlap of 8.52% C1, a δ of 19.765˚, a C2 of 42.471% C1, and a α of -4˚. Based on the CL/CD ratio, the selected multi-element design had a better performance than that achieved by the Eppler 420 (i.e., the traditional one).    The pressure ( Figure 8) and velocity ( Figure 9) contours of the traditional and optimized multi-element hydrofoils were compared. In the pressure contours, the multi-element hydrofoil was observed to provide lower pressures than the traditional one, resulting in a |min C pre | higher than in the traditional hydrofoil. Additionally, an increase in the pressure was observed at the stagnation point. In the traditional and multi-element hydrofoil (main element), due to the increase in velocity, pressure diminished in the upper region and increased in the lower region. In the multi-element hydrofoil, the flow is trapped beneath the hydrofoil, leading to a decrease in the flow velocities and a build-up of the pressure below the hydrofoil. Therefore, there is a higher pressure difference that is available for the lift generation. In the velocity contours, it can be observed how the fluid leaves the trailing edge of the main element and is deflected and accelerated, adhering to the flap and, subsequently, postponing the occurrence of the boundary layer separation, which only occurs in a small area near the trailing edge of the flap. In the traditional hydrofoil, the boundary layer separation is much larger than in the multi-element hydrofoil. The delay in the separation of the boundary layer is one of the factors that contributes to the multi-element hydrofoil, achieving a higher C L than that observed in the traditional one [59].
trailing edge of the main element and is deflected and accelerated, adhering to the flap and, subsequently, postponing the occurrence of the boundary layer separation, which only occurs in a small area near the trailing edge of the flap. In the traditional hydrofoil, the boundary layer separation is much larger than in the multi-element hydrofoil. The delay in the separation of the boundary layer is one of the factors that contributes to the multi-element hydrofoil, achieving a higher CL than that observed in the traditional one [59].  In order to obtain a more-detailed knowledge of the selected multi-element hydrofoil performance, in Figure 10, the CL/CD ratio with respect to α for the Eppler 420 one and the multielement optimized hydrofoil is represented. As observed, the CL/CD ratio is higher in the multielement hydrofoil in comparison with the traditional one for several points. The flap produces the highest maximum lift for a lower α value compared to the traditional section [59]. Regarding the selected multi-element hydrofoil, a difference of 7° was observed for α when maximizing the CL/CD ratio compared to the Eppler 420 hydrofoil. In order to obtain a more-detailed knowledge of the selected multi-element hydrofoil performance, in Figure 10, the C L /C D ratio with respect to α for the Eppler 420 one and the multi-element optimized hydrofoil is represented. As observed, the C L /C D ratio is higher in the multi-element hydrofoil in comparison with the traditional one for several points. The flap produces the highest maximum lift for a lower α value compared to the traditional section [59]. Regarding the selected multi-element hydrofoil, a difference of 7 • was observed for α when maximizing the C L /C D ratio compared to the Eppler 420 hydrofoil. As a performance measurement of the surrogate model, the hypervolume (Hv) values are reported [60]. The hypervolume criterion is based on the theory of the hypervolume indicator, which is defined as a metric exhibited by the non-dominated solutions. This metric consists of the size of the hypervolume fronted by the non-dominated set bounded by reference maximum points [61,62].
This performance measurement estimates the non-overlapping volume of the set of hypercubes formed by the reference point (CL = 4 and CD = 0.14) and every vector in the Pareto set approximation [59]. This performance measurement is congruent with the Pareto front [63,64], and it is used to assess both convergence and maximum dispersion of the solutions from the approximation of the Pareto front obtained throughout the surrogate model. Figure 11 illustrates Hv convergence history, which indicates that the convergence tends to 5.58 × 10 Hv value. High Hv values of this measurement indicate that the solutions are closer to the true Pareto front and cover a wider extension of it. As a performance measurement of the surrogate model, the hypervolume (Hv) values are reported [60]. The hypervolume criterion is based on the theory of the hypervolume indicator, which is defined as a metric exhibited by the non-dominated solutions. This metric consists of the size of the hypervolume fronted by the non-dominated set bounded by reference maximum points [61,62].
This performance measurement estimates the non-overlapping volume of the set of hypercubes formed by the reference point (C L = 4 and C D = 0.14) and every vector in the Pareto set approximation [59]. This performance measurement is congruent with the Pareto front [63,64], and it is used to assess both convergence and maximum dispersion of the solutions from the approximation of the Pareto front obtained throughout the surrogate model. Figure 11 illustrates Hv convergence history, which indicates that the convergence tends to 5.58 × 10 Hv value. High Hv values of this measurement indicate that the solutions are closer to the true Pareto front and cover a wider extension of it.
It is widely known that, during the design optimization of the hydrokinetic turbine rotor, maximizing the power production is a crucial aspect. In this regard, the selection of the hydrofoil has an important effect on the performance of the system. This study shows that the multi-element hydrofoil can be applied for the design of a hydrokinetic turbine due to multi-element hydrofoil configuration has shown to be a promising positive effect on improving the hydrodynamic characteristics in comparison with the traditional hydrofoil configuration. This fact can be reflected in the overall performance enhancement of the hydrokinetic turbine, as it has been previously reported for wind turbines [21][22][23][24]. Therefore, multi-element hydrofoils are highly recommended for the blade design despite the fact that the manufacturing process costs associated with it could be higher than those generated during the manufacture of standard hydrofoils. It is widely known that, during the design optimization of the hydrokinetic turbine rotor, maximizing the power production is a crucial aspect. In this regard, the selection of the hydrofoil has an important effect on the performance of the system. This study shows that the multi-element hydrofoil can be applied for the design of a hydrokinetic turbine due to multi-element hydrofoil configuration has shown to be a promising positive effect on improving the hydrodynamic characteristics in comparison with the traditional hydrofoil configuration. This fact can be reflected in the overall performance enhancement of the hydrokinetic turbine, as it has been previously reported for wind turbines [21][22][23][24]. Therefore, multi-element hydrofoils are highly recommended for the blade design despite the fact that the manufacturing process costs associated with it could be higher than those generated during the manufacture of standard hydrofoils.

Conclusions
From the initial sampling obtained throughout LHS, multi-element hydrofoils with a higher CL/CD ratio than those ones provided by the traditional starting hydrofoil were obtained, which proves the hypothesis that multi-element profiles can provide better properties than the traditional ones. Nevertheless, the multi-element hydrofoil parameters must be carefully selected, because several combinations of these parameters can result in a multi-element hydrofoil that does not exhibit a performance higher than that exhibited by the traditional one. Indeed, in the Pareto front, several geometric configurations of the multi-element hydrofoil providing a lower performance than that obtained by the traditional Eppler hydrofoil were evidenced.
The simulation results showed that the optimized multi-element had a higher CL/CD ratio at a low value. In the multi-element hydrofoil, a larger degree of flow separation on the upper surface of the main element and a lower pressure difference for higher α values were observed. This could lead to a lower lift generation in the multi-element hydrofoil.
By obtaining better multi-element hydrofoils than in the initial sampling plan collection by using the proposed surrogate model, it was demonstrated that the use of a surrogate model is a viable and Figure 11. Hv convergence history.

Conclusions
From the initial sampling obtained throughout LHS, multi-element hydrofoils with a higher C L /C D ratio than those ones provided by the traditional starting hydrofoil were obtained, which proves the hypothesis that multi-element profiles can provide better properties than the traditional ones. Nevertheless, the multi-element hydrofoil parameters must be carefully selected, because several combinations of these parameters can result in a multi-element hydrofoil that does not exhibit a performance higher than that exhibited by the traditional one. Indeed, in the Pareto front, several geometric configurations of the multi-element hydrofoil providing a lower performance than that obtained by the traditional Eppler hydrofoil were evidenced.
The simulation results showed that the optimized multi-element had a higher C L /C D ratio at a low α value. In the multi-element hydrofoil, a larger degree of flow separation on the upper surface of the main element and a lower pressure difference for higher α values were observed. This could lead to a lower lift generation in the multi-element hydrofoil.
By obtaining better multi-element hydrofoils than in the initial sampling plan collection by using the proposed surrogate model, it was demonstrated that the use of a surrogate model is a viable and effective alternative for the design of multi-element hydrofoils for horizontal-axis hydrokinetic turbines.
Properly designed multi-element profiles have higher C L , owing to the fact that the flap increases the contact surface and the deviation of the flow caused by the main element improving the flap, decreasing the flow separation and resulting in an insignificant C D decrease. By increasing C L and C L /C D ratio values, the amount of torque is increased, which is of special interest for the hydrokinetic turbine performance.