Comparison of wind turbine wake properties in non-uniform inflow predicted by different rotor models

The wake of the 2MW NM80 wind turbine subject to non-uniform and laminar inflow conditions is simulated using CFD with a fully resolved rotor geometry, an actuator line method and actuator disc method, respectively and in all simulations the wake properties are compared. Based on the comparison the strengths and limitations of the models are pointed out.


Introduction
Most CFD studies of wind turbine wakes use actuator disc (AD) [1] or actuator line (AL) [2,3,4] models for representing the rotor because these models do not require the rotor boundary layer to be resolved, and thus demand much lower spatial and temporal resolution than simulations where the rotor geometry is fully resolved (FR) [5,6]. The use of AD and AL models are usually justified by the assumption that after a few rotor diameters the details of the wind turbine rotor geometry become irrelevant to the flow behaviour. However, it is important to quantify how much physics is lost in the wake when using AD and AL models in comparison with FR computation. To address this question, a systematic comparison of the three different types of wind turbine models has been initiated. In the initial part of this comparison the three methods were used to simulate the wake of the NREL 5MW wind turbine operating in a non-sheared inflow and the predicted wake properties were compared [7]. In laminar inflow the three methods showed good agreement in terms of predicted mean velocity deficit up to 2 diameters downstream of the turbine. Further downstream the FR simulation predicted a faster smearing of the mean gradients in the wake than the AD and AL, which was caused by higher turbulence levels in the shear layer generated by the tip and root vortices. Generally, the wake properties in uniform inflow were found to be quite different between the FR simulation on one hand and the AL and AD methods on the other hand. However, when imposing a turbulent inflow, it was shown that the wake velocity and turbulent kinetic energy predicted by the FR and AL were in very close agreement. Thus, the study indicated that the differences observed in uniform inflow does not play an important role if the inflow is turbulent. In the present work the comparison of the methods is extended by using them for simulating the wake of the 2MW NM80 wind turbine subject to non-uniform and laminar inflow conditions. A comparison of different rotor models subject to atmospheric boundary layer flow have been carried out previously by Wu and Porté-Agel [8] and Porté-Agel et al. [9] using LES. They simulated both a stand alone turbine and a small operational wind farm using both AD and AL representations of the rotor. Model predictions were compared with both a wind tunnel experiment and full scale measurements and generally good agreement was shown. However, it was shown that using a constant loaded drag AD which do not model wake rotation is less accurate than using actuator discs or lines with airfoil data. The present work goes a step further by comparing the AD and AL representations based on airfoil data with FR simulations. Furthermore, the AD model used in the present work is more advanced than those used in most previous work since it can simulate a fully non-uniform loading and do not assume axi-symmetric loading (see section 2.1.3). Two cases are considered in the present work: a vertical inflow shear and a yaw misaligned turbine, respectively. Thereby, it can be investigated whether the findings from the comparison in uniform inflow still hold when the inflow and blade loading is non-uniform.

EllipSys3D
EllipSys3D is an incompressible finite volume Navier-Stokes flow solver developed by Sørensen [10] and Michelsen [11,12]. The code is formulated in general curvilinear coordinates and the flow variables are collocated in the mesh to facilitate complex mesh geometries. The SIMPLE algorithm [13] has been used to solve the pressure equations and pressure-velocity decoupling is avoided by using the Rhie-Chow algorithm [14]. The solution of the pressure equation is accelerated by a multi-grid technique [10]. The system is parallelized in a multi-block structure, where the problem can be distributed across multiple processor. The communication of data between each processor is done through the MPI libraries. The present simulations are carried out using the DES version k-ω SST model by Strelets [15] and the convective terms are discretized using a hybrid scheme combining the QUICK scheme and the fourth order CDS scheme.

Fully Resolved Rotor Simulations Using an Overset Grid Method
The simulations with the fully resolved (FR) rotor geometry use an overset grid method [16] where a combination of curvilinear body-fitted and off-body meshes are used to resolve the domain. The boundary conditions on the overlapping interfaces are based on interpolated values of velocity from neighbouring grids using trilinear interpolation. An explicit correction of the conservation error associated with the non-conservative interpolation is implemented to ensure a divergence free field. Along the internal overlapping boundaries the solution of the pressure is based on the mass fluxes calculated from the momentum equations.

Actuator line model
In the actuator line (AL) model [2,3] the blades of the rotor are represented as lines along which body forces, determined from the local angle of attack and a look-up table of aerofoil data, are distributed using a suitable smearing function. In the present work the lines are represented with 32 elements and the loading is smeared in a 3D Gaussian manner using a smearing parameter corresponding to approximately three times the grid spacing of the background mesh.

