Optimization of a Small Wind Turbine for a Rural Area: A Case Study of Deniliquin, New South Wales, Australia

The performance of a wind turbine is affected by wind conditions and blade shape. This study aimed to optimize the performance of a 20 kW horizontal-axis wind turbine (HAWT) under local wind conditions at Deniliquin, New South Wales, Australia. Ansys Fluent (version 18.2, Canonsburg, PA, USA) was used to investigate the aerodynamic performance of the HAWT. The effects of four Reynolds-averaged Navier–Stokes turbulence models on predicting the flows under separation condition were examined. The transition SST model had the best agreement with the NREL CER data. Then, the aerodynamic shape of the rotor was optimized to maximize the annual energy production (AEP) in the Deniliquin region. Statistical wind analysis was applied to define the Weibull function and scale parameters which were 2.096 and 5.042 m/s, respectively. The HARP_Opt (National Renewable Energy Laboratory, Golden, CO, USA) was enhanced with design variables concerning the shape of the blade, rated rotational speed, and pitch angle. The pitch angle remained at 0° while the rising wind speed improved rotor speed to 148.4482 rpm at rated speed. This optimization improved the AEP rate by 9.068% when compared to the original NREL design.


Introduction
Energy demands are increasing worldwide and exponentially, given the rising population and needs for economic growth [1]. Consumption of energy is expected to increase by 56% from 553 quadrillion kJ to 855 quadrillion kJ for the period 2010 to 2040 [2]. The extensive consumption of fossil fuels is the primary source of CO 2 emissions into the atmosphere, which is estimated to increase from 31 to 36 billion metric tons during 2010-2020 and may reach 45 billion metric tons by 2040 [3]. The demand for what is termed "clean energy" has increased enormously in recent years due to the fact of people's environmental awareness, desire for energy security, and governments enacting increasingly strict environmental policies [4]. Of all the renewable energy sources, wind energy seems to be favored due to the fact of its low price and rapid global development [5]. The global power installation from wind energy rose from 296,581 MW in 2013 to 539,291 MW in 2017, and it is predicted to reach 817 GW by 2021 [3].
Location has a significant effect on the power output of wind turbines. Accounting for environmental conditions when designing the HWAT for a specific area or region could improve the power output. Many researchers have optimized the rotor shape of wind turbines to maximize the annual energy output. Optimization of the wind turbine blade shape was undertaken at Gökçeada [6] HWAT for a specific area or region could improve the power output. Research has focused on the aerodynamic optimization of the shape of the wind turbine; this is critical in the manufacturing and design of the wind turbine. Australia has widespread and plentifully distributed wind resources. There is a lack of research about the design of wind turbines that takes into account the prevailing environmental conditions in Australia. This study aims to optimize a 20 kW wind turbine for the rural region at Deniliquin, New South Wales, using the horizontal axis rotor performance optimization (HARP_Opt) code. This paper optimized the wind turbine shape to maximize AEP, depending on the Weibull distribution function. The present study is structured as follows. Sections 2 and 3 present the numerical modelling and optimization methods, respectively. Section 4 discusses the simulation results, and Section 5 summarizes and concludes the critical findings of this study.

Governing Equations of Selected Turbulence Models
The principle of the mathematical concept of the RANS equations is based on the calculation method of the Navier-Stokes equation which is divided into the instantaneous fluctuating part and the average flow part. In the present work, the aerodynamics of the wind turbine is predicted with a commercial CFD code known as Ansys Fluent 18.2 [20]. The flow around a wind turbine blade is considered to be incompressible and is modelled utilizing the RANS method. The software uses the finite volume method for solving the mass and momentum equations in addition to equations of turbulence for each control volume cell. The mass and momentum conservation equations are: whereú i is the fluctuating velocity, u i the mean velocity, and v the fluid kinematic viscosity.
The realizable k-ε [21], k-ω SST [22], Spalart-Allmaras [23], and transition SST [24] models are the four RANS models which were investigated. The governing equations of the four RANS models are described in Appendix A.
This section of the paper aims to investigate the effect of those four turbulence models on predicting the aerodynamic characteristics of the twisted wind turbine where the mechanical torque and blade pressure distribution are used for model validation when compared with the NREL test results. Secondly, the differences among the turbulence models under different wind speeds that include stall conditions are documented through simulation of the wind turbine.

