Abstract
In the present work the actuator line (AL) technique is proposed to model a complete wind turbine (rotor and tower) in downwind configuration. The case study used is the Unsteady Aerodynamics Experiment (UAE) Phase VI turbine to analyze its suitability in capturing the tower shadow effect. Computational Fluid Dynamics (CFD) simulations were carried out using the open-source code OpenFOAM and the class of horizontal wind turbine AL from the toolbox SOWFA developed by the National Renewable Energy Laboratory (NREL). The objective of the present work is to validate the use of AL model for tower simulations and to establish whether such simplification is valid to capture the effect that tower shadow exerts over the rotor aerodynamics and to verify the accuracy of the proposed model regarding prediction of loads in the rotor. Computations were done first on an isolated rotor modelled with the AL technique; second on a full downwind turbine with AL in both rotor and tower; and, finally, on a 3D rotor with AL tower. Numerical results show good agreement with the experimental data concluding that the AL model is capable of a fast and accurate load prediction under certain conditions.
1. Introduction
The use of renewable energies to supply an increasing demand of energy over the world seems to be the trend in the near future. Alternative sources of energy such as wind and solar energy are the key to a sustainable way of progress. The rapid increment in the wind turbines size and the complexities associated with it are the reason of the growing need for highly reliable turbines to transform the energy from the source to everyday electricity. The size that wind turbines have achieved nowadays demands more accurate design models in order to avoid expensive maintenance cost, which in some locations and seasons is impossible to carry out.
Furthermore, the current design models are still subject of large uncertainties due to several phenomena produced by three-dimensional effects, unsteady effects, de-attached flow effects (stall), and tower effects, among others [1]. For this reason, the wind energy industry needs are clear: reliable, accurate, and affordable design models that improve the prediction of loads during the operation of the wind turbine on the field. This is just possible by understanding the physical effects behind the operation and states of these machines.
Efforts have been made to provide better and explanatory data about the effects a wind turbine is bearing on the field. In 2001 NREL invited several experts from different institutes to participate in a blind test of loads and performance prediction of an instrumented wind turbine that was tested under controlled conditions in the NASA-Ames wind tunnel [2]. The results computed by the participants did not agree and the predictions showed an important lack of accuracy. Moreover, variations among the different code applications of the participants were found [3].
Currently, the aerodynamic modelling of wind turbines ranged from the basic Blade Element and Momentum theory (BEM) to vortex wake methods to numerical simulations devoted to solving the Navier-Stokes equations. Although engineering models are currently the most common methods for load prediction, numerical simulations have become a very useful tool from which it is possible to get a deep insight and a physical realistic and explanatory visualization of the turbine flow field.
On the other hand, numerical simulations need a considerable amount of computational resources (memory and simulation time); the numerical issues associated and the complexity of the phenomena to simulate have led to developing simplified models based on actuators that capture the essence of the wind turbine behavior reducing the computational costs. Although these models do not describe in detail the underlying physics behind the operation; they are accurate enough when modelling a specific phenomenon. The actuator line, for example, can be used as a tool to predict operational parameters of the turbine, becoming a fast and easier way of design.
The actuator line (AL) model is a technique in which body forces are computed and distributed along a line that represents the actual blade geometry, capable of resolving tip and root vortices and of simulating the dynamics of the flow structure [4]; hence information about the wake structure and fluid dynamics can be obtained using simple geometries and grids.
Studies using the AL model have been often directed to analyze wakes and tip vortices for yawed and nonyawed turbines loads by using several turbulence models, RANS, URANS, and LES [4–8], for instance. This technique has been used to predict wind turbine power production and compared to actuator disk it has deployed good agreement; additionally, good practices of its use were attempted to be achieved [9]; numerical tests to determine mesh density and force distribution (important parameters in the method) guidelines using an infinite and finite span wing have been done [10]. All those studies have proven the accuracy of the AL model for wake analysis and its satisfactory performance [11, 12] as it resolves the wake velocity profiles.
However, guidelines regarding the use and the influence of the parameters associated with the model have not been found and a more detailed analysis is necessary. Furthermore, it is needed to analyze the tower and nacelle implementation and their effect on the turbine’s predicted power, forces, and wake profiles [12].
In the field of wind turbines, the possibilities of designing more flexible and lighter blades and of operating the turbine free yawed are some of the advantages that downwind turbines present [13]. However, while upwind has become the predominant configuration downwind turbines have been led apart due to, among others, the low frequency noise associated with the interaction of the blade while passing through the tower wake and the so-called tower shadow effect.
The current engineering models for tower shadow are based on semiempirical approach as wake width and amplitude must be known in advance; moreover, they consider the tower wake in an averaged sense that does not include the unsteady effects produced by the tower vortex shedding which is a possible cause of large fluctuation of the loads [13].
The presence of the tower causes an unsteady behavior of the flow that models have to take into account; for this reason, studies of tower interaction have been done using different techniques to handle the relative movement between rotor and tower, such as sliding mesh [14], dynamic mesh [15], and overset or chimera grids. Moreover, potential flow models have attempted to model the tower and the flow around it [16] and to couple the solution with Navier-Stokes solvers to simplify the simulations [17].
In the present paper two main contributions are considered: first is the implementation of AL to model the tower of a wind turbine; second is a numerical study of the influence of two AL model parameters (cell size in the rotor region and ε) on the prediction of the global performance of the turbine. The implementation of the AL model to the tower not only includes the influence of the drag force but also the fluctuation of this force in time. Including this effect in low order simulation of wind turbines in downwind configurations is extremely important in order to correctly capture the complex interaction between the tower wake and the rotor dynamics (tower shadow). For this purpose, unsteady computations of the flow over the NREL Phase VI turbine model with the open-source code OpenFOAM are presented and results of the full turbine (rotor and tower) modelled with AL are shown. First, simulations of an isolated rotor and the parametric study of the AL model to validate the predictive capability and the effects of the numerical parameters of AL in the numerical results are presented. Second, the implementation of the AL model to simulate a tower is performed and compared with the flow around a cylinder simulation. Lastly, the interaction between rotor and tower is presented as AL represents both, tower and rotor, in downwind configuration. Additionally, three-dimensional rotor (including all the geometric details) simulations are performed to validate the AL tower performance and comparison of the aerodynamic forces and power variation with experimental data is done aiming to conclude whether the unsteady effect of tower shadow is successfully captured through the AL technique and its implications on the turbine behavior.
2. Material and Methods
2.1. Experimental Data (Base Case)
The turbine of the Unsteady Aerodynamics Experiment (UAE) Phase VI conducted by NREL was selected as base case. Table 1 shows the most important dimensions of this wind turbine model.
The purpose of this experiment was to acquire accurate quantitative aerodynamic and structural data representing a wind turbine in field operation and generate a database [2] devoted to studying different phenomena to get a better understanding of the turbine behavior and to improve the models to predict the aerodynamic forces.
Several measurement campaigns were carried out to test the wind turbine under different operational conditions. The turbine used was a two bladed one of 5 m span with an airfoil S809, in both upwind (sequence H) and downwind configurations (sequence B). For those sequences, the turbine was yaw locked at 0° and the blade tip pitch was 3°.
Data was collected for wind speeds ranging from 5 m/s to 25 m/s and at a rotational speed of 72 RPM. In this study, simulations were performed at wind speeds between 5 and 8 m/s; for these speeds the turbines are under laminar regime, while from 9 m/s the turbine started to present stall delay, an effect whose analysis is not within the scope of this work. In the present work, results and comparisons with experimental data were performed at 7 m/s because at this velocity the turbine presented the maximum power before stall.
Data of pressure along the blade was collected for 30-second campaigns. Pressure taps were placed along the suction and pressure side of one blade at five spanwise locations at = 0.30, 0.47, 0.63, 0.80, and 0.95 was presented in Table 2.
The sequences used in this study were the B and H, downwind and upwind, respectively. For the sequence B the turbine has a cone angle of 3.4° and yaw angles of ±180° and blade tip pitch was 3°. In this work, the yaw angle was kept 0° and the turbine was modelled without teeter and coning since the combination of them results in blades being parallel to the tower axis during the tower passage [15]. For the sequence H, the turbine was rigid with a 0° cone angle.
Blade pressure measures were collected for both sequences and aerodynamic forces computed. For a comprehensive description of the turbine and experiment see [2].
2.2. Computational Set-Up
The open-source code OpenFOAM was employed for all computations as well as the class implemented for horizontal wind turbines with AL model included in the Simulator fOr Wind Farm Applications (SOWFA) developed by NREL. The turbine was simplified using the AL model as tower, the nacelle was omitted, and the rotor was assumed to be completely rigid and fluid-structure interaction was not implemented.
Grids were generated with the meshing tool provided by OpenFOAM that creates structured mesh of hexahedra cells. The domain for the simulations of the isolated AL rotor is shown in Figure 1, large enough to avoid boundary interference with a rectangular shape of 10D (D: rotor diameters) length and 5D width and height; the rotor was placed at 1/3 of the domain length from the inlet. For the case of three-dimensional rotor, the mesh was generated with the in-house automatic Blade Block Mesher BBM (Developed by Fraunhofer IWES and ForWind, Oldenburg), a tool based on OpenFOAM’s blockMesh and snappyHexMesh. The resulting domain consisted of 11.6 million cells, was of cylindrical shape, and had 10D length and a radius of 3D. The rotor was also placed at 1/3 of the domain length from the inlet.
(a)
(b)
The Unsteady Reynolds Averaged-URANS computations were carried out using the PIMPLE algorithm to solve the coupled velocity pressure equations and the turbulence was modelled with the k-ω Shear Stress Transport model [18]. The time step was fixed keeping the condition that the tip of the actuator line does not cross more than one cell per time step which according to [9, 10] is a more restrictive condition than CFL = 1. This led to obtaining CFL number lower than 0.2 with a time step of s; this agrees with computations reported previously [10]. For the three-dimensional rotor simulations, a time step of s is proven to keep stable simulations. The time discretization was done by the second-order implicit backward scheme.
In the case of the AL model the spatial discretization was carried out by Gauss linear interpolation schemes while for three-dimensional case Gauss second-order upwind was employed.
Boundary conditions are listed in Table 3, where is the turbulence intensity, is turbulence kinetic energy, ω is turbulent frequency, is turbulence length scale (l = 0.07D), and is an empirical nondimensional constant equal to 0.09. The turbulence intensity at the inflow was set to be 0.7% as reported for the NASA-Ames wind tunnel [2].
2.3. Actuator Line Model
The AL model, proposed by Sorensen and Shen [19] and implemented in the OpenFOAM-based class for the study of wind farms by NREL [20], is a technique where the actual geometry of the blade is replaced by a virtual line over which the body forces are included as a source term to the Navier-Stokes equation.
The body forces acting on the rotor blades are computed on an element of differential size by using the Blade Element approach as it is shown in (2), where and are the lift and drag coefficients respectively, is the chord length, and and denote the unit vectors in the directions of lift and drag.
The local velocity relative to the rotating blade is given by (3), where and are the tangential and axial velocity in the inertial frame of reference.
To avoid singular behavior is formed by taking the convolution of the computed force and a regularization kernel η which takes the form of a Gaussian distribution where is a constant that adjusts the strength of the regularization function and is the distance between the measured point within the computational domain and the initial force points distributed on the rotor.
The main advantage of this model is the need of few grid points to capture the aerodynamics and fluid dynamics of the blades. Moreover, compared with other simplified methods (actuator disk), actuator line provides information about the dynamics of the flow and allows a more detailed study of wake structures and tip and root vortices [9, 12]. However, boundary layer phenomena like flow separation and dynamic stall cannot be properly captured by the model since there is no surface and the boundary layer over the surface of the blades is not resolved.
The parameters of the actuator line model were set according to SOWFA developers in order to ensure a smooth force projection over the actuator line [20]. There are two main parameters in this model, the grid resolution and the kernel regulator ε. Regarding grid resolution, the cell sizes through the actuator line-tower were set such that there were at least 50 grid points; this value ensures that the forces are distributed uniformly over the mesh. Regarding the regularization parameter ε that depends on cell size it was established that for a smooth projection of the forces the blade/tower width must be related and the ratio between blade’s width and cell sizes must be approximate to 0.75 as suggested. The good performance of this set-up was proven in the results as they showed good agreement with experimental data [9].
3. Results and Discussion
3.1. Isolated Rotor Simulations with AL
The base case at a wind speed of 7 m/s was selected and implemented following the guidelines from SOWFA developers [20]. The positive direction of the axes , , and was oriented according to reference [21]. The cell size was such that across the rotor diameter there are at least 50 grid points giving a cell size = 0.2; then the regularization parameter ; hence . The actuator points should be enough to predict a smooth force projection; therefore they are set as [9], where is the width of the blade.
To evaluate grid independence, results from the base case and three more refined meshes were compared. Characteristics of the meshes are presented in (Table 4) and the refinement boxes are shown in the schematic view of Figure 1 (from this point on, experimental data will be represented on figures by dotted line).
Comparison of the output power with experimental results at different wind speed can be seen in Figure 2. Convergence is achieved from a mesh of 7 million cells and a cell size of 0.15 m.
The regularization parameter was set such that the forces are smeared over an area that represents the actual blade geometry. Therefore, it should be related to a blade characteristic length; this has been suggested in [10, 22]. In the present work, it is proposed to choose the regularization parameter equal to the chord of the blade; ideally, it should be equal to the local chord; however, this is not possible since the model and the computational tool consider it constant along the blade span. The influence of has been studied for an actuator disk [23] and for actuator line [4, 8] with special attention to the flow and wake dynamics. In this study, this parameter was chosen equal to the blade chord on the 5 sections in which the experiment took place.
Figure 3 shows good agreement between the normal force coefficients and the experimental data varying ; moreover the results fit the standard deviation.
(a)
(b)
Something different is shown for the tangential force coefficient (see Figure 3(a)). Differences are found especially on the mid-section of the blade for higher values of . For both cases, best results are deployed by = 0.38 which is the chord of the tip.
3.2. AL Model for the Tower and Full Wind Turbine Simulation (Rotor and Tower) with AL Model
Regarding simulation with AL tower, the velocity field around the actual tower was obtained by numerical simulations of flow around a 2D cylinder. A cylinder of 0.4 m diameter at wind speed between 5 and 7 m/s was simulated. For these conditions, the Reynolds number ranges between and which is within the so-called critical regimen characterized by a sudden drop of the drag coefficient and a disorganized wake behind the tower. Significant differences were found in literature, as the drag coefficient value ranges from ~0.3 to 1.2 and the starting dropping point is in the range of Reynolds number between and .
CFD simulations of the flow around a cylinder at wind speed of 7 m/s were used in order to estimate forces and propose an AL model for the tower. Based on the averaged drag and lift coefficient an AL model for the tower was implemented. Numerical results were compared with those reported in [24] and a disagreement was found in the velocity field and lift and drag forces. This discrepancy was related to the fact that the proposed model used constant values of lift and drag coefficients for computing forces and this assumption was not producing the unsteady behavior of the flow around the actual cylinder.
To improve the proposed AL model for the tower, the source code was modified and the lift and drag coefficients were included as function of time as shown in
The result of this modification showed the same oscillation of lift and drag forces obtained with the CFD simulation of the flow around the cylinder.
Comparing the deficit of velocity behind the AL model proposed for the tower with the expected deficit (CFD of the cylinder), the modified AL tower improves the results regarding width and depth as it is presented in Figure 4. This improvement in those parameters is important as the width of the deficit is the space the blade will pass through and will influence the rotor aerodynamics.
Finally, the turbine was completely simulated using the AL model aiming to simplify the simulation and to verify the capability of the model to represent the effect that the tower exerts over the rotor aerodynamics. An FFT analysis was performed for the evolution of the predicted output power of the turbine for both AL and experimental data. Experimental measurements were sampled at 520.83 Hz over 30 seconds.
The power spectrum in Figure 5 shows important content up to 30 Hz and spectral peaks of significant magnitude in the range 0 Hz to 10 Hz. Comparing the actuator line model and the experimental data at wind speed of 7 m/s, it is clear that the effect of the presence of the AL tower as the highest spectral peaks is found for 2P, 4P, and 6P meaning that the power drops as the blades encounter the tower shadow once per revolution ( = 1.2 Hz); for the experimental data, although the highest peak is at 2P, there are spectral peaks every 1P due to unbalance of the turbine during the experiment.
Normal and tangential force coefficients were compared with experimental data; Figure 6 shows the averaged coefficients for 5 spanwise positions along the blade corresponding to the pressure tap locations in the experiment. Good agreement is found proving that the rotor aerodynamics is well represented by the AL turbine. For the root and the tip of the blade bigger discrepancies were found due to three-dimensional effects present in those locations and due to the fact that they are not properly captured by the model. Although at the root and tip the Glauert correction was applied, it is not enough to accurately predict the forces.
(a)
(b)
Figure 7 shows the evolution of the normal force coefficient for one blade at different spanwise positions. It is clear that good agreement is found up to 80% of the span length of the blade. Close to the tip of the blade (95% of the span) the model does not represent properly the tower shadow effect. In general, it is clear that the AL model for the tower captures the average behavior of the wake and its interaction with the rotor.
(a)
(b)
(c)
(d)
(e)
3.3. Detailed Rotor Simulation with AL Model of the Tower
Rotor tower interaction was also studied with simulations performed by using a three-dimensional model of the rotor and keeping the tower as an AL. Two approaches to handle the relative movement between rotor and tower while keeping a low computational cost were tested: Multiple Reference Frame (MRF) and dynamic mesh (Arbitrary Mesh Interface-AMI).
Multiple Reference Frame computes the flow using both the rotational and stationary reference frame without moving the mesh. In the rotating zone the Coriolis force term is added to the governing equations (see (1)) as a source and the flux is calculated from the relative velocity. In this way, a virtual rotation of the flux is imposed without moving the mesh which is an advantage in terms of computational cost. Since the rotor is not rotating but the flux is, the AL tower is rotating at the same rotational speed but in opposite direction to the MRF Zone modelling the relative movement. To do this it was necessary to modify the solver to compute the AL body forces from the velocity in the inertial frame of reference and the flux around the rotor from the relative velocity.
The results presented in Figure 8 show the presence of the tower by a drop of power every 2P as it is the case for two bladed turbines. However, a high increment of power is shown which follows the drop that is not related to the tower shadow and requires further investigation of the implementation of this technique.
In the dynamic mesh approach, the mesh is rotating; however the AL remains in the same position; therefore, the tower can be placed within and outside the rotating domain and the relative movement between rotor and tower is well represented. Such rotating domain can be as big as the computational domain leading to reduced computational cost avoiding the communication with an AMI interface between the rotor and the rest of the domain.
Unsteady RANS simulations were performed both with and without tower with the AL tower placed in and out of the rotating part. The AL tower presence is proved by the drop of power in both cases. The difference between the isolated three-dimensional rotor and the full downwind turbine is presented in Figure 9 showing the effect of the presence of the tower as the power drops every 0.41 seconds when the blades pass through the tower wake.
4. Conclusions
The use of the AL technique to simulate the full wind turbine and analyze the capability of such model to capture the tower effect on the rotor aerodynamics was the motivation of this work. To do this, numerical simulations of the NREL Phase VI turbine were performed. First, we presented the isolated rotor where a parametric study was performed. Second, a full turbine simulation was presented, as AL was employed to model the tower in downwind configuration.
The simulations were performed with OpenFOAM using the actuator line technique implemented in SOWFA. Computations of the rotor were done for a wind speed of 7 m/s at a fixed yaw and pitch angle of 0° and 3°, respectively, with a constant rotational speed of 72 RPM.
Three aspects were studied: the accuracy of the actuator line model in predicting loads, the effect that the parameters of the model have on the results, and the performance of the model to be used as tower when simulating a full turbine.
Results from the AL model are strongly dependent on the parameters and the mesh resolution. The guidelines found in the literature regarding mesh characteristics are not completely accurate and grid independence studies have to be performed for each case. Other works have a set related to the grid resolution. However, this parameter has a physical meaning because it establishes the fact that the distances over the forces are projected on the flow.
This work proposed to choose such factor from the chord length of the tip of the blade. Although for the output power the results have shown the solution start to converge after a regularization parameter equal to the chord at 63% regarding forces prediction, this value has shown slight discrepancies with experimental results. And as the tip speed is mostly responsible for the power of the turbine it is recommendable to use a regularization parameter related to the chord of the outermost part of the blade at 80% or 95%.
Regarding the use of the actuator line model as tower, the capability of the technique to influence the flow and affect the aerodynamic of the rotor was demonstrated. For full AL turbine simulations, good agreement between numerical simulations and experimental data was found, especially in mid sections (63% and 80%). The reason for this is that the flow on those sections is mostly two-dimensional while in the inner and outer part of the blade the forces are subject of three-dimensional effects that the AL model does not capture completely.
Since URANS simulations were performed, fluctuations of the vortex structures from the tower were not evident; hence the unsteady nature of the interaction between the vortex shedding from the tower and the blade were not observed. To obtain such effect, Large Eddy Simulations with the use of the actuator line technique could deploy better performance.
The capability of the AL to model the tower was proved as well with three-dimensional rotor simulations, as the power drops every time the blades pass behind the tower. This was done with two different techniques to handle the relative movement between the tower and the rotor. The dynamic mesh was the more suitable technique as consistent results were obtained. Although simulations with AL full turbine needed less computational resources, it is not possible to obtain information about the effect of the rotor on the tower structure, as well as effects related to the blade boundary layer, and the limitations regarding stall delay and tip loss make the simplification useful just under certain conditions.
Furthermore, the use of the AL as tower with three-dimensional rotor is a good compromise to get insight of the rotor aerodynamics while simulating full turbines and the effects associated with the tower interference.
In order to compare the proposed method with other computational/engineering low order numerical techniques, two research works, found in the literature, are presented and discussed. Rahimi et al. [25] modelled the UAE phase IV turbine in downwind configuration using the Blade Element Momentum (BEM) model for the rotor and a potential flow solution to model the tower. The simulations were performed with the FAST code also developed by NREL. A full detailed CFD simulation was also performed for comparison purposes. Numerical results show that FAST predictions in velocity deficit are 20% larger than the mean values. Larger discrepancies (between BEM and CFD) were also observed closer to the tip of the blade. The increment in the normal force close to the azimuth angle of 170° was not captured by neither of both approaches, BEM and CFD. From Figure 7, it is clear that, in the actuator line model, this increase in the force coefficients is slightly captured close to the root of the blade but it is not captured in those sections closer to the tip.
The second research work found in the literature deals with the simulation of the UAE phase IV turbine in upwind configuration using the vortex lattice method (VLM). Hogeon et al. [26] show numerical results that are in good agreement with experimental data, but the effect of the tower is not captured since they report no significant disparity whether the tower is present or not. It is important to mention that a model such as VLM requires a special treatment of the wake as it interacts with the tower. In the methodology proposed in the present paper, no treatment of the wake is required since this is resolved in the computational domain by the Navier-Stokes solver.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
Funding for this work was obtained from the Turbine Simulation, Software Development and Aerodynamics Group of Fraunhofer IWES Institute. Many thanks are due to the Turbulence, Stochastic and Wind Energy Group of ForWind Institute, especially the Computational Fluid Dynamics Simulations team, for providing all resources needed and for the productive discussions and support from colleagues. Thanks are due to the Facility for Large-scale Computations in Wind Energy Research (FLOW) at University of Oldenburg for providing computer time. Many thanks are due to Dr. Scott Schreck from the National Renewable Energy Laboratory who provided the experimental data.