Effect of Blade Leading and Trailing Edge Conﬁgurations on the Performance of a Micro Tubular Propeller Turbine Using Response Surface Methodology

: With the recent rise in importance of environmental issues, research on micro hydropower, a kind of renewable energy source, is being actively conducted. In this study, a micro tubular propeller turbine was selected for study of micro hydropower in pipes. Numerical analysis was conducted to evaluate the performance. Response surface methodology using design of experiments was performed to efﬁciently investigate the effect of the blade leading and trailing edge elliptic aspect ratios on the performance. The trailing edge conﬁguration was found to be more related to the performance, because of the drastic pressure variation due to the stagnation point formed, regardless of the leading edge conﬁguration. To improve the performance, a NACA airfoil was introduced. The results show that the ﬂow became more stable than the reference model, and the efﬁciency was increased by 2.44%.


Introduction
In 2015, the Paris Agreement was adopted to reduce greenhouse gas emissions as climate change was becoming a more serious problem worldwide [1]. However, even after the adoption of the treaty, global greenhouse gas emissions have been increasing year by year, and are projected to increase in the future [2]. Greenhouse gases emitted worldwide mostly come from the energy sector (73%), and in particular, the energy use in industry accounts for 24%, which is the largest part [3,4]. Therefore, it is urgent to replace fossil fuel energy with clean energy that does not produce carbon dioxide; but in 2019, renewable energy accounted for only 27.3% of the global energy production [5]. Hence, the importance of renewable energy is expected to increase further.
Renewable energy, which is often referred to as clean energy, is defined to be energy collected from renewable sources. It includes hydropower, solar energy, wind energy, biomass energy and geothermal energy [6,7]. In 2019, hydropower accounted for 50% of the total renewable energy capacity across the globe [8]. This indicates that hydropower has steadily gained ground in the renewable energy market. However, only a fraction of the hydropower energy potential is currently being developed and used [9]. Accordingly, it is necessary to further research hydropower.
Hydropower generates electricity by rotating a hydro turbine using a head of the water, and it has the advantage that it does not produce carbon dioxide at all, while offering high energy density [10]. Micro hydropower, which is typically classified as hydropower with less than 100 kW of installed capacity [11,12], has the advantage of longer life, lower operation costs, and fewer restrictions on space, relative to macro hydropower [13]. In particular, the in-line type hydro turbine for micro hydropower is highly safe and sustainable, because regardless of the external environment, it can supply water at a constant flow rate [13,14]. In addition, it can be applied to wherever there is a pipeline, such as water and sewage, agricultural water, plant cooling water, and high-rise building drainage.
Research on the in-line type hydro turbine for micro hydropower has been conducted for many years, and various types of hydro turbine have been proposed. Since the inline type turbine is a reaction turbine that utilizes pressure energy remaining in the pipe, classical turbines such as Francis and Kaplan turbines for small hydropower have been researched [15][16][17][18]. A pump as turbine (PAT) is also known to be suitable for small hydropower applications [19] and has been investigated in a variety of forms such as mixed flow or centrifugal types [20,21]. Furthermore, in addition to classic turbines, many new types of turbines have been developed. Chen et al. [22] developed a vertical axis turbine to generate power from water pipelines, and conducted numerical analysis and experimental tests for over 20 models through design variables, such as drag-type or lift-type turbine, and deflector configurations. Payambarpour et al. [23] used a Savonius turbine as an in-line type hydro turbine, and improved the performance by adjusting design variables, such as deflector configurations and aspect ratio, through numerical analysis and experiments. A helical turbine, known as the Gorlov vertical axis turbine, was introduced into the LucidPipeTM system to generate energy [24]. Yeo et al. [25] studied a vertical axis turbine based on a helical turbine numerically, and investigated the effect of tip-speed ratios on the performance. A hydrocoil turbine, which is an in-line type turbine version of the Archimedes screw turbine, was developed by HydroCoil Power, Inc. [26]. Sinagra et al. [27] adopted an in-line type crossflow turbine for pressure regulation and energy production, and verified the performance through numerical and experimental analysis. Research on contra-rotating turbines for micro hydropower was conducted, and the performance was evaluated through numerical analysis and experiments [28,29].
Among the various types of in-line type turbines, a micro tubular propeller turbine was selected for this study. A tubular turbine, often referred to as a bulb turbine, is a kind of axial turbine that is usually used in low head conditions. Unlike a propeller or Kaplan turbine, the tubular turbine does not have a spiral casing, and so it is easy to design an in-line hydro turbine [30]. Studies on hydropower using the tubular turbine have been actively conducted before [31][32][33]. However, typical tubular turbines have the disadvantage of complicated structures with generators installed inside the pipes and maintenance problems due to clogging [34]. Samora et al. [35] designed a micro tubular propeller turbine that can be installed in pipes with a diameter of less than 200 mm. It is a miniaturized tubular turbine whose generator is located outside the pipe like an S-turbine. In addition, it does not include a guide vane, so it has the advantage of a simpler configuration and reduction of clogging issues, enabling easier maintenance. In this study, the effect of blade configuration variation on performance was investigated through numerical analysis, and it was expected that modifying the blade configuration would improve the performance.

