Computational Fluid Dynamics Analysis of Multi-Bladed Horizontal Axis Wind Turbine Rotor

Nasim A. Mamaghani, Peter E. Jenkins

Department of Mechanical Engineering, University of Colorado, Denver, CO, USA.

DOI: 10.4236/wjm.2020.109009                

Abstract

The principal objective of this work was to investigate the 3D flow field around a multi-bladed horizontal axis wind turbine (HAWT) rotor and to investigate its performance characteristics. The aerodynamic performance of this novel rotor design was evaluated by means of a Computational Fluid Dynamics commercial package. The Reynolds Averaged Navier-Stokes (RANS) equations were selected to model the physics of the incompressible Newtonian fluid around the blades. The Shear Stress Transport (SST) k-ω turbulence model was chosen for the assessment of the 3D flow behavior as it had widely used in other HAWT studies. The pressure-based simulation was done on a model representing one-ninth of the rotor using a 40-degree periodicity in a single moving reference frame system. Analyzing the wake flow behavior over a wide range of wind speeds provided a clear vision of this novel rotor configuration. From the analysis, it was determined that the flow becomes accelerated in outer wake region downstream of the rotor and by placing a multi-bladed rotor with a larger diameter behind the forward rotor resulted in an acceleration of this wake flow which resulted in an increase the overall power output of the wind machine.

Keywords

Computational Fluid Dynamics, Horizontal Axis Wind Turbine, Multi-Bladed Rotor, Aerodynamic Torque

1. Introduction

The escalation of the global energy demand has led countries to consider renewable energy sources more seriously. In the 2019 released version of the International Energy Outlook, the U.S. Energy Information Administration projects that the world energy consumption will grow by nearly 50% between 2018 and 2050 [1]. The 2019 report of the International Renewable Energy Agency claims that the cost of renewable energy technologies will continue to decline throughout the decade [2], and that solar and wind power are the two most affordable energy solutions for markets worldwide. Following the decline of the weighted-average cost of on-shore wind power during the past years, the report also anticipates the on-shore wind power technologies will be considerably less expensive than the existing conventional coal-fueled plants. The Annual Energy Outlook 2020 predicts that wind and solar energy would increase to 80 percent of the total renewable energy produced by 2050 and the growth of renewable energy would take place in all parts of the U.S. with the West and Mid-Continent areas experiencing the biggest increase in energy production from wind [3]. The steady advancements in wind energy technologies, such as rotor designs and manufacturing processes, have caused a reduction of the levelized cost of wind produced electricity which has resulted in its growth in the global energy market [2].

Over the past two decades, the traditional experimental methods of studying wind turbine fluid dynamics have decreased [4] due to lack of prototypes for the designs [5] and the increase in the large scale of the modern wind rotors. The experimental studies are labor-intensive, whereas the computational fluid dynamics and semi-empirical methods are cost- and time-effective and spatial friendly. These methods enable designers to apply geometrical modifications to their existing models and these methods are used in the wind turbine industry.

The work presented in this paper is a numerical study of the flow characteristics around a multi-bladed rotor wind turbine, a novel concept patented by the Thunderbird Power Corp [6]. The invention introduces a horizontal wind machine comprised of multiple sails that are supported by several arms which span outward from the hub. Also, a shroud is attached around the rear side of rotor, enabling the machine to capture energy of the incoming wind. The patent also discloses the addition of a second rotor placed in series with the first one. Each rotor is coupled with an individual shaft, where both shafts are coupled through a clutch mechanism. The patent claims that the multi-rotor concept was capable of delivering more power compared to a single stand-alone rotor of the same diameter. This turbine can be used for electric power generation and water pumping [6]. Figure 1 and Figure 2 illustrate the multi-bladed and multi-rotor design configurations.

A high-solidity windmill generates more torque at low-tip speeds compared to a modern horizontal axis wind turbine [7]. This means that a wind machine with a higher solidity can operate in wider range of wind speeds, leading to an increase in overall power output of the turbine. The model which was studied throughout this paper was a simplified version of the above-described patented wind machine design and consists of a single multi-bladed rotor. The work conducted analyzed the effects of the high-solidity rotor on the wake flow and studied

Figure 1. Multi-bladed wind machine concept with a shroud, isometric view [6].

Figure 2. Multi-rotor multi-bladed wind machine concept, side view [6].

