Abstract
The hybrid vortex method, in which vortex panel method is combined with the viscous-vortex particle method (HPVP), was established to model the wind turbine aerodynamic and relevant numerical procedure program was developed to solve flow equations. The panel method was used to calculate the blade surface vortex sheets and the vortex particle method was employed to simulate the blade wake vortices. As a result of numerical calculations on the flow over a wind turbine, the HPVP method shows significant advantages in accuracy and less computation resource consuming. The validation of the aerodynamic parameters against Phase VI wind turbine experimental data is performed, which shows reasonable agreement.
1. Introduction
Among various aerodynamic theories of the wind turbine rotor aerodynamic models [1, 2], vortex theories have received wide attention due to their advantages of affordable computational costs and good calculation results. Since the vortex methods were discovered to be suitable for the calculation of the unsteady and incompressible flows with a viscous boundary, relevant studies [3–7] have been conducted in the past decades. Referring to the introduction in [8], the aerodynamic model of wind turbine blade constructed by vortex method in lift line, lift surface, or panel was used to describe the flow around the blade surface. To calculate the trailing edge vortices induced by the wake, the free-wake models [9–11] were employed to describe the wake geometry and the wake. The free vortex wake method based on the Lagrangian approach was used to simulate the free distorted vortex wake structure, which has superior capability to predict aerodynamic performance of wind turbine compared with the BEM [12, 13] method. It was proved especially that the free vortex wake model was more applicable to simulate a wider range of working conditions such as yaw and tower-shadow [14]. However, considering the computational cost, many previous studies [15, 16] on free vortex wake methods applied the simplified wake models neglecting the intimal vortex sheet and the root vortices. Actually, the vortex sheets downstream shed by the wind turbine blades have a full spanwise distribution with complex structure. The neglect of the intimal vortex sheet and its development will result in large errors in the calculation of the blade aerodynamic performance.
The flow field around wind turbine blade is characterized with high Reynolds number and low Mach number. The viscous diffusions of vortices are mainly reflected in the blade surface and wake region. According to the principle of the vortex method, an accurate solution to Navier-Stokes equations will be obtained in the condition of the full-scale vortex particle simulation combined with the vortices field of the blades and wakes. It is hard to perform this simulation because the blade surface eddy scale is very small. The simulation of this surface eddy requires huge of vortex particles and the calculation equals the direct numerical simulation. Actually, when the flow separation is small, the viscous diffusion of the blade surface is concentrated in a thin viscosity layer, which can be regarded as a potential flow vortex sheet and the vortex sheet strength can be quickly obtained by applying the theory of the potential flow. When vortex sheet releases from the blade trailing edge to nonboundary wake zone, the vortex diffusion becomes clear and the vortex scale becomes bigger and bigger. Therefore, the discrete vortex particles can be used to simulate the flow in viscosity eddy which is far from the wake region. The idea of hybrid vortex method was derived from the vortex sheet method by Chorin [17]. In this hybrid vortex method, the sheet vortex motions in the boundary layer are described by the Prandtl approximation of the Navier-Stokes equations. The computational field is decomposed into an interior domain where the conventional vortex particle method is applied and a user-specified numerical boundary layer where the vortex method is adopted [7].
The work presented here targets establishing a hybrid panel and viscous-vortex particle method (HPVP) and developing corresponding computation procedure for the aerodynamic analysis of wind turbine. The vortex surface based on potential flow theory is employed to solve the flow around the blade surface and the region near trailing edge. Meanwhile, the turbulent boundary layer surface friction formula is used to correct the blade surface viscous force. In the zones far from the blade wake, the viscous-vortex particle method is adopted to solve the flow equations. The coupling transfer between the vortex particle and vortex panel is built based on the equivalent principle of the wake panel and the vortices strength. This calculation program is used to conduct a computation case of NREL Phase VI rotor and verified by experimental measure data.
2. Numerical Formulation
In this section, the governing equations for the hybrid vortex method and its formulations for numerical calculations are presented.
2.1. Panel Method
2.1.1. Governing Equations and Boundaries
In the rotating coordinate system (), the rotation direction of the wind wheel toward the -axis, the flow velocity is . The boundary surface of the flow field of the wind turbine is composed of the wheel surface , the trailing vortex surface , and the surface of the infinity incoming flow , as shown in Figure 1. The velocity field is potential due to the induction of wind wheel, and the disturbance potential velocity meets the Laplace equation:
According to the Green function, any point of the perturbation potential in the flow field is represented by the following equation [18]:where, when included in the , ; at the boundary of , ; and at outside of , . is the length of the field point and the point . represents the normal derivative of the point , which satisfies the various boundary conditions as follows [19]. (1)On infinity of the outside interface, the disturbance velocity of the wind turbine tends to zero:(2)On the wheel surface , when the normal velocity is zero, the boundary conditions of the movement are as follow:(3)It assumes that the wind turbine wake is zero thickness plane and parallel to the flow line. On both sides of the trailing vortex surface, the normal velocity difference and pressure difference should be zero: where represents the point located on wake vortex surface and superscripts + and −, respectively, represent the wake of the upper surface and the lower surface. The above three boundary conditions are substituted into (2), written by where is the velocity potential difference on both sides of the surface of the wake vortex:
2.1.2. Equation Discretion
The calculation grids of the blade surfaces and the wake regions are generated, respectively; each mesh surface is called the surface element. It is assumed that the perturbation potential , the source strength , and the wake vortex strength on the blade surface element are evenly distributed. Therefore, the blade surfaces are discreted by surface elements, including unknown velocity aswhere and are the blade surface panel and the trailing vortex surface element, respectively. When , , or, else, . , , and are defined as
2.1.3. Numerical Calculation
Firstly, the body-fitted coordinate of the blade surface is established. Then, the approximate basic functions are selected which are the three adjacent velocity potentials on each direction of the blade surface. Then, each coefficient of the model function is solved by the method of undetermined coefficients and the derivative of each panel control point can be obtained, namely, the tangential disturbance velocity. Finally, the general tangential velocity is obtained by adding tangential disturbance velocity with the tangential component of relative inflow velocity .
When the surface velocity of the blade is determined by the velocity potential, the pressure distribution of the blade surface can be calculated according to the unsteady Bernoulli equations. The dimensionless pressure coefficient can be expressed as follows:
The rotor thrust and torque can be calculated in numerical integration method based on the solved pressure :
2.2. Viscous-Vortex Particle Method
It is assumed that the flow field of wind impeller is incompressible; thus the blade vorticity field is described by the three-dimensional incompressible viscous flow. According to the vortex principle, the Navier-Stokes equations can be expressed as velocity-vorticity ():wherein is the dynamic viscosity coefficient, is the vorticity, and is the velocity.
In order to avoid the solution error caused by the numerical dissipation in Euler coordinate system, the Lagrangian system is used to describe the vorticity field, and the vorticity field is discreted into vortex particles carrying vorticity, volume, and quality:where is vortex particle position, is the particle vortex vector, is a smoothing parameter, and is Gaussian distribution function:where represents the dimensionless distance.
The equations can be divided intowhere the first equation describes the convection effect of vorticity and the second equation represents the vorticity dissipation and stretch effect.
According to the Biot-Savart law, the particle velocity is expressed aswhere is the Kernel function of the Biot-Savart law.
For the vorticity dissipation term, many methods can be used to simulate it, for example, random-walk method [20], vortex core diffusion method [21], Fishelov method [22], and particle strength exchange (PSE) [23]. In this paper, the PSE method is adopted to simulate vorticity dissipation term: represents the Gaussian function; that is, .
2.3. Coupling Method of Panel and Viscous-Vortex Particle (HPVP)
To accurately calculate the viscous diffusion effect of the wake, the viscous discrete vortex is adopted to replace the dipole wake surface element, and the relationship between the unsteady panel method and viscous-vortex wake method is established, as shown in Figure 2; the trailing edge dipole surface element is called Kutta surface panel which meets the Kutta conditions. A group of dipole surface elements behind the Kutta surface panel are called the transition surface element, which is used to be equal to the vortex particle. It is the dipole surface element obtained from the last time step of Kutta surface panel in time-marching; the intensity and location are known. The wake vorticity fields after transition surface elements are discretized by the vortex particle, and these vortex particles are obtained by equating the transition surface element to the vortex particle in time-marching. Therefore, the vorticity field is constituted by the blade dipole surface elements, Kutta surface panels, transition surface elements, and discrete vortex particles. All of the dipole surface elements are solved in unsteady panel method and the vortex particles are solved by viscous-vortex particle method. The transition surface elements enable the coupling of the two methods transferring the data effectively.
According to Hess equivalence principle, the induced velocity of dipole bin with the strength at any point is equivalent to vortex strength and boundary line vortex ; namely,wherein is the location vector; is the line vortex vector; and represents the boundary line of wake surface element.
At time step , the fixed Kutta surface elements are generated and shed from the blade trailing edge. According to the unsteady panel method, the dipole distributions of blade surface and the Kutta surface elements can be solved, and then the move velocity of Kutta panel is calculated. The location for the next time step is solved in the fourth-order Runge-Kutta method, which constitutes the transition bins for the next time step. At time step , the dipole surface elements are converted into vortex particle according to the equivalent principle, and then the dipole intensity distributions of the blade surface and Kutta panel are calculated; this process is repeated step by step.
When the upstream wake flows through downstream wind turbine blades, the interference of wake vortices and blade vortex occurs; therefore the impact of wakes on the blade surface element should be considered. In this paper, the surface element distributions are determined by the Neumann object plane boundary conditions and calculated by the induced velocity of vortex particle on blade surface element. The intensity of the panel element is defined as follows:where is the induced velocity of vortex particle, which is calculated by Biot-Savart law.
3. Calculation Process
The calculation programs are developed based on the hybrid panel and viscous-vortex particle method (HPVP); the process diagram is shown in Figure 3. Firstly, the blade surfaces and Kutta panels are generated according to the geometric model of wind turbine blade. Then, the panel method is employed to calculate the dipole intensity of blade surfaces and the Kutta surface elements. And then, the Kutta bins are propelled by time-marching method, and the viscous-vortex particles are generated by equivalent principle; meanwhile, the downstream existing particle vorticity field is updated. Then, the dipole intensity of blade surfaces and Kutta surface elements in new time step are calculated. Repeat the above process until it meets the convergence criteria.
4. Numerical Results
The NREL Phase VI wind turbine is selected as a validation numerical example. The aerodynamic performance results of the NREL Phase VI wind turbine is gotten by the measurement in NASA Ames wind tunnel of 24.4 m 36.6 m. The experiments are detailed and complete; the measuring is accurate, which has become the verification basis of the academic research and engineering application [24]. The experiment data will be used for comparisons with the calculation value in the following sections. The NREL Phase VI rotor has two blades with the radius of 5.03 m; each section airfoil is NREL S809.
In calculation, the 7 m/s inflow wind speed condition is adopted. In this condition, most of the flow area is attached to the blade surface apart from a small area flow separation near the blade root, which accords with the flow field characteristics of variable pitch wind turbine. The numerical validation shows that 50 bins arranged along blade spanwise, 70 bins arranged along blade circumferential, and 120 vortex particles divided at wake region are enough for a grid-independent solution. The calculation program can be running in normal quad-core processor and finished with less than 300 seconds computer CPU time.
4.1. Wake Vortex and Streamline
Figure 4 shows wake vorticity calculated in hybrid vortex method and CFD method. As can be seen from Figure 4(a), the blade wake vortex is smooth and complete due to the superior vorticity conservation and small numerical dissipation in vortex particle method. As a reference, Figure 4(b) presents the wake vorticity field calculated in CFD method by Sezer-Uzol and Long [25]; in their work, almost 9,600,000 grids are generated and the vortex core zones are refined to conduct the CFD computation, but the wake vortex still shows obvious dissipation. It is mainly due to the inherent numerical dissipation which causes the rapid decay of the vorticity.
(a) HPVP method
(b) CFD method
Figure 5 shows the limit streamline distributions on blade suction surface (SS) calculated by the CFD and HPVP method. As can be seen from Figure 5, the simulation of the flow on blade surface in HPVP method is consistent with the CFD method in most blade regions except for the zone at blade root, some small flow separation is not accurately simulated by HPVP method. The essential reason is that the panel method is based on the assumption that the body-flow principle is used to solve the potential flow through the blade surface.
(a) HPVP method
(b) CFD method
4.2. Pressure Coefficients
For this comparison, the conditions of the measurement and the calculation are set by the same conditions which ensure the comparability. Figure 6 presents the experimental verification of the pressure coefficients distribution at various blade cross-sections. The sections, respectively, located at , 47%, 63%, 80%, and 95%, are chosen. The test and simulation in HPVP method are under the axial uniform wind speed is 7 m/s.
(a)
(b)
(c)
(d)
(e)
By comparing the pressure distributions at various blade cross-sections in Figure 6, it can be found that the numerical simulation of blade pressure in HPVP method match well with that in experimental test. The errors between the test value and the calculated one fall within a reasonable range. The average error gradually decreases from the blade root to the blade tip because the degree of separation flow from the root of the blade to the tip of the blade is decreasing. At the cross section next to the blade root (), the average error is about 10%, while at the blade tip position (), the average error drops to about 4%. In addition, it can be found that the errors of the leading edge and the trailing cross-sections are slightly larger than those at the middle chord zones. This is mainly because the panel method is based on potential flow principle for solving, while the actual fluid contains viscous term which results in the friction on blade surface. Meanwhile, this is also because the varying degrees of the flow separation do exist at the trailing edge of the blade. Therefore, the assumption of the inviscid potential flow in panel method will underestimate this effect of viscous drag, which causes the calculation errors at the leading edge and trailing edge of the blade to be slightly larger.
4.3. Tangential and Normal Force Coefficients
Figure 7 shows the comparisons of the tangential and normal force coefficients between the calculated value and the experimental data. The wind turbine is operating at 7 m/s wind speed condition. It can be found that the tangential and normal forces calculated by HPVP method are consistent with the test values in most areas along the blade span. But at the blade root region , the calculation error of the tangential and normal force coefficients is large; the normal force coefficient especially is larger. This is mainly due to the inaccuracy simulation of the separated flow in HPVP method at blade root area.
(a) Tangential direction
(b) Normal direction
Table 1 shows the comparisons of the torque and axial thrust of the wind turbine. The error of the torque and axial thrust between the calculated values and the experimental data is 4% and 3.4%, respectively, which indicates the good accuracy of the HPVP method for calculating the aerodynamic performance of the wind turbine. However there are some calculation errors compared with the experimental data when simulating the pressure coefficient and aerodynamic loads at blade root zones in HPVP method, thus contributing little to the overall turbine performance. This is mainly because of the mutual influence of integral compensation effects in the HPVP method and the little contribution of the blade root to the overall aerodynamic of the wind turbine.
5. Conclusions
In this paper, a hybrid panel and viscous-vortex particle method (HPVP) has been established and corresponding computation procedure for the aerodynamic analysis of wind turbine has been developed. Phase VI wind turbine was selected as a simulation example to validate the accuracy of the HPVP method with the measured data. It showed the good capability of the HPVP method in the performance prediction of the wind turbine.
The HPVP method has the advantages of less empirical correction, less computing resource consuming, and accuracy. Compared to the CFD method, the calculation in HPVP method can quickly obtain accurate results and take up only one-thousandth of CFD computing time. The simulation in HPVP method considers the viscous effect of wake vortex, which can be closer to the exact solutions to the equations compared to the potential flow vortex method.
Some other detailed flow parameters of wind turbine are also solved in HPVP method besides the surface pressure distributions and the limit streamlines of the blade. It is helpful for conducting the investigation of the aerodynamic and aeroelastic issue of the wind turbine complex, which has great significance on the optimization design process of the wind turbine.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was funded by Henan Provincial Key Science and Technology Projects (142102210059). This support is most gratefully acknowledged.