Medium fidelity modelling of loads in wind farms under non-neutral ABL stability conditions – a full-scale validation study

The aim of the present paper is to demonstrate the capability of medium fidelity modelling of wind turbine component fatigue loading, when the wind turbines are subjected to wake affected non-stationary flow fields under non-neutral atmospheric stability conditions. To accomplish this we combine the classical Dynamic Wake Meandering model with a fundamental conjecture stating: Atmospheric boundary layer stability affects primary wake meandering dynamics driven by large turbulent scales, whereas wake expansion in the meandering frame of reference is hardly affected. Inclusion of stability (i.e. buoyancy) in description of both large- and small scale atmospheric boundary layer turbulence is facilitated by a generalization of the classical Mann spectral tensor, which consistently includes buoyancy effects. With non-stationary wind turbine inflow fields modelled as described above, fatigue loads are obtained using the state-of-the art aeroelastic model HAWC2. The Lillgrund offshore wind farm (WF) constitute an interesting case study for wind farm model validation, because the WT interspacing is small, which in turn means that wake effects are significant. A huge data set, comprising 5 years of blade and tower load recordings, is available for model validation. For a multitude of wake situations this data set displays a considerable scatter, which to a large degree seems to be caused by atmospheric boundary layer stability effects. Notable is also that rotating wind turbine components predominantly experience high fatigue loading for stable stratification with significant shear, whereas high fatigue loading of non-rotating wind turbine components are associated with unstable atmospheric boundary layer stratification.