the torque and power outputs of the wind turbine with different rotational speeds and varying wind velocities. In order to evaluate the proper location of the subsequent rotors, velocity profiles of the rotor wake were thoroughly investigated.

2. Methodology

Fluid motion can be described by a set of differential equations, namely the mass and momentum conservation equations, known widely as Navier-Stokes Equations (NSE). In the absence of gravity, for an incompressible Newtonian fluid, the NSE reads as Equations (1) and (2), respectively.

$\frac{\partial U_i}{\partial x_i}=0,$ (1)

$\frac{\partial U_i}{\partial t}+U_j\frac{\partial U_i}{\partial x_j}=\frac{-1}{\rho }\frac{\partial P}{\partial x_i}+\nu \frac{{\partial }^2{U}_i}{\partial x_j\partial x_j}+\frac{F_i}{\rho },$ (2)

where $x_i$ denotes the position vector, ρ stands for the fluid density, P is the pressure, U is the fluid velocity and ν is the kinematic viscosity. External body forces, denoted by F, act on a section of the fluid and were evaluated from the fluid rotational forces, such as the Coriolis and centrifugal forces. Osborne Reynolds developed a proposal for eliminating the turbulence unsteady fluctuations by averaging the flow quantities which lead to the so-called Reynolds Averaged Navier-Stokes (RANS) equations for the mean flow, Equations (3) and (4).

$\nabla \cdot u=0$ (3)

$\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)u=\frac{-1}{\rho }\nabla p_{eff}+\nu {\nabla }^{2}u-2\Omega \times u-\nabla \cdot \overline{u^{\prime }{u}^{\prime }}$ (4)

The stress tensor term in the momentum equation (Equation (4)), $\overline{u^{\prime }{u}^{\prime }}$, corresponds to the effects of turbulence on the mean flow and results in a closure problem in the RANS equations. The semi-empirical standard k-ε model introduces two equations (Equations (5) and (6)), involving turbulent kinetic energy, k, and turbulent dissipation, ε, to RANS which closes the system.

$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+P-\rho \epsilon$ (5)

$\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\frac{\partial }{\partial x_i}\left(\rho \epsilon u_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\frac{C_{1\epsilon }\epsilon }{k}P-\frac{C_{2\epsilon }\rho {\epsilon }^2}{k}$ (6)

P is the production of k, ${\mu }_t$ is the eddy-viscosity and are defined respectively as in Equations (7) and (8),

$P=-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}\frac{\partial u_j}{\partial x_i},$ (7)

${\mu }_t=\frac{\rho C_{\mu }{k}^2}{\epsilon }$ (8)

The k-ε model assumes the flow to be fully turbulent and is reliable for high Reynolds regions only. The k-ω model, on the other hand, utilizes the transport equations for k (Equation (9)) and turbulence frequency, ω (Equation (10)), to solve the turbulent viscosity.

$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k-Y_k+S_k,$ (9)

$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_i}\left(\rho \omega u_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right]+G_{\omega }-Y_{\omega }+S_{\omega },$ (10)

where G is for the production terms, Y the dissipation and S represents the user-defined source expressions. The k-ω model showed better agreements with the real flow behavior for the viscous sublayer regions and hence was chosen over the k-ε model. Coupling the k-ε and k-ω models with a blending function introduces the k-ω Shear Stress Transport (SST) model [8]. This model shifts between the k-ω model for the near wall regions and the k-ε for the far field regions through the domain depending on cell distance from the closest wall boundary. The k-ω SST can also predict the flow separation in regions close to the walls by using its embedded viscosity limiter expression which damps out the shear stress in the vicinity of walls.

For the rotating fluid domains, such as flow around rotating blades and impellers, it was possible to simulate the unsteady nature of the problem in a steady-state manner, where the mesh motion was such that it followed the motion of the geometry. This approach was more time efficient compared to the unsteady analysis and reduced the amount of computational power needed and was a good method for conducting preliminary studies of turbomachinery problems. In this study, a steady-state Single Moving Reference Frame (SRF) method was utilized to analyze the flow domain where the RANS equations, (Equations (3) and (4)) were solved in a rotating frame of reference to simulate the steady-state condition and an additional source term was applied to the flow equations throughout the domain. The equations of energy, momentum, continuity and transport of species were solved in a segregated manner based on pressure-based solver algorithm which was ideal for incompressible flows and resulted in memory and time efficiencies.

3. Physical Model and Boundary Conditions

