High-Order Numerical Simulations of Wind Turbine Wakes

Previous attempts to describe the structure of wind turbine wakes and their mutual interaction were mostly limited to large-eddy and Reynolds-averaged Navier–Stokes simulations using finite-volume solvers. We employ the higher-order spectral-element code Nek5000 to study the influence of numerical aspects on the prediction of the wind turbine wake structure and the wake interaction between two turbines. The spectral-element method enables an accurate representation of the vortical structures, with lower numerical dissipation than the more commonly used finite-volume codes. The wind-turbine blades are modeled as body forces using the actuator-line method (ACL) in the incompressible Navier–Stokes equations. Both tower and nacelle are represented with appropriate body forces. An inflow boundary condition is used which emulates homogeneous isotropic turbulence of wind-tunnel flows. We validate the implementation with results from experimental campaigns undertaken at the Norwegian University of Science and Technology (NTNU Blind Tests), investigate parametric influences and compare computational aspects with existing numerical simulations. In general the results show good agreement between the experiments and the numerical simulations both for a single-turbine setup as well as a two-turbine setup where the turbines are offset in the spanwise direction. A shift in the wake center caused by the tower wake is detected similar to experiments. The additional velocity deficit caused by the tower agrees well with the experimental data. The wake is captured well by Nek5000 in comparison with experiments both for the single wind turbine and in the two-turbine setup. The blade loading however shows large discrepancies for the high-turbulence, two-turbine case. While the experiments predicted higher thrust for the downstream turbine than for the upstream turbine, the opposite case was observed in Nek5000.


Introduction
The wind-energy sector has been subject to continuous progress throughout the last decades due to the attempt of many countries to emancipate themselves from fossil fuels. Landscape restrictions and judicial regulations lead to wind turbines being placed closely together in small clusters or wind farms. Thus, apart from being subjected to the changing characteristics of the atmospheric boundary layer, individual turbines are also exposed to the wakes of other turbines. The wake from an upstream turbine usually leads to an increase in turbulent intensity and fatigue loading of downstream turbines. Due to the decreased energy content of the flow in the wake of a turbine the downstream turbines extract less energy from the atmosphere than the first row of turbines. Using large-eddy simulations Nilsson et al. [1] show that the power production of downstream turbines in the Lillgrund wind farm drops to less than 50% of the upstream turbines. If the flow approaching the upstream turbines is already turbulent, then the added momentum exchange between the wake and the freestream causes an increase in the mean streamwise velocity of the flow approaching the downstream turbines [2]. To construct and operate closely clustered wind turbines it is imperative to be able to estimate the influence of wake interaction on wind turbine loading. When using numerical simulations the limiting factor is the Reynolds number Re c = W ∞ c/ν, where W ∞ is the free-stream velocity, c is the chord, and ν the viscosity. The length scales of a fully-resolved wind-turbine simulation range from a fraction of the boundary-layer thickness of the rotor blades to the size of the largest eddies in the atmospheric boundary layer. To reduce the range of length scales to be resolved, modeling assumptions need to be introduced. Crespo et al. [3] provide an overview of wake models ranging from the assumption of turbines in wind farms as distributed roughness elements to individual wake models with superposition of interacting wakes. At the other end of the complexity spectrum lies the possibility to resolve the blades, as done in Zahle et al. [4] with a single wind turbine based on the Reynolds-averaged Navier-Stokes equations. Widely used modeling approaches of intermediate complexity are the actuator methods. The actuatordisk method (ACD) e.g. uses the blade element method to compute the blade forces on a disk, coupled with the Navier-Stokes equations [5]. The ACD simplification has been shown to provide accurate predictions of the blade loading and the averaged wake properties. However, the spiral tip and root vortices that are in reality released by wind turbines are replaced by a continuous vortex sheet due to the fact that the blade forcing is applied along the entire rotor disk plane. As a more advanced technique, the actuator-line method (ACL) was proposed in Sørensen and Shen [6] which models each rotor blade as a line force, thus enabling the release of distinct tip and root vortices. A main source of uncertainty in the actuator methods is the airfoil data with which, using two-dimensional airfoil theory, the sectional lift and drag forces are computed. Both ACL and ACD are currently being used in the context of large-eddy simulations (LES) to analyze single wind turbines and wind-turbine interaction. The solvers mostly employ finitevolume techniques of varying orders of accuracy for spatial discretization, e.g. EllipSys3D [7] or OpenFOAM [8]. In Ivanell et al. [9] and Sarmast et al. [10] the dominant modes responsible for the onset of instability of a single wind turbine wake are investigated while in Nilsson et al. [1] the ADM method is used to model the complete Lillgrund wind farm and analyze its performance. Other approaches such as free vortex wake methods have also been employed to investigate wind turbine wakes.
This paper studies the wake of a single wind-turbine setup and the interaction between two consecutive wind turbines by means of numerical simulation and compares them to experimental results ("Blind Tests") by Krogstad and Eriksen [11] and Krogstad et al. [12] (refered to as NTNU in the figures) and numerical simulations using the finite-volume solver EllipSys3D [7].

