Simplified wake modelling for wind farm load prediction

This paper presents a simple numerical wind farm model, where pragmatic choices are made in the modelling of underlying physical processes, with the aim of making useful power production and wind turbine load estimates. The numerical model decomposes the wind farm, inspired by the approach of the dynamic wake meandering model (DWM), into simple sub-models for a single wake deficit (1D Gaussian), wake meandering (statistical), and wake added turbulence (eddy viscosity based). Particular attention is given to selecting a momentum conserving wake summation method, because of its critical role in coupling the influence of individual wakes. Results are presented to illustrate the influence that wake summation methods have on equilibrium velocity and power production in a row of turbines, for different inter-turbine spacing and inflow velocities. Comparisons against published data from the Lillgrund wind farm illustrate that the suggested modelling approach reproduces important trends observed in the field data.


Introduction
Making optimal use of available wind resources for renewable energy generation, implies extracting the highest possible power output, at lowest possible cost, thus minimizing the levelized cost of energy (LCOE). In a wind farm, the wakes behind individual turbines influence both the power production and structural loads of wake affected turbines. Hence, the LCOE is directly influenced by wakes, and therefore, wake modelling plays an important role in wind farm design.
Various fast analytical models for predicting the velocity distribution downwind of turbines have been suggested [1,2,3], and have in common that they prescribe a certain velocity deficit profile and wake expansion rate. The intention of some of these methods is only to give accurate estimates of power production, and as a consequence, the velocity fields they predict do not necessarily lead to accurate structural loads.
For accurate prediction of structural load, aeroelastic simulations need to be performed, where the dynamic wake meandering (DWM) approach has proven to be very successful. Originally suggested [4] nearly two decades ago, the model has been continuously improved and adapted for aeroelastic simulations [5], and is now included in design standards [6].
Validation studies have indicated that the manner in which the influence of individual wakes are combined can have a significant influence on the accuracy of predicted fatigue loads, and different summation approaches are recommended for below [7] and above rated [8] wind speeds. In the current work, a modelling approach is suggested that intends to have computational efficiency closer to that of the fast analytical wake models, yet which incorporates more of the physics of advanced and computationally demanding DWM type models. This approach results in a tool that determines time averaged estimates for velocity and turbulence intensity at the locations of wake affected turbines in a wind farm, for low order power production and structural load estimates at early wind farm design stages. Since it is based on the same building blocks as a DWM type model, it can be used in setting up higher fidelity aeroelastic simulations for eventual refined loads analysis.
This paper is divided into three further parts. In the next section a more detailed description of the numerical model's sub components is given. The section after that shows results of applying the model to some generic wind farm setups, as well comparisons to published measurements at the Lillgrund wind farm. Conclusions are drawn and presented in the last section.
2. Method 2.1. Single wake model 2.1.1. Velocity deficit Consider a single turbine wake, i.e. the area behind the turbine with reduced velocity compared to ambient conditions. The proposed velocity model is based on axis-symmetric Reynolds averaged Navier-Stokes equations for a stationary wake with thin shear layer approximation, as originally proposed by [9]. The symmetry axis, x with velocity component u, is aligned with the downwind direction, with radial axis, r and velocity component v, normal to this. An eddy viscosity assumption −u v = ε ∂ u ∂r relates the mean velocities and their fluctuations. The corresponding momentum and mass conservation equations are: Assuming a Gaussian self-similar profile for the velocity deficit [10,11], the axial velocity, u, in the wake is defined as the difference between reference velocity u 0 and the deficit ∆u, i.e.: The deficit is defined by its mean value at the wake centre-line, ∆u c , and spatial variance ω 2 . This standard deviation ω becomes the wake (half-)width parameter, consistent with other wake models, e.g. [3]. Momentum and mass conservation (ignoring viscous and pressure terms) relate the velocity deficit in the wake to the thrust over the turbine, i.e., T = ρ A u(u 0 − u) dA, resulting (see Appendix A for a detailed derivation) in an expression for the variance: Substituting equations (3) and (4) into (1) and (2), simplifying and then evaluating at the centreline (i.e. r = 0), yields an evolution equation for centre-line velocity, that is solved numerically:  Figure 1. Parameters defining the Gaussian velocity deficit: ∆u(x, r) = ∆u c (x) exp −r 2 2w 2 2.1.2. Eddy viscosity model An axially varying eddy viscosity model is used, which assumes that the total eddy viscosity is made up of two contributions [9,11], a spatially constant contribution related to the ambient wind shear, and a spatially varying contribution related to the shear introduced by the wake. For neutral atmospheric conditions, the ambient contribution is [12]: with κ (≈ 0.4) the Von Karman constant and I a the turbulence intensity at hub height z. While the wake contribution is: with k = 0.015 √ 7.12 from [9] (the √ 7.12 term accounts for the use of different wake width definitions). Since the governing equations are not applicable in the developing flow of the near wake, [9] suggests scaling the wake related contribution by an empirically derived axially varying correction factor: 2.1.3. Initial conditions An initial centre-line velocity deficit is applied at a distance, x 0 /D = 2, downstream of the wake generating turbine. Its magnitude is given in terms of ambient and operating conditions, as suggested by [9] based on wind tunnel studies of rotors: 2.1.4. Statistical wake meandering In this work, the time averaged effects of dynamic wake meandering are taken into account in the steady wind farm model using a statistical approach [13] based on Taylor's dispersion theory [14]. Wake meandering can be described as the displacement of the wake deficit due to large scale turbulence in the ambient wind. The variance σ m of the wake centre offset δm is related to the variance of the ambient wind, with time and length scales that can be related to ambient conditions in the wind farm. According to [13]: where the lateral component of turbulence intensity is found from the mean velocity and turbulence intensity at hub height [6]: and the corresponding integral length scale is given by Furthermore, a Gaussian probability density function for wake offsets is assumed: This allows a convolution to be done of the wake velocity deficit (3) and wake centre displacement (15) functions, giving an expression for the average velocity deficit under meandering.
Comparing (16) to the original velocity deficit given by (3), shows that the effect of meandering is to reduce the magnitude of the velocity deficit, while increasing the effective wake width.

