The consideration of atmospheric stability within wind farm AEP calculations

The annual energy production of an existing wind farm including thermal stratification is calculated with two different methods and compared to the average of three years of SCADA data. The first method is based on steady state computational fluid dynamics simulations and the assumption of Reynolds-similarity at hub height. The second method is a wake modelling calculation, where a new stratification transformation model was imposed on the Jensen an Ainslie wake models. The inflow states for both approaches were obtained from one year WRF simulation data of the site. Although all models underestimate the mean wind speed and wake effects, the results from the phenomenological wake transformation are compatible with high-fidelity simulation results.


Introduction
The atmospheric stratification largely affects the flow conditions in the atmospheric boundary layer (ABL) relevant for both off-shore [1] and on-shore [2] wind energy sites. The atmospheric stability is driven by temperature gradients and buoyant forces and has strong impact on the generation of turbulence and thus the wind turbine wake recovery. It also determines the amount of mixing between different layers of air, and hence it is crucial for the transport of energy into the wind farm from above (and below).
Stable conditions often imply longer and more pronounced wind turbine wakes, due to the above mentioned reasons, unstable conditions have the reverse characteristics. In general, wakes affect both the energy yield and the loads of the wind turbines, hence also their lifetime and the operational costs of the wind farm. Thermal stratification is therefore linked to wind farm planning, and, more generally, to the calculation of the annual energy production (AEP).
However, stratification effects are often ignored for such calculations, and sometimes even the average of different stability situations is considered, due to the lack of better methods. This is the point that we aim to address in this work, we demonstrate that the effects can be represented in computational fluid dynamics (CFD) calculations as well as in the engineering model wake superposition approach.
In the following, we include the effect of atmospheric stability into the AEP calculation of an existing onshore wind farm in quasi-flat terrain for the year 2015. The time series of inflow wind and stability conditions are obtained from a mesoscale simulation with the Weather Research and Forecasting (WRF) model [3]. Afterwards, a classification based on wind speed, direction and the Obhukov length (atmospheric stability) is carried out.

The site
The on-shore site under investigation is located in the USA and consists of 17 Nordex wind turbines with hub height H = 100 m, rotor diameter D = 100 m and nominal power P = 2.5 MW. The layout is shown in Fig. 1. The terrain is quasi-flat and the site is dominated by agricultural land, which is modelled by a uniform roughness length of z 0 = 0.15 m.

Method
In this paper the results from one-year of mesoscale model simulations are used as boundary condtions for subsequent microscale and engineering model simulations to calculate the AEP for an onshore site. Therefore, this section subdivides into the description of the mesoscale, the microscale and the engineering model setup.

Mesoscale model simulations
The mesoscale model simulations presented in this publication were conducted with the WRF model [3]. The model is well established in research and industrial wind energy applications (e.g. [4], [5]). We used a setup similar to the one in [1] which was specifically chosen for wind energy applications and thus vertically well resolves the heights most relevant for wind turbines. The setup consisted of three domains (D1-D3) nested within each other. Fig. 2 shows the distribution of the domains centered around the wind farm location. The location of interest is located in eastern Iowa in mildly complex terrain. Tab. 1 gives an overview of the WRF version as well as the parametrisations and boundary conditions used. An entire year (2015) was simulated separated into four simulations each with a duration of three months using spectral 3D analysis nudging [3]. The first two days (48 hours) of each quarterly simulation were disregarded in the subsequent analysis to allow for enough spin-up time for the model. The results from the model calculations were put out every 10 min. Afterwards, the original wind data (on eta-levels) were interpolated to hub height and together with the dimensionless Obhukov length L 90 = H hub /M OL extracted from the models result files. Based on these three parameters a binning into bins of 1.0 ms −1 in wind speed, 10 • in wind direction and five stability classes in ∆L 90 (see Tab. 2) was performed.  . The roses generally demonstrated prevailing south south-westerly winds with a secondary maximum for north-westerly winds. Generally all wind direction sectors are of importance as none of the twelve sectors has a probability of occurrence of less than 5 %. Stably stratified conditions are dominating at the wind farm location.
The second postprocessing step for the WRF simulation data was a selection of the relevant states based on the distributions of these three parameters. These states then served as boundary conditions for the subsequent set of microscale simulations to determine the AEP at the wind farm location.

