Global Trends of Large Wind Farm Performance based on High Fidelity Simulations

A total of 18 high fidelity simulations of large wind farms have been performed by three different institutions using various inflow conditions and simulation setups. The setups differ in how the atmospheric turbulence, wind shear and wind turbine rotors are modelled, encompassing a wide range of commonly used modelling methods within the LES framework. Various turbine spacings, atmospheric turbulence intensity levels and incoming wind velocities are considered. The work performed is part of the International Energy Agency (IEA) wind task Wakebench, and is a continuation of previously published 5 results on the subject. This work aims at providing a methodology for studying the general flow behavior in large wind farms in a systematic way. It seeks to investigate and further understand the global trends of wind farm performance, with a focus on variability. Parametric studies first map the effect of various parameters on large aligned wind farms, including wind turbine spacing, wind shear and atmospheric turbulence intensity. The results are then aggregated and compared to engineering models as well 10 as LES results from other investigations to provide an overall picture of how much power can be extracted from large wind farms operating below rated level. The simple engineering models, although they cannot capture the variability features, capture the general trends well. Response surfaces are constructed based on the large amount of aggregated LES data corresponding to a wide range of large wind farm layouts. The response surfaces form a basis for mapping the inherently varying power characteristics inside very large wind farm, including how much the turbines are able to exploit the turbulent fluctuations 15 within the wind farms and estimating the associated uncertainty, which is valuable information useful for risk mitigation.

both horizontal and vertical staggering. The power output in the farms was further found by Stevens et al. (2015b) to be well correlated with the vertical kinetic flux, in accordance to what is obtained by Calaf et al. (2010), who used LES of large wind farms to quantify the vertical transport of momentum and kinetic energy across the boundary layer. Stevens et al. (2015a), in a different work, also developed a so-called coupled wake boundary layer model of wind farms to predict the power output 60 in large wind farms. The model coupled a wake model for the turbines with a top-down boundary layer model. The coupled model is much simpler and faster than LES, and was shown to compare well with averaged LES results. LES was further used by Stevens (2016) to study how the optimal wind turbine spacing depends on the wind farm length to find that it is remarkably larger for large wind farms when compared to smaller, conventional wind farms. LES was also used by Wu and Porté-Agel (2013) to study turbulent flow inside and above large wind farms considering a neutral boundary layer, with both an aligned and

EllipSys3D
EllipSys3D is a 3D flow solver that was developed at DTU, Michelsen (1992), and the former Risø, Sørensen (1995). It solves the discretized incompressible Navier-Stokes equations in general curvilinear coordinates using a block structured finite volume 125 approach. It is formulated in primitive variables (pressure-velocity) in a collocated grid arrangement. Additional details about this code can be found in Mikkelsen (2003) and Troldborg (2008).

PALM
PALM (Parallelized LES Model) was developed at Leibniz Unversity Hannover and has been applied several years for the simulation of a variety of atmospheric and oceanic boundary layers. Recently, it has been enhanced by a wind turbine model, 130 see Witha et al. (2014). It is an open source, highly parallelized LES model which solves the filtered, incompressible, nonhydrostatic Navier-Stokes equations under the Boussinesq approximation on an equidistant Cartesian grid. The sub-grid scale turbulence is parameterized by a 1.5th order closure after Deardorff (1980). Further details about this code can be found in Maronga et al. (2015).

EllipSys3D
The wind turbines are modelled by DTU and Uppsala University (UU) by using the actuator line (AL) and actuator disc (AD), respectively. In the former, body forces are distributed along rotating lines, while they are distributed along a rotating disc in the latter. Details about the implementation of the AD and AL in EllipSys3D can be found in Mikkelsen (2003) and Sørensen and Shen (2002), respectively. Local blade forces are determined using tabulated airfoil data and the local inflow conditions. 140 In the DTU-AL model, the 2D airfoil data are corrected for 3D effects, see e.g. Hansen et al. (2006). The body forces in the DTU implementation of the AL are further calculated through a coupling with Flex5, which is a full aeroelastic code used for calculating deflections and loads on wind turbines, see Øye (1996) for details on Flex5. The body forces are determined using local velocities along the rotating and potentially deflecting lines, which are transferred from EllipSys3D to Flex5, see Sørensen et al. (2015) for additional details.