Actuator disc model
The actuator disc (AD) model used here is described in [17] and was originally developed by Mikkelsen as an extension of the model presented in [3]. The actuator disc is discretized with a polar grid and the loading applied to each differential area dA = rdrd is determined from the local flow conditions at the disc using aerofoil data and a blade element approach in the same manner as described by Réthoré et al. [18]. In order to avoid pressurevelocity decoupling the body forces on the disc are smeared radially and perpendicularly to the disc in a 1D Gaussian manner as described by Troldborg et al. [19]. It should be emphasized that the AD model used here differs from most other AD methods since it uses local flow conditions for determining the loading everywhere on the disc and thus is fully three-dimensional. In contrast, other AD models typically use the velocity at the position of the blades to determine the loads and subsequently spread these loads in the azimuthal direction, which means that these models essentially assume axisymmetric flow conditions. In the present work the disc is discretized with 32 radial and 80 angular elements and the smearing parameter in the normal and spanwise direction, respectively, is set to a value corresponding to approximately three times the local grid spacing as suggested by Troldborg et al. [19].

Numerical Setup 2.2.1. Rotor Properties
The simulations presented in this work has been conducted using the 2MW NM80 turbine. This turbine has has a rotor diameter of 80 m and a hub height of 57 m. At 8 m/s the rotational speed of the rotor is 17.25RP M . The rotor is modelled without tower, nacelle and spinner and is assumed completely stiff in the simulations. The NM80 turbine is tilted 5 o but this is not included in the simulations. Further details about the turbine and the airfoil data used for the AD and AL simulations are presented in [20].

Computational Meshes
The simulation with a fully resolved geometry consisted of four overlapping block groups, one body-fitted grid around the rotor and three background grids to resolve the wake and farfield. The surface mesh around each blade of the rotor had 256 cells in the chord-wise direction and 128 cells in the span-wise direction and using 96 cells in the normal direction. The first cell in the boundary layer had a height of 1×10 −6 m corresponding to a y+ of less than 2. The mesh was grown outwards 7 m using HypGrid3D [21]  . This block group consisted of 308 blocks of 32 3 cells. The total grid assembly contained 26.7×10 6 cells. Figure 1 shows a front-and sideview of the mesh. The boundary condition on the rotor surface was set as no-slip, and slip-wall (symmetry) at the ground surface. The simulations with the actuator line and actuator disc methods used the same background mesh as the fully resolved simulations. The mesh resolving the rotor area and near wake (within 1D downstream) was thus coarser in the AD and AL simulations since the refined near-wake grid used in the fully resolved simulations was not included.

Simulation Setup
As mentioned in the introduction, simulations have been carried out both in sheared and yawed inflow, respectively. In the simulation with vertical shear the velocity at the inflow boundary was specified according to a power law profile:  Here α = 0.55 is the power law exponent, H = 57m is the hub height and y is the height above ground. In the yaw misalignment simulations a uniform inflow was used and the yaw error was set to 20 o . In all the simulations the mean free-stream velocity at hub height was V ∞ = 8m/s. The time step was ∆t = 0.0014s, 0.017s and 0.1s for the FR, AL and AD simulations, respectively.

Sheared inflow
This section presents the comparison of the three methods in sheared inflow. The spanwise distribution of normal and tangential loads predicted by the three methods at four azimuth angles are compared in Figure 2. Note that an azimuth angle of θ = 0 o corresponds to the blade pointing upwards and θ = 90 o corresponds to the blade pointing horizontally to the left when seen from downstream of the turbine. In all cases the loads have been normalized as follows: where index i refer to the given load component, i = n, θ and F max i,AL is the maximum value of load component i of the AL simulation. The normal forces predicted by the three methods is seen to be in quite good agreement for all azimuth position. However, the AD and AL methods generally under predict this load component in comparison to the FR simulation. The tangential loads predicted by the AD and AL methods show good resemblance with each other but are significantly higher than what is obtained in the FR simulation. The reason for this difference could be due to the used aerofoil data or that the tangential velocity gradient at the position of the disc is very steep and thus the loading determined from aerofoil data is very sensitive to the position where the local velocity triangle is established. Using higher grid resolution near the rotor or a more advanced method for extracting the local velocity at the disc could alleviate this problem. It should also be mentioned that the loads predicted by the AL and AD methods show some dependency on the choice of smearing parameter. This will be a subject for future investigations. Figure 3 and 4 show vorticity magnitude and streamwise velocity contours, respectively, plotted in a horizontal cross-section at hub height. The wake behaviour predicted by the three methods is seen to be quite different from each other. In the FR simulation the tip vortices are far more distinct and the associated vorticity much stronger than in the AD and AL simulations. In effect this causes the wake of the FR to become unstable and break up much earlier than predicted when using AD and AL. The same observation was made in uniform inflow for the NREL 5MW turbine [7] and is mainly caused by an effectively higher grid resolution in the near wake of the FR. However, the models also predict many of the same wake features: In all cases the wake undergoes an asymmetric development, where the deepest deficit initially occurs on the left side of the wake    (when seen from downwind). A similar behaviour has been reported in previous work [22] and is caused by the interaction of wind shear and the rotation of the wake. Furthermore, the root vortex is in all cases clearly skewed towards one side of the wake and the wake becomes unstable earlier on the left side (for an observer in the wake looking upstream). The agreement between the wake behaviour predicted by the AD and AL appears to be quite good, although the far wake of the AL is somewhat more unstable. Figure 5 presents the contours of the streamwise velocity at various downstream positions. Qualitatively, all methods predicts a similar behaviour where the interaction of the shear and the rotation of the wake causes a redistribution of the wake deficit. However, the wake of the FR clearly approach a fully turbulent state before the wake of the two other methods.  show rather good resemblance. Further downstream the predictions of the FR simulation start to deviate significantly from the AD and AL simulations. This is caused by the higher build-up of turbulence in the near wake of the FR computation, which is clearly visible in the contour plots of vorticity in Figure 3. The wake velocities predicted by the AD and AL are generally in good agreement within the first 3D. However, at 5D rather larger discrepancies are observed in both the streamwise and tangential velocity at an azimuth of 90 o , i.e. at the left side of the wake when seen from downwind. The same is also apparent in Figures 3 and 4.

