Evaluating wind farm wakes in large eddy simulations and engineering models

We study wind farm wakes with large eddy simulations (LES) and use these results for the evaluation of engineering models such as the Jensen model, the coupled wake boundary layer model (CWBL), the Turbulence Optimized Park model (TurbOPark), and the wind farm model developed by Niayifar and Porté-Agel (Energies 9, 741 (2016)). We study how well these models capture the wake effects between two aligned wind farms with 72 turbines separated by 10 kilometers in a neutral boundary layer. We find that all considered models over-predict the wind farm wake recovery compared to what is observed in LES. The TurbOPark model predictions on the wind farm wake effect are closest to the LES results for the scenario considered here.


Introduction
Due to the increasing number of offshore wind farms, studies of low-velocity zones far downstream of wind farms, also known as wind farm wakes, are of utmost importance [23]. Measurements of these wakes are performed with diverse methods such as Synthetic Aperture Radar (SAR) [7,18,10,9,2], Doppler radar [28,2], Laser Imaging Detection And Ranging (LIDAR) [33], research air-crafts [31,35,30,21], and supervisory control and data acquisition (SCADA) power data [17]. Long-distance wind farm wakes have been observed up to 55 km in stable atmospheric conditions, up to 35 km in neutral conditions, and up to 10 km in unstable conditions [33]. Additionally, Reynolds-averaged Navier-Stokes solver (RANS) [27], large eddy simulation (LES) [11], and Weather Research and Forecasting (WRF) simulations [33,30] have shown that wind farms can influence each other. As a consequence, long-distance wakes behind wind farms have to be considered when planning new farms in the vicinity of existing ones. Engineering models are widely used in this planning process. Their main advantage is that they are computationally efficient and can evaluate a wide range of possible scenarios. Due to these advantages, various wake models are continuously developed [39,32]. The most simple wake model goes back to Jensen [19] who assumed a linear wake growth. More recent wake models consider more details, e.g. the model of Bastankhah and Porté-Agel [5] relies on the conservation of mass and momentum and predicts a more realistic Gaussian wake shape. Furthermore, different superposition methods have been proposed to account for wake interactions. These interactions are modeled by considering either a linear superposition of the velocity deficits [22,26] or a linear superposition of the energy deficits [20,41]. Besides, one can consider the wake deficit with respect to the incoming upstream wind speed [22,20], or with respect to the incoming flow speed for that turbine [41,26]. While these models are extensively tested for wind farms, the ability of these models to accurately capture the interaction between different wind farms is still relatively unknown. Hansen et al. [17] compared the wind farm wake effects predicted by engineering models with SCADA data, the WRF mesoscale model, and computational fluid dynamics models and highlighted the necessity of further investigations to achieve more robust predictions. Recently, Nygaard et al. [29] introduced the Turbulence Optimized Park (TurbOPark) model that is designed to model the interaction between wind farms. The TurbOPark model combines the Jensen and the Katić et al. [20] superposition method, often referred to as the Park model, and improves on this modeling approach by accounting for the turbulence intensity in the wind farm as described by the model by Frandsen et al. [13]. Nygaard et al. [29] found that the wind farm power output predicted by TurbOPark agrees favorably with SCADA data for the neighboring offshore wind farms Humber Gateway and Westernmost Rough. In this work, we compare the predictions by the Jensen model, the TurbOPark model, the coupled wake boundary layer (CWBL) model, and the model developed by Niayifar and Porté-Agel [26] for the wind farm wake development against results from a reference LES. In the LES we consider two aligned wind farms with 72 turbines that are separated by 10 km. Due to the high computational costs of the LES we only consider one wind farm layout in this study. In the remainder of the paper, we first summarize the different engineering models and the LES. Subsequently, we present the comparison between the LES results and model predictions before discussing the conclusions.