Numerical Setup
We integrate the incompressible three-dimensional Navier-Stokes equations with the spectralelement code Nek5000 [13], which has been shown to provide accurate results of wind turbine wakes while requiring lower resolution than comparable finite-volume codes [14]. This is due to the use of high-order Legendre polynomials on Gauss-Lobatto-Legendre (GLL) quadrature points for spatial discretization. However, due to non-equidistant meshing the required time step in Nek5000 was smaller than in the finite volume code.

Domain and boundary conditions
In the following we non-dimensionalize all quantities with the wind-turbine radius R and the free-stream velocity W ∞ . The numerical domain is constructed such that it matches the wind- tunnel dimensions adjusted for the displacement boundary-layer thickness of the tunnel walls. This enables the use of a slip condition on the lateral boundaries without the need to resolve the boundary layers forming on the wind-tunnel walls. The distance between the inlet and the single turbine is l in = 7. In the two-turbine cases with low turbulence the distance of the upstream turbine is also l in = 7 while in the high-turbulence case the upstream turbine is placed at l in = 4. We employ a convective outflow boundary condition at the outlet. When comparing with wind-tunnel experiments, the turbulence in the freestream must be taken into account. This study includes the turbulence as random sinusoidal perturbations superposed with uniform plug flow as an inlet boundary condition [15] as opposed to the EllipSys3D simulations [16,17], where the Mann turbulence model [18] is used to generate turbulence which is then included as a force field upstream of the turbines. The turbulence properties in the domain are set at the inlet to match the decay of the turbulence intensity throughout the experimental wind tunnel [11,12] and the turbulence intensity at each wind turbine. In the high-turbulence case the experiments report a turbulence intensity of T i up = 10% at the upstream turbine and T i down ≈ 5% at the downstream turbine. The grid is Cartesian and contains approximately 80000 equidistantly spaced spectral elements with 9 th order polynomials. The actuator lines are discretized with N ACL = 70. The radius-based Reynolds number used in the simulations is Re R = 50000. The blade data is extracted using the chord-based blade tip Reynolds number is Re c,tip ≈ 100000 which is the Reynolds number reported in the experiments. Sørensen and Shen [6] noted that the Reynolds number does not affect the overall wake behaviour above a certain minimum therefore justifying the chosen Reynolds number in Nek5000.

Wind turbine model
The rotor blades are modeled as actuator lines which compute the blade forces at each time step based on the inflow velocity, the angle of attack α and the tabulated, chord-based Reynolds number dependent lift and drag coefficients C L,D obtained by Sarmast and Mikkelsen [19] and corrected for solid, wake and streamline blockage. After the lift and drag forces F 2D = (F L , F D ) have been computed at the ACLs using two-dimensional airfoil theory the blade forces are spread out using the convolution of the force with a three-dimensional Gaussian kernel [6].
The wind turbine towers are also modeled using a body-force approach similar to the ACL method. The experimental data taken from Sumer and Fredsøe [20] is dependent on the Reynolds number and the incoming turbulence intensity. The tower forces are computed Both nacelle and tower forces are distributed on to the GLL points using the three-dimensional Gaussian kernel. Figure 2 shows the wake development for the optimal tip-speed ratio (TSR) λ = 6 of a single wind turbine. The blade-tower interaction creates an additional velocity deficit.