Model Description
A reference model of the micro tubular propeller turbine was designed referring to the precedent study [35]. Figures 1 and 2 show a system schematic and a runner model, respectively, and Table 1 describes the specification. The design of an in-line type hydro turbine while locating the generator outside of the flow path necessitates the use of bent pipes. The bend angle of the pipe is set to 45 • to minimize the head loss inevitably generated.

Evaluation Indicator
Among the output values obtained from numerical analysis, the value adopted as an indicator to evaluate the performance is efficiency. The hydraulic efficiency of the hydro turbine can be calculated as the ratio of the mechanical power generated by the turbine shaft to the maximum power that can be available at the inlet of the runner, as indicated in Equation (3) [36]. Here, the mechanical power can be obtained by multiplying the axial torque and rotational speed, as expressed in Equation (1). In addition, the maximum power available at the runner inlet can be acquired by the product of the density of the working fluid, the gravitational acceleration, the flow rate, and the head. In this case, the head can be obtained using the difference of static pressure, ∆ at the inlet and outlet of the runner, as described in Equation (2).

Evaluation Indicator
Among the output values obtained from numerical analysis, the value adopted as an indicator to evaluate the performance is efficiency. The hydraulic efficiency of the hydro turbine can be calculated as the ratio of the mechanical power generated by the turbine shaft to the maximum power that can be available at the inlet of the runner, as indicated in Equation (3) [36]. Here, the mechanical power can be obtained by multiplying the axial torque and rotational speed, as expressed in Equation (1). In addition, the maximum power available at the runner inlet can be acquired by the product of the density of the working fluid, the gravitational acceleration, the flow rate, and the head. In this case, the head can be obtained using the difference of static pressure, ∆p at the inlet and outlet of the runner, as described in Equation (2).
where the terms are as defined in the nomenclature below.

Grid Generation
For numerical analysis, the model representing the fluid zone of the micro tubular propeller turbine system was constructed with a grid system. The pipe part was composed of a mixture of structured and unstructured mesh using ANSYS Mesh, and the runner part was composed of structured mesh using ANSYS TurboGrid. To resolve the problem of node mismatch at the interface between pipe and runner, the grid at the interface was densely generated. Y-plus is the dimensionless wall distance, which can be calculated as y + = ρyu τ µ where ρ is the working fluid density, µ is the working fluid dynamic viscosity, y is the wall distance, u τ = τ w ρ is the wall shear velocity, and τ w is the wall shear stress. In this study, y-plus was limited below 2 to consider the shear stress on the wall, such as blades.
A grid dependency test was performed to minimize the effect of the grid size on the results, and consider the economics of the numerical analysis time. Figure 3 shows the results of the variation in the power, head, and efficiency represented in Equations (1)-(3) by the number of elements. The results show that those values do not change greatly from more than about 5 million grid elements. Therefore, the number of grid elements used in this study was set to approximately 5 million.

=
(3) where the terms are as defined in the nomenclature below.

