Abstract
Analyses of the unsteady flow behaviour of a 5 MW horizontal-axis wind turbine (HAWT) rotor (Case I) and a rotor with tower (Case II) are carried out using a panel method and a RANSE method. The panel method calculations are obtained by applying the in-house boundary element method (BEM) panMARE code, which is based on the potential flow theory. The BEM is a three-dimensional first-order panel method which can be used for investigating various steady and unsteady flow problems. Viscous flow simulations are carried out by using the RANSE solver ANSYS CFX 14.5. The results of Case I allow for the calculation of the global integral values of the torque and the thrust and include detailed information on the local flow field, such as the pressure distribution on the blade sections and the streamlines. The calculated pressure distribution by the BEM is compared with the corresponding values obtained by the RANSE solver. The tower geometry is considered in the simulation in Case II, so the unsteady forces due to the interaction between the tower and the rotor blades can be calculated. The application of viscous and inviscid flow methods to predict the forces on the HAWT allows for the evaluation of the viscous effects on the calculated HAWT flows.
1. Introduction
One of the best alternative renewable energy resources today is wind energy: it is free and clean, meaning that designing and constructing wind turbines (WT) is a rapidly growing field of technology. However, maximising power output and minimising operation costs require a continuous improvement of their aerodynamic performance. An accurate prediction of the aerodynamic properties is still a challenge, especially since the dimensions of system have become larger. In particular, wind turbine blades can experience large changes in angles of attack associated with sudden large gusts, changes in wind direction, atmospheric boundary layer influence, and strong interaction with tower.
There are several methods of varying levels of complexity that can be used to predict the aerodynamic loads on the wind turbine (WT) aerodynamic parts. The blade element momentum method has been very popular for WT design and analysis as shown in Ingram [1]. This method is highly efficient and inexpensive, but it is incapable of modelling three-dimensional flow effects, the interaction between rotor and tower, and tip relief effects. Empirical corrections are used to consider tip vortex losses. Many researchers, such as Shen et al. [2] and Ceyhan [3], have attempted to increase the accuracy of the blade element momentum method by improving the model for considering the tip vortex effects on the flow in the outer radii region.
Potential flow based methods have been introduced to simulate the WT aerodynamics with high computational efficiency. Different formulations of the potential flow method can be applied, such as lifting line, vortex lattice, and panel method. Generally, in these methods, the blade can be modelled as a lifting line, lifting surface, or panels and the wake can be modelled by either trailing vortices or vortex ring elements. Compared with the blade element momentum method, potential flow methods can be applied to study more complex flow phenomenon such as the consideration of the nonaxial inflow velocity distribution due to the tower influence on the inflow as well as the interaction between the rotor and the tower flows. But such methods suffer from the limitation of the potential flow theory regarding the consideration of the viscous effects. Potential flow methods are not able to predict important viscous phenomena such as stall. Abedi et al. [4] used three different approaches which are lifting line, lifting surface, and panel method to calculate the aerodynamic loads on rotor blades. The results show the ability of the applied panel method to provide detailed information on the pressure and velocity distributions over the blade surface as compared with other applied approaches. A three-dimensional panel method was also used by Bermúdez et al. [5] for simulating the aerodynamic behaviour of horizontal-axis wind turbines, and the comparison between experimental data and the computed results with the panel method shows a good agreement. Gebhardt et al. [6] utilized the vortex-lattice method to simulate the unsteady aerodynamic behaviour of large horizontal-axis wind turbines. The aerodynamic of the blade-tower interaction has been satisfactorily captured as well as the effects of land surface and the boundary layer of the inflow. Kim et al. [7] have used the same method to simulate the blade-tower interaction over the NREL Phase VI and additionally used the nonlinear vortex correction method to investigate the rotor turbine while considering wind shear, yaw influence, distance from blade to tower, and different tower dimensions. The lifting lines model was used by Dumitrescu and Cardos [8] to simulate HAWTs, and the results obtained by using a nonlinear iterative prescribed wake analysis were more accurate than the results calculated by the free wake model.
More accurate simulations can be achieved by employing flow simulation methods such as Reynolds Averaged Navier-Stokes Equation (RANSE) solvers with an accompanying turbulence model. RANSE methods are commonly used for this purpose and are more accurate than potential flow methods but also more expensive in terms of computational time. Lee et al. [9] used RANS Equations solver in combination with the Spalart-Allmaras turbulence model to evaluate the performance of a blade with blunt aerofoil which is adapted at the blade’s root by increasing the blunt trailing-edge thickness to 1%, 5%, and 10% of the chord. The blunt trailing-edge blade helps to improve both the structure and the aerodynamic performance of the blades. Further, the SST turbulence model, see Menter [10], is widely used for wind turbine simulations due to its ability to simulate attached and lightly separated airfoil flows. This turbulence model is also used in Keerthana et al. [11] to calculate the aerodynamic characteristics of a 3 kW HAWT. Derakhshan et al. [12] compared the results obtained by Spalart-Allmaras, , and SST turbulence models for estimating the aerodynamic performance of wind turbine blades. The results show that the three turbulence models predict nearly the same power at low wind speeds, but the results obtained by the are more accurate at higher wind speeds.
The aim of the present work is to analyse the unsteady flow behaviour of a 5 MW HAWT utilizing two different numerical methods: a three-dimensional first-order panel method, which is the in-house panMARE code, and a RANSE method, which is the commercial ANSYS CFX solver. A general comparison between the results of both methods will highlight the viscous effects on HAWT flow that are not considered in panel code. Based on the computed results, the capabilities and limitations of the applied inviscid flow method to predict the complex WT flow will be illustrated. The panMARE code is successfully applied for simulating the flow of ship propellers and return reliable results for complicated flow conditions, such as in Bauer and Abdel-Maksoud [13] and Greve [14].
2. Wind Turbine Geometry
Due to the availability of its geometry and reference data [15], the NREL 5 MW offshore wind turbine is considered in the numerical study. It is a modern wind turbine with a complex geometry shape. While only the rotor is investigated in the first simulation Case I, the full geometry is considered in Case II, including the rotor, tower, and nacelle. The wind turbine rotor has three blades with a blade span of 61 m and attached to a hub with a radius of 2 m, which gives a total rotor radius () of 63 m. The blade is composed of several aerofoil shapes which are described at different radial blade sections in Table in [16]. The NREL 5 MW offshore wind turbine tower is typically circular with a height of 90 m and has a diameter of 6 m at the tower base and 3.87 m at the tower top.
3. Numerical Model for Panel Method
The in-house boundary element solver panMARE (panel code for MAritime Applications and REsearch) is applied in the simulations. The method is able to simulate potential flows for arbitrary steady and unsteady applications. The code is based on a three-dimensional first-order panel method, where the body’s surfaces are discretized by means of quadrilateral panels, where each panel has a constant-strength singularity distribution of sources and dipole. For each body’s surface panel, the solver’s governing equation and the boundary conditions are satisfied at a control point (in the centre of the panel). Potential theory drives the governing equation, where the flow is considered to be irrotational, incompressible, and inviscid. The equations for the conservation of mass and momentum are then simplified to Laplace’s equation for the total potential:The solution of the Laplace equation is a linear combination of several sources and dipoles superimposed with the inflow velocity potential (). The general solution of (1) on the body’s surface and the wake ( is defined in Katz and Ploklin [17] aswhere the vector () is the local normal to the surface (SB) and () is the distance between the panel control point and any other control point; () denotes the dipole strength and () source strength, which are defined asSince at an inner point of the boundary the inner potential can be chosen arbitrary, it will be an advantage to consider ; in this case the dipole strength and source strength can be defined as follows: , .
The velocity potential can be calculated by applying the following boundary conditions: (i)Direct boundary condition (Dirichlet): when (2) is considered in an inner point of the body, the overall potential will be defined as Because the inner potential is arbitrary at an inner point of the boundary, is chosen and the following equation for the velocity potential is obtained: (ii)Neumann boundary condition: the total velocity normal to the surface is required to be zero: and by applying the Neumann boundary condition on the lifting panel, the value of the source can be determined by where is the inflow velocity vector; so in the lifting body case, only the strength of the dipoles is unknown, and in the non-lifting body case, there are no dipoles and the strength of the source can be determined by (2).
The boundary conditions can be transformed into a set of linear equations, the solutions of which provide the required strength of each dipole () or source (). Moreover, the calculated and values are used to compute the induced velocities on each body surface panel. The pressure distribution ( on each panel of the body is calculated by applying the Bernoulli equation:where is the induced velocity, is the time, and and are the water density and the gravity constant, respectively.
The wake surface is a sheet or thin layer of dipoles. On the wake surface, it is necessary that the oscillation of the normal velocity component as well as the force difference between the lower and upper sides vanish.
The Kutta boundary condition is applied at the blade trailing edge in order to obtain the dipole strength of each first wake panel connected to the trailing edge which can be calculated by usingwhere and are the dipole strengths on the upper and lower sides of the aerofoil’s trailing edge. The geometry of the wake surfaces is iteratively adapted to be parallel to the local velocity vector, so that the normal velocity component to the wake surface becomes zero. In the steady case, the strength of the dipoles varies in the radial direction in order to meet the Kutta boundary condition, which remains constant in the axial direction. In contrast, the strength of the wake surface dipoles changes in radial and in axial direction in the unsteady case. In each time step, the Kutta boundary condition is applied and the dipole strength of the panel connected to the trailing edge is calculated. In the next time step, a new panel is generated behind the trailing edge, while the old panel, which was connected to the trailing edge, will keep its strength but will be transported further downstream according to the local flow velocity.
The BEM uses an iterative procedure to solve the rotor-tower interaction as shown in Figure 1(b). First, it solves the rotor and tower problems separately, after which the time-dependent influence of each component on the other components is considered by including the induced velocities of the adjacent structure, respectively.
(a)
(b)
The unsteady solution is initialized based on a steady calculation. The steady calculation is very important in predicting the initial wake panel’s primary location. This prior calculation is followed by unsteady calculation, in which at each time step the free wake moves with the induced velocity that is calculated at each panel. Each unsteady time step, or “explicit step,” involves few implicit steps (i.e., inner iterations).
In order to avoid unstable fluctuations on the tower panel’s properties that intersect with the rotor wake panels, all the wake panels inside the blade-tower interaction region are excluded from the calculations by using the split technique presented in Figure 1(a). Accordingly, all recognized wake panels’ dipole strengths will have a value of zero, so it will not have any effect on the tower panels. It should be mentioned that these dipoles which have already been saved before the collision will return to set again after passing the tower region.
The in-house CAD-code is used to build all geometry details and to generate the surface grid of the panels. The numerical grid on the blade is refined in the leading edge and trailing edge areas in order to capture the high pressure gradients. A smooth change of the panel size on the blade into the wake and on the blade tip is achieved to avoid panels with high aspect ratios. The blade grid in the radial direction is limited to to avoid high distorted panel shapes.
Two numerical studies have been conducted to investigate the dependency of the result on the blade grid and the simulated wake length. The results of the sensitivity studies are shown in Figure 2. The torque coefficient value is used to show the monotonic convergence toward a constant value, which is obtained by the finest mesh and longest wake length, where the numbers of panels (NOPs) in the blade’s radial and in circumferential directions are verified with a constant ratio between them. According to the results of the sensitivity study, the lowest panel number that can be applied is 1113 per blade ( panels in circumferential and radial directions, resp.), while a wake length of 1.5 revolutions is sufficient.
The total number of used panels in Case I for the rotor is 4000. In Case II 4050 panels are applied for the rotor and 450 for the tower, as shown in Figure 3. The investigated operation condition is a uniform wind speed of 11.4 m/sec with zero yaw angles and a rotor rotational speed of 1.267 rad/sec, which gives a tip speed ratio of 7. Therefore, one revolution is completed in 5 seconds. The Reynolds number at the radius is . The -axis is the rotation axis and the plane is the rotation plane. The time step in the simulation is chosen such that the rotor blades advance 6.118 deg. within a single time step.
4. Numerical Model for RANSE Simulation
The RANSE solver ANSYS CFX 14.5 is employed for calculating the viscous flow on the wind turbine. The governing equations are discretized using the Finite Volume Method, while the Shear Stress Transport (SST) turbulence model is applied. The SST model is a combination of for the outer region and for the near-wall region, which makes the best use of a good performance of near the wall while avoiding the sensitivity of this model to the free stream omega value; see Menter [10]. An improved near-wall formulation used in this model leads to a reduction of the near-wall grid resolution requirements, meaning that the solution is insensitive to the value. Thus, employing SST in complicated cases like wind turbines is convenient because of the difficulties in achieving a low value.
A tetrahedral mesh is generated using the ANSYS ICEM CFD. The applied hybrid tetrahedral mesh includes layers of prism elements near the wall boundary surfaces and tetrahedral elements in the interior region. The prism-meshing layer has elements perpendicular to the surface and makes it possible to keep the value nearly constant over the wall regions. The surface mesh type is a triangle with a maximum size element of 0.5 m2 and a maximum deviation of 0.01, where the use of such low deviation values allows for the subdivision of the element surface on the centroid of a triangle element in order to enhance mesh deformation around the surfaces. Figure 4 shows the surface mesh around the blade with a 2D section of the blade boundary-layer mesh. The resulting boundary-layer mesh for each blade consists of 15 layers with a first-layer thickness of about m and a growth ratio of 1.2. Tetrahedral elements are chosen for the flow volume domain so that the mesh refined in the inner region allows for a better flow resolution near the rotor.
Due to the different boundary conditions it is appropriate to divide the computation domain into two domains, which are stationary and rotating domains. The simulations domains considered in this study are shown in Figure 5; in Case I the rotor is housed in a cylindrical rotating domain that is enclosed by another cylindrical stationary flow domain. Interfaces have been defined between the stationary domain and the rotating domain to weakly impose the continuity of the kinematics and traction. Case II has a similar setup as Case I, but the stationary domain contains the nacelle and tower. Due to the near-wall refinement, the maximum calculated value is 200. Finally, the total numbers of elements in Case I and Case II are about 10.218 and 11.326 million cells, respectively.
The next step is to define the boundary conditions for the entire domain. First, a uniform wind speed and an air temperature of 25°C are set at the inlet boundary. At the outlet, an opening boundary condition is applied, which means the flow can exit or enter through this boundary surface. A static pressure value is specified at zero so that the pressure at the outlet is equal to the atmospheric pressure (standard pressure at sea level is 101,325 Pa). The opening boundary condition is also applied on the side boundaries as well as on the top and bottom boundaries of the computational domain. A no-slip boundary condition is applied on the wall of the blades, tower, and ground.
In order to solve the transient problem, the initial conditions must be specified in the stationary and rotating subdomains. For the stationary subdomain, a parallel inflow velocity is applied. For the rotating subdomain, the rotor's angular velocity is considered. The selected size of the time step is sec. The time step allows for an appropriate time resolution of the flow problem and for maintaining the courant number below 1.
The computation using the RANSE solver is carried out by employing 32 processors in an in-house cluster, with 2.2 GHz for each processor and 121.7 GB of memory. The CPU time is thus 96 and 104 hours for Case I and Case II, respectively. The BEM uses 4 processors with 3.3 GHz for each processor and 7 GB of memory. The CPU time is thus 12 hours for Case I and 24 hours for Case II.
5. Results
As mentioned above the flow on a 5 MW HAWT has been analysed using two different methods, a BEM and a RANSE method. The aerodynamics loads are calculated in the RANSE simulation by considering both pressure and viscous forces, where in the BEM simulation only the pressure force is considered. The comparison between the results of both applied numerical methods is given for the pressure distribution on the rotor blade, the thrust, and torque coefficients.
Figure 6 gives an overview on the calculated pressure on the blade surface. In general, a good agreement between the results of both methods can be noted, with the exception of results obtained for the region close to the blade root.
(a)
(b)
(c)
(d)
The first portion of the blade has a cylindrical similar shape which follows by a series of DU (Delft University) airfoils which twists with high angle of attack. The BEM underpredicts the drag in this region, whereas the results of the RANSE solver show a thick boundary layer and a severe flow separation due to the strong adverse pressure gradient, which leads to a high drag force. This region starts at the blade root and continues until . Compared with BEM, the applied RANSE solver is able to capture the most important viscous effects in this region of the rotor blade such as the high thickness of the boundary layer displacement and its influence on the pressure distribution.
In order to carry out a more detailed comparison, different results will be presented for four blade sections (, 0.63, 0.80, and 0.95). The blade radius is considered to be representative to aerodynamic rotor at the middle of the blade, while at section strong influence of the tip vortex must be expected; and the radius of is considered to be representative to the flow between the middle and the tip. The blade section at is close to blade root and has a pronounced separation region as mentioned before. Figure 7 shows pressure coefficient distributions for the four 2D blade sections against the normalized chord. The nondimensional pressure coefficient is calculated using the following formula:where is the inflow velocity, is the rotor angular velocity, and is the radius. The BEM results include the pressure distribution on the blade, which can be integrated to calculate the force vector acting on the blade surface. As is well known, the component of this force vector normal to the inflow velocity is the lift force, and in the direction of the inflow, it is the pressure-induced drag force. The BEM is able to calculate the pressure-induced drag force; however, the viscous component of the drag force cannot be estimated because important viscous effects such as boundary layer and flow separated flows are neglected. It should be mentioned that the viscous drag force can be added as an external force which can be calculated using an empirical formula for the friction coefficient. The friction coefficient can be calculated based on the Reynolds number built by the local velocity on the body surface panel. In the presented results of the study, no external drag force is taken into account.
(a)
(b)
(c)
(d)
As shown in Figure 7, the stagnation pressure value at the aerofoil tip is higher than the middle and root sections on both the suction side (bottom curve) and pressure side (upper curve), and its value equals 1. An acceptable agreement can be seen on the pressure side (upper curve) where the pressure is almost positive, while there are also differences on the suction side (bottom curve) and the trailing edge regions. At the suction side, the pressure until has the lowest value with the highest Reynolds number; after this region, the pressure increases from its minimum value to the value at the trailing edge. The value of at is almost zero at the trailing edge. The pressure at the trailing edge is influenced by the aerofoil thickness and shape near the trailing edge; if the aerofoil is thick, the pressure is slightly positive (the velocity is a bit less than the free stream velocity). But with sharp trailing-edge aerofoils, the value is close to zero, as in the presented aerofoils. Large positive values of at the trailing edge imply more severe adverse pressure gradients. Determining the pressure coefficient is a harder test for the simulation since it is a local quantity.
Figure 8 shows the streamlines over the same four blade sections. Only the flow at the blade section shows separation region near the leading edge, which may be a result of the high angle of attack (about 10°) at this blade section. The effect of boundary layer separation increases with an increase in the angle of attack. The results of the other blade sections show tiny separation areas near at the trailing edge.
(a)
(b)
(c)
(d)
It is evident that the presence of the tower has a strong effect on the inflow velocity and on the wake behaviour downstream, with a considerable reduction in wind speed in front of the tower. The axial distance between the rotor plane and the tower is generally kept as small as possible to limit the length of the nacelle. On one hand, a high nacelle length will increase the moments induced by the rotor blades on the tower. On the other hand, a small distance between the rotor and the tower increases the aerodynamic interaction between flow around the tower and the rotor blades. The reduced flow velocity to the rotor blades as they pass through the tower region has a significant impact on the performance of the whole unit. The reduced velocity in front of the tower can lead to a strong change in the effective aerodynamic angle of attack of the rotor blade. Both effects lead to a sudden change of the lift force on the rotor blade, which strongly affects the aerodynamic loading of the wind turbine. The actual 5 MW wind turbine uses blades with increasing tower clearance without a large rotor overhang. Therefore, precone and tilt angle are considered relative to the baseline wind turbine to reduce these negative effects [15].
Figures 9(a) and 9(b) show pressure coefficient distribution against a normalized chord predicted by both the BEM and CFX at span wise cross section 80% from the blade length on two rotor blades, which have an instantaneous angular position of 0° and 120° degrees relative to the tower. The maximum value of the pressure coefficient at 0° underlies a pressure drop of 10–15%.
(a)
(b)
The influence of the rotor on the pressure distribution of the tower is illustrated in Figures 10 and 11. Figures 10(a)–10(c) include the results for 0° relative angle between the tower and the blade. The corresponding results for 60° degrees are presented in Figures 11(a)–11(c). The calculated pressure distribution obtained by the CFX method can be seen in Figures 10(a) and 11(a) and the corresponding values by the BEM are presented in Figures 10(b) and 11(b). Moreover, Figures 10(c) and 11(c) contain a comparison of the stagnation pressure values on the tower leading edge. The results show that the pressure decreases on the tower leading edge when the blade passes the tower, where this sudden drop creates a definite impact on the tower and the blade loading during every blade passage, and has a great effect on the fatigue life of the blade and the tower. The differences between the BEM and RANSE method results in the region close to the tower top may be because the nacelle geometry is neglected and the wake treatment in the BEM is simplified.
(a)
(b)
(c)
(a)
(b)
(c)
The unsteady nature of the flow, due to unsteady interaction between the rotor and tower, is introduced for a period of 35 sec (7 revolutions) in order to ensure that the converged periodical behaviour of the results has been achieved.
Figure 12 shows the time history of the torque and thrust for Case I and Case II (with and without presence of the tower). The results of both solvers show that the rotor develops an aerodynamic torque, which corresponds to 5 MW, similar to the results presented in Bazilevs and Hsu [18]. As the calculated values of both methods give similar results for the torque, it can be assumed that the applied numerical grid and the employed boundary conditions as well as the used setup of the computations are accurate enough to characterise the temporal and spatial natures of the flow of the HAWT. The results presented for Case II show that the calculated thrust and torque values have a periodical behaviour due to the blade-tower interaction where every peak refers to one blade passage in front of tower. The mean thrust value and the amplitude of the thrust oscillation obtained by panMARE are slightly higher than the results calculated by CFX; however, the amplitude of the torque oscillation obtained by panMARE is lower. It is clear from this result that RANSE solver needs more than two revolutions for an accurate results predication.
Figure 13 shows the force time history on the tower in the flow direction (-direction). Both computational methods predict nearly the same amplitude and force fluctuation. It is also observed that the time history curve shows multipeaks for the force in both results. The wake effects on the tower force predicted by the BEM are stronger and the peak takes place directly after the blade passing. The magnitude of the force is reduced after passing the tower. The calculated tower force characteristics in the BEM may be strongly influenced by the applied method to treat the wake singularity. The tower force obtained by the RANSE solver shows different behaviour but the calculated amplitude of the force fluctuation by both methods is nearly the same.
Furthermore, the cylindrical shape of the tower means a large separated flow region takes place; the separated flow region and its influence on the tower forces cannot be predicted by the panel method so the devolution of the tower forces cannot be expected to show a similar tendency.
6. Conclusions
The flow on a 5 MW HAWT is simulated by using two different methods: BEM and a RANSE solver. The BEM is based on the potential theory, where the flow is assumed to be incompressible, irrotational, and inviscid. The RANSE method is the ANSYS CFX, which solves the RANS equations. In the presented study, the SST turbulence model is used. The BEM and RANSE methods have shown a strong potential for simulating the aerodynamics of wind turbines. The RANSE method is able to capture complex flow phenomenon such as separation and tip vortex flows as compared with the BEM; however, the application of the RANSE method is more costly. The BEM predicts the attached flow around the blade surfaces without the consideration of the boundary layer region, separation, and stall flow effects, which can be calculated by a RANSE solver. The comparison of the calculated local pressure distribution and the velocity streamlines around four blade sections by the two solvers shows an acceptable agreement on blade pressure side but small deviations on the blade suction side, especially near the leading edge area in the root region. The flow over the leading edge of the blade does not show any noticeable separation with an exception in the root region until , due to the cylindrical shape of the blade sections in the root region and the high angle of attack of its aerofoils. The results show a small separation region near the trailing edge of the sections near the blade tip. The calculated torque and thrust values by the BEM method are nearly similar to the results obtained by the RANSE solver. The calculated time histories of torque and thrust due to the blade-tower interaction by the two solvers are compared. The time histories obtained by the two solvers show different characteristics but the amplitude of the fluctuation is nearly the same. The computation using the RANSE solver is carried out in an in-house cluster and the CPU time is thus 96 and 104 hours for Case I and Case II, respectively. The BEM uses the personal computer and the CPU time is thus 12 hours for Case I and 24 hours for Case II.
The simulation results confirm that a BEM such as panMARE is able to capture the main flow properties around wind turbine blades, including the blade-tower interaction, with less cost and time than RANSE methods and that the BEM can be a powerful tool to support the design process of HAWTs.
Competing Interests
The authors declare that they have no competing interests.