Computational Domain
The NREL (CER) extends the NREL (VI) experiment rotor [25] which was performed in the National Renewable Energy Laboratory. The geometry of the 20 kW wind turbine blade was optimized to maximize the annual energy [26]. The results explained the prediction complexity of the aerodynamic performance of the tapered-twisted HAWT turbine when compared with an untwisted blade. As such, the NREL CER performs excellently compared to the commercial blades. The NREL CER was used as a reference for validating the aerodynamic performance of a three-blade wind turbine with variable speed operations. Table 1 summarizes the main characteristics of the wind turbine blade with S809 airfoil applied from a 25% span at the root to the tip. The blade geometry was created using SolidWorks [27] and disregarded the effect of the tower and nacelle to reduce computational time and enhance numerical stability [28]. The wind blade consisted of 18 sections as shown in Figure 1a  The blade geometry was created using SolidWorks [27] and disregarded the effect of the tower and nacelle to reduce computational time and enhance numerical stability [28]. The wind blade consisted of 18 sections as shown in Figure 1a. These sections had different twist angles and chord lengths along the blade as shown in Figure 1b. The geometry of the cylindrical part ended at 0.66 m from the axis of rotation and then started to change until the transition area ended at 1.25 m. From Figure 1b, the maximum twist angle reached 20.040° at 0.25 span and later became zero at a 0.75 span and a negative value at a 1.00 span. The fluid domain was divided into two zones as shown in Figure 2a. The first zone with a 120° radial stream tube was generated with periodic faces to decrease computational time due to the symmetrical flow around the wind turbine model. The upstream velocity was specified The fluid domain was divided into two zones as shown in Figure 2a. The first zone with a 120 • radial stream tube was generated with periodic faces to decrease computational time due to the symmetrical flow around the wind turbine model. The upstream velocity was specified with an 18 m Energies 2020, 13, 2292 5 of 26 radius, offset 12 m in front of the blade. Meanwhile, the downstream outlet was defined with a 36 m radius, offset 24 m behind the blade, and specified as an atmospheric pressure outlet. The second zone was near the blade with a 10 m radius and 5 m from the center of the root. This zone was created to separate the rotating and stationary zones and to increase the number of mesh cells near the blade. The upstream velocity and upper surface of the domain were specified as free stream wind speed. This conical shape of the domain was used to permit wake conical expansion on the back of the blade. An important step is choosing the most suitable computational domain. First of all, a proper computational domain that permits rotation of the wind turbine blade with a no-slip wall effect is essential when considering the optimization of mesh quality. Consequently, a large computational domain will not permit enough grid generation around the wind turbine. On the other hand, a vast domain would increase the computational time corresponding with an increase in the calculated number of cells.

Computational Mesh Generation
The S809 airfoil has a very sharp trailing edge along the blade which produces nonorthogonal face cells [29]. These cells lead to low-quality mesh and inaccurate or unstable CFD solutions. Generating a sharp trailing edge is not possible in the experimental work, so rounding this sharp edge through a radius of 1 mm will improve the quality of the mesh with an insignificant effect on the CFD results.
The present study used different turbulence models to predict the output power and pressure distribution under a variety of wind speed conditions which was crucial to refining the mesh around the blade rotor. As shown in Figure 2b, the cells were refined gradually away from the blade to reduce the computational time. The ANSYS meshing was used to generate unstructured mesh, and a local and global sizing was used to produce a high-quality grid. The minimum mesh size was 0.008 m 3 around the blade. After that, inflation was used to refine the prismatic cells at the blade surface, and it generated fifteen prismatic layers with a growth rate of 1.2 as shown in Figure 2c.
Also, he proximity and curvature with the fine relevance center were defined as the specification of the local mesh size function which helped to refine the mesh grid. Mesh quality plays a crucial role in the accuracy of CFD results. However, a smaller mesh size requires a longer computation time and more computer memory. For this reason, it is essential to compromise between accuracy and computational time. On the other hand, it is crucial to achieve mesh independence. Nine meshes were tested at a 7.2 m/s wind speed to achieve grid independence by monitoring the mechanical torque. As shown in Figure 3, the CFD results converged when the number of mesh cells was 3,559,082, which had an error rate of 3.82% when compared with the measurement value. Any further increase in the number of mesh cells will significantly raise the computational time but will lead to no improvement in the accuracy. For this reason, 3,559,082 mesh cells is deemed to be the most suitable mesh configuration in terms of computational time and efficiency.

Computational Mesh Generation
The S809 airfoil has a very sharp trailing edge along the blade which produces non-orthogonal face cells [29]. These cells lead to low-quality mesh and inaccurate or unstable CFD solutions. Generating a sharp trailing edge is not possible in the experimental work, so rounding this sharp edge through a radius of 1 mm will improve the quality of the mesh with an insignificant effect on the CFD results.
The present study used different turbulence models to predict the output power and pressure distribution under a variety of wind speed conditions which was crucial to refining the mesh around the blade rotor. As shown in Figure 2b, the cells were refined gradually away from the blade to reduce the computational time. The ANSYS meshing was used to generate unstructured mesh, and a local and global sizing was used to produce a high-quality grid. The minimum mesh size was 0.008 m 3 around the blade. After that, inflation was used to refine the prismatic cells at the blade surface, and it generated fifteen prismatic layers with a growth rate of 1.2 as shown in Figure 2c.
Also, he proximity and curvature with the fine relevance center were defined as the specification of the local mesh size function which helped to refine the mesh grid. Mesh quality plays a crucial role in the accuracy of CFD results. However, a smaller mesh size requires a longer computation time and more computer memory. For this reason, it is essential to compromise between accuracy and computational time. On the other hand, it is crucial to achieve mesh independence. Nine meshes were tested at a 7.2 m/s wind speed to achieve grid independence by monitoring the mechanical torque. As shown in Figure 3, the CFD results converged when the number of mesh cells was 3,559,082, which had an error rate of 3.82% when compared with the measurement value. Any further increase in the number of mesh cells will significantly raise the computational time but will lead to no improvement in the accuracy. For this reason, 3,559,082 mesh cells is deemed to be the most suitable mesh configuration in terms of computational time and efficiency.
Three models are used in ANSYS Fluent for handling the rotational effect, namely, dynamic mesh, sliding mesh, and moving reference frame (MRF) models. The sliding mesh model is appropriate for the transient flow problem but requires a full-scale model. Both dynamic mesh and sliding mesh models require high computational resources. The MRF model is the simplest way for modelling the flow of steady-state rotating objects without using rotating mesh to reduce the computational time [30]. Thus, MRF is applied in the small zone near the blade to apply the rotational speed of the wind turbine with periodic boundary conditions. The transformation from a stationary to a moving frame in terms of relative fluid particle velocities, Coriolis, and centripetal accelerations have been used to calculate the MRF model. The MRF was set-up for computation domain by applying a rotational speed of 72 rpm with the absolute reference frame. The blade was assumed to be a non-slip stationary wall that had a zero-relative velocity with other adjacent cells. Three models are used in ANSYS Fluent for handling the rotational effect, namely, dynamic mesh, sliding mesh, and moving reference frame (MRF) models. The sliding mesh model is appropriate for the transient flow problem but requires a full-scale model. Both dynamic mesh and sliding mesh models require high computational resources. The MRF model is the simplest way for modelling the flow of steady-state rotating objects without using rotating mesh to reduce the computational time [30]. Thus, MRF is applied in the small zone near the blade to apply the rotational speed of the wind turbine with periodic boundary conditions. The transformation from a stationary to a moving frame in terms of relative fluid particle velocities, Coriolis, and centripetal accelerations have been used to calculate the MRF model. The MRF was set-up for computation domain by applying a rotational speed of 72 rpm with the absolute reference frame. The blade was assumed to be a non-slip stationary wall that had a zero-relative velocity with other adjacent cells.

