Power Generation Enhancement of Horizontal Axis Wind Turbines Using Bioinspired Airfoils: A CFD Study

Abstract

This research investigates the performance implications of employing a bioinspired airfoil (seagull’s wing cross-section) in horizontal-axis wind turbines. Specifically, we replaced the S809 airfoil from NREL Phase VI with an airfoil modeled after seagull wings. Initially, we calibrated four coefficients of the GEKO turbulence model for both the S809 and the bioinspired airfoil, utilizing experimental data. Subsequently, using the calibrated generalized k-ω (GEKO) model, we conducted a comparative analysis between the S809 and the seagull airfoils, revealing the considerable superiority of the seagull airfoil in terms of lift and drag coefficients. Furthermore, we numerically simulated the original NREL Phase VI turbine and a modified version where the S809 airfoil was replaced with the seagull airfoil using 3D computational fluid dynamics (CFD) with the airfoil-based-calibrated GEKO turbulence model. This investigation spanned a wide range of air speeds, including 7 m/s, 13 m/s, and 25 m/s. At these wind speeds, we observed a substantial increase in turbine power generation, with enhancements of 47.2%, 204.4%, and 103.9%, respectively. This study underscores the significant influence of nature’s designs in advancing energy extraction within industries, particularly within the wind energy sector.

1. Introduction

Until a few decades ago, horizontal-axis wind turbines (HAWTs) primarily employed airfoils originally developed for aviation applications, including the well-known NACA and NASA series [1]. However, the escalating interest in wind energy has prompted research institutions worldwide, such as the National Renewable Energy Laboratory (NREL) in the United States [2], the Delft Wind Energy Institute (DU-WIND) [3], the Risø DTU National Laboratory in Denmark [4], and the Energy Centre of the Netherlands (ECN) in the Netherlands [5], to engineer an array of specialized airfoils tailored to diverse wind turbine blade operating conditions. These airfoils, along with their refined iterations, have become the cornerstone for the design of modern industrial wind turbines.

Optimizing airfoil shapes serves multiple objectives, with the most critical one for wind turbines being the improvement of energy extraction from the wind [6,7,8]. However, it is essential to recognize that optimization has its inherent constraints, including:

(a)

Biasing toward the initial values, or the base airfoil shape [9].

(b)

Constraints imposed by the optimization algorithm, which encompass limitations associated with the dimensions of the search space [10], the constraints and unpredictability of mutations [11], and the probability for convergence at a locally optimal point [12].

(c)

Inherent limitations in geometry parameterization methods, which are employed to reduce the number of variables and expedite the optimization process [13].

Geometry parameterization methods play a crucial role in the optimization process by reducing the number of variables involved. To achieve this, mathematical functions, often of degree 6 or 7, are employed to represent the intricate airfoil shape. This results in a significant reduction in optimization variables, typically down to seven coefficients for a 6th-degree polynomial, from the multitude of points required to define an airfoil, which can exceed 60. Nonetheless, it is important to acknowledge that using polynomials to precisely represent airfoil shapes comes with five inherent limitations: completeness, flawlessness, orthogonality, parsimony, and intuitiveness [14]. These limitations manifest in practical applications. For instance, attempts to parameterize the cross-section of a seagull wing airfoil using the CST method, a well-established technique in this field, have proven to be challenging, yielding errors that exceed expectations. Similarly, alternative methods for parameterizing geometry also come with their own set of constraints. Thus, the limitations inherent in these approaches underscore the value of nature as a rich source of inspiration for the ongoing development of highly efficient airfoils.

In computational fluid dynamics (CFD) simulations, particularly concerning airfoils operating under stall conditions, the primary challenges revolve around accurately predicting the separation point and determining the size of the separated bubble. Given the significant cost associated with direct numerical simulation (DNS) and large eddy simulation (LES) methods [15], most engineering applications resort to utilizing Reynolds-averaged Navier–Stokes (RANS) equations, often coupled with Boussinesq’s hypothesis. It is worth noting that the predictions of RANS turbulence models exhibit variations in terms of the onset of flow separation [16]. These differing approaches have a substantial impact on the calculated aerodynamic coefficients, including drag and lift forces [17]. Additionally, there are notable disparities in modeling free shear flows [18]. While these distinctions may not significantly affect the modeling of simple flows, their consequences become apparent in complex flow scenarios [17]. Furthermore, certain turbulence models, such as SST-k-ω, demonstrate sensitivity to specific parameters, such as near-wall spacing [19].

