Computational Fluid Dynamics Prediction of a Modified Savonius Wind Turbine with Novel Blade Shapes

The Savonius wind turbine is a type of vertical axis wind turbine (VAWTs) that is simply composed of two or three arc-type blades which can generate power even under poor wind conditions. A modified Savonius wind turbine with novel blade shapes is introduced with the aim of increasing the power coefficient of the turbine. The effect of blade fullness, which is a main shape parameter of the blade, on the power production of a two-bladed Savonius wind turbine is investigated using transient computational fluid dynamics (CFD). Simulations are based on the Reynolds Averaged Navier-Stokes (RANS) equations with a renormalization group turbulent model. This numerical method is validated with existing experimental data and then utilized to quantify the performance of design variants. Results quantify the relationship between blade fullness and turbine performance with a blade fullness of 1 resulting in the highest coefficient of power, 0.2573. This power coefficient is 10.98% higher than a conventional Savonius turbine.


Introduction
Wind turbines are generally divided into two categories, horizontal axis wind turbines (HAWTs) and vertical axis wind turbines (VAWTs), based on the relative alignment between their axis of rotation and the wind direction. VAWTs are often favored for small-scale power generation because their performance is not depended on the relative wind direction, allowing generation equipment to be located on the ground shaft [1], resulting in reduced maintenance costs.
The Savonius wind turbine is a common type of VAWT that generates torque through the combined effects of drag and side forces. These turbines are typically composed of two or three arc-type blades. Savonius turbines have the following advantages over other types of wind turbines: (1) ability to operate under complex turbulent flows [2]; (2) low rotation speed and noise emission [3]; (3) high starting torque, good starting performance [4]; and (4) simple structure and low cost.
However, the traditional Savonius turbine has relatively lower power efficiency, with a maximum efficiency of 0.25 [5]. This power coefficient is significantly lower than most lift-type wind turbines, such as HAWTs and Darrieus [6] wind turbines. To increase the efficiency of Savonius turbines comprehensive studies have been conducted to quantify the effects of various design parameters, such as rotor aspect ratio, overlap, number of blades, and endplates, on performance using both experimental and numerical methods [5,[7][8][9][10].
Researchers worked to enhance the Savonius turbine performance by changing the structure of the turbine. Kamoji et al. [11] and Gupta et al. [12] studied the aerodynamic characteristics of a modified Savonius turbine with helical blades, while Golecha et al. [13] and Altan and Atilgan [14,15] placed a guide vane in front of the turbine to deflect flow for the returning blade. McTavish et al. [16] proposed a modified blade shape and carried out both steady and transient CFD simulations. Kamoji et al. [17] and Kacprzak et al. [18] studied the performance of modified turbines with spline-type and Bach-type blades, with an increase in efficiency 16% found in the case of using spline-type blades. Tian et al. [19] carried out CFD simulations of Savonius with elliptical blades. This paper concerns a Savonius wind turbine that has novel blade shapes developed from the Myring Equation [20]. The effect of the blade fullness on the turbine performance is studied over a range of tip speed ratios. The analysis has been performed by solving numerically the incompressible Navier-Stokes equations with the aid of the CFD code Fluent 13.0 [21].

The Modified Savonius Wind Turbine and Parameters Definition
The two-dimensional schematic view and geometrical parameters of a two-bladed Savonius wind turbine with zero gap distance are presented in Figure 1, where U is the wind velocity,  is the rotational velocity of the turbine, c is the chord length of the blades, D is the diameter of the turbine, and  is the azimuth angle of the blade defined as the angle between the flow direction and the line from the blade tip to the rotation center. Unlike traditional Savonius turbine with semi-circular blades, this modified Savonius turbine has its blade shape generated from the Myring Equation [20]. The blade of a traditional Savonius turbine generates drag during the returning half cycle,  between π and 2π, which limits the efficiency of the turbine. Reducing this drag force could enhance the performance of the turbine. The Myring Equations are commonly used in the design of the nose and tail of autonomous underwater vehicles (AUVs) and have been proven to give an ideal coefficient of drag [22]. Therefore the new blade shape could possibly improve the turbine's performance. The Myring Equation for designing the nose of an AUV is: When n = 2 this equation defines a half-ellipse. Additionally, when n = 2 and a = b this equation denotes a half-circle, the shape of traditional Savonius blades. In [19], the authors studied the performance of elliptical blades by holding n = 2 and varying a and b. This paper builds upon this research by focusing on the effect that n has on the performance of the turbine. In the simulations, constant a = b= 0.25 m is used, while n varies from 0.5 to 3. Figure 2 shows the blade shapes for a range of n values. It can be seen that as n increases, the blade becomes steeper at both ends and has a larger integral area along the chord direction. For illustration purpose, n will be defined as the fullness of the blade in the later pages.  The specific values of the turbine geometric parameters for each case are shown in Table 1. The case 5 turbine has the same structural parameters with that tested in the wind tunnel by Blackwell et al. [5], and will be used as the validation case later.

Numerical Method
Though two-dimensional simulations fail to consider the three-dimensional effect, previous studies have shown that two-dimensional simulations give acceptable results for Savonius [1,6,14,16,18,19]. Therefore the three-dimensional effects are ignored and two-dimensional transient simulations are carried out to reduce the time cost in this study. A sliding mesh model was applied to realize the rotation motion of the rotor.

Computation Domains and Boundary Conditions
To allow for the full development of the upstream flow, as well as to decrease the blockage effect, the computational domain was given the dimension of 18D × 12D. The rotor was placed midway between the top and bottom boundaries, at a distance of 6D from the inlet boundary ( Figure 3).
The overall domain is split into three subdomains to facilitate simulating the rotating of the rotor. The highest resolution subdomain contains the grid elements surrounding the rotor (referred to as the rotor domain). The second mid resolution subdomain contains the cells in the near wake region (referred to as the wake domain), and has a dimension of 6D × 2D. The lowest resolution subdomain is the outer domain containing the cells in the outer region. A uniform and steady velocity profile was applied at the inlet of the computation domain. The velocity was set as 7 m/s, the same as that in the experiment and corresponding to a nominal Reynolds number of 4.32 × 10 5 [5]. As the turbulence intensity was not given in [5], a simplification was made that the turbulence intensity was 1% and uniform at the inlet. A pressure outlet boundary was applied at the outlet of the domain. To improve the stability of the numerical simulations, symmetry boundary conditions were applied at both the top and the bottom edges of the domain. The symmetry boundary condition is useful because it allows the solver to consider the wall as part of a larger domain, avoiding the wall effects [6,19,23]. No-slip boundary conditions were imposed at the surface of the blades. Two interfaces were imposed at the overlap edges between the adjacent domains, allowing the transport of the flow properties.

Grid Generation
The computation grid was generated using the MESH tool in ANSYS 13.0. In two-dimensional grid generation, quadrilateral elements are more desirable because they are less memory-occupying than triangular elements (the number of quadrilateral elements is only a half of the triangular elements', considering the same grid nodes number) and provide highly accurate numerical solutions. Therefore, the three subdomains were discretized with quadrilateral elements (Figure 4a). The grid node density utilized in the rotor domain was higher than in the other domains. The grid in the wake domain, where high velocity gradients were expected, was set between the grid resolutions of the rotor and outer domains.
Prism layer grid elements were extruded from the edges of the blade to improve the grid quality and describe with sufficient precision the boundary layer flow. The height of the first prism layer above the surface was set such that the y + value for the first elements from the wall was between 30 and 100, depending on the rotation velocity of the rotor and the position of the elements on the blade (Figure 4b). This y + value made it deal for the use of wall functions. There were a total of eight prism layers, with a layer growth rate of 1.2.

Turbulence Model and Solution Sets
The RNG k-ԑ turbulence model was employed. This model is known to give highly accurate predictions of rotating machinery because the effect of swirl on turbulence is included [21]. Moreover, the RNG k-ԑ model provides an analytically-derived differential formula for effective viscosity, which accounts for low-Reynolds-number effect. Additionally, standard wall functions were utilized to model turbulent effects below the first grid element, increasing the solution accuracy at low Reynolds numbers. For each case listed in Table 1, several simulations were carried out with the tip speed ratio varying from 0.4 to 1.2. Tip speed ratio represents the ratio of the blade tip speed to the wind speed, and has the following expression: A sliding mesh model was used to simulate the rotation of the rotor. Each simulation lasted for five revolutions to allow for convergence. The results of the last revolution were used for the analysis. The time step for each simulation was set so that 1° of rotor rotation was achieved each time step. This time step size was shown to be small enough to lead good results compared to experimental data [19].
Convergence was determined by the order of magnitude of the residuals. The drop of all scaled residuals below 1 × 10 −5 was utilized as the convergence criterion. The maximum number of iterations in a time step was set as 100, which, as we observed, enabled all the residuals to meet convergence in the simulation. The SIMPLE (Semi-Implicit Method for Pressure Linked Equations) pressure-velocity coupling method was used in all simulations. A Second Order Upwind spatial discretization algorithm was used for all the equations, including pressure, momentum, and turbulence, and a Least Squares Cell Based algorithm was used for the gradients. The second order algorithms were shown give more accurate results than first order ones since they reduce interpolation errors and false numerical diffusion [21], and therefore these were utilized.

Verification
A grid convergence study was performed to evaluate the influence of grid resolution on the calculated rotor torque. Simulations were conducted on a conventional Savonius turbine with n = 2 (the case 5 turbine in Table 1) at a tip speed ratio of 1. Figure 5 shows the coefficients of dynamic torque on the turbine for a full rotation of the rotor with three different grid densities, with approximately 60,000, 80,000, and 120,000 elements, respectively. The grid density was controlled by changing the size of elements near the blade and that in the rotor domain and wake domain.  The coefficient of torque is defined as follows: where M is the generated torque and S is cross section area, given by the relationship S = DH with H being the height of the blade. Since only two-dimensional simulations were performed, the unit height H = 1 m was used. It can be seen that grids with approximately 80,000 and 120,000 elements gave approximately the same results. Considering the economy of time in the simulation, the grid with 80,000 elements was chosen for the following simulations.

Validation
Validation studies were conducted using a conventional two-bladed Savonius turbine, i.e., the case 5 turbine. The simulation results were then compared with the experimental data from [5]. Figure 6 shows the comparison between the averaged coefficients of toque as a function of tip speed ratio, with torque averaged over a complete rotor rotation.   Figure 6 suggests that the simulation results are in good agreement with the experimental data, especially at higher tip speed ratios. The maximum relative error occurs at λ = 0.6, where the CFD result is about 5% lower than the experimental data. Therefore, it is acceptable to use the RNG k-ԑ turbulence together with a standard wall treatment to predict the performance of a low-Reynolds-number Savonius wind rotor.

Torque and Power Characteristics
Rotor performance for the Savonius turbine design variants presented in Table 1 are presented to quantify their power and torque characteristics. These simulations were performed using a constant wind speed constant, at tip speed ratios from λ = 0.4-1.2. Figure 7 shows the coefficients of the averaged torque on the turbines with respect to λ. It can be seen that the averaged torques decrease with the increase in λ. Except for the case n = 0.5, torque decreases approximately linearly for λ > 0.6. The curve corresponding to n = 0.5 shows a different trend, dropping more slowly with λ before λ reaches 1. The turbine design with n = 1.5 has the highest torque coefficient at lower tip speed ratios (λ ≤ 0.6). While at higher tip speed ratios (λ > 0.6) turbine with n = 1 highest torque coefficient.
where P is the shaft power. It can be seen from the figure that there is a peak value for each curve and the coefficient of power increases with λ up to a certain point after which it drops down as λ further increases. Maximum power coefficients are achieved for each design at tip speed rations near λ = 0.8, except for turbine with a fullness factor of n = 0.5 which reach has a maximum value near λ = 1.
Taking the case 5 turbine as reference and comparing the coefficients of power at λ = 0.8, we find that as n increases, the power coefficient becomes smaller; while as n decreases, the power coefficient increases until n = 1 and then decreases.  The peak power coefficients produced by these designs are compared with that achieved by the highest performing design presented by (n = 2) in Table 2. Table 2 shows the benefit of quantifying the effect of blade fullness with a peak power coefficient of 0.2573 obtained at a tip speed ratio of 0.8, which is 10.98% higher than a conventional Savonius turbine.  Figure 9 shows the maximum coefficients of power with respect to blade fullness. By using polynomial curve fitting, the peak coefficient of power can be expressed as: The peak value of Equation (4) is 0.2592, obtained at n = 1.099, which is near to the case 3 turbine.

Mechanism of the Effect of Blade Fullness on the Turbine Performance
Blade fullness not only affects the average power and torque values presented in the previous section, but also the instantaneous torque values that are a function of blade rotation angle. Figure 10 shows the dynamic torque on a single blade for different blade fullness as a function of azimuth angle for λ = 0.8. The torque generated by a single blade becomes positive as azimuth angle increases above θ = 320° and reaches maximum values at between θ = 30°. The blade torque becomes negative between θ = 130° and reaches minimum values between θ = 270°. Compared to the traditional Savonius turbine (case 5/n = 2), we find that as n decreases the maximum torque increases, with n = 0.5 achieving the highest maximum torque. However, as n increases the maximum torque values only minimally change, but the torque curves are lower over most of the curve. It can be seen that the blade designs with the higher maximum torque also generates lower minimum torque. These two impacts have opposite contributions to the averaged torque: higher positive torques increase the averaged torque while lower negative torques work decreases averaged torque. Therefore, there must be an optimal blade fullness, at which the turbine generates the maximum averaged torque, and eventually averaged coefficient of power. In order to further investigate the reasons that caused the difference of torque in Figure 10, two typical transient azimuth angles were chosen for study: θ = 60° (one typical point of the positive torque region) and θ = 270° (one typical point of the negative torque region). Figure 11 shows the pressure distributions on all seven blade designs at λ = 0.8 and θ = 60°. The curves are plotted from the rotation center and along the chord direction. It can be seen that the pressure remains nearly constant along the pressure sides (the concave sides) of the blades, especially at higher n values. The blade with n = 0.5 has the lowest pressure on the pressure side, while the pressure on other blades remains near atmospheric pressure (0 Pa).  Figure 11. Pressure distributions on the blades for λ = 0.8 and θ = 60°. Figure 12 shows the velocity contours and streamlines near all seven blade designs for λ = 0.8 and θ = 60°. In the case n = 0.5 (Figure 12a), a vortex that rotates counterclockwise can be observed above the pressure side of the blade. As n increases this vortex becomes weaker (n = 0.75 and n = 1), disappearing for n > 1. This vortex causes the low pressure on the concave side. Another vortex occurs near the center of the rotor when n > 1, however, this vortex is far from the pressure side of the blade and has little impact on the pressure. In Figure 11, the pressure on the suction sides (the convex sides) of the blades is negative across the whole blade, with a much greater dependence on location than on the pressure side. Overall, the pressure on the suction side increases as n increases, but the difference is very small for n > 2. Three negative peaks can be observed in pressure on the suction side ( Figure 11). The first peak occurs at the tip of the blade with a value of about −250 Pa on all the curves. This peak is caused by the shedding of vortex at the tip of the blade, as shown in Figure 12. It should be noted that though the peak values have little difference, the curves near this peak become wider as n decreases, contributing to an increase in the torque of the turbine. The second peak of pressure lies at between 60% and 90% of the blade chord, depending on the blade fullness. Except for the case of n = 0.5, this peak of pressure has a value of about −110 Pa. This peak is caused by the high-speed stream which flows over the convex surface of the blade (as indicated by the streamlines in Figure 12). The position of the second pressure peak moves outwards along the chord direction as n increases, due to the increased curvature of the blade near the tip. In Figure 12, a vortex on the convex side of the blade and near the rotation center can be seen. This vortex is caused by the rotation of the rotor and leads to the third negative peak of pressure (near the root) on the suction side of the blade. Figure 13 shows the pressure distributions on the seven blade designs for λ = 0.8 and θ = 270°. On the suction side (the concave side) of the blade, the pressure is negative and nearly constant for each blade. The case n = 0.5 has the lowest pressure, followed by n = 0.75. The pressure on other blades are much the same, all around −30 Pa. Figure 14 shows the velocity contours and streamlines near the turbine for λ = 0.8 and θ = 270°. For n = 0.5 (Figure 14a), a counterclockwise rotating vortex can be observed above the suction side of the blade, which causes the low pressure on the concave side. This vortex becomes weaker at n = 0.75 and disappears as n increases further. Therefore, a pressure recovery is seen at the suction side of the blade. Additionally, a second clockwise rotating vortex was observed closely behind the above mentioned vortex. This vortex has a diameter of approximate 0.5D and has approximately the same speed for n > 1, therefore the pressure at the blade changes little when n > 1. On the pressure side (the convex side) of the blade, a peak of pressure (about 55 Pa) is observed at about 50% chord. This peak corresponds to the stagnation point on the convex side of the blade ( Figure 14). In general, a blade with smaller n generates higher pressure at the tip and lower pressure near the root (Figure 13). It can be seen in Figure 14 that a blade with a smaller n changes smoother in geometry near the root and the pressure is reduced. Blades with smaller n have a flatter shape at the tip and fully attached flow at the tip is observed for cases n ≤ 1. However, when n further increases separation occur at about 0.9c and the pressure is reduced to negative. The higher positive pressure at the convex side and lower negative pressure at the concave side of a blade with smaller n leads to a larger drag and reduces the total torque at this azimuth position. Figure 15 shows the velocity contours of the wake for λ = 0.8 and θ = 60°. For all rotor designs the wake expands as the turbine is approached from upstream. The turbines with n ≤ 1 have wider wakes than those of n > 1. The wake of turbines with n > 1 enlarges to approximately 3D in width at 5D downstream of the turbine and 5D in width at 10D downstream of the turbine. For the cases of n ≤ 1, the wake enlarges to approximately 4D in width at 5D downstream of the turbine and 6D at in width at 10D downstream of the turbine. The upstream fluid is slowed when it passes the turbine and forms a low-velocity zone just behind the turbine. The maximum velocity deficit is 90% of the free stream. This low-velocity zone extends up to about 5D behind the turbine and is well characterized in the wake domain where a finer grid is used. The design with n = 0.5 has the smallest low-velocity zone but the strongest far wake. Discrete high-velocity zones lie on both up and down sides of the low-velocity zones. These high-velocity zones are caused by the vortexes periodically shedding from the tips of the blade.

Conclusions
In this paper we introduced a design technique for defining the geometry of Savonius wind turbines that is based on the Myring Equation. Two-dimensional CFD simulations were performed on this type of turbine to quantify rotor aerodynamic performance and the wake generated by design variants. More specifically, the effect of the blade fullness on the turbine's performance was studied. Based on the results in the paper the following conclusions were drawn: 1) The rotor with a blade fullness of n = 1 has the highest coefficient of power, 0.2573, which is 10.98% higher than a conventional Savonius turbine. 2) During a rotation period, the blade with a smaller fullness generates both higher positive torque values and lower negative torque values. 3) Turbines with n ≤ 1 have wider wakes than those of n > 1.