Numerical Method and Boundary Conditions
Four turbulence models were used to predict wind turbine aerodynamics. All regions in the domain were applied to the boundary conditions. Domain outlet was applied to the pressure boundary conditions of zero gauges pressure with 101,325 Pa of atmospheric conditions. The current investigation focused on 15 inlet velocities ranging from 4.9m/s to 15.2 m/s which was enough to investigate the prediction of different turbulence models for stall delay phenomena. The air density waas approximately constant due to the assumption of incompressible fluid [11]. In this study, the standard air properties were used, where air density and dynamic viscosity were 1.225 kg/s and 1.7894 × 10 −5 kg/ms −1 , respectively.
The steady-state, pressure-based method was used to solve the incompressible RANS models. A semi-implicit method was used to model the velocity and pressure in momentum and continuity equations. The convergence rate is improved when using a combined algorithm rather than a simple algorithm.
The first-order upwind was used for solving turbulent kinetic energy and the turbulent dissipation rate, while a second-order upwind scheme was used for solving the momentum equations. The least squares cell-based method was used in a gradient spatial discretization scheme. A standard discretization scheme was used in the pressure values interpolation.
In the simulation process ( Figure 4), it is essential to monitor the convergence of the simulation analysis. In this study, two methods were used to assess the convergence of the fluent analysis. Firstly, the residual values method is a popular method for evaluating the convergence

Numerical Method and Boundary Conditions
Four turbulence models were used to predict wind turbine aerodynamics. All regions in the domain were applied to the boundary conditions. Domain outlet was applied to the pressure boundary conditions of zero gauges pressure with 101,325 Pa of atmospheric conditions. The current investigation focused on 15 inlet velocities ranging from 4.9m/s to 15.2 m/s which was enough to investigate the prediction of different turbulence models for stall delay phenomena. The air density waas approximately constant due to the assumption of incompressible fluid [11]. In this study, the standard air properties were used, where air density and dynamic viscosity were 1.225 kg/s and 1.7894 × 10 −5 kg/ms −1 , respectively.
The steady-state, pressure-based method was used to solve the incompressible RANS models. A semi-implicit method was used to model the velocity and pressure in momentum and continuity equations. The convergence rate is improved when using a combined algorithm rather than a simple algorithm.
The first-order upwind was used for solving turbulent kinetic energy and the turbulent dissipation rate, while a second-order upwind scheme was used for solving the momentum equations. The least squares cell-based method was used in a gradient spatial discretization scheme. A standard discretization scheme was used in the pressure values interpolation.
In the simulation process ( Figure 4), it is essential to monitor the convergence of the simulation analysis. In this study, two methods were used to assess the convergence of the fluent analysis. Firstly, the residual values method is a popular method for evaluating the convergence of the CFD solution. In this study, during the calculation process, different variables of the residual values were monitored such as continuity, x-velocity, y-velocity, z-velocity, specific dissipation rate, and turbulent kinetic energy. The solution was considered to be converged when these residual values were below 10 −5 . Secondly, a net mass imbalance was used to check the convergence of the solution. This method is the difference between the inlet and outlet mass flows. It was considered to be converged when the net mass imbalance was less than 0.001 kg/s [31].
Due to the non-linear nature of the fluid flow, the solution should be calculated iteratively. In this study, the solution was achieved after the 1500 iterations. Also, the study used the standard initialization method, where the inlet boundary layer was used for calculating the initial values. After the solution was converged, the aerodynamic power output results were validated against experimental data.
Energies 2020, 13, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies [31]. Due to the non-linear nature of the fluid flow, the solution should be calculated iteratively. In this study, the solution was achieved after the 1500 iterations. Also, the study used the standard initialization method, where the inlet boundary layer was used for calculating the initial values. After the solution was converged, the aerodynamic power output results were validated against experimental data. The computing system details are discussed in Table 2a. The corresponding computation time is illustrated in Table 2b. As can be seen in the table, the computation times for the four models exhibit no noticeable differences. Consequently, the selected model will depend on the accuracy rather than computational time. The computing system details are discussed in Table 2a. The corresponding computation time is illustrated in Table 2b. As can be seen in the table, the computation times for the four models exhibit no noticeable differences. Consequently, the selected model will depend on the accuracy rather than computational time.