RANS turbulence models underwent evaluation through research conducted by the AIAA Fluid Dynamics Technical Committee (FDTC) in a coordinated program. This study investigated the efficacy of nine models in predicting the aerodynamic coefficients of the NACA-0012 airfoil under stall conditions, which included models like SA, SST, and LRR/SSG. Notably, all models exhibited close predictions of aerodynamic coefficients at low angles of attack. The combined models, like the SST-k-ω and LRR/SSG, showed a slight advantage in lift coefficient prediction. However, when the angles of attack exceeded 17 degrees, all models faced limitations and failed to provide satisfactory results [20].

To address these challenges, Menter et al. [21] developed a novel model based on the k-ω formulation, referred to as the generalized k-ω (GEKO) turbulence model. They demonstrated that, through the calibration of specific coefficients referred to as “free coefficients,” GEKO offers a flexible and robust approach for turbulence modeling in a broad range of problems, particularly excelling at high angles of attack [21].

This paper aims to investigate enhancements in power generation using a bioinspired airfoil for a HAWT. To achieve this objective, the seagull’s wing airfoil has been substituted for the S809 airfoil in the NREL Phase VI HAWT, and the power generation of these turbines is compared. The information presented in the remainder of this article is structured as follows: Section 2 provides details about the GEKO turbulence model and its coefficients. It also includes information related to the calibration of GEKO for both the S809 and the bioinspired (seagull) airfoil, utilizing CFD simulations and experimental validation of the results. Section 3 describes the 3D CFD domain and the mesh arrangement used for CFD simulations of both the original 3D NREL Phase VI turbine and a new turbine constructed using the bioinspired airfoil. Section 4 is dedicated to the results of 3D CFD simulations and discussions about how the seagull airfoil enhances the power generation of this HAWT. Finally, conclusions and references are presented in Section 5.

It is important to note that the S1223 airfoil is inherently thin and may not be the default choice for multi-megawatt turbines. However, in the context of the present research, which primarily focuses on enhancing the aerodynamic performance of the wind turbine, the strategy employed here is rational and well-justified.

2. Generalized K-Omega (GEKO)

In a recent study, Menter et al. [21] introduced a revolutionary turbulence model known as the generalized k-ω model, abbreviated as GEKO. This model incorporates six parameters, commonly referred to as “free parameters.” These free parameters play a pivotal role in shaping the behavior of the model to match the flow characteristics in various regions, such as turbulent boundary layers, shear flows, and flow in corners, among others. The inclusion of these calibration coefficients, or free parameters, renders the GEKO turbulence model an adaptable and robust approach for RANS-based turbulence modeling, suitable for a broad spectrum of problems [21].

In the shear stress transfer model, turbulence kinetic energy (k) and specific dissipation rate (ω) are obtained from the following transport equations:

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_j}\left(\rho k{u}_j\right)=\frac{\partial }{\partial x_j}\left((\mu +\frac{1}{{\sigma }_k}{\mu }_t)\frac{\partial k}{\partial x_j}\right)-{\beta }^{\mathrm{*}}\rho k\omega +P_k$$

(1)

$$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_j}\left(\rho \omega u_j\right)=\frac{\partial }{\partial x_j}\left((\mu +\frac{1}{{\sigma }_{\omega }}{\mu }_t)\frac{\partial \omega }{\partial x_j}\right)+C_{\omega 1}{F}_1\frac{\omega }{k}{P}_k-C_{\omega 2}{F}_2\rho {\omega }^2+\rho F_3\frac{2}{{\sigma }_{\omega }}\frac{1}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}$$

(2)

where β*, ${\sigma }_k$, and ${\sigma }_{\omega }$ represent empirical coefficients. Additionally, $u_j$denotes velocity components, ρ stands for density, t represents time, μ represents molecular viscosity, and S represents the added source term. The turbulent eddy viscosity (${\mu }_t$) is defined as:

$${\mu }_t=\rho v_t=\frac{\rho k}{max\left(\omega ,S/C_{Realized}\right)}$$

(3)

where S is the invariant measure of the strain rate, and $C_{Realized}$is an empirical coefficient.

The limiting factor for turbulence production in stagnated areas concerning turbulence kinetic energy (Equation (1)) is defined as:

$$P_k=-{\tau }_{ij}\frac{\partial U_i}{\partial x_j}$$

(4)

and

$$\begin{gathered} {\tau }_{ij}={\tau }_{ij}^{EV}-C_{Corner}\frac{12{\mu }_t}{\mathit{max}\left(0.3\omega \sqrt{0.5\left(S^2-{\Omega }^2\right)}\right)} \\ \times (S_{ik}{\Omega }_{kj}-{\Omega }_{ik}{S}_{kj}) \end{gathered}$$

(5)

where $C_{Corner}$ is another empirical coefficient. The term ${\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}^{\mathrm{E}\mathrm{V}}$can be calculated via:

$${\tau }_{ij}^{EV}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_{t}2S_{ij}-\frac{2}{3}\rho k{\delta }_{ij}$$

(6)

In the above equations, other variables are defined as:

$$S_{ij}=0.5\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_j}\right), {\Omega }_{ij}=0.5\left(\frac{\partial U_i}{\partial x_j}-\frac{\partial U_j}{\partial x_{ji}}\right), S=\sqrt{2S_{ij}{S}_{ij}} , \Omega =\sqrt{2{\Omega }_{ij}{\Omega }_{ij}}$$

(7)

The six free coefficients of the GEKO model are implemented through functions ($F_1, F_2,F_3$). The default values used for these free coefficients are presented in

Table 1. Default values of free coefficients of the GEKO model.

${\mathit{C}}_{\mathit{C}\mathit{u}\mathit{r}\mathit{v}}$	${\mathit{C}}_{\mathit{C}\mathit{o}\mathit{r}\mathit{n}\mathit{e}\mathit{r}}$	${\mathit{C}}_{\mathit{J}\mathit{e}\mathit{t}}$	${\mathit{C}}_{\mathit{M}\mathit{i}\mathit{x}}$	${\mathit{C}}_{\mathit{N}\mathit{W}}$	${\mathit{C}}_{\mathit{S}\mathit{e}\mathit{p}}$
1.0	1.0	0.9	1.0	0.5	1.75

.

It is worth noting that the separation parameter ($C_{Sep}$) stands out as the most critical among the free parameters that significantly impact the behavior of the GEKO model, particularly when dealing with an airfoil in stall conditions. To ensure the utmost precision in our current research, we not only calibrated this pivotal parameter but also compared the accuracy of the GEKO method in simulating deep stall conditions with some other established turbulence models. Subsequently, we harnessed the GEKO method, incorporating the calibrated ($C_{Sep}$) coefficient, to conduct flow calculations around the NREL Phase VI HAWT.

3. Case Studies Using S809 and S1223 Airfoils

In this research, the simulations of the NREL Phase VI turbine are conducted using both S809 and S1223 airfoils. The S809 airfoil serves as the primary airfoil for the NREL Phase VI HAWT. In contrast, the second airfoil, S1223, bears a close resemblance to a seagull’s wing cross-section and hence is employed for calibration of the seagull’s wing cross-section. The latter is then used throughout the 3D investigation in this study to illustrate its potential for enhancing the power generation of the NREL Phase VI turbine.

The S809 airfoil is notably characterized by its relatively thick profile and has been proven to exhibit remarkable performance at high angles of attack (AOAs), as well as in the near-stall and post-stall regimes. Furthermore, it maintains a commendably low drag coefficient, falling within the range of 0.44 to 0.77. This wide operational range makes it well-suited for power generation across various attack angles [22]. In contrast, the S1223 airfoil stands out as a thin airfoil with significant camber and high-lift characteristics. It has been subjected to testing in wind tunnels at both the University of Illinois [22] and the University of Santa Catarina [23]. The pronounced curvature of its camber reduces pressure on the suction surface while concurrently increasing pressure on the lower surface, leading to the generation of a high lift coefficient. Figure 1 illustrates the geometries of the three airfoils: S809, S1223, and the seagull airfoil.