Grid Generation
For numerical analysis, the model representing the fluid zone of the micro tubular propeller turbine system was constructed with a grid system. The pipe part was composed of a mixture of structured and unstructured mesh using ANSYS Mesh, and the runner part was composed of structured mesh using ANSYS TurboGrid. To resolve the problem of node mismatch at the interface between pipe and runner, the grid at the interface was densely generated. Y-plus is the dimensionless wall distance, which can be calculated as = where is the working fluid density, is the working fluid dynamic viscosity, is the wall distance, = is the wall shear velocity, and is the wall shear stress. In this study, y-plus was limited below 2 to consider the shear stress on the wall, such as blades.
A grid dependency test was performed to minimize the effect of the grid size on the results, and consider the economics of the numerical analysis time. Figure 3 shows the results of the variation in the power, head, and efficiency represented in Equations (1)-(3) by the number of elements. The results show that those values do not change greatly from more than about 5 million grid elements. Therefore, the number of grid elements used in this study was set to approximately 5 million.

Boundary Conditions
A numerical analysis was conducted using ANSYS CFX ver. 2020 R2, and Table 2 describes the boundary conditions applied in this study. Here, the inlet velocity, according to the inlet condition is 0.78 m/s and the Reynolds number, = ⁄ is 7.43 × 10 4 , which indicates high turbulence. Due to the absence of a guide vane in this model, the water flows straight until just before the runner entrance. For this reason, the steady state with frozen rotor interface was set to an analysis state, because the flow is uniformly distributed to the blades, and there is no wake at the runner inlet. With respect to the turbulence model, shear stress transport (SST), which is known to be able to simulate the secondary flow at the blade [37,38], was selected to consider the effects of adverse pressure gradient on the blade.

Boundary Conditions
A numerical analysis was conducted using ANSYS CFX ver. 2020 R2, and Table  2 describes the boundary conditions applied in this study. Here, the inlet velocity, V ∞ according to the inlet condition is 0.78 m/s and the Reynolds number, Re D = ρV ∞ D p /µ is 7.43 × 10 4 , which indicates high turbulence. Due to the absence of a guide vane in this model, the water flows straight until just before the runner entrance. For this reason, the steady state with frozen rotor interface was set to an analysis state, because the flow is uniformly distributed to the blades, and there is no wake at the runner inlet. With respect to the turbulence model, shear stress transport (SST), which is known to be able to simulate the secondary flow at the blade [37,38], was selected to consider the effects of adverse pressure gradient on the blade.

Governing Equations
The governing equations for the incompressible Reynolds-Averaged Navier-Stokes (RANS) model to express the motion of fluid are the following, including continuity and momentum equations [39,40]. Here, the Einstein notation was used to denote the coordinate axis.
where the mean strain-rate tensor S ij = 1 , and the Reynolds stresses In this case, µ is the dynamic viscosity, µ t is the turbulent dynamic viscosity, k is the turbulence kinetic energy, and δ ij is the Kronecker delta.

Validation
In order to verify the validity of the steady state numerical analysis, the results were compared to those of the unsteady state numerical analysis with unsteady rotor-stator interface, which is known to more realistically simulate the flow. In this case, the result values were obtained by calculating the average values during one revolution of the turbine. Table 3 states the comparative results, and the relative errors of the power, head, and efficiency are 0.82%, 0.53%, and 0.29%, respectively, showing very small values of less than 1%. In addition, the results of the steady state numerical analysis were compared with the experimental results of the precedent study [35], and the relative errors of the power, head, and efficiency are 0.32%, 2.94%, and 2.54%, respectively, indicating small relative errors of less than 3%. Therefore, it can be determined that the numerical analysis method applied in this study has validity.