Jensen Model
The Jensen model [19] is a simple, classical wake model. The fundamental assumption of the model is that the width of the wake behind a wind turbine increases linearly with the downstream distance. Following Jensen [19] the velocity u w in the wake of a turbine is expressed by: Here u ∞ is the incoming free stream velocity, and x is the downstream distance with respect to the turbine. The turbines are assumed to be actuator disks with a rotor diameter D and the thrust coefficient C T = 4a(1 − a) with the induction factor a. The wake diameter growth rate D w is assumed to be linear: where is the wake expansion coefficient estimated based on the logarithmic wind profile. Here z 0 is the surface roughness, z hub the turbine hub-height and κ the von Kármán constant [12].
To account for the wake interactions the Jensen model [19] is combined with the Katić et al. [20] superposition model, which sums up the squared velocity deficits [20], i.e.
Here the summation i is over all turbine wakes at that location. The power of each turbine P T , normalized by the power of the turbines in first row P 1 , is estimated as where the brackets stand for an average over the turbine disk and x T is the turbine position.

Coupled wake boundary layer model
The Jensen model is one of the so-called bottom-up models in which wake deficits are combined using some superposition model to account for the wake interactions. The CWBL model [36,37] suggests to improve the Jensen model by coupling it to the Calaf et al. [6] top-down model such that the predictions from the Jensen and the top-down model are consistent in the fully developed regime of the wind farm. The top-down model parameterizes the wind farm, instead of the individual turbines, using an increased surface roughness z 0,hi . This model is based on the assumption of two constant momentum flux layers, one above the turbine hub-height and one below. Each has a characteristic friction velocity and surface roughness, such that the velocity at hub-height is given as where δ IBL is the height of the internal boundary layer in the fully developed regime of the wind farm, . Here w f indicates the effective wake area coverage in the fully developed regime of the wind farm. The roughness length of the wind farm is defined as: The value of w f and the value of the wake coefficient in the fully developed regime of the wind farm k w,∞ are determined and updated iteratively until the turbine velocity in the fully developed regime of the wind farm is the same (up to a tolerance of 0.1%) in the Jensen (equation (4)) and the top-down model (equation (6)). The effects of the entrance region of the wind-farm are considered by assigning a wake coefficient to the wake originating from each individual turbine: where m is the number of turbine wakes that overlaps with the turbine of interest and ζ = 1 is determined empirically. k w is the wake expansion coefficient in the entrance region of the wind farm, see equation (3). For further details we refer to Ref. [37].

Turbulence Optimized Park Model
Nygaard et al. [29] extended the Park model based on the idea that the wake expansion rate is dependent on turbulence intensity. To estimate the turbulence intensity behind a wind turbine, the sum of the ambient turbulence intensity I ∞ and the wake added turbulence I w is considered: The wake added turbulence is described empirically with the constants c 1 = 1.5 and c 2 = 0.8 [14].
Assuming that the wake diameter growth rate increases linearly with the turbulence intensity dD w (x )/dx = AI(x ), results in the following wake expansion rate: with α = c 1 I ∞ and β = c 2 I ∞ / √ C T and the model calibration constant A = 0.6 [29]. TurbOPark incorporates equation (11) into the calculation of the wake deficit in equation (2). Further, a prefactor u 0 /u ∞ is added in front of the √ 1 − C T term in equation (1) to take the normalized rotor-averaged inflow wind speed at each turbine position into account. Nygaard et al. [29] set the turbine thrust coefficient based on the velocity at the turbine, which is calculated iteratively based on the wakes induced by upstream turbines and the induced wind farm blockage. This effect is neglected here as the reference LES considers a constant value for C T . We note that Nygaard et al. [29] couple TurbOPark with a wind farm blockage model [29], which is not considered here.

Niayifar and Porté-Agel (2016) Model
Based on the assumption of a wake velocity deficit with an axisymmetric, self-similar Gaussian shape and using conservation of mass and momentum, Bastankhah and Porté-Agel [5] describe the velocity in a wake as Here u 0 is the local average inflow velocity in front of each turbine [43], which is different from u ∞ in equations (1) and (4). y is the spanwise distance with respect to the turbine center. The standard deviation of the Gaussian velocity deficit is: The growth rate k * w (x ) is linked to the turbulence intensity by the empirical expression [26]: k * w (x ) = 0.3837 · I(x ) + 0.003678.
The turbulence intensity is determined using equation (9) in combination with the empirical model for the added turbulence intensity proposed by Crespo and Hernández [8]: Note that equation (15) is only valid in the range 0.065 < I ∞ < 0.14, 5 < x /D < 15 and 0.1 < a < 0.4 [8]. Further, Niayifar and Porté-Agel [26] only consider the turbulence intensity from the closest upstream turbine, i.e. the wake added streamwise turbulence intensity at turbine j is given by the maximum of the added streamwise turbulence intensity induced by turbine k at turbine j: I w j = max I w kj · 4A w /πD 2 . Here A w indicates the intersection between the wake and the rotor area. We note that the wake velocity (equation (12)) is only defined starting from approximately two rotor diameters downstream of each turbine [5]. Considering that the inter-turbine distance in a wind farm is typically much larger this does not affect the model applicability. In contrast to the previous models the wake interaction is modeled as [26]: Different from Katić et al. [20] (equation (4)), Niayifar and Porté-Agel [26] consider the local mean inflow velocity u 0 in front of each turbine instead of the velocity in front of the farm u ∞ . Furthermore, this method considers a linear superposition of velocity deficit instead of a linear superposition of the energy deficit [20].

Large Eddy Simulations
We perform LES of a neutral atmospheric boundary layer (ABL) flow driven by a geostrophic wind. The LES data is used as a reference to evaluate the ability of the various engineering models to capture the development of the wind farm wake. The used LES code is an updated version of the code developed by Albertson and Parlange [3]. The governing equations are the filtered, incompressible continuity and Navier-Stokes equations: Here the tilde represents spatial filtering with a spectral cut-off filter at the LES grid-scale ∆ and u i represents the filtered velocity field components. τ ij = u i u j −ũ iũj is the trace-less part of the sub-grid scale (SGS) stress tensor and it is modeled with the anisotropic minimum dissipation model [1] with a Pointcaré constant of C i = 1/12 in horizontal and C i = 1/3 in vertical directions. The trace of the SGS stress tensor is absorbed into the filtered modified pressurẽ p + =p/ρ − p ∞ /ρ − τ kk /3. The Coriolis parameter is given by f c = (0, 2Ωcos(Φ), 2Ωsin(Φ)) with the Earth's rotation angular speed Ω = 7.3 · 10 −5 rad/s and the latitude Φ = 52 • . G is the geostrophic wind velocity with a magnitude of |G| = 11 m/s. The effect of the wind turbines is added as a body force f i using an actuator disk approach [6]. Effects of resolved viscous stresses are neglected, since a very high Reynolds number limit is assumed. The wall shear stress at the ground is modeled using the Monin-Obukhov similarity theory [24]. The boundary conditions at the top of the domain are zero vertical velocity and zero shear stress. Time integration is performed using a second-order accurate Adams-Bashforth scheme. Derivatives in the vertical direction are calculated using a second-order central finite difference scheme. In the horizontal directions a pseudo-spectral method is applied, and therefore a concurrent precursor inflow method is used to remove periodicity and generate the inlet ABL flow [38]. To ensure that the incoming wind direction at hub-height is aligned with the x-axis, we apply a proportionalintegral controller [4,34]. This simulation approach has recently been validated for various benchmark cases, see Refs. [16,15,40].
The computational domain has a size of L x = 54 km, L y = 7.2 km, and L z = 4 km, in the streamwise, spanwise, and vertical direction, respectively. The domain is resolved on a grid with 1800×480×480 nodes. This results in a resolution of ∆ x = 30 and ∆ y = 15 m in the streamwise and spanwise directions. A stretched grid with a constant ∆ z = 5 m up to z = 1.5 km and larger grid cells above is employed. The fringe region in the concurrent precursor method is ∆x Fringe = 3 km in streamwise direction and ∆y Fringe = 1 km in spanwise direction. We use a surface roughness of z 0 = 0.002 m, which is a typical value for offshore conditions. We consider two aligned wind farms with 12 × 6 turbines. The inter turbine spacing is s x = 7D in streamwise and s y = 5D in spanwise direction. The first wind farm is positioned 7 km downstream of the inflow region, and the distance between the wind farms is 10 km. The turbines have a diameter of D = 120 m and a hub-height z hub = 100m. The thrust coefficient is C T = 3/4 and a = 1/4. The boundary layer in the precursor domain of the LES reaches a quasi-steady state with a fully developed turbulent flow at the end of the 9th hour. Subsequently, the simulation in both domains is continued concurrently for two more hours before the statistics are collected in 3.5 additional hours. From the LES we obtain that the ambient turbulence intensity at hub-height is I ∞ (z hub ) = σ u /u h = 9.02%. The internal boundary layer height (δ IBL ), based on the height where the time-averaged velocity is 99% of the incoming flow speed at that height, is about 700 meters at the end of the second wind farm. This height is used in the CWBL calculations. To conclude, we mention that all engineering model calculations are only performed at hub-height.    Figure 1 shows the time-averaged velocity at hub-height obtained from the LES and the four engineering models under consideration. Figure 1(a) shows that the wakes of the individual turbines are distinguishable up to 6 km downstream of the wind farm. Further downstream, the individual wakes merge into a more homogeneously mixed wind farm wake. Similar observations have been made by synthetic aperture and dual-Doppler radar measurements [33,28]. Further, these observations [33,28] and the LES show that the wind farm wake does not expand much beyond the spanwise extent of the wind farm. We note that the spanwise variations in the wind farm wakes are caused by the streamwise elongated rolls in the ABL [25]. Figure 1(b-e) show that the various models qualitatively capture the trends that individual wakes can be observed up to a certain distance behind the farm, after which a more homogeneous wind farm wake is formed. Additionally, the figure reveals clear differences between the analytical models and the LES. A comparison of the different panels shows that the Jensen, CWBL, and Niayifar and Porté-Agel [26] models strongly overpredict the recovery rate of the wind farm wake. The Jensen model prediction gives a stronger wake deficit directly behind the wind farm, due to which the stronger wind farm wake recovery only becomes apparent further downstream of the wind farm. The TurbOPark model captures the wind farm wake recovery observed in the LES most accurately. In figure 2 we show the velocity at hub-height, averaged over the spanwise locations where the turbines are located. In agreement with the studies by Christiansen and Hasager [7] and Wu and Porté-Agel [42] the velocity deficit at 10 km downwind of the farm, i.e. where the second farm starts, is about 3% of the velocity in front of the first farm. Interestingly, the LES and the analytical models only show slightly lower velocities behind the second farm than in the wake of the first farm. Additionally, this figure shows that all considered models overestimate the recovery of the wind farm wake behind both farms when compared to the LES. In the Jensen model, the wake deficits inside the wind farm are strongly overestimated. As a result, the Jensen model overestimates the wind farm wake up to 2 km downstream of the wind farm, while it underestimates the wind farm wake further downstream. The other models capture the wake deficit inside the wind farm reasonably accurately, but these models also overestimate the wind farm wake recovery. The TurbOPark model predictions for the wind farm wake recovery are closest to the LES results for the scenario considered here. Figure 3 compares the wind turbine power output as a function of downstream position obtained from the models and LES for the first and second wind farm. The LES results in  Figure 3. Comparison of the power output as function of downstream direction obtained from the models for a) the first and b) the second wind farm. The results are normalized by the power production of the first row of the first farm. Error bars are the standard deviation of the power per row over time.