In the present study, a simplified version of the original patented multi-bladed rotor was modeled by a SOILDWORKS CAD package. Since the blades were evenly placed around the rotor at 40 degree intervals, one-ninth of the whole rotor was chosen to be studied and the effects of the remaining blades were taken into account through applying a periodic boundary condition where it was required. Analyzing the effects of other parts of the wind turbine, including the hub, tower and arms was beyond the scope of this work. For the current design the rotor diameter was set to 3.5 ft and it consisted of 18 identical sail with each of them have a length of 6.3 inches. Figure 3 and Figure 4 illustrate several views of the rotor geometry along with the one-ninth periodic unit shown by dashed lines.

Defining a domain to any problem introduces errors to the flow solution as the boundaries should naturally be at an infinite distance from the object. However, this distance was limited in practice. The dimensions of computational domain must be large enough to predict the turbulence phenomenon, pressure and velocity distribution to a reasonable degree. For this study, the domain was a single rotatory region where the free stream flow entered the flow domain at a distance of 1.5 L upstream of the rotor and exited 2 L distance downstream of the rotor, where L stands for length of the sail. The entire domain was one-ninth of a cylinder with a diameter of 1.3 D, with D being the rotor diameter. The geometry of computational fluid domain is illustrated in Figure 5.

Considering the rotation of the rotor in the clockwise direction, the frame motion was activated in the model setup with its axis of rotation being set in the flow direction, rotating as the negative of rotational velocity of the rotor. The inlet velocity boundary condition was used for the upstream surface of the domain by changing the magnitude of the wind velocity. The downstream surface boundary condition was set to be the outlet gauge pressure. As the upper and lower surfaces of the computational domain had no impact on the flow solution, they were adjusted as symmetry boundaries. Periodic boundary conditions were used for the side surfaces as the crossing flow from these faces are identical with each other. Lastly, the surface around the sails were set to no-slip boundaries. The turbulence intensity and viscosity ratio remained as software defaults at 5% and 10%, respectively. Figure 6 illustrates the boundaries of the computational domain.

Figure 3. Multi-bladed rotor geometry views.

Figure 4. Multi-bladed rotor geometry views.

Figure 5. Computational domain geometry.

Figure 6. Boundary conditions on fluid domain.

The governing equations, Equations (3) and (4), were solved in the rotating frame of reference at the same speed as the rotor rpm. The fluid domain was discretized into approximately 440 thousands cells using the meshing software included in the ANSYS Fluent package. The solution to the problem reached to a stable condition for this number of elements and this behavior was shown in the grid independency study (Figure 21). The domain was divided into structured and unstructured mesh regions to speed up the solution process as was shown in Figure 7. The unstructured tetrahedral elements were used to cover the inner section of the fluid domain, due to their ability to adapt to complex geometries. The outer sections of the computational domain were tiled with structured hexahedral elements. Figure 8 demonstrates the hexahedral and tetrahedral elements. To capture the fully turbulent behavior of the flow, several layers of the prismatic wedges with a growth rate of 1.2 were added around the blade surfaces and targeted the y+ values of higher than 30 (Figure 9). Use of the prism layers resulted in a reasonable convergence rate due to the reduction of the numerical diffusion. Global and local meshing techniques, such as edge sizing and body sizing, were engaged to achieve a high overall mesh quality. Average mesh metrics of 0.24 and 0.75 for the skewness and orthogonal qualities were attained, respectively. The match control method was applied on the periodic interface boundary conditions. Fluid properties and the applied numerical scheme for the simulation are listed in Table 1 and Table 2, respectively.

Figure 7. Computational domain was discretized into structured and unstructured regions.

Figure 8. Fluid domain was meshed in hexahedral and tetrahedral elements.

Figure 9. Use of prism layers to capture the fully turbulent flow behavior (near blade view).

Table 1. Applied fluid properties.

Table 2. Applied numerical scheme.

4. Results and Discussion

For this analysis, the incoming free flow was parallel to the ground in the negative Z direction. Flow visualizations were done with four different wind velocities of 5, 10, 15 and 25 MPH which convert into approximately 2.2, 4.5, 6.7 and 11 m/s. The tip speed ratio of the turbine was kept constant at 0.7 which historically results in the highest achievable power coefficient for the American multi-bladed horizontal axis wind turbine. Sectional torque for each case was attained and was multiplied by factor of 9 to get the overall rotor torque, as the simulations were done through a one-ninth portion of the rotor. The rotor power was calculated by multiplying the torque value with the rotational speed of the rotor.