Optimization of Wind Turbine
Optimizing the wind turbine design's operating parameters has significant impacts on the amount of energy output. It is therefore critical when specifying the shape design of wind turbines, taking into Energies 2020, 13, 2292 9 of 26 consideration the wind data available in the specific location. The following section deals with the modification of the blade geometric parameters to meet wind velocity at the Deniliquin site.

Wind Data Modelling
Australia has significant wind resources, and many of the best places are located in New South Wales. Wind speed plays a vital role in the performance of the wind turbine, since it is the primary source of energy. Wind speeds at a specific site vary according to annual, seasonal, and daily changes. It is crucial to conduct a statistical analysis with a probability distribution to study the potential of wind resources. Deniliquin is located in the Riverina region of New South Wales close to the border with Victoria, Australia. This area functions economically as a service center for the surrounding agricultural area. Thus, in this area, small wind turbines could be used for a remote, small electrical generation and agricultural application project. The wind data of the Deniliquin Airport AWS station (station number 074258) were chosen to analyze the wind energy potentials in this study. This station has a coordinate of −35.56, 144.95, and an elevation of 94.0 m. It is necessary to use some statistical analysis for the wind data. The probability distribution describes the occurrence frequency of wind speed [32]. Observations for the period from June 2018 to July 2019 were studied using Weibull probability density function. The Weibull probability density function, as shown in Equation 2, depends on two factors: scale and shape factor [33,34].
where, k is the shape factor (dimensionless), c is the scale factor (m/s), and v is the wind speed (m/s) [35]. The shape parameter decides the curvature of the probability distribution; any variation in the shape parameter is affected by the estimated wind potential. Determination of the Weibull probability density function requires defining the shape and scale parameters, using different estimation methods. Employing the maximum likelihood method is typically done for defining shape and scale parameters as seen in the following equations [36,37]: where v i is the wind speed at (i) time, and n is the number of readings of wind speed data [38,39].
To define the AEP of the wind turbine, the probability distribution f (v) is combined with the power curve of wind turbine P (v) as shown in the following equation: where v cut,in is the cut-in wind speed (m/s), and v cut,out is the cut-out wind speed (m/s). Assessment of the available resource at Deniliquin was defined by the shape and scale of the Weibull probability density function. The shape and scale factors are used as inputs for the optimization process, of which the objective is to maximize the AEP of a 20 kW wind turbine depending on the wind speed data in Deniliquin. The following section describes the optimization methodology.

Optimization Blade Shape Methodology
National Renewable Energy Labs (NREL) developed HARP_Opt in the USA [40]. The HARP_Opt open-source code is used to conduct the horizontal axis rotor performance optimization process. The HARP_Opt gathers a BEM theory code with a genetic algorithm (GA) code to optimize and design Energies 2020, 13, 2292 10 of 26 the wind turbine's rotor shape. Genetic algorithms [41] are evolutionary algorithms, and they are robust and reliable search techniques depending on the mechanism of natural selection [42]. The optimization process of the GA is done by iterating a set of individual solutions, where a set of solutions is called a population. An iteration is carried out from one population to the next to obtain subsequent populations of superior individuals. The WT_Perf software (National Renewable Energy Laboratory, Golden, CO, USA) [43] is used as the essential BEM code in HARP_Opt. The WT_Perf is developed to analyze the aerodynamic performance of the wind turbine using the BEM code. The HARP_Opt uses the WT_Perf code to predict the performance of the rotor wind turbine and the MATLAB GA code to carry out the optimization. The HARP_Opt could be used for a single or multiple objective optimization code for the objective function in the HARP_Opt software, either maximization of the AEP of the wind turbine or wind turbine efficiency.

Design Variables and Objective Function
The blade shape of the wind turbine is defined by the airfoil, chord length, and twist angle of each section along the blade. The validated wind turbine model of 20 kW, which has been discussed before, is used as a baseline for the optimization process and the same family of S809 airfoil was used. The airfoil lift and drag polar was done using spreadsheet AirfoilPrep v2.02 (v2.02.03, Windward Engineering, LLC, USA) This type of Excel formatting was used to generate airfoil data to be imported into WT_Perf. After preparing the airfoil data file, the file was imported into HARP_Opt.
The objective of the optimization process was to maximize the AEP of 20 kW wind turbine depending on the wind speed data in Deniliquin using the output of optimal chord length and twist distributions along the wind turbine blade. In this study, a single objective method was considered where the optimization was focused on the optimization of the aerodynamic shape without including structural optimization studies [44][45][46][47][48]. The Bezier curves were used to define the chord length and twist angle distributions to smooth the span-wise distribution along the blade. There were five control points for each chord length and twist angle parameters; thus, the total overall decision variables amounted to 25 control points. The same rated power, the rotor diameter, and the hub height were used as input in HARP_Opt, with the baseline validated wind turbine model as shown in Figure 5. In this study, the variable rotor speed and variable pitch control were used for the control system to produce more energy output [32]. The allowable rotor speeds range from 25 rpm to 150 rpm. Table 3 below summarizes different parameters of turbine configurations.