Tower influence
The tower introduces a significant asymmetry in the wake shown in figure 3(a) while the nacelle model leads only to a slight decrease in the streamwise velocity in the wake center and a small increase in turbulence intensity. At approximately x = −0.3, in figure 3(a) the tower wake is observed in the averaged streamwise velocity at ∆z = 2 (position shown in figure 1(a)). Further downstream a distinct tower wake cannot be observed in the averaged streamwise profiles. However, the velocity at the centerline is decreased. Due to the fact that the nacelle model only consists of a constant drag force, the unsteady three-dimensional turbulent nature of the flow downstream of the nacelle is not captured and the high streamwise turbulent stress at the rotor center of the experiments in figure 3(b) at ∆z = 2 is not replicated in the simulation. The interaction between the tower and blade wake contributes to an early breakdown of the helical tip and root vortices. The turbulent stress peaks generated by the tip vortices and the wake shear layer in Nek5000 is five times smaller than in the experiment due to the larger vortex core size of the ACL method which is defined by the width of the Gaussian and the finite resolution of the simulation. Additionally, as the blade boundary layer is not resolved in the ACL method any boundary layer turbulence is not represented in the ACL simulation. A simulation conducted at a higher Reynolds number shows that the overall wake behaviour is similar to that of the low Reynolds number. The tower in EllipSys3D, though modeled similarly to Nek5000, influences the wake differently, creating a much more centered "tower dent" as seen in figure 3(a). The wake deficit generated by the blades is well predicted by EllipSys3D.
To the authors knowledge the influence of this type of tower model on the wake center has not yet been investigated numerically. By computing the minimum of the time-averaged azimuthal velocity the wake center position along the streamwise axis is determined. Experimental data provided by Blomhoff [21] and Pierella [22] show that the wake center defined by the minimum azimuthal velocity is shifted downwards and spanwise due to the tower wake: At ∆z = 6 (z = 13 from the inlet) the wake center is shifted to (x c , y c ) = (−0.1, −0.25) for λ = 5.5 [21] and (x c , y c ) = (−0.157, −0.26) for λ = 6 [22]. A similar effect may also be shown through numerical simulation, though at a smaller scale. Figure 4 shows the deviation of the wake center from the rotor centerline at TSR λ = 6. Further upstream the azimuthal velocity minimum is difficult to determine exactly due to blade-tower interaction. Without the tower the wake remains centered. The simple tower body force used in this study leads to a wake center which is shifted downwards to (x c , y c ) = (−0.06, −0.05) at ∆z = 6. A reason for the discrepancy of the spanwise shift of the wake center may be that the tower model does not take into account the variation in the angle of attack with which the flow approaches the tower due to the tangential induction of the blades. If the variation of the tower angle of attack is taken into account, then the deviation in the spanwise direction x at ∆z = 6 is closer to the experiments. The vertical displacement of the wake center is probably underestimated due to the fact that the tower model is based on two-dimensional theory and does not take the radial velocity component of the wind turbine wake into account, when computing the tower lift and drag forces.

Tip-speed ratio variation
This section discusses the TSR influence on the wake structure and loading. Figure 5 shows a comparison of the power and thrust coefficient (P is the power and T the thrust) with the experimental curves [11] and EllipSys3D [17]. The power coefficient is in good agreement with the experimental data at off-design conditions. However, at the optimal TSR λ = 6 the power coefficient of Nek5000 differs significantly (∆C P ≈ 9%), while EllipSys3D agrees well. The numerical thrust coefficients are generally lower than the experimental data, which is in accordance with most simulation results comparing with this set of experimental data. In Krogstad and Eriksen [11] it is shown that a fully resolved turbine simulation achieves a better agreement of the thrust coefficients with the experiments. Figure 6 shows the averaged streamwise velocity and turbulent stress profiles at off-design conditions. At λ = 3, the blades are fully stalled. Due to the larger distance between the vortices the tip and root vortices remain stable for a long distance downstream, causing the two distinct turbulent stress peaks at the tip vortex position to persist in the numerical simulations. In Nek5000 additionally a dip in the turbulent stress peaks are visible at ∆z = 2, where the vortex center is located and the streamwise velocity fluctuations are smaller. Discrepancies with the experiments may be caused by the fact that the turbulence produced by the nacelle geometry, which is a main source of turbulence at λ = 3 is included in neither the Nek5000 nor the EllipSys3D simulations. Additionally, the heightened level of turbulence due to the flow separation along the entire blade is not represented in the simulations. For λ = 10, close to the runaway TSR, the rotational velocity of the turbine is very large. The distance between consecutive tip vortices is decreased, which leads to increased interaction between the vortices and a quicker breakdown to turbulence. vortices are unstable and experience pairing and merging already before ∆z = 2 giving rise to the high turbulent stress at the tip vortices. Where the experiments and Nek5000 predict similar turbulence intensity levels at ∆z = 10 the streamwise turbulent stress in EllipSys3D increases and is significantly higher than the Nek5000 results.