For a better understanding of the flow behavior in the wake region, several planes in parallel and perpendicular positions to the axis of rotation were created. The parallel plane was placed along the periodic surface of the computational domain, covering an axial distance of 3.5 L. Three perpendicular planes were assigned at distances of 0.4 L, 0.75 L and 1.25 L behind the rotor plane, as shown in Figure 10, and were used to monitor the axial velocity of the wake flow.

Figures 11-14 illustrate the axial flow characteristics of the wake behind the rotor at wind speeds of 5, 10, 15 and 25 MPH, respectively. The regions with highest axial velocities are circled with white color for the first set of incoming velocities.

As observed in the above figures, it was evident that there was a consistent flow behavior for all cases where the axial flow velocity increased in the outer regions behind each rotor sail, which was an expected trend as a result of the disk rotation. This increase in axial velocity tended to decrease as the flow traveled further downstream from the rotor. The drop in the axial velocity magnitude in the rotor disk’s shadow regions resulted due to the rotor power extraction from the free flow. The low velocity wake flow also recovered at far distances downstream from the rotor plane. Figure 15 shows the flow behavior in a plane parallel to the axis of rotation and supports the earlier conclusions. The darker blue regions above the blades in the wake region showed an acceleration in the axial velocity which determined the appropriate position for placing the second rotor, as was proposed in the patent [6].

Figure 10. Three planes at distances of 0.4 L, 0.75 L and 1.25 L behind the rotor, perpendicular to the axis of rotation.

Figure 11. Axial velocity contours at a wind velocity of 5 MPH (=2.2 m/s) at a parallel plane to the axis of rotation (top left), perpendicular planes at distances of 0.4 L (top right), 0.75 L (bottom left) and 1.25 L (bottom right) behind the rotor. The flow regions with the highest axial velocities are denoted with white colored circles.

Figure 12. Axial velocity contours with a wind velocity of 10 MPH (=4.5 m/s) at a plane parallel to the axis of rotation (top left), and perpendicular planes at distances of 0.4 L (top right), 0.75 L (bottom left) and 1.25 L (bottom right) behind the rotor.

Figure 13. Axial velocity contours at a wind velocity of 15 MPH (=6.7 m/s) at a plane parallel to the axis of rotation (top left), and perpendicular planes at distances of 0.4 L (top right), 0.75 L (bottom left) and 1.25 L (bottom right) behind the rotor.

Figure 14. Axial velocity contours at a wind velocity of 25 MPH (=11 m/s) at a plane parallel to the axis of rotation (top left), and perpendicular planes at distances of 0.4 L (top right), 0.75 L (bottom left) and 1.25 L (bottom right) behind the rotor.

Figure 15. Axial velocity vectors in a plane parallel to the axis of rotation for a wind speed of 25 MPH (=11 m/s). Black lines show the distance occupied by the rotor sails. The increase of the axial velocity was indicated by the darker blue colors.

By studying the pressure contours on the blade surface it showed that the surfaces which face the wind have higher pressure magnitudes compared to the rear surfaces, as expected. This pressure difference between both sides of the rotor sails caused lift in the direction of the rotor rotation. A closer look at the pressure-side of the blade showed the front edges, which are facing in the direction of rotation, have a higher pressure with respect to rest of the sail surface. However, the opposite holds true for the suction side of sails where the edges that faced the rotation direction are the regions with lower pressure. In Figure 16 the regions of high pressure on each face are denoted with an oval. It should be noted that regions of high pressure move into the whole sail surface as the wind velocity increases, as expected.

As mentioned earlier, the tip speed ratio throughout this analysis was kept at a fixed value of 0.7 and, therefore, the rotational speed increases linearly as the oncoming wind velocity increases. This behavior is shown in Figure 17. Figure 18 shows the increase of rotor torque at increasing values of the wind speed. Figure 19 shows that the power coefficient of the rotor improved at higher wind velocities. Finally, the increase of power extracted from the wind is shown in Figure 20.

Considering the limitations in maximum cell number of the fluid solver’s Academic Teaching license, the grid independency analysis was carried out to study the stability of the main output of the flow solution for the chosen number of mesh elements. It was showed that the fluctuations of the sectional torque were very minor for cell number of around 400,000 (Figure 21).

