Cfd in Wind Energy: the Virtual, Multiscale Wind Tunnel

Over the past two decades, computational fluid dynamics and particularly the finite volume method have been increasingly used to predict the performance of wind turbines within their environment. Increases in available computational power has led to the application of RANS-based models to more and more complex flow problems and permitted the use of LES-based models where previously not possible. The following article reviews the development of CFD as applied by the wind energy community from small to large scale: from the flow around 2D airfoils to the flow through an entire wind farm.


Introduction
Computational fluid dynamics (CFD) consists of solving the differential equations governing fluid flow using approximate numerical means.In principle, with sufficient grid refinement, the solution of the discretized governing equations yields a flow that is a reasonable representation of reality, at least within the context of the underlying assumptions used to derive the model equations (i.e., incompressibility, instationarity, etc.).
The earliest uses of CFD within the context of wind turbine performance analysis was in the prediction of two-dimensional airfoil properties.However, with increases in computing power it has come to be used at all scales: from the airfoil boundary layer to the atmospheric boundary layer.The nature of turbulent flows is such that their exact solution is simply impossible, especially at high Reynolds number.However, in many instances we are satisfied with simply modeling the effects of turbulence on the mean flow and, although use of large eddy simulation for wind energy applications is slowly increasing, the majority of models are based on the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations derived from the principles of conservation of mass and momentum: U represents the mean velocity vector, p is the modified mean pressure and ρ is the fluid density [1].f represents a body force (i.e., Coriolis, buoyancy, etc.).τ is the specific Reynolds stress tensor.It appears as part of the averaging process and represents the effect of turbulent transport of momentum, which is assumed to dominate viscous terms.
To close the RANS equations, the Boussinesq linear isotropic eddy-viscosity hypothesis is often applied where S is the mean strain rate tensor and ν t must be modeled.The choice of turbulence model depends on the problem at hand and is a balance between desired accuracy and computational resources.Within this basic framework, a wide range of theoretical and practical problems can be investigated.As such, CFD has become a virtual, multiscale wind tunnel.
The following article reviews the application of CFD-based approaches in the field of wind energy for aerodynamic analysis.A separate topic not presented involves the aeroelastic analysis of wind turbines, which also often uses a CFD approach.The interested reader is referred to Hansen et al. [2] for a comprehensive review.The following paper is divided into two broad categories: the first involving simulations at length scales smaller than the rotor diameter, the other at larger ones.

CFD at Small Scales: Aerodynamic Analysis of Wind Turbines
The following three sections present a review of the current efforts in the field of CFD applied to the study of wind turbine aerodynamics.First, the smallest scale application of CFD in wind energy is reviewed: the prediction of two-dimensional airfoil characteristics used in the design of wind turbine blades.Although extensive research in this field has been carried out by the aeronautical community for decades, the unique conditions under which wind turbines operate pose new modeling challenges for the aerodynamicist, not the least of which is stall.In this case, the challenges in applying CFD are not so much numerical but rather relate to the physical modeling of the boundary layer; the focus of current research is on turbulence and transition to turbulence.Conversely, three-dimensional CFD modeling of the rotor requires considerable computational resources and much effort in recent years has been invested in the development of numerical techniques to solve the three-dimensional flowfield around an operating turbine.This will be discussed in depth in the second section.As an alternative to full 3D CFD studies, actuator disk or surface techniques appear promising for the study of rotor aerodynamics, especially for modeling wake effects in a wind park.A review of these simplified approaches constitutes the third, and final, section.
Other reviews have been written on the subject of wind turbine aerodynamics [2][3][4][5].The focus of the present review is on work making systematic use of CFD methods to carry out aerodynamic studies.A common and important part of such methods is model validation and calibration based on experimental data.As regards 2D aerodynamics, data is available from the extensive wind tunnel testing that has been performed on a variety of airfoils and under various operating conditions (see [6] for summary).Experiments have also been performed specifically concentrating on wind turbine aerodynamics: The NREL Unsteady Aerodynamics Experiment program has generated multiple datasets for several 10-m-diameter wind turbines, with Phase VI being the experiment with the most detailed measurements [7].Here, a two-bladed twisted and tapered rotor was tested in the NASA Ames Research Center 24.4 × 36.6 m 2 Wind Tunnel.Several flow characteristics were measured (pressure along blade sections, angles of attack, etc.) in addition to blade loading over an extensive matrix of operating conditions, including yawed conditions.A blind comparison of CFD codes based on this data, involving twenty different research groups using various modeling approaches, showed a disturbing range of predictive capabilities [8], indicating that rotor aerodynamics are still not fully understood.Since that study, a number of authors have used the Phase VI data to validate or improve predictions of wind turbine performance.The very recent MEXICO program [9] is a joint European effort to test a smaller 4.5-m diameter turbine at the DNW wind tunnel.In this program, the blades were again instrumented to monitor pressure and loads.In addition, Particle Image Velocimetry (PIV) measurements were recorded in the near wake, which allow for an even more detailed evaluation of model capabilities.