3.1. GEKO Calibration for S809 and S1223 Airfoils

To calibrate the GEKO coefficients, two-dimensional simulations were conducted to evaluate the turbulence models’ accuracy and estimate the aerodynamic coefficients of airfoils in stall conditions. For the S1223 airfoil, the highest available Reynolds (Re) number for the lift coefficient experimental data was obtained from Selig and Guglielmo [22] at Re = 200,000. However, for the S809 airfoil, higher Re numbers were accessible. In the current study, Re = 2,000,000 was employed, which is in the range of the physical Reynolds number related to the NREL Phase VI blade in operation.

To facilitate this, a C-type mesh was employed in the CFD domain surrounding the airfoils. The domain’s height was maintained uniformly at 20C, with the inlet positioned at a distance of 10C from the airfoil and a downstream domain extending 30C. To enhance cell resolution in regions proximate to the airfoils, the CFD domains were subdivided into three distinct subregions: inner, middle, and outer.

A mesh sensitivity study was conducted using 2000 points along the chordwise direction of the airfoils. Meshes comprising 155,000 and 153,000 cells were found to exhibit mesh independence for the S809 and S1223 airfoils, respectively. The distribution of points along the chordwise direction was chosen such that mesh density at the leading and trailing edges was 20 times greater than at the airfoil’s midpoint. A 30-layer mesh was implemented for simulating the flow within the boundary layer, with the first layer’s height set at 5 × 10⁻⁶ (m) and an expansion ratio of 1.2. Cell sizes were adjusted to 0.008 C in the near area, 0.04 C in the middle area, and 0.16 C in the far area. Figure 2 provides a comprehensive view of the mesh arrangement near the airfoil and within the domain.

The inlet featured a constant velocity, while the outlet maintained a constant pressure equivalent to 1 atm. The top and bottom sides, along with the two normal-to-surface planes, were assumed to be symmetric. Spatial discretization methods employed the least squares cell-based method for gradients and the second-order method for pressure equations. For spatial discretization of momentum, intermittency, turbulence kinetic energy (k), and its specific dissipation rate (ω), a third-order MUSCLE approach was chosen.

The simulations were performed using steady-state assumptions. To expedite convergence, high-order term relaxation methods were employed. Convergence was deemed achieved when the results showed minimal fluctuations, with the absolute convergence criterion set to 10⁻⁵ for the residuals.

The GEKO turbulence model was calibrated for the separation parameter ($C_{Sep}$) for application to the S809 and S1223 airfoils. Given the significance of model accuracy at high angles of attack (AOAs), calibration was performed at AOA = 20°, where both airfoils experience substantial separation, and the experimental data are available [22,24,25]. As detailed in

Table 2. Calibration of the separation parameter ($C_{Sep}$) of the GEKO for S809 and comparison with experimental data.

Turbulence Model		Lift Coef.	Error (%)
1	Experiment [24]	0.886	0
2	GEKO, $C_{Sep}$ = 1	1.736	95.95
3	GEKO, $C_{Sep}$ = 1.75	1.283	44.78
4	GEKO, $C_{Sep}$ = 2.5	0.953	7.60
5	GEKO, $C_{Sep}$ = 3.25	0.771	−12.96
6	k-ω SST	1.418	60.06
7	Spalart–Almaras (2D)	1.259	42.07

and

Table 3. Calibration of the separation parameter ($C_{Sep}$) of the GEKO for S1223 and comparison with experimental data.

Turbulence Model		Lift Coef.	Error (%)
1	Experiment [22]	2.04	0
2	GEKO, $C_{Sep}$ = 0.4	2.881	41.23
3	GEKO, $C_{Sep}$ = 0.9	2.081	2.01
4	GEKO, $C_{Sep}$ = 1	1.733	−15.05
5	GEKO, $C_{Sep}$ = 1.75	1.351	−33.77
6	GEKO, $C_{Sep}$ = 2.5	1.339	−34.36
7	k-ω SST	1.413	−30.74
8	Spalart–Almaras (2D)	1.819	−10.83