Considering the novelty of the evaluated rotor design, there was no available experimental data to be compared with the calculated result. However, the reported values for power coefficient in Figure 19 displayed a good consistency with the published statistical diagram of rotor power coefficient as a function of tip-speed ratio. According to Figure 22, for tip-speed ratio of 0.7 one can expect values of about 0.13 to 0.14 for the power coefficient for a multi-bladed American style HAWT. This behavior was well observed for the rotor of interest (Figure 19).

Figure 16. Views of the pressure contours at front surfaces (left) and rear surfaces (right) of the sails at a wind velocity of 5 MPH (= 2.2 m/s) (top) and a velocity of 15 MPH (= 6.7 m/s) (bottom).

Figure 17. Rotational speed of the rotor vs. wind velocity.

Figure 18. Torque output of the rotor vs. wind speed.

Figure 19. Overall power coefficient as a function of wind speed.

Figure 20. Power output of the rotor vs. wind speed.

Figure 21. Sectional torque values of the rotor at different. number of cells.

Figure 22. Statistical diagram of tip-speed ratio (horizontal axis) versus Rotor power coefficient (vertical axis) [9] [10].

5. Conclusions

A high-solidity windmill generates more torque at low-tip speeds compared to modern horizontal axis wind turbines [7]. This means that wind machine with higher solidity can operate in a wider range of wind speeds, leading to an increase in the overall power output of the turbine. The model that was evaluated was a simplified version of a recently patented wind machine [6], consisting of a single multi-bladed rotor. This evaluation analyzed the flow behavior around the rotor and evaluated the torque and power output of the wind turbine with different rotational speeds for varying wind velocities. The axial velocity trend of the rotor wake was thoroughly investigated in order to properly place the added rotors. It was evident that the flow becomes accelerated in the outer wake region downstream of the rotor and, as a result, placing a multi-bladed rotor with a larger diameter behind the front rotor would make use of the accelerated wake flow and therefore would improve the overall power output of the wind machine. For a specified rpm, the increase of solidity increased the torque production and, hence, leads to a higher wind turbine power output [11]. However, higher values of solidity cause undesirable blockage of incoming wind and impact the power extraction of rotor in a negative manner. Therefore, finding the optimum solidity degree which balances the high torque and flow blockage was of a high priority.

A more comprehensive study analyzing the effects of varying the number of blades, varying the sail geometries and evaluating the addition of rotors should be conducted. Considering the deficiencies of applying the steady-state modeling approach for flows with an unsteady nature, a study to apply an unsteady simulation algorithm in order to better evaluate different performance behaviors should be considered.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	EIA (2019) International Energy Outlook 2019 with Projections to 2050. U.S. Department of Energy, Energy Information Administration, Washington.
[2]	IRENA (2019) Renewable Power Generation Costs in 2018, International Renewable Energy Agency, Abu Dhabi.
[3]	EIA (2020) Annual Energy Outlook 2020 with projections to 2050, U.S. Department of Energy, Energy Information Administration, Washington.
[4]	Muiruri, P.I., Motsamai, O.S. and Ndeda, R. (2019) A Comparative Study of RANS-Based Turbulence Models for an Upscale Wind Turbine Blade. SN Applied Sciences, 1, Article number: 237. https://doi.org/10.1007/s42452-019-0254-5
[5]	Gómez-Iradi, S. and Barakos, G. N. (2008) Computational Fluid Dynamics Investigation of Some Wind Turbine Rotor Design Parameters. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 222, 455-470. https://doi.org/10.1243/09576509JPE526
[6]	Sutz, P.E. and Jenkins, R.K. (2018) Multiple-Blade Wind Machine with Shrouded Rotors. Patent No. 10066597.
[7]	Tangler, J.L. (2000) The Evolution of Rotor and Blade Design.
[8]	Menter, F.R. (2011) Turbulence Modeling for Engineering Flows.
[9]	Ragheb, M. and Ragheb, A.M. (2012) Wind Turbines Theory—The Betz Equation and Optimal Rotor Tip Speed Ratio. Fundam. Fundamental and Advanced Topics in Wind Power, 1, 19-38. https://doi.org/10.5772/21398
[10]	Eldridge, F.R. (1980) Wind Machines. 2nd Edition, The MITRE Energy Resources and Environmental Series, Van Nostrand Reinhold Company, New York.
[11]	Kentfield, J. (1996) The Fundamentals of Wind-Driven Water Pumpers. Gordon and Breach, Amsterdam.