Prediction of two-dimensional airfoil aerodynamics
The principle of energy extraction by horizontal-axis wind turbines relies on the creation of lift and drag forces by each of the blade sections whose cumulative effect is a net torque.As a first approximation, the flow around the blades may be analyzed as an assembly of annular streamtubes that cross the rotor surface, whose aerodynamics are essentially independent from one another.In this framework, lift and drag forces are defined with respect to the wind velocity relative to the rotating blade section.Since drag acts to reduce torque, most design methods aim to maximize the lift-to-drag ratio at operating angles of attack.However, this is not the only objective as many other parameters can influence the design process: the type of control (stall or pitch), the desired sensitivity to surface imperfections, structural requirements, etc.
All recent efforts in the design of airfoils for wind turbine blades, essentially pursued by NREL, TUDelft and Risø DTU (in historical order), have been realized using computationally inexpensive tools like integral boundary layer (IBL) methods (i.e., XFoil or Eppler codes) and optimization methods.Presently, CFD methods are almost exclusively used as a validation tool to assess airfoil aerodynamic characteristics and provide a deeper understanding of the flow physics.In this spirit, Yang et al. [10] as well as Wolfe and Ochs [11] have used 2D CFD methods to analyze airfoils designed by NREL, with a focus on the S809 airfoil used for the Phase VI rotor.Fuglsang et al. [12] as well as Bertagnolio et al. [13] have used 2D CFD methods for the aerodynamic analysis of the family of profiles designed by Risø DTU.In general, these methods solve the incompressible RANS equations.Turbulence closure is achieved using algebraic models like the Baldwin-Lomax model, one-equation models like the Baldwin-Barth model or two-equation models like the k − , k − ω or k − ω-SST models (see [1] for descriptions of these models).
Despite all recent efforts, there are still major challenges to overcome in the 2D analysis of airfoils, especially when stall operation has to be considered, as underlined by most workers in the field and by Rumsey and Ying [14] in their review.Among these challenges, modeling the transition to turbulence in the boundary layer is among the most important as it significantly affects airfoil performance.Based on typical blade section chord lengths and relative velocities, blade sections contributing to torque typically operate at a Reynolds number, Re, defined as on the order of 10 6 and, as underlined by Mayda and van Dam [15], airfoils are very sensitive to Re in this range with the possible development of laminar separation bubbles (see Figure 1), and stall regimes changing from trailing-edge to leading-edge stall.Since it is not always possible to experimentally reproduce the same range of Reynolds number, CFD can provide insights into these phenomenon that wind tunnels cannot always offer.An increasing number of authors model the onset of transition and its effect on the airfoil aerodynamics, so that their models are able to capture both laminar and turbulent parts of the flow.To predict the location of transition, different strategies are available.The simplest and least computationally expensive approaches applied in CFD consist in using empirical-correlation methods [16] or database methods [17,18].The e N method, although regularly used in IBL methods, is seldom directly used in CFD methods due to inadequate resolution of the boundary layer characteristics [19].The recent works of Windte et al. [20] in the field of Unmmaned Aerial Vehicle (UAV) aerodynamics and of Freudenreich et al. [21] on the TUDelft airfoil series, indicate however that the e N method can be coupled successfully with a CFD method to study transition for this range of Re.A recent approach consists in using a set of two transport equations modeling the intermittency of turbulence and the evolution of the transition Reynolds number based on the boundary layer momentum thickness Re θ,t as devised by Menter et al. [22].Sørensen [23] has applied this model with success to the study of the transitional flow around the S809 and NACA 63-415 profiles at angles of attack below and beyond stall as well as to the three-dimensional flow around the NREL Phase VI test rotor.
Surface roughness can also play an important role in airfoil aerodynamics as it is known to increase viscous drag, decrease maximum lift and reduce abrupt stall.However, these general features are not necessarily always encountered as pointed out by Freudenreich et al. [21], who studied the importance of transition and roughness, both experimentally and numerically, for the TUDelft airfoil series.In CFD simulations, roughness is usually introduced using the equivalent sand grain approach according to which the roughness height is considered small compared to the boundary layer thickness and the roughness effect on the flow is mimicked by increasing the turbulent eddy viscosity in the wall region to obtain higher skin friction [24].