, the GEKO turbulence model exhibited exceptional precision, with a separation factor of $C_{Sep}$ = 2.5 for the S809 airfoil and $C_{Sep}$ = 0.9 for the S1223 airfoil. It is worth noting that the performance of GEKO surpasses that of the SST-k-ω and Spalart–Almaras (SA) turbulence models significantly.

In Figure 3, a comparison is presented between the experimental data [22,24] and the computed values for the lift, drag, and lift-to-drag coefficients of the S809 and S1223 airfoils. The computed results align well with the experimental data, indicating a high level of agreement between the two.

Figure 4 compares the vortices generated by the GEKO model, employing the previously calibrated separation coefficients, and other turbulence models around the S809 and S1223 airfoils at AOA = 20°. These visual representations vividly showcase the distinctions in flow separation predictions across diverse turbulence models, offering valuable insights into the disparities in the lift and drag values predicted by each turbulence model.

The Q-criterion is employed in this research to visualize vortices. It identifies eddies as regions where the flow’s rotation component exceeds the strain [26] The Q-criterion is calculated as follows:

$$Q=\frac{1}{2}({\Omega }_{ij}{\Omega }_{ij}-S_{ij}{S}_{ij})$$

(8)

where ${\Omega }_{ij}$ is the vorticity tensor and ${\mathrm{S}}_{ij}$ is the strain rate tensor, which is defined as:

$${\Omega }_{ij}=\frac{1}{2}(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i})$$

(9)

$$S_{ij}=\frac{1}{2}(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})$$

(10)

Figure 5 elucidates the vortices generated by diverse turbulence models surrounding the S809 and S1223 airfoils. These vortices are visualized using the Q-criterion and color-coded based on velocity values. Notably, within both airfoils, we observe the presence of a discernible boundary layer characterized by blue to green regions, while the external flow exhibits hues ranging from green to red. The majority of methods have forecasted a heightened turbulence level within the boundary layer of the S1223 airfoil in comparison to the S809 airfoil. This contrast can be attributed to the more pronounced curvature along the centerline (camber line) of the S1223 airfoil, setting it apart from the S809 airfoil.

A conspicuous observation in Figure 5 is the severe reduction in vortex intensity and boundary layer thickness surrounding the S1223 airfoil, particularly at high angles of attack (AOA = 15° and 20°). This phenomenon implies an enhanced pressure recovery on the upper surface of the S1223 airfoil and a decrease in drag attributed to flow separation at elevated attack angles. The configuration of the boundary layer has a direct impact on velocity and pressure profiles, consequently influencing the airfoil’s aerodynamic coefficients. A thinner boundary layer facilitates a more favorable alignment of flow streamlines with the airfoil’s curvature, ultimately augmenting lift coefficients. Furthermore, the curved configuration of airfoils such as the S1223 provides an ideal surface for aerodynamic loading, further contributing to increased forces.

In

Table 4. Increment ratio of lift and lift-to-drag coefficients for S1223 compared to S809 at different AOAs.

AOA (α=)	Lift Increment Ratio (%)	Lift-to-Drag Increment Ratio (%)
−14°	108%	103%
−5°	129%	103%
0°	675%	106%
5°	117%	−47%
10°	97%	1%
15°	109%	75%
20°	118%	100

one can find a tabulation of lift ratios and lift-to-drag coefficients for the S1223 airfoil to the S809 airfoil across different angles of attack (AOAs).

Table 5. Increment ratio of lift and lift-to-drag coefficients for seagull airfoil compared to S1223 at different AOAs.

AOA (α=)	Lift Increment Ratio (%)	Lift-to-Drag Increment Ratio (%)
−14°	50.5	77.6
−5°	127.4	369.7
0°	−0.6	142.6
5°	0.9	150.0
10°	−4.5	119.2
15°	−4.0	68.9
20°	21.7	163.1

presents a similar comparison of the seagull airfoil with respect to the S1223 airfoil.

As demonstrated in Figure 1, the S1223 and seagull airfoils are geometrically close to each other. This can also be seen in the aerodynamic results shown in Figure 6.

3.2. Validation of Simulation and Power Improvement 3D Wind Turbine

Domain Size, Mesh Density, and CFD Setting