Model Validation
This section firstly investigates the effect of those four turbulence models on predicting the aerodynamic characteristics of the twisted wind turbine, where the mechanical torque and blade pressure distribution were used for model validation when compared with the NREL test results. Secondly, the differences among the turbulence models under different wind speeds that included stall conditions were captured by simulating the wind turbine. Thirdly, the main aerodynamic parameters, such as lift coefficient, were extracted at a different span-wise sections along the blade. Fourthly, and finally, the aerodynamic flow of the S809 airfoil was visualized under different angles of attack, and they revealed the variation of the flow from the attached separated flow conditions. 4.1.1. Mechanical Torque Figure 6 compares the modelled shaft torque values using different RANS models and the measured results for the NREL rotor, which operated at a fixed speed of 72 rpm and a pitch angle of 5°. The CFD results demonstrated good agreement with the measurements at low and medium wind speeds between 4.9 and 9.0 m/s. Transition played an insignificant role in the prediction of the flow behavior of the wind turbine at this speed range. All RANS models made an excellent prediction of the flow around the wind turbine at the area where the flow was still attached.
It is important to note here that at a flow velocity between 9.0 and 10.5 m/s, the boundary layer was driven from the laminar to the turbulent transition area where the onset of a stall occurs. The transition SST model made the best prediction for that region due to the fact of its ability to resolve  Trial and error were used to tune the GA's parameters. The values input into this optimization run were the final ones after using trial and error. Table 5 summarizes the GA configuration employed.

Model Validation
This section firstly investigates the effect of those four turbulence models on predicting the aerodynamic characteristics of the twisted wind turbine, where the mechanical torque and blade pressure distribution were used for model validation when compared with the NREL test results. Secondly, the differences among the turbulence models under different wind speeds that included stall conditions were captured by simulating the wind turbine. Thirdly, the main aerodynamic parameters, such as lift coefficient, were extracted at a different span-wise sections along the blade. Fourthly, and finally, the aerodynamic flow of the S809 airfoil was visualized under different angles of attack, and they revealed the variation of the flow from the attached separated flow conditions. the laminar transition that considers the turbulent boundary layer which starts at the stall phenomenon. After 10 m/s, the results showed a marked difference for mechanical torque between the measured and RANS models. The significant separation of the flow and stall phenomena played a role in the difficulty of mechanical torque prediction using the RANS models.  (7): where is the local static pressure, is the free stream pressure, is relative wind speed, Ω is the rotational wind speed (rad/s), is the radius of the section (m), and is the air density (kg/m 3 ). For a 7.2 m/s inlet wind speed, as shown in Figure 7, the modelled pressure distributions matched well with the experimental values, where the flow was almost attached to this wind speed with no stall boundary layer separation having yet started. Thus, all RANS models agreed well with the experimental data except realizable k-ε. In effect, this model had some limitations when the domain included two fluid zones, i.e., stationary and rotating zones [51]. In this case, non-physical turbulent  It is important to note here that at a flow velocity between 9.0 and 10.5 m/s, the boundary layer was driven from the laminar to the turbulent transition area where the onset of a stall occurs. The transition SST model made the best prediction for that region due to the fact of its ability to resolve the laminar transition that considers the turbulent boundary layer which starts at the stall phenomenon. After 10 m/s, the results showed a marked difference for mechanical torque between the measured and RANS models. The significant separation of the flow and stall phenomena played a role in the difficulty of mechanical torque prediction using the RANS models.  (7): where P is the local static pressure, P rel is the free stream pressure, V rel is relative wind speed, Ω is the rotational wind speed (rad/s), r is the radius of the section (m), and ρ is the air density (kg/m 3 ). For an inlet wind speed of 10.2 m/s, which was classified as the onset of stall and shown in Figure  8, the best way to predict the airfoil's pressure coefficient was to use the transition SST model. Any changes in the adverse pressure gradients would be reflected directly on the laminar boundary layer, and hence the early separation for the laminar boundary layer will happen when compared to the turbulent boundary layer. The transition SST model could predict the boundary layer for the flow region that changed from the laminar flow region to the turbulent transition region. Changing between boundary layers enabled the transition SST model to better predict the aerodynamic flow of the wind turbine for a stall wind condition when compared with other RANS models as shown in Figure 8a. At the 47% span location, the stall phenomena begin, and the separation was evident among the RANS models for predicting the wind turbine's aerodynamics.