Results
this figure reveal that the power output of the second farm's first row is about 11% lower than the power production of the first row of the first farm. In agreement with the wind farm wake recovery discussed above, all models significantly overestimate the production of the first row of the second farm. Particularly the reduction in power output of the second farm's first row compared to the the power output of the first row of the first farm is about 9% for the TurbOPark model, 8% for the Niayifar and Porté-Agel [26], 7% for the Jensen and 4% for the CWBL model. Consequently, the TurbOPark model captures the effect of the wind farm wake most accurately. Furthermore, the TurbOPark model most accurately represents the power production as a function of downstream direction. However, in contrast to the CBWL and Niayifar and Porté-Agel [26] model, the TurbOPark overestimates the production in the entrance region of the wind farm. Overall, the TurbOPark, CBWL, and Niayifar and Porté-Agel [26] model predictions are closer to the LES results than the Jensen model.

Conclusions
In this study, we performed LES of a wind farm cluster with two wind farms separated by 10 km and evaluated the performance of four engineering models taking the LES results as a reference. An important finding is that all considered models overestimate the wind farm wake recovery compared to what is observed in LES. From the models considered here, the TurbOPark model provides the best estimate for the recovery of the wind farm wake. However, it is essential to emphasize that more work is required to assess the performance and the capability of various engineering models to capture wind farm wake effects, as only one scenario is considered in this study.