Following the establishment of the seagull airfoil’s superiority over the S809, the subsequent phase entails the replacement of the S809 airfoil in the NREL Phase VI with the seagull airfoil, followed by a comprehensive comparison of the outcomes for the two wind turbines. The NREL Phase VI is characterized by a two-bladed horizontal-axis turbine, with the rotor positioned upstream of the tower. The turbine’s blade radius measures 5.05 m, and specific details regarding the local twist angles and chord lengths of the blades can be found in

Table 6. Local chord and twist angle of the NREL Phase VI at different radial locations.

Radius (m)	Nondimensional Radius (r/R)	Local Chord Length (m)	Local Twist Angle (deg.)
0	0.00	0.218	0
0.508	0.10	0.218	0
0.66	0.13	0.218	0
0.883	0.18	0.183	0
1.008	0.20	0.349	6.7
1.067	0.21	0.441	9.9
1.133	0.23	0.544	13.4
1.257	0.25	0.737	20.04
1.343	0.27	0.728	18.074
1.51	0.30	0.711	14.292
1.952	0.39	0.666	7.979
2.257	0.45	0.636	5.308
2.343	0.47	0.627	4.715
2.562	0.51	0.605	3.425
2.867	0.57	0.574	2.083
3.172	0.63	0.543	1.15
3.185	0.63	0.542	1.115
3.476	0.69	0.512	0.494
4.023	0.80	0.457	−0.381
4.086	0.81	0.451	−0.475
4.696	0.93	0.389	−1.352
4.78	0.95	0.381	−1.469
5.029	1.00	0.355	−1.815

[27].

The blades of both HAWTs were simulated within CFD domains of similar dimensions. In order to mitigate computational expenses, only one half of each turbine, effectively representing a single blade, underwent simulation (Figure 7). As it can be seen from this figure, the CFD model encompasses two distinct sub-domains: an inner subdomain encircling the blade and an outer subdomain. The dimensions of the domain were set based on a sensitivity analysis, not presented here, and the final length of the CFD domain was set to 70 m (=14R), with inner and outer subdomain radii of 7.5 m (=1.5R) and 15 m (=3R), respectively.

Regarding the boundary condition, a periodic boundary condition was imposed on the model. The inlet was defined by a constant velocity ranging from 7 to 25 m/s, while the outlet was treated as a constant pressure outlet (P = 1 atm).

In general, a higher mesh resolution was applied to the inner subdomain, while the outer subdomain featured a relatively lower mesh density. During the investigation for mesh independence, various mesh parameters were explored until it was determined that further increasing mesh density would not substantially influence the power output of the turbines. Specifically, cell counts along the airfoil’s perimeter were tested, ranging from 60 to 90 cells. It was observed that augmenting the number of cells beyond 80 did not result in any significant impact on the power output. To further refine the mesh resolution near the leading and trailing edges, a biasing technique with a factor of 20 was implemented.

In the final mesh arrangement, the first layer of the boundary layer was positioned at a distance of ${5\times 10}^{-6}$ m from the blade surface. Subsequently, 15 layers were generated with an expansion ratio of 1.6. With these settings, the total number of elements in the mesh was approximately ${11\times 10}^6$.

To ensure that the mesh satisfied the requisite near-wall spacing, the Y+ value was assessed on the blade surface at a nominal speed of 12 m/s. As illustrated in Figure 8, the maximum Y+ value on the blade surface was 0.5, which is below the stipulated threshold of one. Furthermore, another criterion involved an examination of the vortices generated around the turbine blade. This was conducted using the Q-criterion contour with a consistent value of 100,000 [s⁻²)] based on velocity values, as visualized in Figure 9. The presence of these vortices signifies the high quality of grid cells and adherence to the selected convergence criteria.

4. Results and Discussions

Initially, to validate the accuracy of the current CFD results, a comparison is made between the output power of the NREL turbine at various wind speeds obtained from this simulation and the results from other researchers [28,29], as well as the NREL experiment [27], as illustrated in Figure 10. The CFD results stemming from the current simulation, utilizing the S809 airfoil and the calibrated GEKO model, exhibit excellent agreement with the findings of other researchers and the experimental data (indicated by dashed lines). This alignment underscores the precision of the ongoing CFD simulation.

Subsequently, to explore the impact of the seagull airfoil and juxtapose it with the S809, the power generated by the two blades is compared in Figure 10 and

Table 7. The power output of NREL phase VI HAWT.

Wind Speed (m/s)	Experiment [28]	Present CFD (S809 Airfoil)	Present CFD (Seagull Airfoil)	Power Improvement (%)
7	6.01	5.92	8.72	47.23
10	10.10	9.88	17.77	79.93
13	9.80	7.54	22.94	204.38
15	8.78	8.56	19.87	132.19
20	8.35	9.30	21.14	127.30
25	11.07	11.66	23.78	103.93

. It is crucial to note that the separation coefficient ($C_{Sep}$) in the GEKO model is set at 0.9 when utilizing the seagull wing airfoil, owing to its resemblance to the S1223 airfoil. In contrast, when employing the S809 airfoil, $C_{Sep}$ is set to 2.5. The results reveal that, in comparison to the original NREL Phase VI, the blade featuring the seagull airfoil demonstrates a substantial increase in the turbine’s output power. The power generated by the seagull airfoil-equipped blade exhibits a maximum enhancement of approximately 204% compared to the original NREL Phase VI, aligning with the superior performance of the seagull airfoil observed in 2D simulations, as detailed in the previous section.

The power coefficients ($Cp$) of these two turbines (implanting the S809 and seagull airfoils) were computed using:

$$Cp=\frac{power }{0.5\rho A{v}^3}$$

(11)

where $\rho$ denotes the air density, $A$ represents the area of the turbine’s rotor disk, and $v$ indicates the wind speed. The maximum power that can be extracted from the wind by a turbine is known as the Betz limit, which has been determined to be 59.3% [30]. While the NREL Phase VI reached 39%, the same turbine achieved a maximum $Cp$ of 54% when employing the seagull airfoil (see Figure 11).

To delve deeper into the aforementioned discussion, a comparison of the pressure contours on the suction surface and the pressure surface of the S809 and seagull airfoils is provided. Figure 12, Figure 13, Figure 14 and Figure 15 visually portray the transformation in pressure distribution brought about by the alteration of the airfoil. This change results in an increase in pressure on the pressure surface and a decrease in pressure on the suction surface. These shifts in pressure distribution significantly contribute to the enhancement of the turbine’s aerodynamic power.

5. Concluding Remarks

In this study, the potential of bioinspired airfoils in enhancing the performance of wind turbines, particularly in the context of the NREL Phase VI wind turbine, was explored. To achieve this, the conventional S809 airfoil was substituted with a bioinspired airfoil modeled after a seagull’s wing. Remarkably, the shape of this seagull airfoil bore a striking resemblance to the S1223 airfoil.

The investigation commenced with a rigorous examination of the precision of using the GEKO turbulence model for the numerical modeling of flow over the S809 and S1223 airfoils, particularly in deep stagnation flow scenarios. For that, a comprehensive comparative analysis was conducted, pitting four turbulence models: Spalart–Almaras, SST-k-ω, IDDES, and the GEKO turbulence model—against each other. The results unveiled the superior performance of the GEKO method, specifically using a separation factor of $C_{Sep}$ = 2.5 for the S809 airfoil and $C_{Sep}$ = 0.9 for the S1223 airfoil.

Building upon these promising findings, a bold step was taken by replacing the S809 airfoil with the seagull airfoil in the NREL Phase VI wind turbine blade. The 3D fluid dynamics simulations, executed with the GEKO turbulence model, yielded results that underscored the transformative potential of this modification. A substantial increase in power generation was observed, with a 2.8 kW (47.2%) boost at a wind speed of 7 m/s, an astonishing 15.4 kW (204.4%) surge at a wind speed of 13 m/s, and an impressive 11.12 kW (103.9%) rise at the maximum operational speed of 25 m/s.

This research serves as a testament to the significant impact bioinspired airfoils can have on improving energy extraction in the wind energy industry. Furthermore, it underscores the reliability and value of the GEKO turbulence model for aerodynamic simulations, particularly for scenarios up to an AOA of 20°.