Microscale model simulations
All CFD results for this work were generated with OpenFOAM [7] (version 2.3.1) and in-house extensions by Fraunhofer IWES (e.g. [8][9][10]). Thirty-six wind directon sectors were simulated for one characteristic wind speed, one neutral condition and four dimensionless Obukhov lengths (see Tab. 2). The weights for adding simulation results are derived from the mesoscale   The terrain is quasi-flat and was generated with the in-house terrainMesher. The applied turbulence model is a buoyant adaption of the k--f P turbulence model [11] with parameters as suggested in the reference. The inflow boundary profiles for wind velocity U , turbulent kinetic energy k, turbulent dissipation and potential temperature Θ were consistently generated by a one dimensional cyclic precursor run. The roughness length is uniformly set to z 0 = 0.15 m, forest canopy can be ignored at the site. The solver is an in-house development for steady state stratified atmospheric boundary layer flow, based on the SIMPLE algorithm [12]. The wind turbines were modelled by uniformly loaded actuator disks with about 500 cells each, switched on after running the background wind field simulation first.
For the AEP calculation the hypothesis of Reynolds-similarity is applied. This means that from a single simulation in a certain wind direction sector with a characteristic wind speed at hub height of 10 ms −1 , the wind speeds in each cell for other inflow wind speeds are obtained by local re-scaling. The error that is made by this simplifying assumption can be expected to be smaller than the error that corresponds to a direct evaluation of the sector based on a single simulation. This assumption reduces the number of simulations for a wind rose with 12 sectors to 36 × 5 × 1 = 180 simulations.

Engineering model simulations
All engineering model results were generated with the wind farm modelling software flapFOAM [9,10,13,14]. flapFOAM is based on the principle of wake superposition, and some of its modelling components are briefly described below. For more details, the reader is referred to the references.
3.3.1. Wake models All calculations in this paper are performed with wake models that are based on the models by Jensen [15] and Ainslie [16]. For reasons that will become apparent later, they are filtered by flapFOAM's wake transformation option 'Gaussianize'. This transformation is based on the assumptions of axial symmetry and neglects radial velocity deficit components. It can be applied to arbitrary wake models and translates the deficit into three functions R(x), A(x) and r 0 (x) on a one dimensional grid, for which we choose a spacing of ∆x = 10 m. The function R(x) represents the wake radius in m at a distance x behind the disk, determined by the condition that the deficit has dropped below 1% of the central value.
The other two functions represent a Gaussian function with respect to the radial coordinate r, and they are interpolated at fourth order accuracy as functions of the continuous variable x. The amplitude A(x) is fixed by the central value (r = 0) of the underlying wake model, and the width of the Gaussian function by the requirement that the total mass flux of the Gaussian function should be identical to the calculated flux F (x) of the original wake model within a disk of radius R(x), yielding r 0 (x) = 2F (x)/A(x).
The increase of turbulence intensity in a single wake is modelled according to Frandsen [17]. For this work we chose quadratic summation of wind velocity and turbulence intensity deficits in multiple wake situations.

Inflow states and AEP calculation
The meteorological variables that are particularly of interest for the wake flow are the wind direction φ, the wind speed v and inverse dimensionless Obukhov length L 90 at hub height. With binning ∆φ = 10 • , ∆v = 1.0 ms −1 and L 90 as in Tab. 2 we obtain 2106 states with non-vanishing probability from the mesoscale results. For each of these states one calculation is performed with flapFOAM, and the AEP is then obtained by taking the weighted sum of the predicted wind farm power. For each state and each turbine the latter is obtained from the power curve of the wind turbine model, evaluated at the effective local wind speed as measured by the turbine. For simplicity we consider the wind speed at the disk centre only and ignore partial wake and wind shear effects.
The ambient wind field is modeled by standard log profiles with uniform roughness length z 0 = 0.15 m and thermal stratification corrections.The ambient turbulence intensity I is obtained from the Richards-Hoxey solution [18] of the k-turbulence model where C µ = 0.033 and κ = 0.4 are constants. Notice that stability corrections of turbulence intensity are ignored in the engineering modelling approach throughout this work.