Design of Experiments
Design of experiments (DOE) was conducted to investigate the effect of the variation of the elliptic aspect ratios of the blade leading and trailing edge (see Figure 4) on the performance of the hydro turbine. The DOE is a method of logically analyzing the results, using statistical theory by efficiently placing experimental conditions, rather than performing experiments on all design variables within the range of the experimental area of interest. Unlike physical experiments, there is no random error in computational experiments; therefore, it is known that sampling methods which uniformly fill the experimental areas, the so-called "space-filling" concept, are appropriate [41]. Latin hypercube sampling (LHS) is a method of dividing each of the dimensions in the experimental design into regions with equal levels and extracting one sample from each region [42]. It is known to be a suitable sampling method for computational experiments, because it satisfies the concept of space-filling. One disadvantage of LHS is that it places samples randomly in the experimental area, hence it is likely that the distributed samples will be uneven. Optimal space-filling (OSF) is a sampling method based on LHS, but by maximizing the minimum distance between the points, it provides a more uniform distribution within the design space [43]. In this study, the design variables were sampled using OSF, and Figure 5b shows the distribution of extracted design variables within the experimental area. However, OSF does not include endpoints in the experimental area, which has the disadvantage of being less accurate at endpoints. Therefore, the samples for DOE added endpoints to the samples from OSF as shown in Figure 5c. The ranges of the elliptic aspect ratios of the blade leading and trailing edge were set to 2- 16. to be a suitable sampling method for computational experiments, because it satisfies the concept of space-filling. One disadvantage of LHS is that it places samples randomly in the experimental area, hence it is likely that the distributed samples will be uneven. Optimal space-filling (OSF) is a sampling method based on LHS, but by maximizing the minimum distance between the points, it provides a more uniform distribution within the design space [43]. In this study, the design variables were sampled using OSF, and Figure  5b shows the distribution of extracted design variables within the experimental area. However, OSF does not include endpoints in the experimental area, which has the disadvantage of being less accurate at endpoints. Therefore, the samples for DOE added endpoints to the samples from OSF as shown in Figure 5c. The ranges of the elliptic aspect ratios of the blade leading and trailing edge were set to 2-16.  concept of space-filling. One disadvantage of LHS is that it places samples randomly in the experimental area, hence it is likely that the distributed samples will be uneven. Optimal space-filling (OSF) is a sampling method based on LHS, but by maximizing the minimum distance between the points, it provides a more uniform distribution within the design space [43]. In this study, the design variables were sampled using OSF, and Figure  5b shows the distribution of extracted design variables within the experimental area. However, OSF does not include endpoints in the experimental area, which has the disadvantage of being less accurate at endpoints. Therefore, the samples for DOE added endpoints to the samples from OSF as shown in Figure 5c. The ranges of the elliptic aspect ratios of the blade leading and trailing edge were set to 2-16.

Response Surface
To analyze the effect of each design variable on the performance, genetic aggregation response surface (GARS) was used to generate the response surface. GARS is a method of adopting the best response surface using the weighted aggregation of several other metamodels. It can be expressed as Equation (6): whereŷ ens is the prediction of the ensemble, N M is the number of metamodels used, w i is the weight factor for the i-th metamodel, andŷ i is the prediction of the i-th metamodel. It corrects the deviation by K-fold cross-validation [44]. The commercial code, ANSYS DesignXplorer, was used to implement GARS.

Optimization
Genetic algorithm (GA) was used to find the optimal point of the elliptic aspect ratios of the blade leading and trailing edge having the highest efficiency in the design range. GA is an algorithm that mimics the evolution of individuals as they adapt to the natural environment. It is suitable for optimization as it is likely to choose the optimal solution Appl. Sci. 2021, 11, 5596 7 of 17 as they evolve. It selects a superior population from a random population and uses the crossover and mutation operators to create the next generation. This process is repeated until the fitness is satisfied [45]. Figure 6 illustrates a general flowchart of GA procedures. To implement the GA, the commercial code, ANSYS DesignXplorer, was utilized.