PALM
The PALM implementation considers an AD model with rotation (FW-AD-R) in which local body forces are derived from airfoil data. The PALM simulations were performed by ForWind (FW). In contrast to the AL method, the forces are distributed across the rotor plane. This model also includes tower and nacelle effects that are modelled by a drag force approach. See Dörenkämper et al. (2015) for details of the PALM implementation.

Turbine Controller
The three models used in this work include a turbine controller. This causes the applied body forces to be governed by the inflow conditions, meaning that the turbines are not constantly loaded, but operate as "real turbines". Larsen and Hanson (2007) or Hansen et al. (2005) provide a general description of such controllers.

155
Two different three-bladed horizontal axis wind turbines have been considered in the simulations, i.e. the NM80 and the NREL 5MW. The NM80 turbine, see e.g. Aagaard Madsen et al. (2010), has a radius R of 40 m, a hub height z hub of 80m, and a rated power of P 0 = 2.75M W at a nominal hub height velocity of 14 ms −1 . The radius of the NREL 5MW turbine is 63 m, its hub height is 90 m, and its rated power is P 0 = 5M W at 11.4 ms −1 , see Jonkman et al. (2009). Figure 1 compares the C P and C T of the two turbines, which are comparable although the C T is higher for the NREL 5MW than for the NM80 for below rated.
In the coordinate system used in this work, x, y and z correspond respectively to the streamwise, crosswise and vertical directions. The grids used for the simulations are cartesian, and equidistant in the horizontal direction in all cases. The grids are usually stretched in the vertical direction from a significant distance above the wind turbines.

Atmospheric Boundary Layer and Turbulence
All participants simulated a neutrally stratified atmospheric boundary layer (ABL). Details about the methods used to model the ABL and associated turbulence in respectively EllipSys3D and PALM are provided below.

EllipSys3D
EllipSys3D uses the prescribed boundary layer (PBL) method, in which body forces are used to impose any arbitrary vertical 170 wind shear profile, see Mikkelsen et al. (2007) and Troldborg et al. (2014). A comparison of PBL with a wall model approach was performed by Sarlak et al. (2015). This study showed that these two approaches yield very comparable vertical profiles of mean streamwise velocity, shear stress, and streamwise velocity fluctuations in the rotor region when large wind farms are modelled. The present simulations use a combined parabolic and power law profile as defined in Andersen et al. (2015) and given by: where ∆ P BL is the height, where profile transitions from the parabolic to power law profile. H hub is hub height of the turbine, c 1 , c 2 are shape parameters determined to ensure a smooth transition between the parabolic and the power law expression. α P BL is the shear exponent of the PBL. Ambient turbulence is modelled by introducing pregenerated synthetic ambient turbulence using the Mann model, see Mann 180 (1998). Turbulence planes are imposed at an axial position of 6 R and 13 R in the DTU and UU simulations, while the first simulated turbine is located at 10 R and 30 R from the inlet, respectively.

PALM
PALM uses a no-slip bottom boundary condition and the Monin-Obukhov similarity theory between the surface and the first grid level to model the atmospheric boundary layer. Random perturbations are initially imposed on the velocity fields until 185 atmospheric turbulence has developed in a precursor simulation. The latter is performed on a smaller domain with periodic boundary conditions in the streamwise and lateral directions. The precursor results are used to initialize the full simulations with non-periodic boundary conditions in the streamwise direction. Turbulence recycling is also applied, see Maronga et al. (2015) for details.

