Abstract
Downwind wind turbines have lower upwind rotor misalignment, and thus lower turning moment and self-steered advantage over the upwind configuration. In this paper, numerical simulation to the downwind turbine is conducted to investigate the interaction between the tower and the blade during the intrinsic passage of the rotor in the wake of the tower. The moving rotor has been accounted for via ALE formulation of the incompressible, unsteady, turbulent Navier-Stokes equations. The localized , , and are computed and compared to undisturbed flow evaluated by Panel method. The time history of the , aerodynamic forces ( and ), as well as moments were evaluated for three cross-sectional tower; asymmetrical airfoil (NACA0012) having four times the rotor's chord length, and two circular cross-sections having four and two chords lengths of the rotor's chord. 5%, 17%, and 57% reductions of the aerodynamic lift forces during the blade passage in the wake of the symmetrical airfoil tower, small circular cross-section tower and large circular cross-section tower were observed, respectively. The pronounced reduction, however, is confined to a short time/distance of three rotor chords. A net forward impulsive force is also observed on the tower due to the high speed rotor motion.
1. Introduction
In the quest of finding a clean, renewable, and sustainable energy resource, wind turbines stand “tall” in comparison to solar technologies, that is, troughs and heliostats. It has been reported that the deployment and implementation of wind turbines has increased between 20 to 30% annually throughout the last decade [1]. It is estimated that for offshore turbines, the carbon dioxide emitted indirectly due to material fabrication is paid back within nine months of operation [2]. Essentially, there are two types of wind turbines, classified by the direction in which the rotor rotates with respect to the ground. The older generation turbines, of vertical axis configuration, (i.e., anemometer cup, Savonius, and Darrieus rotor) are characterized by a lower power output. Horizontal axis turbines, such as Vestas series v15, v27, v39, and v66, Windmaster series 500, 750, and 1000, and Nordex N500, are gaining worldwide commercial deployment due to their higher power output. Horizontal axis turbines can either be of downwind or upwind configuration. Both are subjected to rotor tilt aiming to increase the tower rotor clearance. This creates an additional yaw misalignment to the inclined wind. The extended nacelle required for upwind turbine configurations creates a disadvantage since the rotor needs to be positioned far enough away from the tower to avoid blade-tower strike. Therefore blade tailoring is more restrictive and so a stiffer and more expensive blade needs to be employed. However, this induces high stress at the rotor hub during gusty wind conditions, a problem that is not associated with downwind turbines. However, downwind turbines have lower tower shading. The flexing of the blade is an interplay parameter between fatigue and air tailoring particularly with the engineering of new and more reliable blade materials. The higher negative tilt of the downwind rotor creates a larger tower clearance and reduces its shading. Furthermore, for upflow conditions (flow over a hill) of 1–11 degrees, Yoshida [3] have experimentally reported that the misalignment could be reduced to a level below that of the upwind configuration and that the power coefficient values could be improved by 5%–20%. The same authors conducted a simplified CFD simulation and suggested that the energy production could increase by 7%. Therefore, the downwind rotor could be reevaluated as a promising configuration for the future wind turbine market. There is, at present, limited high fidelity simulation dedicated to the downwind turbine.
The flow through a wind turbine can be analyzed through the use of a number of methods that vary in accuracy and complexity, from the simple Blade Element Momentum method (BEM), to medium complexity methods such as the vortex and Panel methods, through to CFD analysis of the wind turbine rotor. For the BEM method, the rotor is dissected into several radial sections (10 to 12) from which the rotor shaft torque is computed [4]. Although it is a fast and reasonably accurate method, it is limited to the availability of adequate aerodynamic (lift and drag) airfoil data [5] at the corresponding angle of attack () and Reynolds number (Re). The Panel method is an inviscid solution to the velocity potential formed from sources (representing the body surface) and dipoles (simulating the lift) that are distributed to form the shape of the turbine airfoil and its potential [6]. Therefore, the Panel method is an intermediate step that can be used for BEM when airfoil data is absent. Modern wind turbines are dramatically larger than their former equivalents and now have different blade, airfoil twists, and 3D features. Simplified static analysis is only limited to small rotor diameters of less than 25 m (200–250 kW rated power) according to Danish standard DS472 [7] and accurate aerodynamic flow simulation is required to account for their aeroelasticity and their dynamic behavior.
The CFD analysis of a wind turbine involves obtaining a digitized description of the computational domain of the 2D or 3D airfoil geometry in order to carry out a high spatial and temporal resolution simulation. Commercial preprocessors including Gambit, ICEM, Proam, and CAD tools (Catia, Unigraphics, Solidworks, and Autocad) could be used to generate the wind turbine rotor geometry and the computational mesh. It is possible for these preprocessors to perform conceptual analyses and optimal designs with access to only a small amount of preliminary data. This is beneficial since there is limited data that accounts for the twist of airfoil baseline positions with respect to the nacelle or tower. Considering the large wind turbine domain size, parallel implementation of the CFD code is used to relax the mesh/domain size constraint, resulting in a finer geometry description and incorporation of other geometry components. Complete NS solvers with turbulent implementation are available commercially (i.e., Fluent, Star-cd, Numeca, Openfoam) and have been used for wind turbine studies [8]. The rotation of the turbine can be accounted for with features such as the moving reference frame, sliding/moving, or mesh fluxes. These codes are subjected to proper discretization (central for diffusion/viscous, and 2nd or 3rd order for time and advection), and implicit solvers to avoid restrictive explicit time stepping.
Wind turbines are subjected to flows that have a low Mach number, yet are turbulent, with Reynolds number from 1–6 million. Flow near the root of the blade can have Mach numbers as low as 0.01. The incompressible NS equations solvers are the ideal choice of solver, since they treat the pressure as the primary variable. Alternatively, the well-developed aerospace compressible NS solvers can be used with preconditioning or with artificial compressibility for transient applications which require further subiteration to enforce incompressibility [9]. The advantage of the compressible solvers is their flexibility in using overlapping grid and high-order upwind schemes when dealing with large and complex geometries or when incorporating several components, that is, nacelle, rotor, tower, and piedmont.
The flow through a wind turbine is turbulent and typically has a domain that spans 100 m with (Re) greater than 106. Kolmogorov eddies are the smallest turbulent scales and are expressed as , where L is the characteristic length, that is, the tower cross-sectional diameter or rotor chord length which varies between 1-2 m. Therefore, the number of grid points required to directly resolve the three-dimensional domain are enormous and impractical for current computational resources ( for two dimensional). Only the time-averaged Navier-stokes system, denoted as Reynolds-Averaged Navier-Stokes (RANS), reduces this computation and renders it feasible. The Spalart-Allmaras model is a one transport equation ( turbulent viscosity) turbulent model which is widely used for aerodynamic applications. The resulting Reynolds stresses initially are modeled with an eddy viscosity model allowing the summation of the diffusion/molecular viscosity term to an equivalent viscosity () [10]. In the Spalart-Allmaras model the transport equation essentially constitutes the closure of the system. There are also other eddy viscosity models including , Baldwin and Lowmax [11], [12], and Reynolds Stress. Nonetheless, RANS models have shown good implementation for wind turbine studies. Hansen [13] has indicated their ability to model the stalled flow regime at high wind speeds. The turbulence at this point becomes less isotropic, and this forms the basis of such modeling. An alternative to Reynolds averaging is a method which filters the NS with Kernel (typically the mesh size is implicitly kernel). This filtering technique includes both LES and DES which are, however, greatly computationally expensive as they require a finer grid and are inherently transient.
As for the application of CFD to wind turbine flow problems, most of the previous work has focused upon a zero-yaw rotor configuration (wind perpendicular to the rotor plane). Steady flow without shear is assumed and it is only the flow over the rotor that has been the area of concern. The flow over the nacelle and tower has been neglected. The interaction between the rotor and the tower in transient flow has, to date, not been considered. The objective of this work is to numerically calculate the aerodynamic forces exerted upon the rotor and the tower of a downwind horizontal axes turbine in operation.
Three tower cross-sections are considered, one is a symmetrical airfoil, and the other two are circles of different diameters. The moving/rotating rotor is accounted for through the use of a sliding mesh with an ALE formulation. A second-order discretization scheme, unsteady and averaged Navier Stokes flow is employed for the simulations. The moving mesh configuration brings valuable insight to the pertinent design parameters: the variation of the tower cross-sectional shape, wind direction, and the consequence of the tower/rotor blade interaction. The results of these analyses will be of great use to the wind turbine designer in the evaluation of the wake interaction, which is convected downstream, and which causes flow unsteadiness that is considered responsible for low-frequency noise generation. The focus of this work is to conduct a detailed two-dimensional flow analysis to quantify the resulting pressure coefficients, aerodynamic lift and drag forces on the tower and rotor. Virtual wind tunnel data will be used as a substitute for the airfoil data essential for the BEM method and will represent the first CFD step, amongst several, required to conduct a complete CFD simulation.
2. CFD Flow Simulation
A detailed downwind turbine has been developed in a unigraphics CAD environment. This geometry is imported into CFD preprocessor as an IGES file as shown in Figure 1. A cross-section domain is established at that cuts through the tower and the rotor blade cross-sections. A CFD domain and mesh, as depicted in Figure 1, are also constructed.
The prescribed rotor motion allows for the domain to be divided into moving and fixed subdomains. One subdomain, that contains the rotor, lies in between the fixed upstream and downstream subdomains. The rotor domain is modeled as a sliding mesh against the two fixed up- and downstream subdomains. This approach has a lower associated computational penalty and higher accuracy than a remesh of the domain, as in the case of nonprescribed motion where the flow is coupled with three degrees of freedom solver and involves exhaustive solution interpolation. Rotor blades and tower section static pressure and aerodynamic forces will be computed at different incident wind angles and wind speeds simultaneously.
2.1. Conservative Law Equations
The wind flow is considered to be transient, incompressible, and viscous Navier-Stokes as well as turbulent. The one transport eddy viscosity equation model, Spalart Allmaras is utilized for turbulence closure in conjunction with the wall function to the near-wall region. The flow is governed by the two-dimensional, transient, incompressible Navier-Stokes equations, written as (1)–(3) represent statements of mass continuity and momentum conservation and are written after expanding the scalar variable expansion and ensemble (over bar) averaging as The term is the Reynolds stresses and is modeled utilizing the mean () velocity via the common eddy turbulent viscosity. It is expressed as where is the turbulent viscosity. Substituting (5) in (3) conveniently allows for the summation of the Reynolds stresses terms to the diffusion term (2nd right hand term in (3)) with an equivalent viscosity. Therefore, closure of the above system is achieved with the integration of the additional transport equations forthat has the following form: where and are the production and the destruction of turbulence viscosity, respectively. and c are constants and is the user-defined term.
2.2. Computational Domain and ALE Representation
Two configurations that will simulate the passage of the rotor in the wake of a symmetrical airfoil and circular tower sections are considered. Figure 1 depicts the baseline geometry of the two-bladed, 85 m hub, 70 m blades diameter, downwind, 1MW-rated power wind turbine. The two-dimensional geometry encapsulates both the wind turbine tower cross-section and the rotor airfoil, and can evaluate the forces and pressure distributions at different locations of the rotor during its passage in the wake of the tower. To do this, the continuity and momentum equations are written in the Arbitrary Lagrangian Eulerian frame following the work of Bathe et al. [14, 15]. Thus, writing (1) and (2) in their conservative form such that where is the flow velocity vector and is the grid velocity that is equal to the translating rotor velocity. D is the diffusion coefficient and is the source term. A sliding mesh is implemented by creating an interior zone across the sliding interface that consists of an overlapped zone and either a periodic or a wall region. Referring to Figure 2, the b-d, d-e, e-f, f-g, and g-h faces are the overlapping faces.
The remaining two faces, a-b, and h-i, are the wall or periodic regions. The concept of interior region is manifested by allowing for b-d, d-e, e-f, f-g, and g-h faces instead of the cell faces of 1, 2, and 3 to be used for the computation of the flux to cell 4, 5, and 6 from cell 1, 2, and 3, respectively. A pressure-based PISO second-order implicit solver is used to solve the conservative finite volume momentum and continuity equations. To construct the computational domain, a horizontal section is taken at of the rotor as shown in Figure 1. The domain consists of a symmetrical airfoil that represents the tower NACA0012 with 3 m chord length and a rotor airfoil blade NACA6312 with 0.72 m chord length. A multiblocked, wall bounded boundary layer with quad cells is constructed for the tower and rotor surfaces. The overall domain is comprised of 128,000 cells. The computational domain and boundary conditions are presented in Figure 3.
3. Results and Discussion
3.1. Mesh Sensitivity
The accuracy of the simulation as a whole, that is, the setup, discretization, BC dependency, and the domain size, is assessed by the accuracy of the numerically evaluated values of the lift coefficient and Strohaul number for the flow around the circular cross-section that the model obtains. Fluent Inc. solver [16] is used for the analysis and the unsteady, implicit, noniterative solver with PISO algorithm at double precision is employed. A second-order implicit model for the time and a second-order upwind model for the pressure (continuity) and momentum are used at a tight tolerance () and relaxation factor ( where is the provisional solution).
Three meshes at different refinement levels and body-fitted type are constructed for the tower cross-section. They are composed of multiple blocks and quadratic cell types at an initial cell height/wall distance of 0.001D from the cylinder for the finest mesh. These meshes are referred to as coarse, medium (at double resolution of the coarse), and fine (at one and half resolution of the medium). The mesh sensitivity is conducted at a laminar Reynolds number of 100 as this solution is well documented in literature. Time history of the lift and drag coefficients show their anticipated oscillatory behavior due to vortex shedding and are depicted in Figure 4. The lift and drag spectra presented in the figure show that the drag coefficient has an associated frequency twice that is displayed by the lift coefficient time histories.
(a)
(b)
(c)
Figure 5 depicts the vorticity shedding of the cylinder at three instants of time. The mean values for the three refined mesh levels, and the relative errors for the coarse mesh based on reported experiments [17, 18] are summarized in Table 1.
The comparison also includes the Strouhal number (, where f is the shedding frequency, U is the upstream undisturbed inlet velocity, and D is the cylinder diameter). Similarly, the relative error of the medium mesh prediction for the drag coefficient, lift coefficient, and Strouhal number are 7.5%, 0%, and 1.8%, respectively. It is 14.3 and 4.3 percentage points below the coarse prediction and 1.7 and 0.6 percentage points above the fine prediction for the drag coefficient and Strouhal number, respectively. Accordingly, the medium mesh produced values with reasonable errors and can therefore be used for subsequent analyses. It is also noted that the extent of the domain is considered appropriate since the velocity gradients are vanishing upon reaching the end of the domain.
3.2. Choice of Rotor Airfoil
Figure 6(a) depicts the lift-to-drag ratio estimation, based on the Panel method and the Schlichting boundary layer drag [19], for four NACA airfoils rotor candidates. The NACA 6312 airfoil generates the highest lift-to-drag ratio when subjected to the widest range of the angle of attack of the incoming wind. Other factors such as the material reliability ease of production and the cost, will influence the airfoil selection. As for the tower, circular cross-sections are currently employed, however, a symmetrical airfoil cross-section may be a better choice, since it will produce a more confined and smaller wake than a circular cross-section, and will therefore cause less instability to the rotor.
(a)
(b)
(c)
3.3. Passage of Rotor
The passage of the rotor is simulated through sliding of the rotor subdomain mesh at a velocity of 50 m/s. This linear velocity corresponds to the rotor’s rotational velocity at a radial distance ratio of . The wind turbine is subjected to a prevailing wind speed of 12.5 m/s. The rotor’s subdomain slides against the fixed upstream domain that contains the tower and the fixed downstream domain that contains the wake. The unsteady term of the Navier-Stokes is discretized through the use of a second-order implicit scheme to allow for the largest possible time step. The stability of the time step is limited by the local dilatation wave speed. It is expressed as (minimum (, )/u), where is the local flow velocity and is not known as a priori. is initially iteratively determined ensuring CFL (Courant number)1.
Figures 6(b) and 6(c) depicts the vorticity contours produced at four different instances of the sliding of the rotor past the tower for a number of tower cross-sections. The unsteadiness of the flow can be observed by the asymmetric vorticity contours, and in particular those behind the circular tower cross-section. The tower with an airfoil cross-section produces a confined wake whereas the circular tower produces a free shear that influences the rotor stability and reduces the rotor aerodynamic forces.
Figure 7(a) shows color contours for the rotor region for the three different tower cross-sections. Figures 7(b), 7(c), and 7(d) show the instantaneous pressure coefficients for the airfoil and cylinder at three positions of the rotor transitions, 2C, 0C, and 2C where C is the rotor’s airfoil chord length. Differences between the computed values and those obtained by the inviscid Panel method solution and those obtained with the presence of upstream tower appear to be in favor of Panel method, that is, higher aerodynamic forces. Figure 7(e) displays the values for the three geometries at the instant the rotor is located behind the tower and displays the distribution calculated by the Panel method. The largest aerodynamic forces exerted upon the rotor have been calculated by the Panel method (undisturbed rotor flow), followed by the Airfoil-, then the small circular cross-section, and finally those exerted upon the large circular cross-section tower. The suction side near the LE shows larger difference than the pressure side during the rotor motion. This is because these regions are subjected to separated flow due to the acute angle of attack and the inviscid flow is limited to unseparated flow. The overall angle of attack due to the incoming wind, the rotor translation, and the rotor blade tilting is 9.5 degrees. During the rotor motion there is a fall in the positive value as soon as the rotor enters the wake of the tower zone by the airfoil tower at the LE. This observation may sound trivial, the consequence of this, however, upon the overall rotor force has not been emphasized in literature. Referring back to Figure 7(c), it illustrates the magnification of this effect to a level that forces the stagnation point to move towards the suction side and with a drastic reduction in the aerodynamic forces of over 57% over a period that elapses the translation of few chord lengths as depicted in Figures 8(a), 8(b), and 8(c).
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
These figures also illustrate the evolution of the aerodynamic forces exerted upon the rotor’s airfoil for the large and small circular cross-section towers, respectively. Furthermore, the effect of the tower upon the values extends several chord lengths (3C, 7.5C, and 5C for the airfoil, large circular, and small circular tower cross-sections, resp.). The tower wake is stretched along the direction of the moving rotor. The beauty of these plots is hidden in their quantitative estimation of both the undisturbed (far field tower presence) and the disturbed (close field tower presence) flow field. Since the nacelle is elongated ahead of the tower and the rotor plane is tilted several degrees, the undisturbed rotor may be considered a representation of an upwind wind turbine configuration.
In the deployment of a circular cross-section tower, both the diameter and the tower-rotor clearance have a large impact upon the rotor aerodynamic forces. This is clearly illustrated in Figures 8(b) and 8(c) that show that the impact upon these aerodynamic forces is reduced when the tower has a smaller circular cross-section. However, the effect of the circular cross-section upon the rotor is much more pronounced than that of an airfoil cross-sectional tower.
It is clear that the airfoil-shaped tower has the least impact upon the rotor instability and reduction of aerodynamic forces. This is also clearly captured in the moment around the tower leading edge point. The drop in these coefficients is considerable for both circular cross-sections, and in particular the larger one, both in value and duration leading to over 17% and 57% drop in the aerodynamic forces. Due to the low computed drag value of the rotor in the airfoil shaped tower configuration, one can infer that similar results (with slight variations) will be obtained for the upwind wind turbine.
3.4. Tower Aerodynamic Load
The tower is subjected to a low fatigue cycle for every passage of the rotor through the tower wake. Nonetheless, based on the relative rotor wind speed, the tower is subjected to a much smaller stagnation pressure than the moving rotor blade. Figure 9 depicts the distribution across the chord at three instances in the rotor motion (3C below, directly behind the tower, and 3C above). The computed stagnation value is as low as (0.12) at the leading edge. This increases progressively as the rotor blade approaches the tower, and a maximum value of (0.3) is reached in the trailing edge region at the instant the rotor passes through the airfoil wake. This is followed by an asymptotic drop to the initial value as the rotor moves away from the tower. The asymmetry of the value of the symmetrical tower is observed when the rotor blade is several chord lengths away from the tower as deduced from Figure 9.
The displacement corresponding to the time history of the lift and drag coefficients of the airfoil and cylindrical cross-section towers are depicted in Figure 10. The drag shows the generation of a net appreciable forward force for the cylindrical tower, due to the considerably high rotor velocity. This is a very important finding which facilitates accounting for accurate tower loading. This forward force/impulse appears to be more pronounced for the tower of smaller cross-section and persists several chord lengths whereas it is negligible for the airfoil cross-section. The lift force goes through a sinusoidal pattern as the rotor approaches the tower section.
4. Conclusion
A numerical investigation of the flow field during the intrinsic rotor blade motion for a downstream wind turbine was conducted. The flow was governed by the incompressible turbulent unsteady time averaged Navier-Stokes equations with adjustment of the convective term to account for moving rotor via ALE formulation. The commercial fluent code with user defined function (UDF) is utilized to simulate the flow for the tower and rotor combination. This work shows a step-sized approach to downwind turbine analysis. The initial step involves the building of a conceptual downwind turbine in Unigraphics CAD environment to establish and import two-dimensional domain within the Fluent preprocessor, Gambit. A multiblock quad mesh is then constructed to the fixed and moving domains, and then the velocity inlet, pressure outlet, sliding interfaces, and surface wall boundary conditions are properly assigned. A second-order implicit model in time and a special discretization noniterative PISO scheme solver is used. A Sensitivity study has been conducted on three tower cross-sections, an airfoil with a chord length four times the length of the rotor chord and two circular cross-sections, one with a diameter four times the length of the rotor chord and another two times the length of the rotor chord. Pressure coefficients at the rotor and tower surfaces were generated and compared to Panel inviscid method results. The following conclusions can be drawn:(1)A net reduction in the aerodynamic load for the large cross-section tower can be substantial and up to 57%, while 17% resulted for the smaller circular and less than 5% for the airfoil cross-section.(2)The reduction in the aerodynamic forces persist over a time period corresponding to the translation to several rotor chord lengths, 7.5C, 5.0C, and 2.5C for the large, small, and airfoil tower cross-sections, respectively.(3)The tower suffers from a net forward force due to the rotor motion. The smaller the tower section, the higher is this force coefficient. (4)Despite the low force exerted on the symmetrical airfoil cross-section tower and the low resulted rotor wake, under highly varying wind direction the advantages of airfoil-shaped tower is reduced.
Nomenclature
| : | Area unit vector |
| C: | Rotor chord length |
| : | Drag coefficient |
| : | Lift coefficient |
| c: | Constant |
| D: | Diffusion coefficient |
| D: | Cylinder diameter/aerofoil cord length (m) |
| F: | Strouhal frequency () |
| : | Gravitational acceleration in direction (m/) |
| : | Destruction of turbulent viscosity |
| n: | Integer |
| : | Static pressure (Pa) |
| R: | Distance along blade length (m) |
| R: | Rotor radius/blade length (m) |
| Re: | Reynolds number |
| : | Source term |
| : | User defined term |
| t: | Time (s) |
| : | Flow velocity vector (m/s) |
| : | Grid velocity that is equal to the translating rotor velocity (m/s) |
| U: | Local flow velocity (m/s) |
| : | Mean velocity (m/s) |
| : | Mean velocity in xi direction (m/s) |
| : | Mean velocity in direction (m/s) |
| : | Reynolds stresses |
| (i = 1, 2, 3): | Velocity components in xi direction (m/s) |
| : | Velocity fluctuations about ensemble average velocity (m/s) |
| U: | Upstream undisturbed inlet velocity (m/s) |
| V: | Volume () |
| : | X-direction unit vector |
| : | Cartesian coordinate. |
Greek Letters
| : | Turbulent/eddy viscosity (Pa.s) |
| : | Equivalent viscosity (Pa.s) |
| : | Molecular viscosity (Pa.s) |
| : | Stress tensor components |
| : | Angle of attack |
| : | constant |
| : | Relaxation factor |
| : | Provisional solution |
| : | Time averaged static pressure (Pa) |
| : | Density (kg) |
| : | Time averaged density (kg). |
Acknowledgment
The sponsorship given from Masdar Initiative and Masdar Institute is highly acknowledged.