Optimization
Genetic algorithm (GA) was used to find the optimal point of the elliptic aspect rat of the blade leading and trailing edge having the highest efficiency in the design ran GA is an algorithm that mimics the evolution of individuals as they adapt to the natu environment. It is suitable for optimization as it is likely to choose the optimal solution they evolve. It selects a superior population from a random population and uses the cro over and mutation operators to create the next generation. This process is repeated un the fitness is satisfied [45]. Figure 6 illustrates a general flowchart of GA procedures. implement the GA, the commercial code, ANSYS DesignXplorer, was utilized.  Figure 7 shows the generated response surface for sensitivity analysis of the desi variables. The x-axis and y-axis represent the elliptic aspect ratios of the leading and tra ing edge of the blade, respectively, while the z-axis indicates the power, head, and ef ciency from left to right. Between the elliptic aspect ratios of the leading and trailing ed of the blade, the more sensitive design variable in terms of performance is the latter, a there is a tendency that as the blade edges sharpen, the power and head decrease. In a dition, the response surface in the region of 2-4 for ALE is rather steep, especially for t efficiency. This can be considered for several reasons. Here, the figures used in the exp nation are represented at 50% span which is the middle of the hub and shroud (see Figu 8). First, for the blunt leading edge (particularly ALE = 2 or 3), the flow separation occu  Figure 7 shows the generated response surface for sensitivity analysis of the design variables. The x-axis and y-axis represent the elliptic aspect ratios of the leading and trailing edge of the blade, respectively, while the z-axis indicates the power, head, and efficiency from left to right. Between the elliptic aspect ratios of the leading and trailing edge of the blade, the more sensitive design variable in terms of performance is the latter, and there is a tendency that as the blade edges sharpen, the power and head decrease. In addition, the response surface in the region of 2-4 for A LE is rather steep, especially for the efficiency. This can be considered for several reasons. Here, the figures used in the explanation are represented at 50% span which is the middle of the hub and shroud (see Figure 8). First, for the blunt leading edge (particularly A LE = 2 or 3), the flow separation occurs due to the strong adverse pressure (see Figure 9). This is because of the drastic change of curvature radius. Furthermore, the sharper the trailing edge, the more the flow separation is delayed (see Figure 10), and the lower the pressure on the pressure side (see Figures 11 and 12). In addition, in the case of the leading edge, drastic variation of the pressure is inevitably formed due to the stagnation point, regardless of the change of the leading edge shape of the blade (see Figure 13). The sharper the leading edge, the greater the tendency of the drastic pressure fluctuation on the suction side to be resolved, but the more insignificant the overall pressure difference between the pressure and suction sides becomes (see Figure 14). Moreover, the results can be considered using the velocity triangle at the inlet and outlet of the blade. The water flows in a straight direction until just before the runner entrance, due to the absence of a guide vane. Figure 15 shows that the absolute velocity and blade linear velocity are perpendicular at the leading edge. Equation (7) indicates Euler's turbine equation, which enables calculation of the theoretical turbine output based on Newton's second law of motion [46]:

Sensitivity Analysis
dicates Euler's turbine equation, which enables calculation of the theoretical turbine output based on Newton's second law of motion [46]: In the equation, ⋅ , corresponding to the velocity component at the leading edge, has the value of 0, due to the orthogonality. As a result, the variation of the curvature of the leading edge configuration does not significantly affect the performance.    edge, has the value of 0, due to the orthogonality. As a result, the variation of the curvature of the leading edge configuration does not significantly affect the performance.       (a) ATE = 2 (b) ATE = 16 Figure 11. Pressure distribution at 50% span of the different elliptic aspect ratios of the trailing edge (a,b). Figure 11. Pressure distribution at 50% span of the different elliptic aspect ratios of the trailing edge (a,b).      For these reasons, there are many studies focused on blade trailing edge rather than leading edge to improve blade performance. For a case in point, there is a study that modified the blunt trailing edge thickness of the airfoil to improve the aerodynamic performance of a wind turbine [47]. In addition, the blade trailing edge profile also has a significant effect on pump performance, and there is a case where pressure pulsations in a centrifugal pump were reduced through modification of the blade trailing edge profile [48].

Results of the Modified Model
The elliptic aspect ratios of the blade leading and trailing edge of the reference model are two. The most efficient points of the elliptic aspect ratios of the blade leading and trailing edge derived through the GA were calculated as 5.5 and 16, respectively. To verify the predicted point, numerical analysis was conducted with the calculated elliptic aspect ratios of the blade edges. Table 4 specifies the comparison results, and the relative errors In the equation, V 1 · U 1 , corresponding to the velocity component at the leading edge, has the value of 0, due to the orthogonality. As a result, the variation of the curvature of the leading edge configuration does not significantly affect the performance.
For these reasons, there are many studies focused on blade trailing edge rather than leading edge to improve blade performance. For a case in point, there is a study that modified the blunt trailing edge thickness of the airfoil to improve the aerodynamic performance of a wind turbine [47]. In addition, the blade trailing edge profile also has a significant effect on pump performance, and there is a case where pressure pulsations in a centrifugal pump were reduced through modification of the blade trailing edge profile [48].