Prediction of three-dimensional rotor aerodynamics
In essence, 3D CFD methods are not fundamentally different from 2D methods with regards to modeling flow physics (especially turbulence) or the discretization of the governing equations.The mesh still needs to be refined in the regions of boundary layers and turbulent momentum exchange must be taken into account.The extra difficulty associated with a 3D analysis is essentially linked to the extension of the mesh to incorporate a third spatial direction and the associated increase in the number of nodes required to adequately discretize the domain.In general, grid independent solutions are much harder to achieve in the 3D case.For example, while total number of nodes on the order of 10 4 are typical for 2D airfoil studies, 3D rotor studies require a number of nodes on the order of 10 6 and the solution of an additional momentum equation.The computational effort of a 3D analysis is therefore much greater and parallel computation methods must be employed to avoid unreasonable simulation times.As such, block-structured meshes are generally preferred since parallelization is easier.When modeling tower and/or ground effects, it becomes necessary to introduce methods such as overset grids (also known as the Chimera method) or sliding planes to decompose the mesh into rotating and non-rotating parts.
A little over a decade ago, the first complete CFD solutions of flow around a wind turbine rotor were produced by Duque et al. [25], Sørensen et al. [26] and Varela [27].Since then, a number of authors have applied CFD methods to the 3D analysis of rotors for various situations, with a focus on the reproduction of the NREL Phase VI rotor experiments.In [28], the results of the CFD analysis are compared with experiments and a vortex method.In [29], the 3D incompressible code Ellipsys3D, which is the cornerstone of Risø CFD methods, as well as another CFD method are used to analyze the NREL rotor under parked conditions.The influence of RANS turbulence closure is discussed and modern approaches for modeling turbulence are employed (LES or DES [30]) with success to model the flow in heavily stalled parts of the blades.Le Pape and Lecanu [31] have applied a compressible CFD formulation (developed at ONERA to study rotating wind turbine aerodynamics) to the case of the NREL rotor.They show the relative advantages of the k −ω-SST turbulence model over the k −ω model.Furthermore, they show the importance of the preconditioning step for solving the system of discretized equations in compressible form, which can significantly alter the flow solution.Tachos et al. [32] as well as Mandas et al. [33] and Hartwanger and Horvat [34] have used the commercial CFD softwares Fluent and Ansys to study wind turbine aerodynamics.Gomez et al. [35] have investigated the wind tunnel wall effects, as well as tower effects, on the flow aerodynamics of the NREL experiments using a compressible CFD formulation.Finally, Bechmann and Sørensen [36] have recently presented results for the MEXICO rotor that exhibit good agreement with experiments.
Although important limitations on the use of full 3D CFD analysis in the field of wind turbine aerodynamics still exist, many authors consider CFD methods to be mature tools for the study of turbine configurations and for some of the particular aspects of blade aerodynamics.Indeed, the study of blade aerodynamics both experimentally and numerically show that the inboard sections of the blade produce lift and drag forces significantly different from the forces that may be calculated using a local blade-element analysis and 2D airfoil characteristics.This phenomenon is referred to as rotational augmentation in the literature [37] and is generally attributed to the effect of local inertial forces and 3D aerodynamics.Three-dimensional CFD methods can naturally reproduce this effect and be applied to generate lift and drag look-up tables for later use in methods that extrapolate local lift and drag characteristics using blade-element analysis, like Blade-Element Momentum (BEM) or actuator disk/surface/line methods, as presented in [28,38,39].Another interesting use of 3D CFD methods is aimed at the comparative and qualitative analysis of rotor designs: Ferrer and Munduate [40] study different blade tip geometries, Mac and Johansen [41] study the effects of adding winglets to the blade tips, Gomez and Barakos [42] study different tip and root configurations and, finally, Chao and van Dam [43] study the aerodynamic effects of modifying the inboard section of the NREL Phase VI rotor for structural reasons.
An important weakness of 3D CFD methods, despite advanced discretization techniques using high-order schemes to handle viscous and momentum fluxes, is the modeling of the wake; its viscous and turbulent diffusion is a particularly difficult problem due to numerical diffusion and to the difficulty in identifying appropriate turbulence models [44].To alleviate this problem, some authors [45,46] have designed hybrid methods where the region of the flow close to the blades is solved by a CFD formulation and outside this region the flow is modeled as being inviscid as in panel, or lifting-line, vortex-based methods.The following section presents actuator surface-based methods that represent an interesting alternative for the study of wake aerodynamics.