190
An overview of the numerical methods described in the previous sections are summarized in Table 1 for each of the three contributions. The performed simulations require substantial computational resources. The AL methodology does typically use O(10 5 ) CPU hours for each case. The AD methodology does in general require a factor 10-20 less due to the increased time step and decreased resolution.   The present analysis is an extension of the previous work on the inherent variability of the flow statistics in LES as presented by Andersen et al. (2015). The long term average velocity within large wind farms is expected to converge towards a constant level deep inside the wind farm, where a balance between the extracted energy and the entrained energy is reached. However, as shown by Andersen et al. (2015) the distributions of instantaneous and even 10 min average velocities show significant variability within the same simulation. Here, the focus of the present study is on mechanical power, as opposed to the electrical 210 power, which requires estimation of the electrical losses in for instance the generator. Hence, the power production calculated as: where T is the torque and ω is the angular velocity.
First, the inherent variability of LES are described, before the effect of free stream turbulence intensity, of turbulence and 215 shear combined, as well as of turbine spacing is investigated using the different numerical setups. Finally, the large amount of data is aggregated, and a more generalized analysis is performed on mechanical power production and variability within large wind farms.

Variability of LES
Simulations DT U 3, U U 2, and F W 5 (cf. Tables 2, 3, and 4) are comparable in terms of turbine spacing as well as freestream 220 velocity and turbulence intensity at hub height, although the difference in methods yields difference in the vertical profiles of velocity and turbulence intensity. Box plots based on the 10 min average mechanical power production normalized by rated power P 0 of the first 16 turbines for are given in Figure  spacing, as also discussed by Andersen et al. (2017b). These finding are furthermore corroborated by the recent experimental study by Turner V and Wosnik (2020), which identified resonance related to the turbine spacing. The simulations performed by FW included 50 turbines, so the full spatial extent of the wind farm is given in Figure 3. The "anomaly" appears throughout 240 the wind farm with distinct peaks at turbines 7, 12, 16, 23, 30, 39, 42, and 45. Furthermore, the variability clearly increases towards the end of the wind farm, where the power production ranges from 0.13-0.20 of rated power for the NREL 5MW turbine.
The variability is investigated more specifically in terms of the turbulence intensity, shear and turbulence intensity, as well as turbine spacing in the following.

Effect of Turbulence Intensity
The simulations from DTU and UU utilize body forces to introduce ambient turbulence into the flow. This enables direct investigation of the isolated effects of changing the ambient turbulence by changing the forcing. Here, the distributions of instantaneous power production of the 16 turbines are compared directly in violin plots in Figure   4 for DT U 2 and DT U 3, i.e. identical setup except an approximate free stream turbulence of 3% and 15%, respectively. The 250 differences in the distributions are clear. An increase in freestream turbulence increases the mean level of power production due to increased energy entrainment. Initially, the distributions are also broader for the high turbulent case than for the low turbulent case, which appears Gaussian, in particular for the second turbine. The width of the distributions becomes more similar further into the farm, but the difference in median level is maintained. Similar trends were reported by Andersen et al. (2016). The effect of the controller is also clearly seen as the distributions are capped around P P0 ≈ 0.33, corresponding to the 255 turbine reaching the maximum rotational speed.

Effect of Shear and Turbulence Intensity
Turbulence intensity and shear are inherently linked in the simulations performed by FW, as a change in equivalent roughness yields different shear and turbulence profiles. Figure 5 shows violin plots of the instantaneous mechanical power production in F W 2 (T I = 3%) compared to F W 5 (T I = 10%) for the first 16 turbines normalized by rated power. The distribution is 260 once again significantly broader for the high turbulence case and the distribution for the second turbine in the lower turbulence case is close to Gaussian. However, the median level appears to be very similar for the following turbines (3-6) with infrequent  higher tails. Further into the farm, the distributions become broader for the high turbulent case with a slight increase in the median level. However, the increase in the median level is not as pronounced as in Figure 4, which indicates that high shear decreases the effects of an otherwise high turbulence intensity.

Effect of Spacing
The initial turbulence intensity and shear discussed in the previous sections develops through wind farms, and the flow development is closely related to the turbine spacing. Fig. 6 shows violin plots of the instantaneous mechanical power production in DTU5 compared to DTU7 for the first 16 turbines normalized by rated power. This allows comparing spacings of respectively 12R × 12R and 14R × 14R. This does however also relate to the turbulence intensity and shear inside the farm, see sections of 12R × 12R is more irregular and seems to consist in two distinct parts, presumably due to how the turbine controller reacts to being in the near wake.

Comparison to Simple Engineering Models
The simulation data are aggregated in terms of 10 min statistics for each operating turbine. Aggregating the statistics from The results are generally encompassed by the results of DTU and FW. All results fall within a clear limit showing how much power can be extracted from a wind farm operating below rated wind speed depending on representative turbine spacing. An upper limit is indicated by DT U 4 and DT U 6 (in gray), which have a freestream velocity above rated (15 m/s), but with different turbulence intensities. The power productions deep inside the farm result in below rated conditions for DT U 6 due to 295 no freestream turbulence, while the turbines in DT U 4 also experience above rated velocities deep inside the farm due to the increased entrainment from the large atmospheric turbulence. Hence, it shows the transition from below rated to above rated conditions.
The effect of atmospheric turbulence is also clear, both when comparing the general trends of the plot and when comparing the DTU and FW results for different turbulent intensities. A higher atmospheric turbulence yields a higher production deep 300 inside the farms, while low or even no atmospheric turbulence results in a lower boundary in terms of production.
Finally, the figure includes the resulting power production based on two asymptotic expressions derived by Jensen and by Frandsen (1992), respectively. The Jensen model is widely used, also by industry, although it is less physical as it is not based on a proper momentum analysis. The model yields a velocity ratio given by: Here, r 0 is the turbine radius, x 0 is turbine spacing and hence, the original Jensen model only has a single input parameter, α, which governs the wake decay and expansion. The recommended values of α ≈ 0.04 for offshore wind farms, see e.g.  The model developed by Frandsen is on the other hand more physical, and involves more parameters, which are interlinked.
Here, the expression given in Frandsen and Madsen (2003) is used, which gives the converged velocity at hub height The geostrophic wind (G) and the roughness length (z 0 ) have an impact on the velocity at a given height. Hence, the geostrophic wind has been calibrated to give a mean wind speed of 8 m/s at a hub height of 90 m for two realistic roughness lengths 315 corresponding to turbulence intensities of 3% and 15%. A latitude of 55 • is assumed and a modified parameter of A * = 4 is used for compute f = 1.2 · 10 −4 · exp(A * ). The geostrophic wind and roughness lengths are summarized in Table 5. The converged velocity is then found using C T = 0.8 for various distances. The mean power production ratio is then computed using the cube of the converged velocity ratio and assuming constant C T .    Stevens et al., 2015Jensen, 1983: Frandsen and Madsen, 2003 z0

Response Surfaces
The total amount of aggregated data in Figure 7 comprise 12, 016 different, albeit overlapping, 10 min realizations, which includes the variability, both within a given 10 min realization and between different 10 min realizations as shown previously.
The power per ground area, or power density, compared to the standard deviation of power normalized by the mean power 330 for different relative spacings is shown in Figure 8, where each dot is a 10 min realization. The green dots show results with low atmospheric turbulence(T I = 0 − 3%), while the red dots show results for high turbulence (T I = 10 − 15%). The bin averaged data is shown in black with the standard error plotted as error bars. The standard error is here defined as: i.e. the standard deviation of the power density within a given bin normalized by the square root of the number of observations.

335
The data show significant spread in both power per area and standard deviation of the power although all simulation results generally cluster together. The binned values are generally very consistent except at low standard deviations, in particular  Figure 8a), which show large standard deviations of the binned data. Furthermore, it appears that the power production per ground area is not very influenced by the standard deviation normalized by mean power.