Multiple wake interactions
A major initial simplification in the numerical wind farm model, is assuming that the wind farm can be decomposed into the combined influence of individually calculated wakes. This leads to an eventual complication when the contributions of individual wakes need to be combined into the calculated wind farm velocity influence. Interactions between individual turbines are fully defined by the superposition of individual wake velocity deficits and turbulence intensities. For deficit summation, the superposition of individual wakes is usually derived from principles of momentum or energy conservation [1,15,16,17]. Despite the same underlying principles, differences in defining reference wind speeds and deficits for individual wakes, give rise to many different superposition methods. The choice of method can have a significant impact on predicted velocities in the farm model.  [18] make convincing arguments for and give a detailed description of a momentum-conserving wake summation method. Their main insight is that the total wind farm deficit, ∆U , is a scaled superposition of individual wake deficits, ∆u i , with scaling factors defined as ratios of the mean convection velocities of individual, u i , and combined, U , farm deficits: Conservation of momentum deficit flux over a cross section normal to the flow requires that: where T is a thrust and v the local velocity.
Since v = f (∆v), there is a non-linear relationship between T and ∆v. Mean convection velocities that represent the spatially varying velocities over a cross-section are defined to linearize this relationship [18]: ∆v(x, y, z) dy dz (19) from which the mean convection velocities follow for: • individual wakes (substituting T = T i and v = u i into (19)) u i (x) = u i (x, y, z)) · ∆u i (x, y, z)) dy dz ∆u i (x, y, z)) dy dz (20) • and for the combined farm flow (substituting T = i T i and v = U into (19)) U (x) = U (x, y, z)) · ∆U (x, y, z)) dy dz ∆U (x, y, z)) dy dz (21) The velocity deficit summation rule (17), then follows from the notion that momentum be conserved during wake superposition. The wake summation strategy has certain implications for the order of operations in the wind farm model. Wakes need to be evaluated sequentially as they are encountered along the mean wind direction, which is done by describing turbine positions in a reference system aligned with the mean wind direction. The black dotted lines in Figure 3 indicate the normal plane on which mean convection velocities are evaluated. Determining the global mean convection velocity at any location requires iteration. Once individual wakes have been convected to the required downstream position and their mean convection velocities determined according to equation (20), the farm deficit is estimated using equation (17) with an initial guess for U . This value is then updated according to equation (21), and the process is repeated until the relative change in U is less than a specified tolerance.

Wake added turbulence intensity summation
As is discussed in [12], there are two main mechanisms for increased turbulence in the wake, namely: rotor-generated turbulence (most influential in the near wake, and therefore only indirectly accounted for through equation (10)), and shear-generated turbulence due to the wake deficit profile. Assuming the latter to be the  dominant contribution in the far wake, the same ratio between eddy viscosity and turbulence intensity, equation (7), is re-used to calculate the wake turbulence intensity from eddy viscosity: Crudely considering the wake turbulence intensity as the sum of the initial ambient turbulence intensity and a wake added component, I + , and treating these as uncorrelated normally distributed quantities, the wake added turbulence intensity for a individual wake is found as In the wind farm model, wake added turbulence intensity contributions from different upstream wakes are then superimposed to define new ambient conditions for wake affected turbines [19]:

Results
All sub-models of the simple wind farm model have been described in the previous section, and simulation results for a number of test cases can be presented and discussed.