Simplified approaches to three-dimensional rotor aerodynamics
First proposed by Rajagopalan [47] and Madsen [48] for the aerodynamic analysis of vertical-axis turbines, the actuator concept has been further developed by several researchers in the field of CFD applied to horizontal-axis wind turbine aerodynamics [49][50][51][52][53].The concept of the Actuator Disk (AD) is central in the heavily used, industrial methods based on BEM with a multiple streamtube integral analysis.It consists in representing the rotor with an equivalent, porous surface whose action on the flow is modeled by an associated system of forces, distributed across volumes or surfaces depending on the exact approach adopted.This singular object can be integrated within regular differential Navier-Stokes based CFD methods with the advantage that unsteady, convective or dissipative effects are naturally taken into account by these methods, contrary to BEM methods which need to explicitly account for these effects in situations involving yawed operation, or in tower shadows, for example.However, in AD-based approaches, the same limitations with respect to the wake exist since the AD does not model tip or root effects of the blades.As such, AD methods are more appropriate for modeling far wake effects and are further discussed in Section 3.3.
Figure 2. Isosurfaces of vorticity for the NREL Phase VI rotor using an actuator surface model [54].
Focusing on the near wake, new methods have been recently proposed that aim to represent the trailing system of vortices, referred to as Actuator Line (AL) or Actuator Surface (AS) methods depending on the object used to model the blades.Sørensen and Shen [55] have developed the AL method in which each of the blades is modeled as a rotating line carrying volume forces distributed across space using convolution functions and whose magnitude is found from local blade element analysis.Ivanell et al. [56] have recently published an analysis of numerically-generated wake structures using this approach.Troldborg et al. [57] have used this technique for large eddy simulation of the wake of the Tjaereborg turbine.Sibuet Watters and Masson [58], as well as Dobrev et al. [59] and Shen et al. [60] have developed the AS method in which each of the blades is modeled as a rotating surface that is a simple, porous 2D plane characterized by velocity and pressure discontinuities whose action on the flow is achieved through an attached system of forces.These discontinuities and forces are determined from blade element analysis and the Kutta-Joukowski relation.Sample results using this approach for the NREL Phase VI rotor are shown in Figure 2. As underlined by Leclerc and Masson [61], in all these approaches, it is important to model not only the system of forces acting on the flow associated with the blades, but also the kinematic conditions the blades impose on the flow due to the attached blade bound-vorticity that, in AD or AS-based approaches, result in velocity discontinuities across these singular surfaces.Otherwise, flow solutions may exhibit spurious, undesirable oscillations in flow properties.

CFD at Large Scales: Simulations in the Atmospheric Boundary Layer
A central issue in the development of a wind energy project is the determination of optimal wind turbine placement within a wind farm.The criteria used to define optimal is always some measure of energy capture; ideally cost of energy is minimized, but often total production is simply maximized.
Given that even modest production increases yield significant financial gains, there is a strong impetus to develop modeling tools which accurately predict the mean flow properties within a park which can be used to improve layouts.
The problem of micro-siting turbines is especially challenging given the implicit nature of the task; turbines themselves affect the local resource by removing kinetic energy and increasing downstream turbulence.A good wind farm flow model might be simply defined as one that provides accurate predictions of these two effects.Unfortunately, as topographic influences become more important, the combined effect of multiple wakes and surface conditions on flow properties becomes increasingly difficult to predict.CFD models are starting to be used to deal with this complexity.At the scale of wind farms, their application is largely focused on the prediction of the flowfield over topography and on the evaluation of turbine wakes.The following sections look at the uses of CFD at larger scales and discuss the simulation of the atmospheric boundary layer (ABL) over ideal and real terrain, and the resolution of wake flow.

Simulation of the homogeneous neutral surface layer
Although axisymmetric simulations are common, many numerical analyses of wind turbine performance assume a fully-developed neutral surface layer as the operating environment.This approach greatly simplifies the true complexity of the ABL while retaining the effect of wind shear over the rotor.In cases where the fully-developed assumption is not valid, for example over complex terrain, the equilibrium neutral surface layer profiles of velocity and turbulence properties are often used to define the approach flow [62].
Faithful reproduction of the neutral ABL using CFD is thus an important task and is often achieved using a RANS approach with two-equation k − ε closure.The standard model as proposed by Jones and Launder [63] is where the turbulence production rate is calculated using and the turbulent viscosity is modeled with C ε1 , C ε2 , C µ , σ k and σ ε are all modeling constants.Several sets have been proposed, a selection of which are discussed briefly here.
Launder et al. [64] made early efforts to re-optimize the original constants for a variety of free shear flows.Although these have come to be known as the standard values they are widely recognized as inappropriate for atmospheric flows.Detering and Etling [65] were among the first to address this issue and showed that, in fact, a non-constant C ε1 was necessary to account for a rotating boundary layer.
Their arguments have since been further refined and extended by Apsley and Castro [66].Focusing on the surface layer, Crespo et al. [67] tailored the standard model constants with the atmospheric measurements of Panofsky and Dutton [68].Richards and Hoxey [69] also calibrated the constants using field measurements.In general, the condition must be satisfied for Equations 10-12 to be an exact solution to the model equations.Since observations are generally used to calibrate C µ and C ε2 , either C ε1 or σ ε are adjusted such that Condition 9 is respected.El Kasmi and Masson [70] have tailored the RNG k − ε model constants in an analogous fashion.Some common sets of model constants are listed in Table 1.Until recently, the inflow conditions defined by the above equations were poorly maintained in an empty domain with uniform surface roughness.Richards and Younis [71] argued that improper boundary conditions were at fault.Soon after, Richards and Hoxey [69] specified the appropriate conditions for use with the k − ε model.Due to differences in wall function implementations between the Richards and Hoxey approach and that of most commercial software, it had been difficult to apply their boundary conditions with success [72].However, Hargreaves and Wright [73] have since illustrated how this can be done.Insufficient grid refinement in the near-wall region may still result in streamwise gradients of turbulence properties; Sumner and Masson [74] have discussed possible solutions to this problem.
Recently, Yang et al. [75] have proposed new inflow conditions for the k − ε model that better reflect the observed turbulence structure in the surface layer.Here, the logarithmic velocity profile is used in conjunction with a turbulent kinetic energy profile that decays with height.Gorlé et al. [76] have corrected an inconsistency in Yang's derivation that ensures conditions for a logarithmic velocity profile are satisfied.In this case, both C µ and σ ε become functions of wall distance, which may present challenges if applying these conditions to flow over complex terrain.