340
The high atmospheric turbulence realizations (in red) generally results in larger power density, as expected. Interestingly, the low atmospheric turbulence realizations yields larger larger standard deviations of the power within the realization, see Figure   8a), e), and g).
A multiple linear regression is applied to the full set of bin averaged data with a freestream velocity of 8m/s from Figure   8, i.e. aggregating all data with comparable C T . The regression is fitted in a least squares sense using the Matlab function 345 "regress" 1 . The regression fits the bin averaged power production per area to the normalized standard deviation of the power production and the relative turbine spacing. The fit is performed to second order, i.e. for combinations of S * = √ S X × S Y and σ * = P 10min S X ×S Y of the following matrix: The fit gives the coefficients b, and the combined set yields a response surface of the fit.
350 Figure 9 shows a contour plot of the response surface. The power density found here for a freestream velocity of 8m/s is in the range of 0.5 − 2.0W/m 2 , which is comparable to the general range of 1 − 11W/m 2 reported by Denholm et al. (2009).
The power density clearly decreases when the relative turbine spacing increases, as expected, because although the power production increases for larger spacing, the area increases faster and hence dominates the ratio. However, it is also clear how the power density varies with the standard deviation of the power production, i.e. how much power are the turbines able to 355 exploit and extract from the turbulent fluctuations. For large spacing, the power density is not influenced significantly by the standard deviation of the power production, cf. Figure 8g). For smaller spacing, there is an increased power density for small standard deviations in the power production, albeit related to the aforementioned small clusters of increased power density for small spacing, particularly seen in Figure 8a). Figure 9 also includes circles indicating the binned data used for the fit. The circles are colored according to the difference 360 between the fit and the binned data. The difference is generally (87% of the binned data points) less than ±0.5W/m 2 and alternating between a positive and negative difference for different relative spacings. The fit is particularly good for larger spacings, but it struggles for smaller spacings with large outliers. The sensitivity of the fit is examined by performing a 10-fold cross-validation. The mean squared error (MSE) is estimated between the fitted response surface and the actual input data. The resulting MSE is 0.0732 ± 0.0011 indicating that the fit is consistent across the data, but that an improved response surface 365 could be created. However, this would presumably also require a larger amount of data.
. Contours based on multiple linear regression fit of bin averaged power production per area to the standard deviation of power production normalized by mean power production and relative turbine spacing. Points mark the binned data, where blue and red shades indicate whether the fit underestimates or overestimates compared to the binned averaged values, respectively.
As shown previously, the simple engineering models by Frandsen (2007) is capable of capturing the average trends, similar to the response surface. However, the inherent variability of LES is important for farm performance and for improving risk assessment during the design phase. Hence, a similar response surface can be fitted to the standard deviation and the standard error of the bin averaged values in Figure 8. The corresponding response surfaces are shown in Figures 10 and 11, which can 370 be interpreted as the variability and the uncertainty associated with the response surface of mean power density.
The variability around the mean is up to 0.4W/m 2 for small relative spacings, which is comparable to the difference in the fit and bin averaged data as shown before. The variability is higher for shorter spacings, where the amount of outliers affect the fit. The outliers can be related to the significant non-linearities in the near wake before the wake breaks down into small scale turbulence, see Sørensen et al. (2015). The 10-fold cross validation on the variability yields M SE = 0.0276 ± 0.0014, which 375 indicates a significant uncertainty in the fit, but again very small variation within the given dataset.
The increased variability for smaller spacing also comes with an increased uncertainty as shown in Figure 11. The standard error decreases for increasing spacing, where the fit is very good, while the discrepancy is larger for the very short distance.
The fit tends to overestimate the standard error for the shortest spacings. Again, a 10-fold cross validation is performed on the response surface for the uncertainty, which gives M SE = 0.0085 ± 0.0015. Again, a very small variation within the given 380 dataset, but that the chosen response surface could be improved.
The response surfaces are only fitted to second order, because the aim here is merely to provide general insights into the global trends and hence to avoid overfitting. It should be strongly emphasized that this is a rather crude approach. However, the response surfaces yields a first attempt at constructing global response surfaces of the power density including the inherent variability based on significant amounts of LES data for a wide range of wind farm layouts operating at 8m/s, which show 385 physical trends. Figure 12 shows a heatmap of how the 10 min realizations are distributed for the standard deviation of power production within each 10 min period normalized by the corresponding power production versus the relative turbine spacing. Clearly, the majority of 10 min realizations are in the range of σ(P10min) P 10min The presented response surfaces are directly dependent on the data availability and hence most representative in this range 390 and care should therefore be taking in the ranges with few realizations. However, Figure 12 can be used to guide which scenarios should be computed next essentially to fill the gaps. As such, additional LES computations should be focused around Figure 11. Contours based on multiple linear regression fit of the standard error of the bin averaged power production per area to the standard deviation of power production normalized by mean power production and relative turbine spacing. Points mark the binned data, where blue and red shades indicate whether the fit underestimates and overestimates compared to the binned averaged values, respectively.
where there is no data, and for √ S X × S Y < 14R, where the uncertainty is large. The response surfaces can therefore continuously be improved by adding more data and a more sophisticated response surface or surrogate can be constructed once the entire parameter space is full.