Two turbines with spanwise offset
This section discusses the influence of partial wake interaction on the loading and the wake development of the two-turbine setup shown in figure 1(b). The experimental data is obtained from Krogstad et al. [12] and also compared to numerical data from EllipSys3D extracted from the same paper. The upstream turbine operates at the design TSR of λ 1 = 6 while the downstream turbine operates in stall λ 2 = 3.5, at optimum condition λ 2 = 4.75 and at a runaway TSR λ 2 = 8. The power and thrust coefficients of both low and high turbulence cases are presented in figure 7. The power coefficient for the low turbulence case agrees well at all operating conditions. The thrust coefficient curves of the upstream and downstream turbines are on top of each other in the experiments. In Nek5000 the thrust coefficient of the upstream turbine T 1 is overpredicted by 11%. However, in the so-called "Blind Test 4" published later by Bartl and Saetran [23] a similar case is presented with T i = 0.23% with aligned and centered turbines. There the thrust of the upstream turbine at λ 1 = 6 is significantly higher at C T ≈ 0.81 and matches the numerical results better. The thrust coefficient of the downstream turbine T 2 is underpredicted by up to 25%, which may be due to the underestimation of the streamwise velocity speedup around the wake of the upstream turbine in the numerical simulations (visible in the wake profiles in figure  8(a-c)). In the high turbulence case, the experimental power coefficient curve is smoother in the stall region thus suggesting significant influence of the turbulence on the stall behaviour [12]. The thrust coefficient of the downstream turbine is in this case higher than that of the upstream turbine. For both solvers, Nek5000 and EllipSys3D, the influence of the turbulence on the thrust and the power coefficients is much smaller than in the experiments; the thrust coefficient of the upstream turbine is higher than that of the downstream turbine in both Nek5000 and EllipSys3D. This trend for the upstream turbine at high incoming turbulence is again consistent with the results published by Bartl and Saetran [23]. At similar inflow conditions (T i = 10%) the same upstream turbine (though centered in this case) has a thrust coefficient of C T ≈ 0.83 and a power coefficient of C P ≈ 0.47, which is very similar to the Nek5000 results. Figure 8(a-c) shows the wake development at two streamwise positions behind the downstream turbine. The second turbine wake is superposed on the upstream wake between −0.5 ≤ x ≤ 1.2. In both codes the averaged velocity profiles for λ 2 = 3.5, 4.75 transition quicker than the experiments to a Gaussian-like profile. For λ 2 = 8 the wake deficit is overpredicted, which is counter-intuitive as the thrust coefficient of the downstream turbine was underpredicted by 25% in Nek5000. The streamwise turbulent stress is depicted in figure 8(d-f). The location of the peaks from the tip vortices of T 2 is well predicted by the numerical simulations. Apart from the turbulence peak in figure 8(e) at ∆z = 2 at x = 1.5 the magnitude of the peaks is also captured well by Nek5000. At ∆z = 6 the tip vortices of T 2 are either broken down or in the process of breaking down, resulting in a much smoother distribution of turbulent stress with the peak at x = 1 − 1.5 being the result of the tip vortex of T 2 , which is not subjected to the upstream wake. At high incoming turbulence intensity the averaged wake profiles compare well with the experimental profiles (not shown). At ∆z = 6 the experimental velocity profiles are Gaussian shaped and match the numerical results. The turbulent stress profiles are well predicted.

Conclusions
The aim of this study was i) to validate the actuator line method in the spectral-element code Nek5000 against an extensive experimental dataset, ii) to test the turbulent inflow boundary condition and iii) to analyze the influence of different additional models (tower, nacelle). The ability of Nek5000 to capture the wake interaction is evaluated in comparison to experiments and additional data from the finite-volume code EllipSys3D. The paper first investigates the wake and loading of a single-turbine setup. The power coefficient is in good agreement with the experimental data, while the thrust coefficient is lower for both numerical solvers. The horizontal velocity and turbulent stress profiles show the main features of the experimental dataset for all tip-speed ratios. Difficulties are only observed where the turbulence generated by the nacelle and blade separation is dominant (i.e. for tip-speed ratio λ = 3). The wake center is identified based on the minimum azimuthal velocity and compared to the wake center identified in experiments [21]. While the wake center is shifted in the same direction as in the experimental case [21] the magnitude of the shift is much smaller. Changing the angle of attack of the tower lift and drag forces to take into account the tangential induction of the blades leads to an improved comparison of the spanwise displacement of the wake center with the experiments. The second investigation reproduces the two-turbine experimental setup by Krogstad et al. [12], where the turbines are offset in the spanwise direction and subjected to different turbulence levels. The downstream turbine tip-speed ratio is varied to include stall, optimal and runaway conditions. For low turbulence intensity the power coefficients agree well with the experiments. The thrust coefficient of the downstream turbine is up to 25% different from the experimental results due to the underestimation of the speedup around the wake of the upstream turbine. The thrust coefficient of the upstream turbine differs both in the high and low turbulence case significantly from the experiments but is found to match data from comparable cases by Bartl and Saetran [23]. The wake profiles agree well for both low and high turbulence cases with both wake deficit and turbulent stresses being well predicted. No large differences are observed between the spectral-element code Nek5000 and the finite-volume code EllipSys3D. The only difference that may be observed is that the streamwise turbulent stress in Nek5000 sometimes showed a better agreement with the experimental results (e.g at tip-speed ratio λ = 10). Using a comparable grid size the time-averaged velocity and turbulent stress do not differ significantly.
The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC). The first author acknowledges Dr. Sasan Sarmast for providing the numerical data from EllipSys3D.