Prediction of the flowfield over real terrain
Evaluation of wind resources at a given site generally consists of a short measurement campaign (∼1 year) where wind speed is measured at a few discrete locations and at several heights.Often, this data is correlated to long-term regional observations to improve climatological representivity [77].The end result is a characterization of the wind resource, in terms of direction, frequency and intensity, at the tower locations.The challenge for the wind energy specialist is then to spatially extrapolate the in situ measurements to the entire region of interest.
Historically, the most common approach has been to use a linear model such as WAsP [78] or MS-Micro [79].Although these models are, to some extent, a reflection of the available computing power at the time of their development, they are clearly inspired by early-generation wind farms.These models are valid for neutral flow over gently sloping terrain and low hills and have performed very well when predicting the flowfield for cases that conform to this limited parameter space [80].However, the modern wind energy industry is expanding and looking to exploit both offshore and mountainous sites.As the terrain becomes more complex, nonlinear effects such as recirculation become dominant flow features and linear models are ill-suited.The calculations of Ayotte [80] over smooth and rough two-dimensional hills suggest that linear models yield unacceptably large error for slopes greater than 0.2.
CFD has never been the tool of choice of wind energy specialists for this task, although use of CFD models is expected to improve the accuracy of resource predictions in areas where flow separation and thermal effects are characteristic of the flow [81].At the research level, CFD has long been used to predict the flow over complex terrain (see Bitsuamlak et al. [82] for a concise review).In fact, CFD models are now starting to be applied to the problem of ideal turbine siting [83], not only for the improved flowfield representation but also the ability to estimate turbulence properties.Examples are also available of CFD being used to help design measurement campaigns in the selection of proper measurements sites (for example, see [84]) and in numerical site calibration [85].Palma et al. [86] have provided guidance on the use of CFD in combination with conventional techniques for wind resource assessment and micro-siting in a recent case study.The development of new commercial CFD software marketed specifically to the wind energy sector will contribute to increased industrial use of such methods.
Keeping that in mind, Landberg et al. [81] have outlined some of the current challenges in wind resource estimation which include: the presence of separated flow in complex terrain, atmospheric stability effects, and the presence of forested regions.These problems may be amenable to a CFD solution and are discussed below.
Flow over analytical shapes Most current use of CFD for flow simulations in complex terrain entails the solution of the incompressible RANS equations with two-equation turbulence closure.Usually, thermal effects and the Coriolis force are neglected.Lower-order turbulence models are avoided as they appear to lack the sophistication to handle recirculation whereas higher-order methods require longer computing times.Presently, the k − ε model, and variants thereof, are the most popular.
Many authors have previously reported on the known weaknesses of the standard model which, in the context of flow over complex terrain, tend to manifest as an overestimation of turbulent kinetic energy and an underestimation of mean flow recirculation.To try to remedy these issues, variations on the k − ε theme are very common.Chen and Kim [87] modified the ε equation by adding a new production term in an effort to balance turbulence production for largely strained flows: In their derivation of the RNG k − ε model, Yakhot and Orszag [88] and Yakhot and Smith [89] modified the standard ε equation in a similar manner: Maurizi has tested these two versions of the k − ε model, along with its standard form, for flow over two-dimensional valleys using wind tunnel data from the RUSVAL experiment [90].For gentle slopes with attached flow, all three models yield similar results for the mean velocity field.However, when recirculation is present, the mean flow solution is much more sensitive to the ε transport equation and results between the models vary considerably; the RNG version yields the best agreement with data.Considering the prediction of turbulent kinetic energy, differences are present even for attached flow, and none of the models provide consistently better predictions.For the Reynolds stresses, the RNG model again appears to provide the best results, however Maurizi suggests a transport equation for uw be included to overcome some fundamental problems with the modeling of this quantity under the eddy-viscosity concept.Maurizi suggests that for flows involving recirculation, the RNG model should be used.
Ying et al. [91] have performed a similar analysis over the two-dimensional analytical hill from the RUSHIL experiment by solving the compressible RANS equations again using three closure schemes: the standard k − ε model, an algebraic Reynolds stress model (ARSM) and an extended k − ε − uw model that includes a transport equation for the uw stress component.All the closure schemes provide reasonable and roughly equivalent results for the mean velocity field but large discrepancies are observed in the calculated turbulent shear stress.Focusing on predictions at the hilltop, both the standard model and the ARSM provide poor underestimations of uw while the k − ε − uw model provides satisfactory results.The improvement is attributed to the ability of second-order closure to account for advection of upstream turbulence.
While some researchers have focused on modifying the ε equation, others have taken a closer look at the prescription of the time scale used in the definition of eddy-viscosity.Whereas standard k − ε uses with the relaxation time being defined as Durbin [92] proposed imposing a realizability constraint, with Nagano et al. [93] and Nagano and Hittaro [94] have developed this idea further and proposed various mixed time scale models based on mean strain rate and vorticity tensors (referred to as the S model and Ω model, respectively, and S − Ω for their hybrid).
Despite some success using the revised nonlinear k−ε model proposed by Shih et al. [95] for flow over a curved hill [96], numerical stability problems prompted Lun et al. to evaluate these improved linear k − ε models for wind energy predictions in complex terrain [97].For flow over a single isolated hill, the Durbin model predicts upstream turbulent kinetic energy well, but severely underestimates its magnitude in the wake.Conversely, the Ω model performs well in the wake, but overestimates k upstream.In terms of mean velocity, the Ω model is in good agreement with measurements whereas the Durbin model grossly overestimates the size of the recirculation zone behind the hill.Use of the mixed time scale model S − Ω somewhat corrects the overestimation of k at the hill top and generally improves estimates of separation and reattachment points in the hill wake.From this analysis and others, Lun et al. conclude that the S − Ω version performs best; Murakami et al. [98] have integrated this approach as part of a wind turbine micro-siting scheme.
In his analysis of computational wind energy assessment methods, Ayotte et al. [80] has also simulated flow over symmetric two-dimensional hills of various slope using a full Reynolds stress model (RSM) and compared with wind tunnel data.Mean flow predictions are in excellent agreement with measurements except in the wake region for large slopes where the mean flow recovers too quickly.This points to limitations in the RANS approach that will not likely be overcome by any of the above treatments and suggests the need for more advanced eddy-resolving techniques, especially if accurate turbulence predictions are desired.It was further confirmed that the predicted mean velocity is relatively insensitive to the closure used: tests indicate that two-equation k − ε and full second-order closures yield negligibly different results.
Flow over real topography Considering flow over real terrain, Kim and Patel [99] have investigated the performance of RNG k − ε by simulating neutral flow through the Sirhowy Valley in Wales, over an embankment on the Rhine in Germany, and over Askervein hill in Scotland.The choice of RNG was motivated by case studies involving flow over a triangular ridge and several two-equation closure schemes.In general, the RNG-based model best predicted mean velocity and turbulence characteristics, including the size and shape of recirculation zones.In a separate work, Kim et al. [100] presented further case studies using the RNG model for Cooper's Ridge, Kettles Hill, Askervein hill, and the Sirhowy Valley.For Cooper's Ridge, the simulation results for mean wind speed at 3 m AGL show good agreement with measurements on the windward slope and at the hill top.Similar conclusions can be made for the flow prediction over Kettles Hill.For Askervein, predicted 10-m velocities are in good agreement, even on the leeside, although hill top wind speeds are underestimated.Some problems predicting hill top and leeside turbulence are noted.The Sirhowy Valley simulations further demonstrated the ability of the RNG model to predict separation and reattachment.El Kasmi and Masson have also applied the RNG model for flow over Blashavel hill [70].
Starting with Raithby et al. [101], many RANS models (for example [62,99,[102][103][104][105]) have been evaluated using the Askervein Hill experiment (see [106] for description, [107] for data).Castro et al. [102] have carried out a grid dependence study using the standard k − ε model in addition to unsteady RANS (URANS) calculations to investigate low frequency time-dependent effects in the lee of the hill.Mean velocities at 10 m AGL are well predicted, but k is overestimated in the upstream region.Recently, Eidsvik [103] presented a down-scaling method for wind power estimation in mountainous terrain for near-neutral flows based, at the smallest scale, on a RANS/k − ε approach which is validated using the Askervein data.To account for the anisotropy of turbulence, Eidsvik employs the nonlinear algebraic stress model proposed by Gatski and Speziale [108].As for Castro et al., mean velocity predictions agree well with measurements at 10 m AGL, even in the lee of the hill.Upstream turbulence is correctly predicted, however hill top values underestimate observations.Prospathopoulos and Voutsinas [62] have used the Askervein case to develop guidelines for RANS simulations in complex terrain.
A new and extensive measurement campaign over the small isolated island of Bolund has been recently carried out by Risø DTU to provide a new database for the validation of flow models over real topography [109].A blind comparison based on these measurements has underlined the challenges involved in making flow predictions over complex terrain (see [110] for details).The average overall error in predicted mean velocity of the top ten models (all RANS-based) was on the order of 13-17% for principal wind directions.The measurements used for the blind comparison should be available early in 2010.The entire database is scheduled for release in 2011.
In addition to variable surface roughness and orography, complex terrain also implies the possible presence of forested regions.Given that forest canopies absorb momentum over a finite depth, a distributed drag force is a more appropriate boundary condition than simply incorporating a displacement height within the velocity wall function [80].Lopes da Costa et al. [111] used the extended k − ε model of Svensson and Häggkvist [112] and an additional drag term in the momentum equation to study the wind over two moderately complex sites with forest cover.Comparisons with wind data above the forest highlight the importance of incorporating the distributed effect of canopies when predicting mean wind speed and turbulence properties.Dalpé and Masson [113] have implemented a similar approach with the modified k − ε closure of Katul et al. [114].Results of one-dimensional simulations within and above three different forests are in good agreement with measurements.Ayotte [80] has implemented second-order closure and compared with LES calculations for flow over a forested hill.Here, the problems observed for an unforested hill are somewhat exacerbated: the increased drag causes the flow to be more prone to separation.The influence of eddies with length scales related to the hill and canopy are not adequately modeled with a RANS approach.
Atmospheric stability is another parameter that may require modeling.Its effect has been largely overlooked by the wind energy community but may be important when considering very large rotors.Alinot and Masson [115] have incorporated thermal stratification to model the ABL under various stability conditions over flat terrain.Eidsvik [103] postulates that stability effects in mountainous terrain will lead to large uncertainties for RANS models.

Large eddy simulation
The vast majority of flow modeling over complex terrain to date employs a RANS approach.Considering wind resource assessment, it appears that the exact closure scheme has little impact on the predicted mean flow velocity for locations of interest for simple cases (i.e. an isolated hill top).The RNG variant seems best at dealing with flow recirculation and is recommended where such effects are important.Concerning the prediction of turbulent properties, there is much less agreement between closure schemes.Given the importance of turbulence predictions for the evaluation of turbine loads, the use of an additional transport equation for uw seems prudent.
Although the aforementioned works have improved RANS modeling of complex flows, fundamental limitations exist.The RANS formulation is inherently incapable of capturing unsteady effects, like intermittant separation and eddy generation and transport, and generally has difficulty modeling turbulence in areas of strong separation.Furthermore, as outlined by Wilcox [24], the Boussinesq approximation ties the Reynolds stresses to local mean flow properties and imposes isotropy-two conditions which are unjustified for many flows-and although second-order closure can be used to account for anisotropy of atmospheric turbulence and allow for stress transport, the dissipation equation is still largely modeled and presents a source of error.
As such, research on the use of large eddy simulation (LES) for flow over complex terrain is increasing (see [1,116] for review of method).LES is analogous to DNS for high Reynolds number flows as a large fraction of the turbulent kinetic energy is directly resolved.Sub-grid scale models are used to handle turbulence at scales smaller than the filter, which is often the grid itself.At this scale, the eddy-viscosity concept has more relevance and the assumption of isotropy may be valid.Although many problems associated with RANS closure can be avoided using an LES approach, the computational effort is considerably greater.In his review on the use of LES for flow over complex terrain, Wood [117] concluded that true LES of atmospheric boundary layer flow over a three-dimensional, rough surface of arbitrary shape was still a long way off based on the grid refinement and averaging time required to properly resolve non-linear interactions at all scales and obtain meaningful turbulence statistics.Citing the work of Chow and Street [118] (see [119] for most recent developments) and Chow et al. [120] regarding LES modeling of flow over Askervein and a valley in the Alps, Ayotte [80] concludes that direct use of LES specifically for wind energy is not yet feasible, although concedes that at some point it will likely be used as part of wind farm design.
Despite these predictions, Uchida and Ohya have developed an LES-based model for analyzing neutral flow over variable orography [121,122] and applied it to the problem of proper site selection [123].Model performance was evaluated using data from wind tunnel tests over simple geometries and from a real site [124].
The difficulties in applying LES to wall-bounded flows are largely due to impractical grid requirements in the near-wall region [125,126]-the region of greatest importance for wind energy purposes.Hybrid RANS/LES methods, in which the near-wall flow is modeled using a RANS approach that is coupled to an LES model away from the surface, may offer a way out.Silva Lopes and Palma [127] were the first to analyze Askervein using an LES approach and a later paper by Silva Lopes et al. [128] elaborates on the strengths and weaknesses of using such a hybrid scheme.More recently, Bechmann and Sørensen [129] have also applied a hybrid model to Askervein that uses RANS/k − ε in the near-wall region and LES with k − ε acting as a sub-grid model for the outer layer.Transition between the two regions is improved with a stochastic backscatter model.Validation with the Askervein data shows, as reported by others, that calculations using RANS/k − ε result in an underestimation of hill top wind speeds and leeside turbulent kinetic energy.However, the proposed hybrid RANS/LES approach yielded excellent agreement with these measurements, although wake velocities were underestimated.

Analysis of wind turbine wakes
An extensive review of the modeling of wind turbine wakes has been previously reported by Crespo et al. [130] and more recently by Vermeer et al. [44].Here, we focus on recent developments involving the use of CFD for the modeling of far wakes in the context of wind farms.As mentioned previously, far wake modeling is dominated by AD methods as the action of the blades need only be accounted for in an average sense.Methods for specifying forces applied by the AD on the flow vary; usually constant loading is assumed or blade-element momentum theory is applied (see [53]), although the choice of method appears to have little effect on resolved far wake properties [131].
Single wake analysis Considering a single isolated rotor in a uniform flow, Sørensen et al. [132] have used the actuator disc concept to analyze wind turbine wake states for laminar conditions; however, most current analyses incorporate turbulence effects.Standard k − ε closure typically underestimates the velocity defect as turbulent diffusion is too high in the wake region.El Kasmi and Masson [131] have applied the Chen and Kim modified ε equation to a discrete volume around the rotor to correct this weakness and improve wake predictions for a single turbine.The Chen and Kim modification effectively limits the turbulent kinetic energy (and viscosity) in this region as the new ε source term is a function of the turbulence production rate.
Cabezón et al. [133] have presented a comparison of single-wake RANS simulations using various two-equation closure schemes, as well as RSM, for the Sexbierum experiment.They have shown that while standard k − ε grossly underestimates the velocity defect, the use of the El Kasmi and Masson approach greatly improves predictions.The realizable model proposed by Shih et al. also performs well.In both cases, the improvement is explained by an increase in the dissipation rate in the region of the rotor.In terms of velocity defect, the results are comparable to RSM and in good agreement with data.All models tend to underestimate turbulence intensity, especially in the near wake, except along the axis of rotation where agreement with measurements is better.With respect to wake turbulence intensity, an earlier study by Gómez-Elvira et al. [134] analyzed the anisotropy of wake turbulence using the Sexbierum case with an explicit algebraic stress model.Prospathopoulos et al. [135] have carried out a similar analysis to Cabezón et al. using the Nibe wake data and k − ω closure for both neutral and stable conditions.Here, the El Kasmi and Masson and Durbin corrections both improve velocity defect predictions.
However, as Réthoré [136] argues, non-physical increases in ε to temper overestimations of ν t makes application of these methods somewhat dubious for multiple wakes.More generally, Réthoré has exposed some fundamental problems with the use of the acuator disk/eddy-viscosity concepts for modeling wind turbine wakes that suggest a completely different approach may be needed.One possibility is LES.Although the eddy-viscosity concept may be used for sub-grid models, the context in which it is applied is more appropriate and should not pose the same problems [136].Jimenez et al. [137] have implemented a simplified LES/actuator disk approach and comparisons of calculated turbulence properties are in good agreement with experimental data.RSM may also be an attractive solution.
Multiple wakes The objective of wake modeling, at the scale of wind farms, is to accurately predict the velocity defect and increase in turbulence to better model power variations and fatigue loading.Early approaches to accounting for the velocity defect in micro-siting relied on empirically-derived guidelines outlining minimum distances between turbines in an array [138].Using an actuator disk approach to analyze a two-row array, Ammara and Masson [139] have shown these guidelines to be overly conservative.Barthelmie et al. [140] have carried out a comparison of wind farm models, ranging from engineering to full CFD models, for predicting power losses due to wake effects in the large Horns Rev wind park.Although models are not specifically identified in the presented results, the RANS/k − ε models tend to overpredict wake losses for narrow measurement sectors; wider sectors yield better agreement with data.Barthelmie et al. [141] have also published a good summary of developments in the field of AD applied to the study of wakes within a wind farm.
The development of LES-based wind farm models is on the rise.Ivanell et al. [142] have presented an analysis of the row-to-row power variation at the Horns Rev wind farm using LES and the actuator disk concept.Meyers and Meneveau [143] have also applied LES techniques in the analysis of an infinite wind park.

Conclusions
Although the cost of a CFD analysis may be comparable to that of a wind tunnel experiment, CFD has the advantages of being infinitely scalable and providing field (not point) data.As with all methods of analysis, the CFD approach has limitations.In the case of wind turbine aerodynamics, these are essentially related to turbulence modeling.However, within this context, CFD simulations can reproduce in situ conditions nearly identically and are currently being used by the wind energy community to carry out aerodynamic analyses at all scales.They constitute a general framework of analysis and are expected to soon become an indispensable tool for the assessment of wind turbine aerodynamics.
With respect to airfoil and rotor performance prediction, the current challenges are in the modeling of turbulence and flow separation.Treatment of transition to turbulence and of roughness effects are important issues which are presently being adressed by the research community.Considering ABL flow, whereas in the past the nonlinearity of flow equations was simplified, increases in computing power are making CFD simulation times reasonable such that these nonlinearities can be resolved.Although current state-of-the-art CFD models already perform much better than simpler methods, the complex and diverse nature of turbulent flow means model improvement will be a mainstay of research.In the short term, practical use of CFD by the wind energy sector will be based on a RANS approach and more comparative studies between various closure schemes are needed to inform their use.In the long term, LES-based methods will likely eventually supplant RANS techniques for the evaluation of flow over complex terrain and in the micro-siting of turbines.

Figure 1 .
Figure 1.CFD calculation of the S809 airfoil showing occurence of laminar bubble phenomenon.

Table 1 .
A selection of proposed k − ε model constants.