395
The response surfaces could also be made dependent on more parameters by adding more LES data. Currently, the turbine spacings in the lateral and streamwise direction have been collapsed to a single dimension, but the dependency could be unfolded. Similarly, the dependency on a number of additional parameters could be investigated, for instance: -Free wind speed However, this would obviously require substantial amounts of computing resources.
One way to circumvent the large computational costs would be to utilize SCADA data in combination with the LES. Similar 405 response surfaces could be constructed based on SCADA data from operating wind farms, which would enable a more global verification of LES and the actuator disc/line methods on wind farm scale. Such a verification would be valuable as direct comparison of time series of specific events between LES and actual wind farms is extremely difficult, if not impossible, to achieve given the complexity and amount of information required on the atmospheric conditions to enable such a comparison.
A successful verification would facilitate the direct integration of LES data and SCADA data to construct more certain 410 response surfaces covering a larger range of scenarios and parameters. It could act as a lodestar and inform researchers in which regions of turbine spacing and turbulence intensity to perform the expensive LES in order to fill the gaps and explain physical trends not captured by the simpler models.
Finally, the response surface could be extended to include e.g. fatigue loads for turbines operating in wind farms. Such a surrogate model for fatigue loads on a single wind turbine was developed by Dimitrov et al. (2018), who compared the accuracy This work aimed at providing a general overview of the global trends of power performance for large wind farms, with a focus on variability. This was done through the analysis of Large Eddy Simulations (LES) performed on large wind farms from the three institutions that co-authored this work. LES results of large wind farms obtained from other researchers as well as simu-420 lations performed using simpler engineering models were also included to provide a more complete envelope for the results.
As LES require large amounts of computational resources, emphasis was made on extracting as much information possible from the existing set of simulations performed using different setups and incoming flow conditions. As such, emphasis is not put on comparing the simulations to each other, but rather on using as many results as possible to cover a wide range of possible 425 scenarios that can provide a global picture of the power characteristics within large wind farms.
Parametric studies were first performed to inform about the effect from atmospheric conditions as well as turbine spacing on the production and its variability. An increase in atmospheric turbulence intensity, by increasing energy entrainment, was shown to raise the mean level of power production. It was also associated to wider distributions of the production values. A 430 larger spacing between the turbines was also associated to greater levels of production, as expected.
The analysis was extended further by aggregating the large amount of LES performed under various conditions. This was done in terms of 10-minute statistics for each turbine operating in deep farm conditions. LES works from other researchers as well as simulations performed with simpler engineering models were also included in a first step when looking at the power produced deep inside the farm as a function of a representative spacing. All results were shown to fall within a clear limit 435 showing how much power can be extracted from a wind farm operating below rated wind speed, as a function of representative turbine spacing. Whereas higher turbulence levels lead to larger production levels deep inside the farms, while cases without incoming turbulence were shown to provide a lower power production. While LES provide more information in terms of variability, simple engineering models were shown to produce a reasonable envelope for the results obtained using the high fidelity methods.

440
As a second step, response surfaces encompassing the total amount of aggregated LES data, i.e. 12,016 different albeit overlapping 10-minute realizations, were created. They revealed information regarding various aspects of the power production within large wind farms, among which the amount of power the turbines are able to extract from the turbulent fluctuations, as well as the variability and uncertainty associated with the mean power densities.

445
The work presented in this paper serves to provide valuable information regarding power and its variability deep inside large wind farms. Nonetheless, the response surfaces presented here would gain in being complemented with more LES results to provide an even more complete picture. This could be done by considering further turbine spacings to fill existing gaps. The dependency of response surfaces to more parameters could also be investigated, including individually-considered spanwise 450 and streamwise spacings, the freestream velocity as well as the atmospheric stability. As LES are known to be very computationally demanding, SCADA data could also be used to provide more complete response surfaces. Future work could also go one step further by investigating the behavior of turbine loads in similar terms as what was performed here regarding power production.