Introduction
For wind farm (WF) production estimation stationary WF flow field modeling as provided by e.g. full CFD RANS models or fast linearized CFD RANS models [1] may suffice. However, for load estimation of wind turbines (WT's) exposed to wake affected inflow conditions, a non-stationary WF flow field description is inevitable. Insisting on models reflecting the basic physics of the problem, and considering direct numerical simulation (DNS) precluded due to its excessive computational demand, we are thus left with either high-fidelity CFD LES actuator disc or line a type of approaches or, alternatively, the medium-fidelity Dynamic Wake Meandering (DWM) [2], [3] type of model approach. CFD LES models must -like the DWM model -be linked to an aeroelastic model of each and every WT in a WF to provide a complete picture of the load conditions inside a WF for all design load cases. This is extremely CPU demanding and considered unrealistic for design purposes even with the capacity of nowadays very powerful state-of-the-ark super computer clusters. This challenge is further enhanced when including atmospheric boundary layer (ABL) stability as an additional design load case dimension and/or considering optimization of WF layout.
Being a medium-fidelity approach the DWM model offers significant savings in CPU-demands compared to the above mentioned non-stationary CFD alternatives, and it thereby facilitates detailed WF production/load design simulations to be conducted even when including the extra design load case dimension constituted by ABL stability. The "classical" (i.e. associated with neutral ABL conditions) DWM model is well established and about to be included in the revised IEC standard as a recommended practice. The model has been extensively validated both directly against full-scale wake affected flow measurements [5], [6], [7] and indirectly [8], [9] by using DWM together with the aeroelastic code HAWC2 [10] to compare predictions of WT production and loading with full-scale data. In [9] a huge scatter in the measured load data was observed, and the motivation for the present study is thus to investigate whether at least part of this scatter can be attributed to ABL stability effects as well as to get qualitative and quantitative insight in ABL stability effects on WT loading in WF's.
The core of the DWM model is a split of scales in the wake affected flow field, with large turbulence scales being responsible for stochastic wake meandering and small scales being responsible for wake attenuation and expansion in the meandering frame of reference as caused by turbulent mixing. Thus, essentially the DWM model assumes that the transport of wakes in the ABL can be modeled by considering the wakes to act as passive tracers driven by a combination of large-scale turbulence structures and a mean downstream advection velocity, adopting the Taylor hypotheses [4]. This basic split in scales further facilitates generalization of the DWM model to non-neutral ABL stability conditions in a simplistic way, as ABL stability is known primary to affect the low wave number turbulence regime.
The paper is structured as follows. First the generalization of the DWM model to non-neutral ABL conditions is briefly described in Section 2; then follows in Section 3 a description of the Lillgrunden case study, which previously have been used in an extensive validation study but, however, limited to neutral ABL conditions [9]. In Section 4 results from comparing model simulations with the Lillgrunden full-scale measurements are presented. Finally, conclusions are drawn in Section 5.

DWM for non-neutral conditions
The large scale turbulence structures, used to describe wake dynamics in the DWM model, are traditionally provided by a fast Navier-Stokes (NS) consistent turbulence model [11] which, however, assumes neutral ABL stability and thereby neglect the effects of buoyancy on turbulence production and thus turbulence characteristics. As ABL stability mainly affects the large scale turbulence structures [4], the effect of buoyancy on wake meandering cannot be neglected. In addition, the significance of buoyancy effects on small turbulence structures relevant for wake deficit expansion and attenuation in the meandering frame of reference must be clarified. The crucial question in this respect is whether the effect of buoyancy can be "decomposed" in analogy with the description of wake dynamics.
The starting point is the following fundamental conjecture launched in [12]: In a wake context, ABL stability affects primary wake meandering driven by large (lateral and vertical) turbulent scales, whereas wake expansion in the meandering frame of reference is a second order effect only. This conjecture has been investigated for a variety of stability conditions in [13] and [14], respectively. Based on full-scale LiDAR and sonic measurements as well as detailed numerical simulations the wake affected flow field behind a 500kW turbine was examined with focus on both organized wake deficit flow structures (resolved in the meandering frame of reference) and wake dynamics. These extensive investigations have justified the above stated conjecture, and we will therefore generalize the applicability of the DWM model to non-neutral ABL stability regimes by adopting this conjecture.
For consistency, we need a generalization of the classical spectral tensor to the non-neutral regime. Such a model has been developed and validated against full-scale data in [15] and [16]. Being a generalization of the classical Mann spectral tensor, this model share some of the basic features of the Mann spectral tensor (e.g. linear mean wind speed shear and eddy life time model), but includes a linear vertical temperature profile as well as buoyancy terms in the linearized NS equations used in the Rapid Distorsion approach. The inclusion of buoyancy comes with a cost of two additional input parameters. The generalized spectral tensor is used to model both the large-and small scale ABL turbulence used with the DWM model.

The Lillgrunden case
The Lillgrund WF consists of 48 Siemens SWT-2.3-93 WT's, and one of these (C-8) is instrumented with strain gauges resolving blade, main shaft and tower loads, respectively. These measurements have been made available by Siemens Wind Power, and they constitute probably one of the most comprehensive sets of wake affected wind turbine load measurements ever recorded. The measurement period extends from 2008-06-03 to 2013-03-19 -i.e. over a period of almost 5 years. The layout of the WF is shown in figure 1. Whereas the Egmond aan Zee WF, previously investigated in [8], is characterized by a "conventional" WT inter spacing, the layout of the Lillgrunden WF is characterized by very small WT inter spacing's -i.e. down to 3.3 D. This makes the Lillgrunden WF especially interesting and challenging as a wakeload validation case, because the close spacing magnifies wake generated load effects compared to more traditional spaced WF's.

Measurements
WT C-8 is instrumented with strain gauges resolving respectively the blade root flap bending moment and the tower base for-aft bending moment. In addition to these high-frequency data, WT SCADA data (pitch setting, rotational speed and electrical power) were available for the WF WT's during the measuring period. Unfortunately no meteorological mast data was available within the recording period, but as the WT power and pitch angle are directly correlated with the inflow wind speed, the ambient undisturbed wind speed has been determined as based on power and pitch angle recordings for corner placed WT's located in free inflow conditions. A similar philosophy was used for estimating the ambient undisturbed wind direction, which was determined from nacelle orientations of corner placed WT's. These nacelle orientations were initially calibrated against power deficit polar's constrained to directions, where wake losses were dominated by the closest neighboring WT's. Full polar load cases, associated with normal WT operation, are available for mean wind speeds ranging from 6m/s to 16m/s, and the data have been classified into 2m/s velocity bins. The blade root flap moments and the tower bottom for-aft moments have been post processed to fatigue equivalent moments using the Palmgren-Miner [17] approach and subsequently normalized with the respective fatigue equivalent moments associated with a inflow wind speed of 9m/s -i.e. here represented by mean equivalent moments associated with the velocity bin [8;10]m/s. Wöhler exponents of 5 and 10 were assumed for the tower and blade composite structures, respectively.

Simulations
For the DWM/HAWC2 validation study the load response of WT C-8 is simulated for mean wind speeds reflecting the median of the defined velocity bins. Measured wind speed dependent turbulence intensities (TI's) are used, reflecting the offshore wind speed dependent "surface" roughness. However, no attempt is done to resolve TI as function of upstream fetch (i.e. direction). Thus, in the mean wind speed regime 6m/s-14m/s a TI of 5.8% is used -gradually increasing to 6.2% at 16m/s. Both the deterministic and the stochastic part of the inflow are affected by ABL stability conditions. The deterministic mean wind shear and the stochastic turbulence input are described in the following.

Mean wind shear
The mean wind shear is stability dependent through the buoyancy impact on mixing in the ABL layer -in stable conditions with little mixing the shear profile gradient tend to increase, whereas in unstable conditions with increased mixing the shear profile tend to become more uniform.
In the IEC-61400 code [18] the neutral wind shear is specified in terms of a power law approach as with U being the mean wind speed, z being the height above terrain, and z hub being the hub height. The power exponent α is specified to 0.14 for an offshore site [19]. However, the new (not yet approved) revised version of the IEC code opens alternatively also for use of a log-profile defined as with z 0 being the roughness length. The log-profile is physical consistent and used in Monin-Obukhov (MO) similarity theory [4] and, as seen, it takes the value 0 at z = z 0 and the value U hub at z = z hub . However, a recommended roughness length for offshore sites, consistent with the specifications of the power law exponent in equation (1), is not yet specified. We will therefore use equation (2) to introduce the effect of atmospheric stability in the power law description in an approximate manner.
As a first step we rewrite equation (2) as where the last identity provides the link to the conventional description of the log-profile, with * being the friction velocity, and κ being von Karman's constant (= 0.4), viz.
To establish a connection between the IEC specified power coefficient, α, and the roughness length, z 0 , we will require the best possible agreement of respectively the power law profile and log profile over the vertical extent of the rotor, and we will define this best fit in terms of minimum least square deviation. For this purpose we define the functional Π as and define z 0 accordingly as the particular value, 0 z  that minimize the above functional, viz.
Adopting M-O similarity theory, the mean wind shear in stable conditions is traditionally formulated as [4] where L M is the M-O length, z i is the height of the atmospheric boundary layer, and ψ m (s) is the stability function associated with stable stratification. Based on the reflections concerning the neutral shear profile, we will now "transform" relation (7) into a power law formulation for the stable ABL condition using eq. (6). The resulting expression is where the value of the roughness length, z 0 , will relate to the power law exponent as described in equations (5) and (6). Again adopting M-O similarity theory, the mean wind shear under unstable conditions may be formulated as [4], [20] Taking the same approach as for the stable stratification, we arrive at the following approximate expression when taking a power law approach where, as in the stable case, the value of the roughness length, z 0 , relates to the power law exponent as described in equations (5) and (6).

Turbulence
Using the coupled DWM-HAWC2 platform [3] requires simulation of three different turbulence fields -a large scale turbulence field dictating the wake meandering, a traditional turbulence field accounting for conventional WT turbulence loading, and an isotropic and inhomogeneous small scale turbulence field accounting for WT loading caused by wake self generated turbulence. The wake self generated turbulence field is assumed to be invariant with respect to the ABL stability condition, and it is simulated using the classic Mann turbulence simulator [11]. The two remaining turbulence fields, which are highly dependent on ABL stability conditions, are generated based on the generalized buoyancy dependent spectral tensor [15], [16], which degenerates to the classic spectral tensor for neutral stability conditions. In all cases the turbulence fields are simulated as three dimensional fields resolved in suitable Cartesian grid configurations. The parameters of the buoyancy dependent spectral tensor are, for each stability class, obtained from fitting model auto-and cross spectra to respective spectra obtained from full-scale sonic measurements from the Høvsøre site in Denmark. Fits are performed based on data recorded at 60m and 40m altitudes, respectively. For both altitudes, data associated with the [8;9]m/s mean wind speed bin were used for the parameter fitting. However, the fits at 40m seem to fit slightly better the fullscale data compared to the 60m data in especially the neutral and the unstable regimes. As 40m altitude moreover is more likely to be within the surface layer, where L M is defined, than the 60m altitude (especially for stable stratification), it was decided to base input parameters to the generalized spectral tensor on the 40m data fits.
Except for the turbulence intensity (described by the αε 2/3 parameter of the spectral tensor when the remaining parameters have been fixed), it was further decided to adopt the (neutral) turbulence input specifications from the IEC code to [18] mimic the turbulence conditions at the Lillgrunden site, where no high-frequency meteorological measurements, as mentioned, are available. The IEC spectral tensor input parameters were originally obtained by fitting the spectral tensor spectra to a target Kaimal spectra [18].
Since we have decided to base the neutral turbulence generation partly on IEC specifications of turbulence tensor input (i.e. turbulence length scale L and eddy life time parameter Γ), the spectral tensor input parameters for neutral conditions do not match the fitted parameter values from the Høvsøre full-scale data. As a consequence, it was decided to scale (L, Γ) for the non-neutral (i.e. diabatic) conditions accordingly, although this approximation might not be completely true. However, it is believed to be a fair approximation -alternatively the neutral case could defined by the directly fitted values for L and Γ at the Høvsøre site, which is just another approximation to the Lillgrunden conditions.

Results
The present analysis deviates from the analysis presented in [9] by the inclusion of ABL stability affected turbulence and shear in the load simulations. Two non-neutral stability cases, corresponding to Richardson number Ri = -0,015 (unstable ABL stratification) and Ri = 0,01 (stable ABL stratification) have been simulated. Based on such numerical predictions, conducted for various ABL stability conditions, it is possible to get insight in how ABL stability affects the WT load levels, as well as to what extend ABL stability can explain the rather high scatter level in the measurements. Results are shown in figure 2 for below rated wind speed WT operation and in figure 3 for above rated wind speed WT operation.  Table 1, simulations are conducted for unstable (-2), neutral (0) and stable (2) ABL stability conditions using the linear perturbation deficit wake merging operator.
For all wind speeds it can be seen that the tower bottom fatigue bending moments are largest in unstable and smallest in stable ABL conditions. This is easily explained by the difference in turbulence intensity associated with the stability conditions, where unstable conditions causes a turbulence level of 9.8%, neutral conditions corresponds to a turbulence level of 6.0%, and stable conditions corresponds to a turbulence level of 1.99%.
In the inflow sector 150deg.-180deg., the investigated WT is not affected by upstream wakes, and it is thus possible to evaluate the impact from ABL stability on solitary WT loading. When observing the blade flapwise fatigue loads in this free sector, the results are significantly different than the tower fatigue loading. Contrary to tower loading, the blade loads are sensitive to both shear and turbulence levels, as both cause varying wind speeds over the rotor area. For rotating WT components these are counter acting effects, because stable conditions result in large shear and low turbulence intensity and vice versa for unstable conditions. For most wind speeds, it is seen that both stable and unstable conditions result in a larger fatigue loads than the neutral case. The only exception is at very low wind speed, where neutral and stable conditions result in approximately identical fatigue load levels, whereas larger fatigue load levels are observed for unstable conditions. Similar results were obtained in [21]. In the wake affected conditions the same trend, with increased tower fatigue loads in unstable conditions and decreased tower fatigue loads in stable conditions, are seen. Also for wake affected inflow conditions the neutral situation seems to cause the smallest blade fatigue loads. However, for single wake situations it is interesting to notice, that stable ABL conditions result in both the largest fatigue load level (i.e. in a half-wake situations) and the smallest load level (i.e. in the full wake situation) compared to unstable and neutral ABL conditions. Finally, it is noted that a significant part of the observed measurement scatter potentially can be explained by ABL stability effects, especially remembering that very stable and very unstable stability conditions are not included in the present analysis.

Conclusion
A medium-fidelity DWM/HAWC2 framework for modelling of wind turbine component fatigue loading, when the wind turbines are subjected to wake affected non-stationary flow fields under nonneutral atmospheric stability conditions, has been developed and demonstrated on full-scale load data from the Lillgrunden WF.
A comparative analysis of measured and simulated blade and tower fatigue loads has revealed that ABL stability effects can explain a significant part of the large measurement scatter observed in fullscale data. The study has further gained insight in important load mechanism, and shown that mean wind speed shear is the dominating fatigue load driver for the rotating WT components in stable ABL conditions, whereas turbulence (including its dictating wake meandering mechanism) is the dominating load driver for tower fatigue loading in general.
After conclusion of this study wind speed and temperature recordings, covering the major part of the measurement period, have been made available from the nearby Drogden lighthouse. In a future study, it is planned to use these measurements to classify the Lillgrunden measurements into stability classes, thus facilitating a one-to-one comparison of ABL stability affected measurements and fullscale data.