Results of the Modified Model
The elliptic aspect ratios of the blade leading and trailing edge of the reference model are two. The most efficient points of the elliptic aspect ratios of the blade leading and trailing edge derived through the GA were calculated as 5.5 and 16, respectively. To verify the predicted point, numerical analysis was conducted with the calculated elliptic aspect ratios of the blade edges. Table 4 specifies the comparison results, and the relative errors are less than 1%; thus we can conclude that the prediction is verified. At this time, the power and head decrease by about 8% and 11%, respectively, compared to the reference model, and the efficiency increases by 2.16%.

Application of the NACA Airfoil
An airfoil refers to the cross-sectional shape of an airplane's wings. In 1929, the National Advisory Committee for Aeronautics (NACA) studied a series of airfoils, and systematically standardized them to be suitable for specific conditions [49]. The NACA airfoils play a role in improving the aerodynamic performance of the airplane, and if applied to the hydro turbine blade, may be expected to stabilize the flow passing through the runner, and improve the performance of the hydro turbine. In particular, the maximum thickness of the NACA 4-digit series airfoil models is located at 30% of the chord from the leading edge [50]. This configuration has a similarity to the results previously derived from the optimization of the blade leading and trailing edge elliptic aspect ratios. Taking these into consideration, a NACA 4-digit series airfoil model was applied to the blade profile of the micro tubular propeller turbine, and numerical analysis was conducted to analyze and compare the performance with the reference model. The blade angle was set to be identical to the reference model. Based on the camber line, the NACA 0004 to the hub face and NACA 0006 to the blade tip were applied to keep the maximum thickness of 1.7 mm, which is the same as the reference model. Figure 16 compares the blade cross-sectional shape on the hub face of the reference, modified, and NACA airfoil model, while Table 5 describes the results of the numerical analysis.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 12 of 18 sectional shape on the hub face of the reference, modified, and NACA airfoil model, while Table 5 describes the results of the numerical analysis.
(a) Reference model (b) Modified model (c) NACA airfoil model     Figure 17 shows the pressure variation along the blade surface. The pressure fluctuation at the trailing edge is resolved for the modified model and the NACA airfoil model, compared to the reference model. In addition, the pressure variation line on the blade surface is the smoothest for the NACA airfoil model, so we can assume that the flow becomes stable. Furthermore, the flow becoming stable can also be determined through Figure 18, which shows the pressure distribution in the flow at the runner. In addition, as can be seen in the velocity distribution shown in Figure 19, the wake generated at the trailing edge is the weakest for the NACA airfoil, which results in the least influence from the wake.
The cases of applying a NACA airfoil to the hydro turbine can support the results of this study. For instance, a NACA airfoil is applied to the guide vane in the Francis turbine [51,52], and also to the runner blade [53]. In addition, there are examples of applications to the axis tidal turbine [54,55].  (c) Trailing edge

Conclusions
In this study, a micro tubular propeller turbine with a simple configuration and easy maintenance was selected for micro hydropower in pipes. The effect of the blade leading

Conclusions
In this study, a micro tubular propeller turbine with a simple configuration and easy maintenance was selected for micro hydropower in pipes. The effect of the blade leading and trailing edge configuration on the performance of the hydro turbine was analyzed numerically. The variation of the performance depending on the design variables was efficiently studied through DOE and GARS, and it turned out that the trailing edge configuration is more sensitive to the performance. This is because the flow separation delays at the sharper trailing edge, and the sudden change of pressure does not disappear even if the leading edge configuration has changed due to the inevitably formed stagnation point. To improve the performance of the hydro turbine, a NACA airfoil was applied to the blade of hydro turbine, and it was effective in stabilizing the flow. It results in 2.44% increased efficiency compared to the reference model.