Thermal stratification effects on wakes
A new phenomenological model for the effect of thermal stratification on wakes has been proposed in [19] and implemented in flapFOAM. The model has been obtained by fitting parametrised curves to central wind speed deficit results of CFD RANS actuator disk simulations. For eight different values of the inverse Obukhov length parameter ξ a single uniformly loaded actuator disk was simulated with methods from Section 3.2. The deviation of the central wake deficit in the stratified situation compared to the neutral case for identical inflow wind speed U 0 at hub height is Furthermore the actuator disk simulations did not show a strong variation of the width of the wake [19]. In flapFOAM the model is realized by two steps. First, the chosen base wake model is transformed into a Gaussian function, as described in Section 3.3.1. Then, the value U 0 f (x/D, ξ) is added to the corresponding amplitude and the scale r 0 (x) is not modified. Examples for the Jensen and the Ainslie wake models and an Obukhov length L = ±200 m are shown in Fig. 6.

Results and discussion
The to the simulations corresponding year (2015) of SCADA data are available for the site from Section 2. The data was filtered according to a status flags of the turbines. The resulting average wind speed and the AEP of the wind turbines is compared to simulated results in Figs. 7-8 and 9, respectively. In figures 7 and 9, the results have been normalized to the average of each simulation, while in 8 the normalisation was done with respect to the average of the SCADA data to allow for an estimation of the bias between the simulation methods. For flapFOAM both neutral and stratified results are shown, for Jensen and Ainslie wake models. Generally, the variations in the wind speed at the turbines are much lower in case of the flapFOAM results compared to the SCADA data. One reason for this is that the variation in the wind speed at the measurement location (downstream of the rotor disk) is generally very strong. CFD and flapFOAM seems to underestimate the wind speeds at the locations of the turbines, which is also visible in flapFOAM results (Fig. 8). However, another possibility is that the SCADA wind speed is biased towards high values.
The lowest mean wind speed and AEP results according to the SCADA data are found for turbines 2, 5, 6 and 12. From the wind farm layout in Fig. 1 and the wind rose from Fig. 3 it is apparent that these turbines are subject to wakes from southerly winds, which are well represented in the wind rose. The wind turbines 4, 15, 16 have the highest average wind speeds and AEP, they are neither hit by southern nor by north-easterly winds or (in case of turbine 4) the next upstream turbine is several tenth of rotor diameters upstream. flapFOAM is generally capable in simulating this variation between the turbines. Here, turbines 9 and 16, which are both located at the edges of the wind farm show the largest wind speeds and AEP.
None of the models, CFD or engineering, seems to able to convincingly capture the complexity of the wake situation at the site. This is likely due to the relatively coarse wind direction resolution of 10 • . Since some wake situations are missed by 36 wind directions, the non-wake situations statistically dominate the results for most of the turbines and no clear picture emerges.  Figure 9. AEP calculation results normalized by the average of each data-source.
The difference between the original Jensen and Ainslie wake models and their thermally transformed counter parts is around 5% or less. However, the stability corrected wake models yield better results compared to measurements. The small impact of the stability transformation of the wake models on the resulting AEP is due to the small wake effect that was observed in all simulations. Since for each state the wind speed at hub height was fixed to the WRF result, and since no shear effects have been considered for the calculation of the effective wind speed at the rotors, the wake effect exclusively enters the flapFOAM calculation through the imposed wake deficit transformation in the chosen setup.

Conclusions
In this work we studied an existing wind farm in North America and compared two simulation approaches for the AEP calculation with averaged one year SCADA data. The two methods are steady-state CFD and the engineering model software flapFOAM. The input for both was based on the same results of a one-year mesoscale simulation of the site with WRF. By imposing a binning in the three variables wind direction, wind speed and dimensionless inverse Obukhov length, and by sorting and filtering according to probability, the number of states was reduced to 2106 states.
For CFD this number is still too high, it was reduced to 180 by assuming Reynolds-similarity for the local wind speed at hub height. However, the results show an underestimation of the wind speed and thus the AEP. Also, variation of the SCADA results due to the wake effect are not captured well. Since this is also observed for the engineering model results we conclude that the sectors of the wind rose were most likely not sufficient for the observation of significant wake effects.
For the very same reason we do not observe a large impact of the presented engineering model for modelling thermal effects on wake deficit models. However, the transformed Jensen and Ainslie wake models have better agreement with the SCADA AEP results than CFD, which underpredicts the wind speed. The stratification transformation of arbitrary wake effects was obtained by fitting the variation of the central wake deficit with Obukhov length to corresponding CFD actuator disk results. This process yielded a phenomenological model that can now be used for fast AEP calculations that include thermal stratification. This model was first applied to an existing wind farm in this work and compared to measurements. We conclude that further validation and more detailed simulations are necessary, also the dependency on roughness length and the performance of the model in complex terrain remain to be investigated.