Influence of wake summation
In addition to the wake summation strategy described in section 2.2.1, a number of common wake summation strategies have also been implemented. To illustrate the significance that wake summation method has for a wind farm model, consider a farm consisting of eight in-line turbines, spaced 4.3 diameters apart (this is the same case considered by [20]).
As noted by [20], moving downwind in the row, the mean inflow to individual turbines decreases before reaching a near constant wind speed, which, when non-dimensionalized with the undisturbed free stream velocity, is found to be independent of free stream velocity. This is due to an equilibrium developing between energy being extracted, and recovered by entrainment, from the free stream flow [20]. Figure 5 shows that this behaviour is not captured when individual wake deficits are calculated with respect to the free stream velocity, and linearly summed (giving nonphysical results for large deficits at low wind speeds). The results in Figure 6 confirm that the expected behaviour is reproduced by the momentum-conserving summation strategy.  The results indicate that an upper bound for equilibrium velocity appears to be given by the summation method that considers only the single largest wake influence, U = U ∞ − max (U ∞ − u i ). As the influence of single wake deficits decreases, either due to increased interturbine spacing, or the effect of wake added turbulence, there is closer agreement between the different wake summation strategies.

Comparison to measurements at Lillgrund
In [20], an axis-symmetric boundary layer solver was used to do simulations of a single row at the Lillgrund wind farm. While restricted to axis-symmetric in-line cases, the solver implicitly handles the wake summation problem by directly solving the relevant governing equations. It was shown that wake effects (especially at high wind speeds) are significantly over estimated if meandering and increased wake turbulence are not accounted for. Significantly increasing ambient turbulence intensity was shown [20] to allow more realistic power estimates to be made, though at the same time creating unrealistic conditions for load estimation. Similar results are obtained with the current wind farm model when disabling meandering and wake added turbulence.  Figure 7. Influence of inter-turbine spacing and wake summation strategy on equilibrium velocity, without considering wake added turbulence. Figure 8. Influence of inter-turbine spacing and wake summation strategy on equilibrium velocity, including the effect of wake added turbulence. Figure 9 shows power (non-dimensionalized using that of the first turbine) generated by the different turbines in row B at a simulated Lillgrund wind farm, for a range of wind speeds, only considering the effect of statistical meandering. Lines with markers are simulated results, and opaque lines wihtout markers are measured results taken from [20]. The same is done in Figure 10, only now including the influence of wake added turbulence. Comparing results between the two figures, it is apparent that qualitative agreement with measurements is only obtained if both meandering and wake added turbulence is accounted for. Figure 9. Power output non-dimensionalized with that of the first turbine, at turbines in row B at Lillgrund, for a range of wind speeds, only considering the effect of wake meandering. Figure  10.
Power output nondimensionalized with that of the first turbine, at turbines in row B at Lillgrund, for a range of wind speeds, considering effects of both wake meandering and wake added turbulence. Figures 11 and 12 show (for the same conditions as Figures 9 and 10, respectively) where on the turbine thrust and electric power curves (in coefficient form) individual turbines are operating. Note that at low wind speeds, where the thrust coefficient is fairly constant, there is a gradient in power coefficient, while towards and beyond rated wind speeds there are steep gradients in both. This helps illustrate how small differences in the calculated velocity deficit can have a large influence on predicted power output.

Conclusion
An engineering type wind farm model was presented with the goal of simulating relevant physics as simple as possible, but as complex as necessary to reproduce realistic inflow conditions at wake affected turbines. The proposed model represents individual wakes as very simple 1D Gaussian velocity deficits, accounts for wake meandering in an averaged sense, and reconstructs the farm velocity field by means of a robust momentum-conserving deficit summation method, and wake added turbulence model. By using only readily available quantities of the free stream flow, the farm model was shown to reproduce expected behaviour in terms of equilibrium velocity in a row of in-line turbines at low and high wind speeds, this is not the case for all wake summation methods. Qualitative agreement with measured power produced in a row of turbines at Lillgrund wind farm was also shown. The agreement was achieved without requiring changes in parameters describing ambient conditions (mean wind speed and turbulence intensity), indicating that qualitatively correct thrust loads (the cause of downstream wake deficits) at downstream turbines were achieved. This gives reason to believe that the proposed wind farm model has the potential to be used for wind farm power production and load prediction, at early design stages.
Work remains to be done to better quantify the contributions of different sub-models to changes and improvements in the predicted wind farm velocity field, and directly relating this to wind turbine fatigue loads (e.g. as calculated from aeroelastic simulations). Validation and verification of the current models will continue by considering more available measurements and advanced simulation results, and sub-models refined and improved as required.
Note that u = u 0 − ∆u (as stated in equation (4)) was used. All terms in equation (A.1) are independent of azimuth, so the integral over θ can be evaluated before making further substitutions. Inserting the expression for the wake deficit (as stated in equation (3)), then gives: Evaluating the remaining integrals gives: Next, thrust is expressed in coefficient form: Finally, one can rearrange terms, substitute ∆u c = u 0 − u c , and simplify to find: