## Abstract

Large eddy simulations (LES) of wind farms have the capability to provide valuable and detailed information about the dynamics of wind turbine wakes. For this reason, their use within the wind energy research community is on the rise, spurring the development of new models and methods. This review surveys the most common schemes available to model the rotor, atmospheric conditions and terrain effects within current state-of-the-art LES codes, of which an overview is provided. A summary of the experimental research data available for validation of LES codes within the context of single and multiple wake situations is also supplied. Some typical results for wind turbine and wind farm flows are presented to illustrate best practices for carrying out high-fidelity LES of wind farms under various atmospheric and terrain conditions.

This article is part of the themed issue ‘Wind energy in complex terrains’.

## 1. Introduction

With the current trends towards rotors of ever increasing size and wind farms with ever larger footprints, there is a clear need for reliable simulation methods that can faithfully replicate the dynamics of wind turbine wakes. Of most practical concern is understanding how wake effects alter turbine loading and power extraction and how such effects are influenced by atmospheric conditions and topography. This is essential to fully characterize the ultimate limit states which turbines must withstand and to correctly predict the accumulated fatigue loading. However, as noted by Churchfield *et al.*, ‘the current state of knowledge concerning wind turbine wakes and how they interact with other turbines and the atmospheric boundary layer (ABL) is not complete’ [1]. With these thoughts in mind, this work is an attempt to outline the current state of the art: the most prevalent modelling techniques, how they are implemented in the most commonly used LES codes and the current level of satisfaction in their results.

Although the engineering models available for modelling wind turbine wakes in farms [2–5] are fast and useful for certain purposes, the information they provide about the flow is quite limited. And while computations solving the Reynolds-averaged Navier–Stokes (RANS) equations are widely used for wind farm simulations [6–8], the results are known to depend on the choice of turbulence closure models that limits their general use. Furthermore, the nature of the time-averaged equations restricts the quantity of information they can provide about the inherently unsteady processes underlying wake phenomena. Direct numerical simulations (DNS), where all scales of the flow are resolved, are still far too computationally demanding for use in wind farm simulations.

By contrast, large-eddy simulations (LES), which resolve the largest eddies and model only the smallest scales, are an appealing alternative for carrying out reliable simulations of large wind farms. They can provide detailed and accurate information about the dynamics of the wake and the associated turbulence length scales along with power production and turbine loading. Recent advances in computational power have made LES accessible for that very purpose and this technique has been applied to many wind farm projects to evaluate, among other things, the effects of various atmospheric conditions, terrain, and turbine placement [1,8–17]. In direct contrast to lower-order approaches, ‘by resolving parts of the spectrum of turbulence, LES captures important features of unsteady, anisotropic turbulence in turbine wakes interacting with the ABL’ [18]. Such information is proving to be invaluable.

This article provides an overview on the use of LES within the research community to produce high-fidelity wind farm simulations. The focus is on the methods available for rotor modelling, for reproducing the relevant characteristics of the atmosphere, and for incorporating topographic effects. As such, it is expected to complement the recent works by Sanderse *et al.* [19] and Mehta *et al.* [20], who have provided a thorough review of computational fluid dynamics (CFD) methods for wind turbine aerodynamics and of LES for wind farm aerodynamics, respectively. An overview of available experimental data to test LES methods against is also provided herein. Typical results obtained from LES employing various submodels and under a variety of configurations are presented and discussed along with the current status of these methods and possible research directions.

## 2. Numerical methods

### (a) State-of-the-art large-eddy simulations codes

An overview of the most common LES codes used by the research community to model wind farm aerodynamics is presented in tables 1 and 2. The intent here is to summarize the range of models implemented in each one. The list of references given next to each code does not necessarily refer to the original developers but to a sample of studies applying a subset of the available models to a specific case. Owing to space limitations, it is not possible to present a detailed mathematical description of each of the codes; the reader is referred to the reference list for more information. A short introduction to LES is nonetheless given below along with some discussion of the codes.

The goal of LES is to explicitly resolve the large-scale energy containing motions and to model the effects of small-scale turbulent motions via so-called subgrid-scale (SGS) closures. The smallest eddy that can be represented on a computational grid is a function of the grid resolution and, in LES, the Navier–Stokes equations are modified to contain the dynamics of these resolved scales alone. This is achieved by applying a low-pass filter. As part of this operation, the velocity *u*_{i} is decomposed into the sum of a filtered velocity (the tilde operator denotes this filtering) containing the large scales of motion and an unresolved velocity *u*_{i}′ representing the small scales and calculated using an SGS model:
2.1
The flow field is then approximated by solving the filtered three-dimensional incompressible Navier–Stokes equations:
2.2
and
2.3
where *i*,*j*,*k*∈{1,2,3} represent the three orthogonal coordinate directions, is the filtered pressure, *ρ*_{0} is the (constant) density and *t* is time, *μ* is the molecular viscosity, *τ*_{ij} is the subgrid stress tensor
2.4
where represents various body forces (to model e.g. the presence of a rotor, to impose synthetic turbulence or a shear profile, etc.), is the Coriolis force. The last term of (2.2) accounts for thermal stratification, i.e. for cases where the ABL is not neutral. This term, referred to as the Boussinesq approximation, accounts for the buoyant force acting on a volume of fluid as a result of small changes in density owing to variations in temperature. The gravitational field strength is denoted by *g*, is the filtered potential temperature (where the angle brackets represent a spatial average [20], and *θ*_{0} is a reference potential temperature.

The subgrid stress tensor is usually modelled based on the assumption that the SGS stresses can be accounted for through an effective eddy viscosity. This is an important assumption. Unlike RANS models where the eddy viscosity model is applied to turbulence at all scales and is of questionable validity, LES limits its use to only the smallest scales where one might reasonably hope that turbulence has a more universal character. Mehta *et al.* [20] discuss the physical validity of this assumption which, for better or worse, ties the direction of the unresolved stress tensor to the resolved strain rate.

Filters are similarly applied to the transport equations for scalars. In the case of potential temperature, which accounts for the thermally generated turbulence,
2.5
where *α* is the thermal diffusivity, *q*_{j} is the SGS flux of heat and
2.6
Obviously, equation (2.6) is considered only in cases where thermal stratification is of importance.

The SGS model is responsible for the coupling between resolved and subgrid scales. One of the simplest SGS closures bases this coupling on an equilibrium between the flow of kinetic energy from the resolved scales to its dissipation by the subgrid scales. The simplicity of the standard *Smagorinsky* model makes it attractive from an implementation point of view and it is one of the models available in the SOWFA code [1,29]. More sophisticated SGS models have been developed such as the scale-dependent Lagrangian dynamic model [37] which is employed in both the EPFL [10,38] and JHU [27,28] codes. This development primarily arose from a difficulty in correctly specifying the closure coefficients in eddy-viscosity/diffusivity models (see [20] for a complete description of these coefficients). The standard Smagorinsky model assumes constant values for these coefficients which have been well established for homogeneous and isotropic turbulence [26]. This assumption, however, loses its validity near walls under large wind shear conditions [38] or strong thermal stratification [26]. Dynamic models were introduced to cope with this issue, where the coefficients are allowed to vary based on the resolved flow using a test filter. But the dynamic model still assumes scale invariance which means that values do not depend on the filter size. LES performed under neutral conditions [39] as well as experiments [40] proved this assumption to be inaccurate near the wall which led to the development of scale-dependent dynamic models. For a detailed discussion of the evolution of SGS models, the reader is referred to [20,38].

EllipSys3D, SOWFA, PALM, SnS, WiTTS and UMN employ finite difference or volume methods to discretize the Navier–Stokes equations. EPFL, JHU and KU Leuven are pseudo-spectral codes which make use of pseudo-spectral discretization in the horizontal directions and finite differencing in the vertical direction.

### (b) Rotor modelling

Actuator methods have been developed over the last 20 years as a means of representing the effect of the rotor on the flow in CFD computations without requiring the boundary layer on the wind turbine blades to be resolved. This section provides an overview of such methods as used in LES. Power control and aeroelasticity in the simulations are also briefly discussed.

#### (i) Actuator methods

*Actuator disc:* the actuator disc (AD) has been employed for many years as a means to include body forces in the energy and momentum equations. It is an intrinsic part of the one-dimensional momentum theory and an important ingredient in the blade-element/momentum approach. In CFD computations, it is often impossible to resolve the entire rotor blade geometry and instead it may be replaced by appropriately determined body forces. In its simplest version, the AD consists of a constant loading determined from the average thrust acting on the rotor. Comparisons with experiments have demonstrated that the method works well for axisymmetric flow conditions and can provide useful information regarding basic assumptions underlying the momentum approach [41,42]; turbulent wake states occurring for heavily loaded rotors [43,44]; and rotors subject to coning [45,46]. However, important flow dynamics are missing if no information regarding the actual loading is included in the forcing. A way to include the actual forcing is to employ a blade-element approach combined with tabulated aerofoil data. The first computations of actual wind turbines employing numerical AD models in combination with a blade-element approach were performed by Sørensen and colleagues [41,42] and Masson and colleagues [47,48] in order to study–unsteady flow phenomena. Masson *et al*. [47,48] also devised techniques for employing their AD model to study wake interactions in wind farms and the influence of thermal stratification in the ABL. The first comprehensive and systematic LES investigations of a wind farm in the ABL using an AD approach for the wind turbines are due to Jimenez *et al.* [49], Calaf *et al.* [13] and Porté-Agel *et al.* [50]. More recent contributions include the papers by e.g. Gopalan *et al.* [51] and Gundling *et al.* [52].

Generally, the thrust force is computed as
2.7
where *ρ* is the density of air, *A* is the swept area of the rotor, *u*_{0} is the *reference* velocity and *C*_{T} is the thrust coefficient. A special issue when using body forces for wind turbines operating in wind farms is the determination of this reference wind speed; the ‘undisturbed’ wind velocity for turbines operating in wake flow is not well defined. In [48], *u*_{0} was taken as the unperturbed upstream velocity at hub height and *C*_{T} was determined from the thrust curve of the wind turbine employed for the experiment. In [49], the reference velocity was taken as the unperturbed velocity of the incident flow in the centre of the disc and *C*_{T} was obtained directly from the thrust curve of the considered wind turbine (a 300 kW Nordtank). In [13], the reference velocity was determined indirectly from the average velocity over the rotor disc. Introducing the axial induction factor *a*, the average axial velocity in the rotor plane is given as *u*_{R}=(1−*a*)*u*_{0}. From one-dimensional axial momentum theory *C*_{T}=4*a*(1−*a*); hence
2.8
where, according to the theory, *a*=1/3 for an optimum rotor. In [13], they employ as a typical value *a*=1/4. Furthermore, they employ two different versions of (2.7): one in which the body force is taken as a constant over the rotor disc and one in which the body force varies according to the local axial velocity. The latter is evidently more correct, as local turbulence properties and wind shear will change the value of the body force both in time and space. In [50], the body force is obtained directly from a blade-element approach whereby the problem of determining the reference velocity is avoided. This however demands detailed knowledge of the rotor geometry and its control system. The latter was obtained directly from SCADA data by [50]. However, in the planning phase of a wind farm, there is usually little technical information about the wind turbines to be used. In this case, the best approach may be to exploit the fact that most modern wind turbines are tip-speed regulated to rated wind speed after which they are power regulated. For wind speeds less than rated, (2.8) may be employed with *a* taken as optimal. At wind speeds above rated, the average power coefficient is given by
2.9
where *P*_{G} is the generator power. Combining (2.9) with one-dimensional momentum theory and stating that *C*_{P}=(1−*a*)*C*_{T} results in
2.10
which is valid for inflow velocities between rated and cut-out.

While relatively simple formulae can be used to obtain the axial force component, the azimuthal force component (which is responsible for rotation in the wake) can normally only be determined via a more detailed blade-element approach. In [50], it was shown, however, that the swirl component may be of importance for describing the dynamics of the near wake. The AD representation taking into consideration the swirl component is usually referred to as the actuator disc with rotation (AD + R). In the case where there is no *a priori* knowledge of the rotor, one might assume that the circulation is a constant. Assuming furthermore that the induction is governed by a root vortex of radius *δ*, typically corresponding to the hub radius, the azimuthal force component is given by (see [53] for more details)
2.11
where *λ* is the tip speed ratio, *x*=*r*/*R*, and 〈*u*_{R}〉 is the average axial velocity in the rotor plane.

A tip correction is usually applied in AD models to correct for a finite number of blades. The tip correction was originally introduced by Prandtl & Betz [54]. In this method, the correction factor is introduced to correct the circulation between an infinitely bladed rotor () and an *N*_{b}-bladed rotor (*Γ*). An approximate formula for the Prandtl tip loss function was introduced by Glauert [55]
2.12
where *ϕ*=*ϕ*(*r*) is the angle between the local relative velocity and the rotor plane. In the axisymmetric model by Sørensen & Kock [42], the tip correction was implemented in the computations of the body forces by dividing the lift coefficient by the tip correction factor *F*. As an axisymmetric model solves the flow equations by distributing the force field uniformly in an annulus, the output of the axisymmetric Navier–Stokes solver corresponds to the averaged solution or (*a*_{B}*F*,*a*′_{B}*F*). In order to apply aerofoil data within the blade-element approach, the local induction factors at the blade position are needed. From the definition of Prandtl's tip loss factor, and . Once *a*_{B} and *a*′_{B} are computed, the flow angle and local relative velocity at the blade position are obtained from the blade-element approach. As the tip loss factor depends on the local flow angle, an iterative procedure is required.

*Actuator line* (AL): the main limitation of the AD method is that it is unable to capture near-wake characteristics like tip and root vortices. To overcome this limitation, the AL model was developed by Sørensen & Shen [56]. This model combines a three-dimensional Navier–Stokes solver with a technique in which body forces are distributed radially along each of the rotor blades. The kinematics of the wake are determined by a full three-dimensional simulation whereas the influence of the rotating blades on the flow field is included via tabulated aerofoil data to represent the loading of each blade. As in some implementations of the AD model, aerofoil data and subsequent loading are determined iteratively by computing local angles of attack from the movement of the blades and the local flow field. This concept allows for detailed study of wake dynamics, i.e. of the tip vortices and their influence on the induced velocities in the rotor plane. The AL technique has in particular been employed in many LES of wind turbines in wind farms [1,29,32,57–61].

In the AL method, it is usually necessary to distribute the forces by using a filtering function based on a Gaussian kernel, which is given as
2.13
where *r* denotes the distance to the actuator point and *ϵ* denotes the width of the kernel. Typically, *ϵ* takes a value of two to three mesh point lengths [58]. In [60], a set of guidelines is proposed for determining the most important parameters in the AL method. Among these is a recommendation to let *ϵ* be determined by the actual shape of the rotor blade. It should be noted that this expression was first introduced in conjunction with the AL method [56], but a similar expression may well be employed for an actuator disc.

*Actuator surface* (AS): an extension of the AL technique is the AS method in which vorticity sources are distributed along the surface of the rotor blade [62,63]. Using this technique, the limitations of using a lifting line approach are avoided which leads to a better representation of the loading near the blade tips.

All of the above methods that rely on aerofoil data should have these data corrected for three-dimensional rotational effects [64]. Aerofoil data are usually taken from wind tunnel measurements under two-dimensional conditions.

#### (ii) Power control and aeroelasticity

*Power control:* when rotation of the rotor is taken into consideration in an LES code, a power controller can be implemented such that the turbine will adapt the angular speed of the rotor according to the prevailing wind conditions. This is particularly useful in wind farm simulations as it does not require a reference wind speed to determine the rotor speed of a given turbine, which can be problematic [65].

The basic strategy for determining the speed of the rotor is as follows. Suppose at a given time step, the rotational speed of a given turbine is *Ω*(*t*) and the wind velocity distribution at the rotor is known. The LES code will, with these data, calculate the aerodynamic torque *T*_{aero} on the AD, AL or AS. The torque produced by the generator (*T*_{gen}) is a function of rotor speed and can be determined from the generator torque curve, whose origin is explained in detail in [66]. Any difference between the generator torque and the aerodynamic torque causes the system to accelerate, viz.
2.14
where *I*_{rotor} is the rotor inertia and *I*_{gen} the drivetrain inertia (comprising the generator and any other rotating components). The ODE is solved discretely: the rotational speed after some small increment of time is approximated with *Ω*(*t*+Δ*t*)≈*Ω*(*t*)+d*Ω*/d*t*|_{t}Δ*t*. All variables are re-evaluated at the next time step, and the process is repeated. Above rated power, the blades are actively pitched to limit the aerodynamic torque. The amplitude of pitch movements is usually determined using a feedback control loop through the use of a PI controller.

*Aeroelasticity:* larger wind turbines with increasingly slender blades give rise to important aeroelastic effects. For such rotors, the non-rigid character of the blades and their interaction with the wind play an important role in the determination of loads. Aeroelasticity can be modelled by performing a two-way coupling between the LES code and a structural analysis model that can prescribe the blade deformations resulting from the aerodynamic loading. The deformed blade and its motion in turn affect the LES calculation. A detailed description of such a coupling can be found in [36], where simulations were performed using the structural analysis module of FAST [67] with the research LES code SnS that modelled the rotor with an AD+R.

### (c) Modelling of the appropriate characteristics of the atmosphere

The two most influential boundary layer properties on turbine loading and power production are the ambient turbulence (intensity and structure) and mean vertical wind shear. As such, the focus of many engineering models is on correctly accounting for these factors. However, the interplay between atmospheric turbulence and turbine-induced turbulence is complicated, making it difficult to derive general rules. LES offers the possibility to account for such complexity and even reproduce other secondary phenomena such as wake meandering. Although results in some situations are dependent on the choice of SGS closure, Sarlak *et al.* [68] investigated the effect of the SGS model on flow structures and turbine loading for two aligned wind turbine wakes and conclude that the choice of SGS model has limited influence. Despite this there seems to be a consensus developing within the wind energy community that a scale-dependent Lagrangian dynamic SGS model [69,70] (or very similar models such as [71]) is the preferred approach to simulate the flow around and within large-scale wind farms. This closure sidesteps any debate over model constants as the Smagorinsky coefficient is both variable and determined as part of the solution.

The vertical and lateral boundary conditions used to simulate the ABL are generally a stress- and flux-free condition on the upper boundary [1,18,32,38,72], periodic conditions in the spanwise direction and a law of the wall based on Monin–Obukhov similarity theory [1,18,32,38,72,73] at the ground. Unfortunately, similarity theory is based on the mean statistics of an equilibrium flow, and it is generally recognized for inhomogeneous flows that wall treatments are the least satisfactory aspect of LES where wall effects are modelled [74]. The absence of a shear stress at the upper boundary requires a prescribed pressure gradient to drive the flow.

Inflow and outflow treatments depend on how the horizontal equations are discretized. High-fidelity spectral techniques are often used for horizontal components for better numerical accuracy but these methods require periodic boundaries in all associated directions [27], in principle, limiting their application to infinite wind farms. Recently, Stevens *et al.* have proposed a method for simulating a finite length wind farm with doubly periodic boundary conditions [27,75]. Abkar & Porté-Agel [72] propose a similar approach. If differencing techniques are used instead of spectral methods, the inflow can be simply prescribed and the outflow treated with a non-reflecting or convective condition [1,32].

Ideally, the inflow should represent in every way a real ABL—this implies the need for a *precursor* (or possibly concurrent) simulation to generate a time series. Often, for the sake of efficiency, the inflow is *synthesized* based on the approach by Mann [76,77]. The turbulence can be created with proper spectral characteristics via body forces injected in the flow [78]. This is not necessarily done at the boundary itself, but rather at some upstream distance from the wind turbine rotor to avoid convecting the turbulence structures over large distances which would require (somewhat useless) grid refinement [79]. Wind shear can also be (independently) generated via body forces distributed throughout the domain usually with a slip condition at ground [57,80–82] (there seems to be some debate as to the correct treatment as Nilsson *et al.* [83] opt for no-slip). The Mann technique was recently generalized to take into account the change of turbulent structures for stable and unstable atmospheric conditions [84,85]. This generalization was performed by adding two new parameters to the three parameters associated with the original Mann spectral tensor [76], i.e. the gradient Richardson number and the rate of destruction of temperature variance.

Despite the attractiveness of this approach, there is concern that synthetically generated atmospheres lack ‘proper ABL coherent structures’ [27]. More specifically, Park *et al.* argue that turbine-scale variables such turbulence intensity and wind shear are strongly related and should not be prescribed independently [73]. The alternative is to let the turbulence properties develop in an empty domain through solution of the model equations. The advantages are many, but additional computational time is needed as well as a longer domain. Bechmann [24] notes that the flow field from a precursor simulation already represents a solution of the model equations and is thus adapted to the solver. Churchfield *et al.* [1] identify another plus: the ability to capture wake meandering owing to large-scale atmospheric motions. This is something stochastic models may struggle to reproduce. Unfortunately, mean profiles often demonstrate an ‘overshoot’ near the ground [32] (although Brasseur & Wie [86] offer remedies). The precursor simulation is driven by an imposed gradient on the filtered pressure. For the subsequent simulation where the precursor data is used as a Dirichlet condition at the inflow, the pressure gradient term is sometimes dropped [27] although Churchfield *et al*. [1] retain the time-average value from the precursor and [10,32] keep the imposed forcing term.

It could be argued that the best approach is a combination of both: use synthetic turbulence to reduce the computational time of a typical precursor simulation. Wind shear could be imposed via body forces and the domain filled with turbulence generated with the Mann algorithm as a good first approximation. The precursor simulation could then be run as usual.

### (d) Complex and forested terrain

#### (i) Complex terrain

Wind farms are located on sites of various complexity and it is important to be able to accurately model topographic effects on both wind flow and wake recovery. Modelling the flow over complex terrain, however, remains quite challenging and even more so in the presence of wind turbines. This has been performed with modest success using RANS [65] and has only recently been attempted using LES. To the best of the authors' knowledge, the only work published on this topic is by Yang *et al.* [12] wherein turbine rotors have been modelled using an AL approach.

Although wind farm studies in complex terrain using LES are rare, there has been extensive work on the central microscale problem, i.e. simulating the flow over variable topography without the presence of turbines. The progress has been considerable since Wood, over 15 years ago, concluded that the ‘LES of flow over any surface other than a very idealized one is still a long way off’ [87]. A well-known test case is the Askervein hill in Scotland [88] and more recently the small island of Bolund in Denmark [89]. The latter was the basis for a blind comparison of microscale flow models [90] and was found to be a challenging case. In general, the RANS-based models outperformed the LES ones. Recently, flow over the complex site at Perdigão has been modelled with LES but validations with field measurements [91] have not been published.

There are two different methods to model the flow over complex terrain with LES. In the first one, the computational grid is set to follow the terrain topography, using so-called body-fitted grids. This has been performed with success by many research groups [24,26,30,92–95]; however, generating a body-fitted grid over a heterogeneous surface can be challenging. Furthermore, some LES codes are poorly adapted for their use, especially those based on spectral methods which are limited to relatively simple geometries [96]. In these cases, the immersed boundary method can be used [97,98]. This technique makes it possible to incorporate topography on a regular structured grid. The presence of a flow boundary representing the topography is simulated through the use of body forces added to the Navier–Stokes equations. Cristallo & Verzicco [98] discuss the efficiency of this approach compared to using body-fitted grids. It has been used by many authors and has been shown to produce reliable results [12,33,99].

In both approaches, the flow near the surface is treated with special wall-layer modelling. The wall-layer modelling can be performed using a wall function in the LES simulations, which relates the wall stress to the tangential velocity at the first near-wall node. Alternatively, the wall-layer can be simulated via the RANS equations in the inner layer, whose solution can be used as a boundary condition for the LES simulations performed in the outer flow. Both techniques have been described extensively by Bechmann [24].

#### (ii) Forested terrain

Forests are known to have an important impact on the flow, causing high wind shear and strong turbulence intensity above the forest canopy [100]. Many wind farms are built in forested areas and a thorough understanding of the effects of forests on the flow is necessary for satisfactory prediction of turbine power and loads. To the authors' knowledge, no work has been published where an LES was used to predict the flow in a wind farm located in forested terrain. The authors are familiar with many ongoing activities in this area, for example within the New European Wind Atlas (NEWA) project; however, no results have been published at this time. The literature on modelling the flow above canopies via LES is however fairly extensive.

The method to handle forests in LES is usually to include a canopy model wherein drag forces are added as source terms in the filtered Navier–Stokes equations to reproduce the aerodynamic drag effects from the canopy [101]. Various, but similar, expressions have been proposed for the drag term, generally involving a drag coefficient for the canopy, the local leaf area density and the filtered wind velocity. [101–105] all present examples on the use of a canopy model in LES. The model of Dupont *et al.* [105] also adds a sink term in the equation for the SGS turbulent kinetic energy to account for the increased dissipation of turbulent eddies in the inertial subrange. According to Liu *et al.* [101], LES, combined with a canopy model, is capable of reproducing many of the turbulent flow features observed over vegetation on flat terrain.

There have been some developments in modelling complex forested terrain using a combination of the methods outlined above [101,105,106] although not in the presence of wind turbines.

## 3. Available experimental research data for validation

In an article from 2009, Barthelmie *et al.* [107] cited the lack of available data from large wind farms as a major hurdle in the effort to characterize wake behaviour. Still, a certain number of experiments have been carried out, in the field or in wind or water tunnels, where the emphasis has been either on wind farm behaviour as a whole or the behaviour of the wake behind a single wind turbine. This section identifies a number of data sources which have already been used for model validation (only research data which have been accessed previously are discussed). The research data traffic light indicates three (simplified) states for accessing research data [108].

— Green is ‘go’ (but always double check).

— Yellow is ‘prepare’ (available soon).

— Red is ‘stay’ (‘no go’ with no change expected in the situation in the near term).

‘Green’ research data can be shared extensively, ‘yellow’ research data may be shared with a non-disclosure agreement or a fee or similar, while ‘red’ research data are highly confidential and access is unrealistic in the short term. ‘Red’ data will not be included here.

As already discussed, a fundamental feature of LES is its inherent capability to predict unsteady effects. Results from measurement campaigns are usually provided as 10 min statistics in the form of a standard deviation around a mean value. When comparing LES results with measurements for validation purposes, it is important to take into account the inherent variability of LES. The importance of this has been shown in a recent work by Andersen *et al.* where the variability of different quantities has been investigated with LES along a line of 16 wind turbines modelled as ALs or ADs [16].

### (a) Experiments focusing on wakes behind single wind turbines

The wake properties behind operating, scaled wind turbines have been measured a few times in large wind tunnels under very controlled conditions and in the field under various flow and terrain conditions.

#### (i) Wind tunnel measurements

The wind tunnel measurements are recorded under controlled inflow conditions, where the focus has been to measure the near wake behind a rotating, scaled wind turbine model. Seven major wind tunnel experiments with measurement of the near wake behind single wind turbines have been identified and are shown in table 3.

#### (ii) Field measurements

Field measurements of wakes are performed on full-size wind turbines located in different types of terrain, ranging among offshore, coastal, flat and even complex terrain. The focus has been both to measure the long-term wake distribution with mast mounted instruments and to scan the near wake distribution with LiDARS or Windscanners. Nine different field experiments with free access have been identified in table 4 and many of these have been used in recent international research projects. A 10th one, not publicly available, was also included in table 4. Please note that dataset no. 8 will be enlarged with more measurements in 2017. The recent double wake experiment no. 11 between two smaller turbines has been published in [125]. As part of the SWiFT project no. 12 the wake measurements will be shared on the NREL open data repository [126].

### (b) Experiments and datasets focusing on wind farm wakes

Data from wind farms of different size, layout and operational mode have been used to validate operational park models and CFD calculations during the last 15 years. At least five different wind farm layouts have been tested in wind tunnels and SCADA data from small and large wind farms have been used in several large research projects.

#### (i) Tunnel measurements

Although a large wind tunnel is required, model wind farms have been tested several times over the last three decades. Table 5 includes five recent experiments but other older experiments exist.

#### (ii) Field measurements

The field measurements mainly include access to a supervisory control and data acquisition (SCADA) database containing 10 min statistics for each wind turbine in the wind farm. Unfortunately, offshore wind farm measurements often do not include inflow conditions owing to a lack of masts. Table 6 lists 12 different wind farms which have been used in recent model validations. Only the two oldest wind farms are publicly available, while the remaining 10 wind farms require a signed NDA before use. The Horns Rev dataset has been used for validation in many research projects, while some of the new wind farms (e.g. no. 11 [137]) have a more challenging location and layout.

## 4. Results and discussion

An overview of simulation results that illustrate the effect of the various considerations outlined above is presented and discussed in this section. Unfortunately, owing to space limitations, it is not possible to present results from all the models outlined above.

### (a) Rotor modelling and wake interactions in a neutral atmospheric boundary layer

Figure 1 shows comparisons of the axial velocity profiles at different distances behind the single model rotor considered in the NTNU blind comparison tests [139,140] obtained from EllipSys3D simulations where the rotor is modelled either as an AD+R or by ALs. Comparisons with experimental data are shown where available. The results show that the AD+R and AL models produce almost identical profiles up to 18*R* and can correctly predict the near- and far-wake behind a single model wind turbine. There are some minor differences with respect to measurements close to the wake centre as the presence of the hub and nacelle is not considered in the simulations.

There are however differences in the wake properties predicted by AD+R and AL models. Martínez-Tossas *et al.* [142] performed simulations of three different wind turbines using both AD+R and AL representations of the rotor. While both these approaches provide similar power output, they produced notably different near-field instantaneous wakes. The wake produced by the AL contains tip and root vortices that are convected downstream. By design, the AD cannot generate tip or root vortices. This difference causes subtle changes in the way the wake breaks down into turbulence. However, while the instantaneous profiles are different in the near wake, the mean profiles are nearly indistinguishable. Porté-Agel *et al.* [50] have also compared LES using AD, AD + R and AL behind a single wind turbine and in a wind farm. These authors found that AD+R and AL model results compare favourably with measurements in terms of turbulence intensity and mean wind velocity profiles. On the other hand, the AD model was found to overestimate the average velocity in the centre of the wake and underestimate the turbulence intensity at the outer radius where turbulence levels should be highest owing to the presence of a strong shear layer. They conclude that resolving the near wake requires a rotor model that induces wake rotation and allows for a non-uniform distribution of the turbine-induced forces. All three methods produce similar results in the far wake. In wind farm simulations, the AD + R model is less computationally expensive compared to AL and yet captures the important turbulent scales of the ABL and of the wake. AL simulations are usually expensive owing to the fine grid resolution (a spatial resolution of at least *R*/30 is recommended) and small time step (such that the movement of the rotor tip is confined to a single grid cell at each time step [143,144]). This is the least attractive aspect of the AL simulations but the price to be paid to study near wake flow and near wake dynamics in greater detail [145]. AD+R is a cost-effective alternative to AL if the goal is to study the far wake and wake interactions of the wind turbines [9] or farm-to-farm interactions [17]. AS methods, for their part, have not been used to model large wind farms with LES.

Nilsson *et al.* [146] investigated the influence of aerofoil data and rotor geometry (the details of which are generally not known for a commercial wind farm) on turbulence intensity and mean streamwise velocity in the wake as well as turbine power and thrust. In order to do so, the performance of a line of turbines modelled with ADs using three different rotor configurations was simulated in a pre-generated turbulent and sheared inflow using EllipSys3D in a neutral atmosphere. The distance between the turbines and level of ambient turbulence intensity was varied and the rotors operated at a constant rotational velocity. Figure 2 illustrates the power relative to the first turbine (located at a position of 17*R*) and turbulence intensity level along the line of turbines separated by 6.6*R* resulting from the use of different aerofoil data and rotor geometries. Notably, very small variation is observed between the three configurations. Of the quantities compared, differences were only observed in the computed turbulence intensity level when using a larger turbine spacing (of 14*R*, not shown here). It was concluded that, provided the inflow is sufficiently turbulent, the choice of aerofoil data and rotor geometry is not crucial if the goal of the LES is to compute the mean wake characteristics within the wind farm.

Figure 3 shows the results of an LES of the Lillgrund wind farm performed using EllipSys3D and an AD model, with a pre-generated turbulent and sheared inflow under neutral atmospheric conditions. Lillgrund is composed of 48 2.3 MW wind turbines (Siemens SWT-2.3-93) located unusually close to each other with a spacing as small as 6.6*R* in the most dense direction. As the blade geometry and aerofoil data were not available, they were approximated by downscaling [9] the 5 MW NREL model wind turbine [66]. This is a well-documented method in the community because, as mentioned above, turbine documentation is often a limiting factor. Results are shown for two lines of turbines with and without a power controller below rated wind speed [147]. The agreement with the experimental data is very satisfactory, and it is worth noting that the simulation captures the flow recovery in a gap between turbines (row E). Again, lack of access to the actual turbine geometry and aerofoil data does not seem to be problematic. Use of a power controller yields only very small differences in the power predictions, owing to the shape of the *C*_{P}–*λ* curve which is fairly flat around its maximum. The calculated thrust however (not shown here, see [148]), whose associated thrust coefficient varies importantly as a function of the tip speed ratio, was found to be very dependent on the use of a power controller. The same conclusion applies to the downstream velocity field. A power controller is therefore an important element in LES of wind farms. Storey *et al.* [36] arrive at a similar conclusion after coupling their LES code with FAST. Their work showed among other things that notable changes in loading occur between the different control regions. The contribution of aeroelastic effects in the simulations was however not discussed in this work. To the authors' knowledge, study of the interaction between aeroelasticity and wake development in LES is so far limited. However, there seems to be a general consensus [149] that considering aeroelastic effects is essential when simulating multi-megawatt turbines; the assumption of small blade deformations upon which most current aeroelastic codes are based loses accuracy with increasing blade size (see [150] for a study on this issue). For a review of the state of the art in aeroelasticity research on wind turbines, the reader is referred to [149].

Churchfield *et al.* [1] have also simulated the Lillgrund wind farm under neutral atmospheric conditions using LES, in this case with the open-source platform OpenFOAM and a precursor simulation to generate the ABL. A power controller was again used and the geometry of the modelled turbine was determined using all publicly available information to match as closely as possible the 2.3 MW turbines used in the park. Satisfactory agreement was obtained in terms of power up to the fifth turbine in a wind-aligned row, while over-predictions varying from 25% to 40% were obtained further downstream (results not shown here). It is worth noting that a substantial one million computational hours (CPU hours) were required to simulate 10 min of real-time operation. However, experience acquired with OpenFOAM since that simulation was first run would lead to the same simulation requiring one-half to one-third of this cost. By comparison, the simulations shown in figure 3 required about ten thousand CPU hours.

Another large wind farm that has been an important subject of interest for LES studies is Horns Rev where the turbine spacing is more typical. Results comparing the power production along a line of turbines as predicted by two LES models under neutral atmospheric conditions (see [27] and [82]) were shown in [20]. Better agreement was found by Stevens *et al.* which was claimed to be owing to a more physical description of the ABL based on a precursor simulation, as opposed to the synthetic turbulence used in the work of Ivanell [82]. However, as mentioned in [20], other parameters such as the grid resolution, SGS model and numerical schemes differed between the two models making it difficult to attribute the difference to ABL modelling alone.

It is also worth mentioning that LES has begun to be used for the optimal control of complete wind farms, see e.g. Goit *et al.* [151], who have dynamically controlled in time the energy extraction of the turbines modelled as actuator discs in a large wind farm so as to influence in an optimal way the flow field in the neutral boundary layer. This strategy, whose details can be found in [151], led to an overall increase in the energy extraction of the farm of up to 7%.

### (b) Wake interactions in a thermally stratified atmospheric boundary layer

Atmospheric instability has been shown to have a significant impact on wind turbine wake development [14,23,72,114,124,132,152]. A higher level of instability in the ABL has indeed been observed to contribute to faster wake recoveries and higher turbulence levels [114]. As such, the effect of atmospheric stability cannot be ignored when attempting to predict the wake characteristics, associated power production and wind turbine loads within a wind farm. Dörenkämper *et al.* [14] recently investigated the impact of a stable ABL on the wakes of offshore wind farms. This was done by simulating a 21 turbine offshore farm with PALM using an AD+R representation of the rotors. They found that wake effects under stable conditions were as much as twice as large as in the unstable ABL. The effect was found to be stronger with turbines separated by at least 20*R*. Keck *et al.* [81] used SOWFA to model both a neutral and an unstable ABL with the aim of improving their dynamic wake meandering model for use under different atmospheric stability conditions. They studied the impact of the different stability conditions at constant turbulence intensity and concluded that the stability had an important effect on the turbulent length scales and that the turbulence intensity level alone is not enough to describe the wake dynamics in the ABL. Lu *et al.* [153] used an LES–AL model to simulate a large idealized wind farm operating in a stable ABL. High wind shears associated with stable conditions were found to create significant asymmetric loadings on the turbines. Furthermore, the Coriolis effect was found to contribute a significant horizontal shear and to drive part of the turbulence energy away from the wake centreline. Churchfield [1] studied the effect of a change in atmospheric stability from neutral to unstable conditions on wake evolution, blade loading and structural response. Two turbines separated by 14*R* were considered, whose blades were modelled as flexible ALs coupled to the dynamic model FAST [67] which includes a power controller. They found among other things that the turbine located downstream produced from 15% to 20 % more power in unstable conditions when compared with neutral conditions. Abkar *et al.* [72] very recently performed simulations using an LES–AD+R model of a 2 MW turbine subject to a convective, neutral and stable boundary layer. The simulations reproduced experimental observations [114,124], i.e. that the higher turbulence levels associated with greater atmospheric instability led to a quicker recovery of the wake. It was also determined that the magnitude and spatial distribution of the turbulence production, dissipation and transport terms were affected by the atmospheric stability. Figure 4 shows isocontours of turbulence intensity for the three thermal stratifications considered, where the significant impact of atmospheric stability on wake development is revealed. Meandering was also a subject of interest in this paper. It was predicted to be stronger with increasing instability, with the lateral component being larger. The findings from the LES simulations were used to generalize the analytical wake model developed by Bastankhah and Porté-Agel [154]. It is worth noting that LES combined with the actuator disc has also recently been used to simulate the ABL throughout a full diurnal cycle [155,156]

The studies cited above all model stability effects *directly*. The other possibility is to model the turbulence synthetically as recently proposed by Machefaux *et al.* [23]. They performed a similar study to the one of Abkar *et al.* but they imposed a wind shear using body forces on which they superposed turbulence generated with a generalization of the Mann algorithm extended to take into account atmospheric stability [84]. External forcing terms were also added to the momentum equations to represent thermal effects. The wind turbine rotors were represented with an AD+R model. Comparisons with LIDAR measurements downwind from the Nordtank turbine (see §3) show that this model, although seemingly faster than the one from Abkar *et al.* as it avoids a precursor simulation, is in fair agreement in terms of axial velocity deficit from 6*R* (differences closer to the rotor were suggested to be due to the normalization procedure used for the LIDAR). It is also in fair agreement in terms of vertical wake meandering magnitude, while lateral wake meandering was found to be underestimated under unstable conditions. Machefaux *et al.* comment on some of the limitations of their approach. First, distortions in turbulence properties occur when inserting synthetic turbulence into the domain. Second, the flow near the ground is not well modelled because of the use of a slip condition. Lastly, a scale-dependent dynamic SGS model is better at resolving the statistics of the flow near the surface than e.g. the mixed model from Ta Phuoc *et al*. [157] (which is echoed in the work by Wan *et al.* in [93]). Still, Machefaux *et al.* claim that their results illustrate that adapting the length and velocity scale of the atmospheric turbulence is a satisfactory way to model the effects of thermal stability, arguing further that the large-scale dynamics of the wake can be decoupled from the evolution of the wake deficit. As hinted at in §2c, this is in opposition with an earlier publication from Park *et al.* [73] that affirmed that many turbine-scale variables are strongly interrelated and should not be prescribed as independent degrees of freedom. As no direct comparison of these two methods has been carried out for wind farm flow, the debate as to the best approach is as of yet unresolved.

Yang *et al.* [12] have actually implemented both techniques to generate the inflow, i.e. imposing synthetic turbulence or using a precursor LES, but to simulate different configurations, making a direct comparison between these approaches impossible. The first approach was used to simulate a utility-scale wind farm of 43 turbines operating under stable atmospheric conditions. The synthetic turbulence generated with the Mann algorithm [158] was imposed on a time-averaged inflow generated using the WRF mesoscale model [159]. It seems however that the version of the Mann algorithm used in this case was not the one which was generalized for non-neutral cases [84] and employed by Machefaux *et al.* The synthetic turbulence was shown to contribute to wake recovery as expected, but the statistics of the turbulence were not studied in greater detail. It was however mentioned by the authors that while using synthetic turbulence was found to be efficient, the associated turbulence structures ‘are not necessarily physical’. Further investigations as regards the effects of synthetic turbulence on wake behaviour and interactions between the turbines were suggested by the authors. The second approach was used by this group to model wind turbines sited in complex terrain.

### (c) Complex and forested terrain

As mentioned in §2d(i), to the authors' knowledge, LES results of an array of turbines in complex terrain have been published only once. In this work by Yang *et al.* [12], an array of eight turbines was considered on a site of moderate complexity under a neutral atmosphere (Prairie Island in Minnesota where turbines are not currently installed). Two met masts installed at the site were used to verify that the wind profile without wind turbines was modelled satisfactorily. The immersed boundary method was applied wherein a triangular mesh was immersed into the background grid to discretize the complex terrain. The blades were modelled as ALs rotating at a constant tip-speed ratio based on the rotor-swept average streamwise velocity one diameter upstream, i.e. there was no power controller. A precursor simulation was performed on a flat terrain to compute a fully developed turbulent inflow to use for the subsequent simulation over the real topography. Figure 5*a* shows the topography of the site as well as the location of the hypothetical turbines and the met masts, whereas figure 5*b* shows contours of the time-averaged turbulent kinetic energy at the horizontal plane at hub height. These results show that the turbulence induced by the terrain is comparable to the turbulence induced by the turbines for this specific case which is significant for power extraction and the dynamic loading of the turbines.

Yang *et al.* continued this work in a more recent paper [33] by using similar methods to systematically investigate the effect of a three-dimensional hill of varying height on a turbine located at varying distances downstream. The flow behind the hill (in the absence of a turbine) was compared with wind tunnel measurements and shown to be in good agreement in terms of time-averaged velocity, turbulence intensity and primary Reynolds shear stress. It was found that the turbulence induced by the hill of a certain minimum height allowed the wake behind the turbine to recover faster. A model was also suggested from the outcomes of this work to represent the turbine-added turbulence kinetic energy in complex terrain.

As discussed above, no studies have been presented wherein LES has been used to predict the flow in a wind farm located in forested terrain. There is one recent work by Nebenführ *et al.* [100], who used LES together with a canopy model to study the effect of the added turbulence caused by a forest canopy on the dynamic loading of a single wind turbine under neutral atmospheric conditions. The flowfield calculated using LES over the forest canopy was used as input for load calculations performed in FAST. The results (not shown here) suggest that turbulent fluctuations for the turbine located above a forest canopy are two times larger than without the forest, which hints at the necessity for further research in this area.

## 5. Conclusion

When correctly designed and properly validated, LES can act as a ‘virtual wind tunnel’ and provide unique, reliable and valuable information about wind turbine wake properties everywhere of interest. They can provide a nearly complete characterization of the turbulence in the field, something that would be difficult and costly if attempted experimentally. LES can predict wake length [145,160], meandering [23,72,81], the impact of large wind farms on local meteorology [10] and wind turbine wake physics as a whole. Information obtained from such studies can be used to develop, test and improve less expensive engineering models used among other things to predict wind turbine wakes as well as to parametrize wind farm effects in weather and climate models [13,25,72,81,107,154,161–164]. LES is likely the key to developing that ‘new generation of models […] required to deal with this complex interaction of wakes with each other and the boundary-layer’ that Barthemie *et al.* [107] envisioned almost a decade ago.

In this review, various methods to model important quantities related to LES of wind farms have been surveyed and discussed. These include the SGS stress, the presence of the rotor, the ABL and its stability, as well as terrain and forest effects. Results obtained using different codes under diverse configurations were shown and discussed as a means of providing a set of good practices that can be used to perform high-fidelity LES of wind farms.

In this regard, the use of the AD+R approach to model the turbine rotor in wind farms is found to be a cost-effective solution to study the far wake and wake interactions. Access to the actual aerofoil data and rotor geometry required as input to this model is not essential if the aim of the simulation is to compute the main characteristics of the wind farm. The use of a power controller to actively determine the rotational speed of the rotors is however important, as well as the consideration of aeroelastic effects for large wind turbines. Two different methods to generate the ABL under different atmospheric stability conditions are available, i.e. via a precursor LES or by imposing the shear profile and turbulence structures synthetically. The answer as to which approach is best is not definitive. LES of wind farms sited in complex and forested terrain remain scarce but techniques tested in the absence of wind turbines look promising.

As features such as forested regions and complex terrain are being considered to a greater extent, the complexity of future LES models is only expected to increase. Recent works by Lu & Porté-Agel [153] and Park *et al*. [73] also stress the need to include a more accurate description of atmospheric physics to account for large-scale atmospheric forcing as a means to better capture ‘real-life’ turbulence. Resolving a very large domain with LES under these circumstances to capture, for example, farm-to-farm interactions would imply an extremely large computational cost. Indeed, given the large spectrum of scales to take into account, compromises are needed. One possibility is to ‘nest’ LES within mesoscale weather models such as WRF [159]. Work along these lines has already started [51,165,166] and mesoscale-to-microscale coupling is currently becoming an important area of research. This is seen for example in the U.S. Department of Energy's Atmosphere to Electrons (A2e) programme [167] wherein one of the main research areas is dedicated to the coupling of mesoscale and microscale models.

## Authors' contributions

S.-P.B. coordinated the writing process, contributed to all sections of the paper, and wrote §§1, 2d, 4 and 5. J.S. also contributed to all sections of the paper, and wrote §§1, 2c and 5. J.N.S. wrote §2b and participated in discussions regarding the structure and content of this work. K.S.H. provided the overview of the available data sources addressed in §3 and performed SCADA data analysis of the datasets in table 6 used in previous model benchmarks [6–8,91,107]. S.S. was involved with §2a and with the results and discussions surrounding figure 1. S.I. provided expertise and guidance in the preparation of this article and participated in discussions regarding its structure.

## Competing interests

We authors declare we have no competing interests.

## Funding

The Swedish Energy Agency is acknowledged for providing research funds for this work. The research of J.S. is supported by FRQNT and Dawson College.

## Acknowledgements

K. Nilsson is acknowledged for providing figure 3

## Footnotes

One contribution of 11 to a theme issue ‘Wind energy in complex terrains’.

- Accepted November 3, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.