Investigation of the Airfoil Characteristics
It is essential to improve our understanding of the main aerodynamic characteristics of each section along the blade. Determination of the lift and drag coefficient of each airfoil section along the blade is a critical parameter for calculating the angle of attack. The Reynolds number increases as the wind speed rises, which affects the angle of attack and corresponding lift coefficient [52]. The results of the lift coefficient were obtained from simulation over NREL at a pitch angle equal to 5° and wind speeds varied from 4.9 to 15.2 m/s at a fixed angular speed, i.e., 72 rpm. As shown in Figure 9a, the lift coefficient varied with both radial sections and velocities.
As shown in Figure 9b, the calculated angle of attack had a high twist angle of 12.13° at 7.2 m/s wind speed on a plane through the blade at a distance of 1.51 m. At the same length, the angle of attack at 10.2 m/s was 22.04°. Thus, the angle of attack increased with increasing wind speed and decreased with the radial position. At 7.2 m/s near the hub region, minor transition and separation may occur due to the slightly high angle of attack, but the flow was still below the stall region and almost always attached. When increasing the wind velocity, the angle of attack increased at 10.2 m/s to the onset of stall. At this wind speed, the full separation occurred between the hub and 80% of the tip sections. Thus, the separation on the flow from the leading edge to the trailing edge occurred at 47% section of the blade and remained until 80% section. For a 7.2 m/s inlet wind speed, as shown in Figure 7, the modelled pressure distributions matched well with the experimental values, where the flow was almost attached to this wind speed with no stall boundary layer separation having yet started. Thus, all RANS models agreed well with the experimental data except realizable k-ε. In effect, this model had some limitations when the domain included two fluid zones, i.e., stationary and rotating zones [51]. In this case, non-physical turbulent viscosities were produced, and these affected the value of the turbulent viscosity. Since the present simulation had two domains, the realizable k-ε was not appropriate for predicting the current wind turbine simulation.
For an inlet wind speed of 10.2 m/s, which was classified as the onset of stall and shown in Figure 8, the best way to predict the airfoil's pressure coefficient was to use the transition SST model. Any changes in the adverse pressure gradients would be reflected directly on the laminar boundary layer, and hence the early separation for the laminar boundary layer will happen when compared to the turbulent boundary layer. The transition SST model could predict the boundary layer for the flow region that changed from the laminar flow region to the turbulent transition region. Changing between boundary layers enabled the transition SST model to better predict the aerodynamic flow of the wind turbine for a stall wind condition when compared with other RANS models as shown in Figure 8a. At the 47% span location, the stall phenomena begin, and the separation was evident among the RANS models for predicting the wind turbine's aerodynamics.

Investigation of the Airfoil Characteristics
It is essential to improve our understanding of the main aerodynamic characteristics of each section along the blade. Determination of the lift and drag coefficient of each airfoil section along the blade is a critical parameter for calculating the angle of attack. The Reynolds number increases as the wind speed rises, which affects the angle of attack and corresponding lift coefficient [52]. The results of the lift coefficient were obtained from simulation over NREL at a pitch angle equal to 5 • and wind speeds varied from 4.9 to 15.2 m/s at a fixed angular speed, i.e., 72 rpm. As shown in Figure 9a, the lift coefficient varied with both radial sections and velocities. The flow behavior at different airfoil sections along the blade was visualized around the airfoil with velocity and pressure contours. Transition SST indicated good results amongst the four turbulence models with the measurements shown above. As such, the velocity vectors and pressure contours results were predicted by the transition SST model. The blade cross-section at a 1.51 m radius investigated the flow on the inner region of the blade. In contrast, the blade cross-section at a 4.78 m radius studied the flow on the outer area of the blade. The figures below explain that the pressure fell to a small value as the flow over the airfoil accelerated. As the results in Figure 10a,c,e,g,i illustrate, the pressure coefficient dropped quickly to zero and reached a negative value. As the flow slowed down, the pressure increased, and the magnitude of the pressure coefficient decreased. As a As shown in Figure 9b, the calculated angle of attack had a high twist angle of 12.13 • at 7.2 m/s wind speed on a plane through the blade at a distance of 1.51 m. At the same length, the angle of attack at 10.2 m/s was 22.04 • . Thus, the angle of attack increased with increasing wind speed and decreased with the radial position. At 7.2 m/s near the hub region, minor transition and separation may occur due to the slightly high angle of attack, but the flow was still below the stall region and almost always attached. When increasing the wind velocity, the angle of attack increased at 10.2 m/s to the onset of stall. At this wind speed, the full separation occurred between the hub and 80% of the tip sections. Thus, the separation on the flow from the leading edge to the trailing edge occurred at 47% section of the blade and remained until 80% section.
The flow behavior at different airfoil sections along the blade was visualized around the airfoil with velocity and pressure contours. Transition SST indicated good results amongst the four turbulence models with the measurements shown above. As such, the velocity vectors and pressure contours results were predicted by the transition SST model. The blade cross-section at a 1.51 m radius investigated the flow on the inner region of the blade. In contrast, the blade cross-section at a 4.78 m radius studied the flow on the outer area of the blade. The figures below explain that the pressure fell to a small value as the flow over the airfoil accelerated. As the results in Figure 10a,c,e,g,i illustrate, the pressure coefficient dropped quickly to zero and reached a negative value. As the flow slowed down, the pressure increased, and the magnitude of the pressure coefficient decreased. As a result, the pressure of the lower surface became far more significant than the pressure of the upper surface, causing the blade to rotate. These results confirmed that the values of the pressure coefficient were negative at high air velocity. The separation will start when the flow on the upper surface at the trailing edge of the airfoil decelerates and mixes with the airflow from the lower surface. The point of flow separation happens earlier at larger angles of attack, and as the intensity of the adverse pressure rises, the separation point shifts forward on the airfoil. of flow separation happens earlier at larger angles of attack, and as the intensity of the adverse pressure rises, the separation point shifts forward on the airfoil. As seen in Figure 10b,d,f,h,j, the velocity vectors on the airfoil sections varied following the differences of velocity along the blade. The velocity vectors increased from the root to the tip where the highest velocity vector was achieved by the outer section of the blade which was affected by the boundary layer separation [53]. The axial velocity increased at a uniform pattern along the blade, explaining the effect of the centrifugal force on the rotating blade. The NREL CER wind turbine will use the original blade geometry design. The numerical modelling of this wind turbine will serve to investigate the mechanical output at different rotational speeds and variable pitch angles. Transition SST is the RANS model that was selected for use when the output power of the original blade geometry was compared with the optimized blade design. Transition SST can predict well the mechanical torque and the pressure distribution of different sections along the blade.