Yawed inflow
In this section we compare the wake properties obtained in yawed inflow using AD and AL representations of the rotor. The FR simulation for this case has not been carried out yet but will be included in the comparison in future work. Figure 8 compares AL and AD predictions of the normal and tangential loads in yawed inflow. Again the shown loads have been scaled according to equation 2. Close agreement is obtained except near the tip where the AL generally predicts slightly higher loads than the AD. Figure 9 shows snapshots of vorticity magnitude contours in a horizontal cross-section at hub height. Furthermore, Figure 10 shows contours of vorticity magnitude at various downstream positions. The wake behaviour predicted by the two methods is similar: In both cases the wake is deflected to the left (when seen from downwind) as a consequence of the reaction thrust imposed upon the air by the rotor. Furthermore, both methods predict a non-symmetric development of the wake where the right side of the tip vortex sheet (when seen from downwind) moves towards the wake center and starts to curl up. This  behaviour is caused by streamwise vorticity, which is generated as a consequence of the wake being skewed compared to the flow direction. This in turn induces velocities traverse to the flow direction which explains the behaviour shown in Figure 10. Despite the similarities there are also some differences: For example, it appears that the shear layer generated by the tip vortices is more unstable in the far wake when using the AL than when using the AD. Figures 11  and 12 shows the streamwise and horizontal cross-flow velocity profiles at different downstream positions. Generally, there is a good agreement in the predictions of both velocity components. The largest differences are seen 5D downstream at r/R = 0.5 and θ = 180 o (i.e. downwards) where the AD predicts larger velocity magnitudes than the AL.

Discussion
The quite large differences obtained in the wake when using the three rotor models in sheared inflow may seem critical for all wake models based on the AD/AL technique. It appears that the small structures of the FR have a signicant impact on the build-up of turbulence and eddyviscosity in the close wake region, which the AL and AD are not able to capture accurately. However, as mentioned in the introduction, a similar study on the NREL 5MW turbine in nonsheared inflow revealed that while large discrepancies in the wake properties were observed in laminar inflow, these discrepancies became insignificant when ambient turbulence was present [7]. It is expected that introducing ambient turbulence will also improve the comparison in nonuniform inflow but this remains to be investigated. Nevertheless, it is important to be aware that disregarding the effect of ambient turbulence in AD/AL wake simulations does give significantly different results in the far wake than what is obtained from FR simulations. A possible remedy to alleviate this difference could be to introduce a source of turbulence or distortion of the vortices at the rotor disc position when using AD or AL models. The overall close agreement between AD and AL wake predictions suggest that the AD is sufficiently accurate for wake studies even when the inflow is non-uniform. However, it should be emphasized that the loads on the AD used in the present study are determined from local flow conditions everywhere on the disc. Other AD models, in which the loads are prescribed or which assume axisymmetry may not compare as favourably to AL simulations. Also it might be expected that larger differences would be observed between the AD and AL models if a finer grid resolution had been used. However, the resolution used in the present study is representative of typical AL wake simulations and thus, unless a very high grid resolution is used, it is justified to use a AD model for wake studies.

Conclusion
The wake of the 2 MW NM80 wind turbine subject to non-uniform inflow has been simulated with three different CFD rotor models (fully resolved rotor, actuator line, actuator disc) and their results compared. Both sheared and yawed inflow conditions have been considered, however, the full rotor simulation was only carried out for the sheared inflow case. In the case of a vertically sheared inflow the three models show good agreement in predicted streamwise velocity up to 2D downstream of the turbine. However the full rotor simulation predicts a much more turbulent wake, which is caused by higher turbulence levels in the shear layer generated by the tip and root vortices. As a consequence the full rotor simulation predicts a faster smearing of the mean gradients beyond 2D downstream than the two other methods. The agreement between the actuator disc and line method was generally found to be good both in sheared and yawed inflow. Based on these findings it is concluded that an actuator disc representation of the wind turbine rotor is as accurate for wake studies as using actuator lines.