Optimization of Wind Turbine
The objective of this optimization process is to maximize the AEP at the Deniliquin site. It is essential to obtain the Weibull probability density function (see Figure 11a) to maximize the energy output. The Weibull function shape and scale parameters were defined as 2.096 and 5.042 m/s, As seen in Figure 10b,d,f,h,j, the velocity vectors on the airfoil sections varied following the differences of velocity along the blade. The velocity vectors increased from the root to the tip where the highest velocity vector was achieved by the outer section of the blade which was affected by the boundary layer separation [53]. The axial velocity increased at a uniform pattern along the blade, explaining the effect of the centrifugal force on the rotating blade. The NREL CER wind turbine will use the original blade geometry design. The numerical modelling of this wind turbine will serve to investigate the mechanical output at different rotational speeds and variable pitch angles. Transition SST is the RANS model that was selected for use when the output power of the original blade geometry was compared with the optimized blade design. Transition SST can predict well the mechanical torque and the pressure distribution of different sections along the blade.

Optimization of Wind Turbine
The objective of this optimization process is to maximize the AEP at the Deniliquin site. It is essential to obtain the Weibull probability density function (see Figure 11a) to maximize the energy output. The Weibull function shape and scale parameters were defined as 2.096 and 5.042 m/s, respectively. The root mean square error (RMSE) and the coefficient of determination (R 2 ) were used to evaluate the accuracy of the Weibull probability density function. The values of R 2 and RMSE were evaluated as 0.998587 and 0.013567, respectively. Those values reflected the fact that the Weibull probability distribution was very accurate in representing the recorded wind data. Depending on Weibull function parameters for the selected site and the 20kW rated capacity, the rotor diameter chosen seemed good with AEP as shown in Figure 11b.
In this study, the optimization process modified the shape of the blade design using chord and twist distribution along the blade. The chord and twist distributions of the optimized wind turbine blade are shown in Figure 12. The power coefficient of the optimized rotor with wind speed is depicted in Figure 13a.The power coefficient had the highest value of 0.433108 until the rated wind speed (9.5 m/s). Due to the Weibull probability density function being high in the low wind speed range, it is essential to have a high-power coefficient to increase the AEP. The blade pitch control is the variable pitch to stall angle; the control system will be able to brake the turbine at high winds safely. As shown in Figure 13b, the variation of rotor speed and pitch angle with different wind speeds shows the change of rotor speed and pitch angle with varying speeds of wind. The rotor speed at the cut-in speed (2 m/s) was approximately 32.98848 rpm. The pitch remained at 0 • , while the rising wind speed improved the rotor speed to 148.4482 rpm at the rated speed of 9.5 m/s. The active pitch control began at this point and regulated the rotor speed to 150 rpm until 25 m/s was reached.
As seen below in Figure 14, the optimized rotor had a higher power output when compared to the reference rotor. According to Weibull probability density function for the Deniliquin site, 2 m/s speed was a well-suited value for cut-in speed. For the optimized rotor, a cut in the was is producing only 0.201682 kW which was approximately 1% of rated power capacity. As the wind speed rose to 8 m/s, the turbine output was 18.37828 kW with a power coefficient of 0.433108. This output was approximately 76% of the turbine's rated capacity. At 9.5 m/s, the 20 kW rated power was reached. These differences between the power output of the tested and optimized rotor will be reflected in the AEP for the tested and optimized wind turbine. In this study, the AEP increased from 30,819.3 kW-hr/year in the original turbine to 33,614 kW-hr/year in the optimized wind turbine. Optimization improved the AEP by 9.068% when compared to the initial tested wind turbine design. In this study, the optimization process modified the shape of the blade design using chord and twist distribution along the blade. The chord and twist distributions of the optimized wind turbine blade are shown in Figure 12. The power coefficient of the optimized rotor with wind speed is depicted in Figure 13a.The power coefficient had the highest value of 0.433108 until the rated wind speed (9.5 m/s). Due to the Weibull probability density function being high in the low wind speed range, it is essential to have a high-power coefficient to increase the AEP. The blade pitch control is the variable pitch to stall angle; the control system will be able to brake the turbine at high winds safely. As shown in Figure 13b   In this study, the optimization process modified the shape of the blade design using chord and twist distribution along the blade. The chord and twist distributions of the optimized wind turbine blade are shown in Figure 12. The power coefficient of the optimized rotor with wind speed is depicted in Figure 13a.The power coefficient had the highest value of 0.433108 until the rated wind speed (9.5 m/s). Due to the Weibull probability density function being high in the low wind speed range, it is essential to have a high-power coefficient to increase the AEP. The blade pitch control is the variable pitch to stall angle; the control system will be able to brake the turbine at high winds safely. As shown in Figure 13b  wind speed improved the rotor speed to 148.4482 rpm at the rated speed of 9.5 m/s. The active pitch control began at this point and regulated the rotor speed to 150 rpm until 25 m/s was reached.  wind speed improved the rotor speed to 148.4482 rpm at the rated speed of 9.5 m/s. The active pitch control began at this point and regulated the rotor speed to 150 rpm until 25 m/s was reached. As seen below in Figure 14, the optimized rotor had a higher power output when compared to the reference rotor. According to Weibull probability density function for the Deniliquin site, 2 m/s speed was a well-suited value for cut-in speed. For the optimized rotor, a cut in the was is producing only 0.201682 kW which was approximately 1% of rated power capacity. As the wind speed rose to  approximately 76% of the turbine's rated capacity. At 9.5 m/s, the 20 kW rated power was reached. These differences between the power output of the tested and optimized rotor will be reflected in the AEP for the tested and optimized wind turbine. In this study, the AEP increased from 30,819.3 kWhr/year in the original turbine to 33,614 kW-hr/year in the optimized wind turbine. Optimization improved the AEP by 9.068% when compared to the initial tested wind turbine design.

Conclusions
This study aimed to optimize the NREL (CER) wind turbine by CFD modelling. It investigated the aerodynamic characteristics of a 20 kW wind turbine at three-bladed twisted angles and tapered blades. Grid independence was achieved in terms of mechanical torque. As well, the study investigated the effect of four turbulence models in predicting the mechanical torque and pressure distributions. Simulated wind conditions varied from attached to separated flow conditions. Finally, a HARP_Opt code was used to optimize the wind turbine design using a genetic algorithm, aiming to maximize the AEP at the Deniliquin site. The major findings of this study are summarized as follows: (1) All four RANS models agreed well with experimental data at low wind speed ranges. Differences appeared among the four turbulence models as the wind speed increased; (2) At the onset of a stall condition of 10.2 m/s, the ttransition SST reported the best accuracy for predicting the pressure coefficient of the airfoil. The angle of attack increased with increasing wind speed and decreased with the radial position. The full separation occurred between the hub and 80% of the tip sections; (3) The shape of the rotor was modified by changing the chord and twist distribution along the blade, leading to 9.1% improvement in AEP. k-ε for dissipation rate (ε) and turbulent kinetic energy (k). However, the turbulent viscosity generation and calculation method are different in these models. The other RANS model used for wind turbine simulation is k-ω turbulence model, which was initially proposed by [57]. The Imperial College group has made a new improvement to this model, but the most outstanding development was done by Wilcox [58]. In some applications, the k-ω model is more accurate than standard k-ε; for example, in boundary layers with adverse pressure gradient and sublayer able to be integrated without any extra damping functions required [58]. However, k-ω is still sensitive to some flows with free stream boundary applications. k-ω SST Shear Stress Transport (SST) is an advanced turbulent model proposed by [59]. This model combines the advantages of k-ω and k-ε turbulence models. In the inner part of the boundary layer, it uses the k-ω models, and then it is converted gradually to k-ε in the free shear layer and wake regions' outer layer. The other advantages of this model are the modification of eddy viscosity, which considers the effect of turbulent shear stress transportation. The Spalart-Allmaras version is the simplest RANS turbulence model which uses one transport equation and has the advantage of less computational time. The computation of turbulence quantity is formulated by one transport equation, in which the kinematic eddy turbulent viscosity is the equation's variable [23]. This model was designed and optimized for compressible flow over airfoils and wings in aerospace applications and demonstrated good results. It could simply give stable and converge results in different practical situations that include adverse pressure gradients. However, this model could generate enormous diffusion, especially in regions of 3D vortices flow [60]. Various improvements were added to the model, including the effects of rotation, near wall, and reduction of diffusion [61,62]. Thus, the main advantage of the fast convergence of this model is reflected in the low computation time when compared to other turbulence models.
Another RANS model is the transition SST (γ-Re θ ) model which was extended based on the k-ω SST [63]. This model has four transport equations that combine k-ω (SST) equations with the momentum thickness Reynolds number transition outset method (Re θ ) and the intermittency (γ) transport equations. The transition SST model is more precise than classical fully turbulent models due to the fact of its ability to deal with the laminar-turbulent transition flow model where the separation of flow and stall phenomena occurred. The governing equations of the four RANS models are described in Table A1. Table A1. Governing equations of RANS models.

Spalart-Allmaras
Turbulence eddy viscosity (v t ): v t = v f v1 , f v1 = X 3 C v1 3 +X 3 , X ≡ v/v where f v1 is a damping function ranging from zero value at the wall to 1 at far away from the boundary.

Transition SST
Specific design parameters should be in place to generate acceptable blade geometry [28,49,50]. In References [49], the maximum and minimum values for twist angle and chord length at the control points were decreased along the radial direction as the following equation: where X i min is the lower limit and X i max is the upper limit for the chord and the twist angle [65].