Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Review of computational fluid dynamics for wind turbine wake aerodynamics

Review of computational fluid dynamics for wind turbine wake aerodynamics 1 INTRODUCTION During the last decades, wind turbines have been installed in large wind farms. The grouping of turbines in farms introduces two major issues: reduced power production, because of wake velocity deficits, and increased dynamic loads on the blades, because of higher turbulence levels. Depending on the layout and the wind conditions of a wind farm, the power loss of a downstream turbine can easily reach 40% in full‐wake conditions. When averaged over different wind directions, losses of approximately 8% are observed for onshore farms and 12% for offshore farms (see e.g. Barthelmie et al. ). When studying power losses and blade loading, wind turbine wakes are typically divided into near and far wakes. The near wake is the region from the turbine to approximately one or two rotor diameters downstream, where the turbine geometry directly affects the flow, leading to the presence of distinct tip vortices. Tip and root vortices lead to sharp gradients in the velocity and to peaks in the turbulence intensity. For very high tip speed ratios, the tip vortices form a continuous vorticity sheet: a shear layer. The turbine extracts momentum and energy from the flow, causing a pressure jump and consequently an axial pressure gradient, an expansion of the wake and a decrease of the axial velocity. In the far wake, the actual rotor shape is only felt indirectly, by means of the reduced axial velocity and the increased turbulence intensity. Turbulence is the dominating physical process in the far wake, and three sources can be identified: atmospheric turbulence (from surface roughness and thermal effects), mechanical turbulence (from the blades and the tower) and wake turbulence (from tip vortex breakdown). Turbulence acts as an efficient mixer, leading to the recovery of the velocity deficit and a decrease in the overall turbulence intensity. Far downstream, the velocity deficit becomes approximately Gaussian, axisymmetric and self‐similar. Wake meandering, the large‐scale movement of the entire wake, might further reduce the velocity deficit, although it can considerably increase fatigue and extreme loads on a downwind turbine. It is believed to be driven by the large‐scale turbulent structures in the atmosphere. The distinction between near and far wakes is also apparent when classifying the existing numerical models for wind turbine wake aerodynamics (see Table ). The first and simplest approach is an analytical method that exploits the self‐similar nature of the far wake to obtain expressions for the velocity deficit and the turbulence intensity. The second, Blade Element Momentum (BEM) theory, uses a global momentum balance together with two‐dimensional (2D) blade elements to calculate aerodynamic blade characteristics. The vortex‐lattice and vortex‐particle methods assume inviscid, incompressible flow and describe it with vorticity concentrated in sheets or particles. Panel methods similarly describe an inviscid flow field, but the blade geometry is taken into account more accurately, and viscous effects can be included with a boundary‐layer code; the wake follows as in vortex‐wake methods. These four methods have been extensively discussed in previous reviews, such as Vermeer et al. , Crespo et al. , Snel and Hansen et al. The last two methods, the generalized actuator disk method and the direct method, are relatively new and are commonly called computational fluid dynamics (CFD) methods. In this survey, we will review these methods, discuss their ability to predict wind turbine wakes and give an outlook into possible future developments. Classification of models. Method Blade model Wake model Kinematic Thrust coefficient Self‐similar solutions BEM Actuator disk + blade element Quasi one‐dimensional momentum theory Vortex lattice, vortex particle Lifting line/surface + blade element Free/fixed vorticity sheet, particles Panels Surface mesh Free/fixed vorticity sheet Generalized actuator Actuator disk/line/surface Volume mesh, Euler/RANS/LES Direct Volume mesh Volume mesh, Euler/RANS/LES RANS, Reynolds‐averaged Navier‐Stokes; LES, large eddy simulation. The paper is organized as follows. First, we will consider the Navier–Stokes equations and discuss their use to predict turbulent flows ( Section 2 ). Section 3 then discusses rotor modeling, Section 4 discusses wake modeling and Section 5 deals with the question how to verify and validate wind turbine wake CFD codes. 2 GOVERNING EQUATIONS 2.1 The incompressible Navier–Stokes equations It is reasonable to assume that the flow field in wind turbine wakes is incompressible, since the velocities upstream and downstream of a turbine placed in the atmosphere are typically in the range of 5–25 m s −1 . Only when calculating the aerodynamics at blade tips that compressibility effects may be important. Since in most calculations of wind turbine wakes the rotor is not modeled directly (which will be discussed later), the incompressible Navier–Stokes equations are a suitable model to describe the aerodynamics of wind turbine wakes: ∇ ⋅ u = 0 , ∂ u ∂ t + ( u ⋅ ∇ ) u = − 1 ρ ∇ p + ν ∇ 2 u , supplemented with initial and boundary conditions, which will be discussed in Section 2.3 . In the case of a non‐neutral atmosphere, the Boussinesq approximation is typically employed to account for buoyancy effects, and an extra equation for the temperature has to be solved. The effect of the rotation of the Earth, given by the Coriolis term, is typically neglected in many wake studies but can have an effect when computations involve large wind turbines and wind farms (see e.g. Porté‐Agel et al. ). Although this set of equations provides a complete model for the description of turbulent flows, it is not easily solved. The difficulty associated with turbulent flows is the presence of the non‐linear convective term, which creates a wide range of time and length scales. For example, in the atmospheric boundary layer (ABL), the largest turbulent scales are of the order of 1 km, whereas the smallest scales are of the order of 1 mm. Inside the blade boundary layers, the scales are even smaller. The range of scales depends on the Reynolds number ( Re), the dimensionless parameter that indicates the ratio of convective forces to viscous forces in the flow. Large values of the Reynolds number, encountered in blade and wake calculations, lead to a large range of scales, making computer simulations extremely expensive. Resolving all scales in the flow, so‐called direct numerical simulation (DNS), is therefore not feasible. Turbulence models need to be constructed, modeling the effect of the unresolved small scales based on the behavior of the large scales. However, even with the cost reduction provided by a turbulence model, one cannot resolve both the boundary layers on the turbine blades and the turbulent structures in the wake. This necessitates a simplified representation of the wind turbine in case of wake calculations and a simplified representation of the wake in case of blade calculations. 2.2 Turbulence modeling A large number of turbulence models have been constructed in the last decennia (see e.g. Wilcox, Sagaut and Geurts ). This section will discuss the two most important methodologies in turbulence modeling for wind turbine wakes, RANS and LES, their applicability and their limitations. 2.2.1 RANS Reynolds‐averaged Navier–Stokes methods aim for a statistical description of the flow. Flow quantities such as velocity and pressure are split in an average and a fluctuation, the so‐called Reynolds decomposition: u ( x , t ) = u ¯ ( x ) + u ′ ( x , t ) The averaging procedure, ensemble averaging, is such that u ¯ ¯ ( x ) = u ¯ ( x ) and u ′ ¯ ( x , t ) = 0 . The Reynolds decomposition [equation ] is substituted into the Navier–Stokes equations, which are then averaged, resulting in ∂ u ¯ ∂ t + u ¯ ⋅ ∇ u ¯ = − 1 ρ ∇ p ¯ + ν ∇ 2 u ¯ − ∇ ⋅ u ′ u ′ ¯ The term u ′ u ′ ¯ is called the Reynolds stress tensor, which appears as a consequence of the non‐linearity of the convective term, and represents the averaged momentum transfer because of turbulent fluctuations. The Reynolds stresses can be interpreted as turbulent diffusive forces. In wind turbine wakes, they are much larger than the molecular diffusive forces v ∇ 2 ū , except near solid boundaries. In order to close the system of equations, a model is needed to express the Reynolds stresses in terms of mean flow quantities. A widely adopted approach of modeling the Reynolds stresses exploits the Boussinesq hypothesis (not to be confused with the Boussinesq approximation mentioned earlier). Based on an analogy with laminar flow, it states that the Reynolds stress tensor can be related to the mean velocity gradients via a turbulent ‘eddy’ viscosity ν T : u ′ u ′ ¯ = v T ∇ u ¯ + ∇ u ¯ T so that the RANS equation becomes ∂ u ¯ ∂ t + u ¯ ⋅ ∇ u ¯ = − 1 ρ 1 ∇ p ¯ + ∇ v + v T ∇ u ¯ + ∇ u ¯ T This approach of modeling the effect of turbulence as an added viscosity is widely used for turbulent flow simulations. It is very useful as an engineering method, because the computational time is only weakly dependent on the Reynolds number. However, the validity of the Boussinesq hypothesis is limited. In contrast to ν,ν T is not a property of the fluid but rather a property of the type of flow in question. Since eddies are fundamentally different from molecules, there is no sound physical basis for equation , and DNS calculations have indeed not shown a clear correlation between u ′ u ′ ¯ and ∇ ū . The Boussinesq hypothesis is therefore inadequate in many situations, for example, for flows with sudden changes in mean strain rate (e.g. the shear layer of the wake), anisotropic flows (e.g. the atmosphere) and three‐dimensional (3D) flows. Many different methods have been suggested to calculate ν T , typically called zero‐equation (algebraic closure, such as mixing length), one‐equation and two‐equation models (see e.g. Wilcox and the references therein). The k ‐ϵ modelis an example of a two‐equation model often encountered in wind energy wake applications; the k ‐ω model [with shear stress transport (SST) limiter] is more convenient near blade surfaces. In the k ‐ϵ model, two additional partial differential equations are introduced, one for the turbulent kinetic energy k and one for the turbulent diffusion ϵ.They contain a number of constants that have been determined by applying the model to some very general flow situations (e.g. isotropic turbulence decay or flow over a flat plate). A fundamentally different approach is the Reynolds stress model (RSM), also called differential second‐moment closure model, which does not rely on the Boussinesq hypothesis. In the RSM, all the components of the Reynolds stress tensor are modeled, which makes it suitable for anisotropic flows. However, it leads to six additional partial differential equations (PDEs), making the approach expensive. Moreover, these PDEs contain terms that have to be modeled again, and often, closure relations resembling the Boussinesq hypothesis are still employed. Lastly, the disappearance of the (stabilizing) eddy viscosity term can lead to numerical problems. 2.2.2 LES In recent years, LES is receiving more attention in the wind energy wake community, because of its ability to handle unsteady, anisotropic turbulent flows dominated by large‐scale structures and turbulent mixing. This is a significant advantage over RANS methods, but the drawback is that the computational requirements for LES are much higher than for RANS. In LES, the large eddies of the flow are calculated, whereas the eddies smaller than the grid are modeled with a subgrid‐scale model. This is based on the assumption that the smallest eddies in the flow have a more or less universal character that does not depend on the flow geometry. Mathematically, this scale separation is carried out by spatially filtering the velocity field, splitting it in a resolved (also called large scale, simulated or filtered) velocity and an unresolved (small scale) part. In general, this filtering operation is defined as a convolution integral: ũ ( x , t ) = ∫ u ( ξ , t ) G ( x − ξ , Δ ) dξ where G ( x − ξ , Δ ) is the convolution kernel, depending on the filter width Δ . The subgrid velocity is then defined as the difference between the flow velocity and the filtered velocity: ũ ( x , t ) = ( x − t ) − ũ ( x , t ) This decomposition resembles that of Reynolds averaging but with the difference that in general u ˜ ˜ ≠ u ˜ and u ′ ˜ ≠ 0 . Applying the filtering operation to the Navier‐Stokes equations (and assuming certain properties of the filter) leads to the following equation: ∂ u ˜ ∂ t + u ˜ ⋅ ∇ u ˜ = − 1 ρ ∇ p ˜ + v ∇ 2 u ˜ − ∇ ⋅ u u ˜ − u ˜ u ˜ As in the RANS equation a new term appears, the subgrid‐scale (SGS) stresses. These stresses represent the effect of the small (unresolved) scales on the large scales. A widely used model to calculate these stresses is the Smagorinsky model, which again employs the Boussinesq hypothesis: τ SGS = u u ˜ − u ˜ u ˜ = − v SGS ∇ u ˜ + ∇ u ˜ T Possible ways to calculate the subgrid‐scale eddy viscosity ν SGS are to use an analogy of the mixing‐length formulation and to use one‐equation or two‐equation models involving kinetic energy and turbulent dissipation. Because of the use of the Boussinesq hypothesis, similar limitations as in RANS are encountered. A large number of other subgrid‐scale models have therefore been proposed, for example, dynamic models, regularization models and variational multiscale models (see e.g. Sagaut and Geurts ). In contrast to RANS, where the computational cost is only weakly dependent on Re, the computational cost of LES scales roughly with Re 2 . Near solid boundaries, where boundary layers are present, LES is extremely expensive because it requires refinement in three directions, whereas RANS only requires refinement in the direction normal to the wall. A possibility is to employ a hybrid approach: RANS to resolve the attached boundary layers and LES outside the wall region, so‐called detached eddy simulation (DES). Since equations and both have a similar form, this switch between RANS and LES can be made by changing the eddy viscosity based on a wall‐distance function. 2.3 Boundary conditions One of the ongoing challenges in CFD simulations of wind turbine wakes, especially when comparing with experimental data, is the prescription of inflow conditions that mimic all relevant characteristics of the atmosphere, such as the sheared velocity profile, the anisotropy of the turbulence and the instationary nature of the inflow. In early CFD simulations, uniform, laminar inflow profiles were used, but LES simulations showed later that both the presence of the shear inflow profile (the ABL) as well as the turbulence in the incoming flow have a pronounced effect on the flow field behind the rotor. For RANS simulations, Monin‐Obukhov similarity theory (see e.g. Panofsky and Dutton ) can be used to prescribe profiles for the velocity components and turbulence quantities, which are independent of time. For LES, unsteady inflow data are necessary, and basically, two different types of methods exist: synthesized inlet methods and precursor simulation methods. The synthesis technique was employed in the models of Mann, Veers and Winkelaar. They are based on the construction of spectral tensors, which model the frequency content of the wind field, and have the advantage that certain parameters like turbulent length or time scales can be directly specified. Mann's model combines the rapid distortion theory (based on a linearization of the Navier‐Stokes equations) and the von Kármán spectral tensor to generate a 3D incompressible atmospheric turbulence field, which is homogeneous, stationary, Gaussian and anisotropic but does not include the effect of the ground boundary. Troldborg used this model to create a turbulent inflow profile. This profile is introduced close to the rotor (one rotor radius upstream) to prevent the decay of turbulent fluctuations before they reach the rotor. In order to ensure a divergence‐free velocity field, the turbulent fluctuations are introduced via body forces in the momentum equation, following Mikkelsen et al. In the second category, a separate precursor simulation is made as an inflow model for a successor simulation. The advantages of this method are that the ground boundary is taken into account and that the resulting turbulent velocity field is a solution to the non‐linear Navier‐Stokes equations. However, it cannot be easily manipulated to obtain certain desired turbulence characteristics, and it is computationally more expensive. In the work of Bechmann, Meyers et al. and Stovall et al. , such a method was used. In Stovall et al. , roughness blocks were introduced to generate the turbulence. When a desired turbulence level is achieved, the blocks are removed, and the turbines are introduced in the domain. The model of Bechmann can handle the anisotropy of the atmosphere, except near the surface, where it switches to a RANS solver. As was stated earlier by Vermeer, Bechmann observed again that ‘imposing realistic and even experimentally obtained inflow conditions was found extremely difficult’. At other boundaries of the computational domain, the prescription of boundary conditions is not straightforward either. The most physical model for the presence of the ground is a body‐fitted mesh with a no‐slip boundary condition, but such a condition does not allow the prescription of a ground roughness and would require very fine grids near the surface. Therefore, a common approach in both RANS and LES models is to use wall functions to prescribe wall friction. Porté‐Agel et al. prescribed the surface shear stress and heat flux by adopting a time‐dependent variant of Monin‐Obukhov similarity theory. An alternative to no‐slip conditions are slip conditions, so that a laminar shear profile can be imposed without depending on the building‐up of a boundary layer associated with the Reynolds number of the atmosphere. A turbulent inflow is difficult to prescribe, because the turbulent viscosity of the ABL is much (10 6 times) larger than the molecular viscosity, making the effective Reynolds number of the blades much too low. At the upper boundary, symmetry conditions are often prescribed, but care should be taken when the height of the ABL is of the same order as the turbine height, something which is not unusual for modern‐size turbines, especially at night time. In streamwise direction, periodicity conditions can be used to simulate infinite wind farms or to generate a turbulence field with the precursor method discussed before, so that an ABL can form without the need for a very large domain. Periodicity conditions are also convenient from a numerical point of view because they allow the application of fast Poisson solvers like the fast Fourier transform. 3 ROTOR MODELING To solve the RANS equation or LES equation in the near and far wakes of a wind turbine, a representation of the blades is necessary. Basically, two approaches exist: the generalized actuator disk approach, in which the blades are represented by a body force ( Section 3.1 ), and the direct approach, in which the presence of the blades is taken into account by discretizing the actual blades on a computational mesh ( Section 3.2 ). 3.1 Generalized actuator disk modeling In many near and far wake calculations, the rotor is represented by an actuator disk or an actuator line. Such a representation circumvents the explicit calculation of the blade boundary layers, reducing computational cost and easing mesh generation. The actuator disk exerts a force on the flow, acting as a momentum sink. This force is explicitly added to the momentum equation : ∫ Ω ∂ u ∂ t d Ω + ∫ ∂ Ω u u ⋅ n d S = − ∫ ∂ Ω 1 ρ n d S + ∫ ∂ Ω v ∇ u ⋅ n d S + ∫ A ∩ Ω f d A which are written in weak form, because the force leads to a discontinuity in pressure. Apart from this momentum sink, one should also introduce sources of turbulence corresponding to the mechanical turbulence generated by the blades. Currently, three different approaches for prescribing the force term f exist: the actuator disk, the actuator line and the actuator surface models (see Figure ). Illustration of the actuator disk (AD), line (AL) and surface (AS) concept. 3.1.1 Actuator disk In case of a uniformly loaded actuator disk, f acts on the rotor‐disk surface A and is usually expressed in terms of the thrust coefficient C T only: ρ f = 1 2 ρ V ref 2 C T The determination of the reference velocity V ref in order to calculate C T is not obvious. For a turbine facing the undisturbed flow, V ref is evidently V ∞ , but for a turbine in the wake of an upstream turbine or in complex terrain, this is not the case. Prospathopoulos et al. proposed an iterative procedure to obtain the reference velocity and the thrust coefficient for downstream turbines modeled as actuator disks; starting with a certain V ref , one determines the thrust coefficient, from which the axial induction a follows, and then a new reference velocity based on the local flow field is computed: V ref = V local /(1 − a ). This procedure is repeated until convergence is achieved. Meyers and Meneveau and Calaf et al. used a similar approach by taking the local velocity and axial induction factor to determine the reference wind speed, but not in an iterative manner. The local velocity was obtained by disk averaging and time filtering (a LES model was used). For non‐uniformly loaded disks, the force is depending on radial position but constant over an annulus (see Figure ). Sectional lift and drag coefficients ( c l and c d ) are then used to find the local forces on the blades, like in BEM: ρ f 2 D = 1 2 ρ V rel 2 c c 1 e L + c d e D where e L and e D are unit vectors in the direction of lift and drag; c l and c d are functions of the Reynolds number and the angle of attack α. The force integral in equation then becomes a line integral; the disk is recovered as a time average of line forces. The relative velocity V rel at the radial positions is found by interpolating the velocity field in the surrounding computational cells. This is different from BEM, where V rel is found from an iterative procedure that employs a global momentum balance. Another difference with BEM is the application of a tip‐loss correction. The assumption of an infinite number of blades is corrected in BEM by locally changing the induced velocity. In Navier‐Stokes computations, this is not necessary because the flow field will notice the presence of the disk so that the induction changes automatically. However, the use of 2D airfoil data still requires a correction to obtain the right flow angle and flow velocity. A related problem is the determination of the local α to find c l and c d . Shen et al. developed a technique with which α can be determined based on information slightly upstream of the rotor. Rajagopalan et al. were one of the first to use the actuator‐type approach in a CFD code, for the calculation of vertical‐axis turbines. Time‐averaged forces were prescribed, and a finite‐difference laminar flow solver was used to solve the steady Navier‐Stokes equations. Ammara et al. and Masson et al. followed the time‐averaging approach of Rajagopalan in a control‐volume finite‐element setting. They included a second grid around the turbine in order to evaluate the surface force integral [see equation ]accurately and linearized the force term in an iterative procedure to a steady‐state solution. In Masson et al. , the lift and drag coefficients were obtained with a dynamic‐stall model. In this work, the tower was, like the rotor, also modeled as a porous surface. The forces on this surface were obtained from experimental data on the drag of a cylinder and were enforced by imposing a pressure discontinuity instead of adding the surface force in the equations. The presence of oscillations as a result of the use of collocated methods was mentioned (i.e. storing pressure and velocity variables at the same location), which was resolved by storing two different pressure values for the points located on the disk surface. Unsteady computations with the actuator disk approach were made by Sørensen et al. by using cylindrical coordinates in a rotor‐fixed reference frame. In Sørensen and Myken, a finite‐difference method was employed to solve the unsteady Euler equations in a vorticity‐stream function formulation (for advantages and disadvantages of such a formulation, see Vermeer et al. . A constant rotor loading was specified, but because of numerical difficulties at the disk edge, it was replaced by an elliptic distribution with equivalent total loading. In Sørensen and Kock, viscous terms were included, but only to stabilize the solution. The force on the non‐uniformly loaded actuator disk was obtained from tabulated airfoil data, and a non‐linear filter was applied to suppress oscillations in the vicinity of the disk. In Sørensen et al. , the vorticity‐velocity formulation was used to study different wake states, including windmill brake, turbulent wake, vortex ring and hover state. This formulation has the advantage of avoiding pressure‐velocity coupling issues, although a regularization kernel is still necessary to smoothly distribute the loading from the disk to the surrounding mesh points: f mesh = f disk * η ∈ The regularization function η ϵ is, for example, a Gaussian. The choice for ϵ, indicating the amount of smearing, is often a trade‐off between stability and accuracy. With this approach, the wiggles present in Sørensen and Kock disappeared. Madsen investigated both uniformly and non‐uniformly loaded actuator disks and the effect of turbulent mixing to show the validity of the BEM theory. It was found that BEM, with the application of a tip correction, gives a good correlation with the CFD results. In later work, Madsen et al. proposed a correction to BEM based on the comparison with actuator disk simulations, which showed that BEM overestimates the induction at the inboard part of the rotor and underestimates it at the outboard part. Apart from the uniform or non‐uniform axial loading described above, one can also introduce tangential forces on the disk surface to account for rotational effects. Meyers and Meneveau applied this in a LES context and showed that the effect of the tangential forces on the wake and extracted power appears to be negligible in the case of moderate power coefficient and high tip speed ratio. However, Porté‐Agel et al. and Wu and Porté‐Agel showed that the inclusion of rotation and non‐uniform loading leads to significant improvement in the prediction of the mean velocity and the turbulence intensity with respect to the uniformly loaded disk (see Figure ). This is especially apparent in the center of the near wake, where the uniformly loaded disk leads to an underestimation of the wake deficit and the turbulence intensity. Further downstream, the effect of rotation and non‐uniform loading disappears. Comparison of actuator disk without (dashed line) and with rotation (solid line), actuator line (solid black dots) and experiments (open red dots). Average horizontal velocity as function of height at different downstream positions. Reproduced from Porté‐Agel et al. Another approach to describe an actuator disk is the actuator shape model by Réthoré. The common cells of two different grids, one for the computational domain and one for the actuator (one dimension lower), were determined, and based on the intersecting polygons, forcing was applied to a cell. A comparison between this actuator model and a full‐rotor computation shows that modeling the wake by using forces is a good approximation for the mean flow quantities at distances larger than a rotor diameter from the wind turbine. It was observed that 10 cells per rotor diameter are sufficient; a similar number is typically found in other studies as well. However, the forces fail to represent the mechanical turbulence generated at the blade location. This turbulence can therefore be added at the disk location, independently of the actuator force, but its effect on the far wake in comparison with atmospheric and wake turbulence is small. As was mentioned, in actuator disk simulations, the boundary layers are not explicitly simulated, but their effect is taken into account via the lift and drag coefficients. Correctly simulating the blade Reynolds number is then less important, reducing the required computer resources considerably. It was shown in Sørensen et al. that the velocity field does not change noticeably when the Reynolds number is larger than 1000 (see also Snel, Ivanell et al. and Sørensen and Shen ). This corresponds roughly to results obtained by research on the nature of the interface of turbulent wakes and the outer flow, which changes around Re ∼ 10 4 . Whale et al. also found that the behavior of the wake is possibly rather insensitive to the blade Reynolds number. 3.1.2 Actuator line As an extension of the non‐uniformly loaded actuator disk approach, Sørensen and Shen introduced the actuator line approach (see Figure ). The line forces are not averaged over the disk but depend on time. Whereas in the actuator disk model vorticity is shed into the wake as a continuous sheet, in the actuator line model distinct tip vortices can be calculated. As for the non‐uniformly loaded actuator disk, the actuator line method requires knowledge of the lift and the drag on the blades. Corrections for Coriolis, centrifugal and tip effects are necessary when 2D airfoil data are used. In order to transfer the rotating line forces to the stationary mesh, a regularization kernel similar to equation is used. Mikkelsen investigated the actuator line method in detail and implemented it in EllipSys3D, a finite volume code for the solution of the incompressible Navier‐Stokes equations in pressure‐velocity formulation in general curvilinear coordinates. Howard and Pereira used point forces to represent the blades but did not take into account the distance between the force location and the cell centers, leading to high‐frequency noise in the power output. They modeled the tower as a square cylinder and found that it generates a wake that partially destroys the blade tip vortices. They recommend for future work to model the tower by point forces as well. 3.1.3 Actuator surface Shen et al. extended the actuator line method to an actuator surface method and applied it to vertical‐axis wind turbines. Whereas in the actuator line model the blade is represented by a line, in the actuator surface model it is represented by a planar surface (see Figure ). This requires more accurate airfoil data; instead of c l and c d , knowledge of the pressure and skin‐friction distribution on the airfoil surface is needed: f 2D AS ξ = f 2D F dist ξ where ξ is directed along the chord and F dist (ξ) is determined by fitting empirical functions to chordwise pressure distributions. These were obtained in Shen et al. with Xfoil, a highly accurate tool to compute pressure‐friction and skin‐friction profiles on airfoils. A comparison with 2D RANS calculations on a body‐fitted mesh around an airfoil shows that the pressure field can be accurately computed approximately one chord away from the airfoil. Furthermore, it was shown that the use of a dynamic‐stall model in vertical‐axis wind turbine calculations significantly improves the agreement with experiments. In Shen et al. , the method was applied to 3D turbine calculations and compared with the actuator line technique. Some improvements were seen in the representation of tip vortices and the flow behavior near the airfoil surface. In the actuator surface approach of Dobrev et al. , the pressure distribution was represented by a piecewise linear function that is directly calculated from the lift and drag coefficients. As in Shen et al. , the flow field was found to be more realistic compared to the actuator line approach, but the method was not yet able to model the wake of an airfoil since the shear forces on the blade were not accounted for. Leclerc, Watters and Masson used a slightly different approach with their actuator surface concept. Inviscid aerodynamic theory was used to relate vorticity distributions to both pressure and velocity discontinuities across a porous surface. This bears resemblance with vortex methods and lifting‐line theory, which describe an inviscid flow field by concentrated vortex sheets or lines. Given the lift coefficient, the circulation was computed, and subsequently, the velocity discontinuity over the surface was obtained by assuming a parabolic or cubic distribution. Viscous drag was not taken into account, but in three dimensions, the actuator surface could still extract energy from the flow because of induced drag. In two dimensions, an actuator surface could not extract energy, and a rotor disk was then modeled by prescribing the vorticity distribution of the slipstream surface, instead of the momentum loss through the disk area. The use of surface forces instead of volume forces was found to be the reason that the solution did not exhibit spurious oscillations. 3.1.4 Conclusions To conclude, a hierarchy of actuator models has emerged, from ‘simple’ actuator disks to more accurate actuator line and surface models. Gaining accuracy is accompanied by higher computational cost and the need for more detailed airfoil data: from C T (uniform actuator disk) to c l and c d (non‐uniform actuator disk, actuator line) to C p and C f (actuator surface). It would seem that the unsteady nature of actuator line and surface methods makes them most suitable for LES simulations and that the steady nature of actuator disk methods limits their application to RANS simulations. However, the computational costs limit the use of the actuator line technique to single‐wake investigations, and most LES simulations of wind farms are being performed with actuator disks. 3.2 Direct modeling The complete or direct modeling of the rotor by constructing a body‐fitted grid is physically the most sound method to compute the flow around a turbine. The focus is here on work that has contributed to the knowledge of wake aerodynamics. Compared to the generalized actuator disk approach, the blade is represented ‘exactly’, instead of a disk/line/surface approximation. However, accurately simulating the boundary layer on the blades, including possible transition, separation and stall, is computationally very expensive. Additionally, compressibility effects at blade tips can require the solution of the compressible Navier‐Stokes equations, whereas the wake remains essentially incompressible. Furthermore, the generation of a high‐quality moving mesh is not trivial. Mesh generation is therefore commonly done with so‐called ‘overset’ or ‘chimera’ grids: different overlapping grids, often structured, that communicate with each other. An example of such a grid is shown in Figure . Layout of overset grids around tower, blade and far field. Reproduced from Zahle et al. The first simulations (for wind turbine applications) with direct modeling were done by Sørensen and Hansen, employing a rotating reference frame and the k −ω SST model. The rotor power is predicted well for wind speeds below 10 m s −1 , but at higher speeds, the power is underpredicted. The strongly separated flow on the blade is not correctly captured at these speeds, which the authors attribute to insufficient mesh resolution and limitations of the turbulence model. Duque et al. used the compressible thin‐layer Navier‐Stokes equations with the Baldwin‐Lomax (algebraic) turbulence model. A comparison with the NREL phase II rotor showed good agreement of the pressure distributions, but the rotor‐tower interaction was not predicted well. In later work, the full equations were solved with the Baldwin‐Barth (one‐equation) turbulence model and compared with the NREL phase VI rotor, showing good agreement with experiments, including stalled flow and cross flow. Different transition models were investigated by Xu and Sankar with a hybrid Navier‐Stokes potential model. In this approach, the viscous compressible flow equations are solved only in a small region surrounding the rotor, whereas the rest of the flow field is modeled using an inviscid free‐wake method. The position of the transition line on the blades was shown to clearly depend on the transition model. Benjanirat and Sankar used an overset grid model to study the effect of different turbulence models. All the models predicted the out‐of‐plane forces and associated bending moments well, but difficulties were found in predicting the in‐plane forces and therefore in predicting power generation. This can probably be explained by the fact that lift prediction is less sensitive to turbulence modeling than drag prediction. Sørensen (with the k −ω SST model) and Johansen (DES) performed simulations for the NREL phase VI rotor with a rotor‐fixed reference frame. No evidence was found that the DES computations improved the RANS results in predicting blade characteristics. Madsen et al. compared direct modeling with a generalized actuator approach and concluded that the local flow angle is generally better predicted by the direct model. In computations by Johansen and Sørensen, the full 3D solution was used to extract airfoil characteristics. A method was developed to find the local angle of attack, based on earlier work of Hansen et al. The resulting c l −α and c d −α curves can be used as an input for BEM calculations and can improve actuator line or surface computations. Shen et al. devised a method to extract the angle of attack for more general flow conditions. For further improvement, BEM calculations in shear were performed, and it was found that although the power variation over a revolution was small, the variation of individual blade forces was much higher. The effect of the rotor on the velocity profile was investigated, and it was shown that disturbances in the velocity profile were still noticeable at five rotor diameters (5 D ) upstream. In later work, the same computations were repeated with an overset grid that indicated that the induction was not felt anymore more than 2 d upstream, showing that the rotor‐fixed reference frame approach is less appropriate when dealing with sheared inflow. Zahle and Sørensen investigated the influence of the tip vortices on the velocity at the rotor plane and the influence of the tower in downwind and upwind configurations, using the k −ω SST model and a fully turbulent flow. A reduction in thrust and torque of 1 − 2% because of the tower shadow was reported, as well as signs of dynamic‐stall behavior on the blades when exiting the shadow. Their results showed the breakdown of the root vortices after 3 d and the tip vortex breakdown after 4–6 d. Furthermore, it turned out that fine grids are needed to track the tip vortices and that the rotor thrust and torque are hardly affected by the far wake resolution. Investigations on transition modeling were done by Sørensen using the k −ω model and the Langtry‐Menter transition model. It was found that transitional computations lead to better agreement with experiments than fully turbulent simulations. A similar comparison with the same turbulence and transition model was made by Laursen et al. This study also showed that the application of the transition model leads to a more realistic performance of the blade: increased lift and lowered drag. Near‐wake calculations without transition model were made by Bechmann and Sørensen on the MEXICO rotor (see Section 5 ). The axial velocity in the near wake corresponds well with measurements, but in the far wake (from 2.5 d), coarsening of the mesh leads to an unphysical velocity increase. The radial velocity at 1 d downstream compares reasonably well with measurements. The blade pressure distributions are used to extract 3D airfoil characteristics (i.e. sectional c l and c d when 3D flow effects are present) that can, as mentioned before, be used for different wake models. Calculations by Wußow et al. were aimed at wake characteristics and resolved the entire wake but used a coarse resolution on the blade to limit computational efforts. An exception to the structured overset grid approach is the work of Sezer‐Uzol et al. who used a moving unstructured grid. This circumvents the interpolation problem inherent to the (structured) overset grid approach and allows for grid coarsening or refinement, but the continuous re‐meshing can be expensive. 4 WAKE MODELING By having discussed the possible ways to model the rotor in the previous section, this section discusses the different ways of modeling the wake. As mentioned before, the simulation of turbulence is of major importance in CFD of wind turbine wakes, and therefore, the work done in this area will be grouped according to the two turbulence techniques encountered in wake modeling, being RANS and LES. 4.1 RANS Although the application of the averaging procedure to the Navier‐Stokes equations, resulting in the RANS equations, and the use of the generalized actuator approach lead to a significant reduction in computational effort, solving equation is still a formidable job for wake calculations. One of the main reasons for this is the divergence‐free constraint in equation , which requires that the pressure is calculated implicitly via an elliptic equation, the pressure Poisson equation. Two simplifications to this approach, parabolization and linearization, and the original elliptic models are discussed in this section. 4.1.1 Parabolic models The first wake calculations based on RANS used axial‐symmetric and parabolic forms of the Navier‐Stokes equations. Since most parabolic models were already reviewed quite extensively in Vermeer et al. and Crespo et al. , only the main developments will be mentioned here. The parabolic model, also called boundary‐layer wake model, was obtained by neglecting both the diffusion and the pressure gradient in streamwise direction. This allowed for a fast solution by space marching, but the (pressure driven) expansion of the wake cannot be predicted properly. Ainslie used a parabolic eddy viscosity model with axial symmetry, steady flow and neglected the pressure gradients in the outer flow. Since the parabolicity assumption is not valid in the near wake, the calculations started at 2 d downstream of the rotor, where the wake is assumed to be fully developed. At that location, a Gaussian velocity profile was used as a starting condition for the far wake, and the solution was computed by advancing in the wind direction. Improvements have been made by using a variable‐length near wake and by refining the initial condition for the far wake computation by linking it to a BEM method. Reasonable results were obtained when compared to wind tunnel experiments, but because of the assumption of axial symmetry, the downward movement of the wake centerline and the upwind movement of the location of maximum turbulence intensity cannot be predicted. A parabolic model that does not assume axial symmetry is UPMWAKE, by Crespo et al. The k −ϵ turbulence model was employed, with the constants adapted for neutral ABL flow. The presence of the turbine was taken into account by imposing perturbation values of velocity, temperature, kinetic energy and dissipation at the turbine location. The UPMWAKE model was the foundation for the WAKEFARM model of the Energy research Centre of the Netherlands (ECN). This model was subsequently improved by including the axial pressure gradient via a boundary‐layer analogy and coupling it to the near wake via a vortex‐wake method, so that the expansion of the wake is taken into account. A version of UPMWAKE that accounts for the anisotropic nature of the turbulence in the atmosphere and the wake, UPMANIWAKE, was made by Gómez‐Elvira et al. An explicit algebraic model for the components of the turbulent stress tensor was used. The near wake was assumed to be fully turbulent, and the analysis was carried out for a neutral atmosphere. Improved results with respect to isotropic models were obtained, and analytical expressions for the added turbulence intensity in the wake as function of thrust coefficient were derived. They suggested that future improvements can be made by using LES. The UPMWAKE model was extended to UPMPARK to calculate parks with many turbines. In a parabolic setting, this is a sequence of single‐wake computations, where the outflow of one wake forms the inflow condition of the next. 4.1.2 Elliptic models Crespo and Hernández extended the parabolic code to a fully elliptic version to resolve discrepancies with experimental results that are mainly found in the near‐wake region, but no significant improvements were observed in single‐wake computations. Rajagopalan et al. used an elliptic approach to study interaction of vertical‐axis wind turbines in a farm with a laminar solver, from which it is concluded that the mutual influence of turbines is significant. Ammara et al. used an algebraic closure (Prandtl' mixing length) with values from Panofsky and Dutton. In single‐wake computations, a good prediction of the velocity deficit is found. In farm simulations, the existence of a ‘Venturi’ effect is demonstrated, positively influencing the power output of a farm. Different researchers found that the standard k −ϵ and k −ω models result in too diffusive wakes: the velocity deficits are too small, and the turbulence intensity does not show the distinct peaks observed in experiments. The reasons for the failure of these standard turbulence models in wind turbine wakes were explained by Réthoré and are basically caused by the limited validity of the Boussinesq hypothesis, as mentioned before. To improve correspondence with experimental data, several adaptations of the k −ϵ model have been suggested. These adaptations all strive to reduce (directly or indirectly) the eddy viscosity and hence the diffusion, in the near wake. Firstly, El Kasmi and Masson added an extra term to the transport equation for the turbulent energy dissipation (ϵ) in a region around the rotor, based on the work of Chen and Kim. This method leads to the introduction of two additional parameters: a model constant and the size of the region where the model is applied. Compared with experimental data, significant improvements over the original k −ϵ model were observed. Cabez‐n et al. and Rados et al. have also applied this approach and found acceptable agreement with other experimental data. Réthoré noticed however that increasing the dissipation proportionally to the production of turbulence contradicts the LES results. A second approach is the realizability model. In this model, the eddy viscosity is reduced to enforce that Reynolds stresses respect the so‐called realizability conditions (see e.g. Schumann ). According to Réthoré, these conditions are not satisfied in the near wake by the eddy viscosity‐based k −ϵ model, because of the large strain rate at the edge of the wake, where the Boussinesq hypothesis is inadequate. Both Réthoré and Cabezón used a realizable model based on the work of Shih. Compared with the standard k −ϵ model, better results are obtained with the realizable model, although the prediction of both wake deficit and wake spreading and turbulence intensity remains unsatisfactory. More adaptations are under investigation. Prospathopoulos et al. adapted the constants of the k −ϵ and k −ω model for atmospheric flows and investigated the effect of complex terrain. Rados et al. used Freedman's model to change the constants of the ϵ equation in case of stable atmospheric conditions. Prospathopoulos et al. obtained good results with Durbin's correction, in which the turbulent length scale is bounded. Another possibility comes from an analogy with forest and urban canopy simulations, where obstacles are also modeled as body forces. Réthoré adapted the k and ϵ equations to take into account the extraction of turbulent kinetic energy by the actuator and found significant improvement with respect to the standard k −ϵ model. The RSM approach, not relying on the Boussinesq hypothesis, was employed by Cabezón et al. and gave more accurate results than both the El Kasmi‐Masson model and the realizability model. However, in predicting the velocity deficit in the near wake, the parabolic UPMPARK code outperformed the elliptic models because in the parabolic code, the streamwise diffusive effects are totally neglected. In the far wake, the velocity deficit is predicted similarly to the elliptic models, but the turbulence intensity is overestimated and does not show the distinct peaks. In case of stratification, the production of turbulence because of buoyancy has to be taken into account in turbulence closure models. Most researchers use the constants of the k −ϵ model proposed by Crespo et al. for neutral atmospheric flows. For non‐neutral stratification, Alinot and Masson changed the closure coefficient of the k −ϵ model related to buoyancy production to match the basic flow described by Monin–Obukhov similarity theory; no wind turbines are present. With respect to Crespo and Hernández, this gives improved results for stable atmospheric conditions, especially when considering profiles of turbulent kinetic energy and dissipation. Sanz et al. found that CFD modeling with Monin–Obukhov theory is not valid for stable conditions in offshore situations. 4.1.3 Linearized models Ott et al. used the k −ω model employing the canopy approach mentioned before, in combination with the fast linearized model of Corbett et al. The linearization of the Navier‐Stokes equations, combined with a spectral method and lookup tables (containing a.o. velocity profiles), allows for a very fast solution of entire wind farms. In such a spectral method, the governing equations are transformed to Fourier space, which has the advantage that no computational grid and therefore no artificial diffusion is present, although restrictions on the boundary conditions apply (typically periodicity in lateral direction). The results for the velocity field are most accurate in the far wake, where the assumption of small perturbations (linearization) is justified. However, compared to a full‐model prediction, the turbulent kinetic energy shows serious mismatches in both the near and far wakes. 4.2 LES As mentioned before, LES has the advantage over most RANS models in that it is better able to predict the unsteady, anisotropic turbulent atmosphere. Jimenez et al. used a dynamic subgrid‐scale model with the rotor represented by a uniformly loaded actuator disk. A good comparison with experiments was found. The same model was used in a paper by Jimenez et al. to study the spectral coherence in the wake and in another paper by Jimenez et al. to study the wake deflection because of yaw; in both cases, reasonable agreement with experiments and an analytical model was observed. Wußow et al. used a Smagorinsky‐Lilly subgrid‐scale model with a wind field generated with a von Kármán spectrum to study loading of wind turbines in a farm. It was shown that when calculating the mean velocity field, the averaging time must be long enough to obtain a good comparison with 10 min averaged experimental data. Better meshing strategies led to improvements of velocity deficit and turbulence intensity. Troldborg et al. used a vorticity‐based mixed‐scale subgrid‐scale model with the actuator line technique to study the influence of shear and inflow turbulence on wake behavior. The ABL was imposed with the force‐field technique of Mikkelsen, and atmospheric turbulence was generated with Mann's method (see Section 2.3 ). It was found that shear and ground effects cause an asymmetric development of the wake with larger expansions upward and sideward than downward. Inflow turbulence, when compared with laminar inflow, destabilized the vortex system in the wake, resulting in an improved recovery of the velocity deficit but also an increased turbulence intensity. Furthermore, it was found that the wake is also unstable and becomes turbulent for uniform, laminar inflow, especially at low tip speed ratios. In the work of Ivanell, the same LES model was used for farm simulations. Twenty non‐uniformly loaded actuator disks, in combination with periodic boundary conditions, were used to simulate the 80 turbines that comprise the Horns Rev wind farm (Figure (a)). The power production of downstream turbines agreed reasonably well with experimental data (Figure (b)). LES was also used to study the stability of tip vortices. By inserting harmonic perturbations at blade tips, an exponential relationship between the perturbation frequency (related to turbulence intensity) and the length of the stable part of the wake was observed. LES of the flow through 20 actuator disks, representing the Horns Rev wind farm. Reproduced from Ivanell. Meyers and Meneveau performed LES with a (non‐dynamic) Smagorinsky model to study infinite arrays of wind turbines in staggered and non‐staggered arrangement, and it was observed that the staggered arrangement had a higher power production. Porté‐Agel et al. and Wu and Porté‐Agel et al. employed a Lagrangian scale‐dependent SGS model, a parameter‐free model that performs better in ABL simulations than the traditional Smagorinsky and standard dynamic models. Similar to what Réthoré showed for C μ in the k −ϵ model, Wu and Porté‐Agel showed that the Smagorinsky coefficient C S is not constant in the wake but increases in the center of the wake and decreases near the ground and in the shear layer. With an actuator‐type approximation for the turbine, this model provides a very good agreement with atmospheric wind tunnel data. Furthermore, the effect of a wind turbine on a stable ABL was investigated, demonstrating that the momentum and buoyancy flux at the surface are reduced and as a result possibly influence the local meteorology. A comparison between RANS (standard k −ϵ) and LES (method of Bechmann ) was made by Réthoré in order to explain the aforementioned discrepancies between the k −ϵ model and the experimental data. The LES results are clearly superior to RANS when comparing both mean velocity profiles and turbulence quantities, but the computational time increases from hours (RANS) to days (LES) for a single‐wake case. In the study of Stovall et al. , the standard k −ϵ model and LES with a one‐equation SGS model were compared. Again, it was observed that the LES results are closer to experimental data in a single‐wake calculation and that the wake recovery because of turbulent mixing with the outer flow is much better captured by LES in a multiple‐wake situation. The difference in computational time was a factor 60, when the same grid for both RANS and LES was used. However, a well‐resolved LES should require a much finer grid than RANS, making it more expensive, as will be explained in the following. An estimate of the necessary grid resolution can be obtained based on the Reynolds number. For a sufficiently resolved LES, one needs cell sizes far enough in the inertial subrange, i.e. of the order of the Taylor microscale. The fact that the Taylor microscale scales with Re −1/2 and the Reynolds number based on the diameter of a modern‐size wind turbine is O 10 8 leads to the requirement of cell sizes around 1 cm 3 . This is somewhat conservative; for atmospheric flows (based on an integral length scale of 1 km), the Taylor microscale was estimated at 10 cm. Most of the LES computations to date use cell sizes of approximately 1–10 m, which might be too coarse to resolve enough scales. Computations require typically 10 7 –10 8 grid points and run on supercomputers for several days or weeks, even for single‐wake situations. With the continuous advance in computer power, finer meshes will be possible, but the associated increase in data analysis will leave this approach unattractive for many engineering purposes. 4.3 Numerical issues Most RANS codes use second‐order accurate finite volume schemes on structured meshes, with upwind discretization of the convective terms and central discretization of the diffusive terms, leading to stable and robust schemes. Implicit methods are normally used to find a steady‐state solution. In LES, the temporal and spatial discretization of the convective terms should be done more carefully. The use of upwind schemes for the spatial discretization can influence the energy cascade from large to small scales because of the introduction of numerical dissipation. Several authors have indeed mentioned premature turbulence decay and related it to possible numerical dissipation (see e.g. Bechmann, Stovall et al. and Troldborg et al. . This forms the rationale for the application of central or spectral schemes (see e.g. Bechmann and Meyers et al. ). Apart from the spatial discretization, the temporal discretization can also introduce numerical dissipation. In Bechmann, this limits the time step, partly removing the advantages of an implicit method. High‐order low‐dissipative explicit schemes, such as the standard four‐stage Runge‐Kutta method, used in Calaf et al. can then be an attractive alternative. Lately, high‐order (typically fourth order) schemes have also been employed for the spatial discretization, reducing interference between the subgrid‐scale model and the discretization. A blend with third‐order upwind schemes is sometimes made to ensure numerical stability. However, it should be realized that the formal order of accuracy of a discretization scheme is only obtained in the limit of very high spatial resolution, something typically not encountered in LES of wind turbine wakes. Mesh refinement studies are often computationally too expensive, and one has to rely on the energy spectrum to check if the inertial subrange is captured well. It should be noted that partial cancelation of subgrid modeling and numerical errors may occur, leading to the counterintuitive conclusion that high‐order accurate schemes, improved subgrid models or finer meshes, may lead to worse results. 5 VERIFICATION AND VALIDATION The process of testing a CFD model should consist of two steps: verification (‘solving the equations right’) and validation (‘solving the right equations’). The first step, verification, is necessary to evaluate the errors that are made by approximating the continuous equations by a discrete set. This gives an indication of the mesh size necessary to reach a required accuracy. Highly accurate numerical simulations (benchmark tests) or analytical solutions may serve as verification data. An example of such tests is the actuator disk simulation performed by Réthoré: the flow through lightly loaded (2D and 3D ) and heavily loaded actuator disks. The second step, validation, requires a comparison with experimental data to investigate the validity of the model and to assess modeling errors. Generally, the turbulence model, the blade approximation and the inflow conditions are the main sources. Wind tunnel tests on rotating wind turbine models yield very suitable experimental data because they take place in a constant and controllable environment. However, the wind tunnel conditions are typically different from atmospheric conditions (scaling and blocking effects, turbulence intensity). One can therefore only conclude how well the CFD model simulates the wind tunnel conditions, not how well it simulates the ‘real’ flow through turbines standing in the field. An elaborate overview of the near and far wake wind tunnel experiments is given in Vermeer. Since that review, three important wind tunnel tests have been performed. First is the tests in the German‐Dutch wind tunnel DNW in 2006, on a 4.5 m diameter rotor, the so‐called MEXICO project. Near‐wake measurements with particle image velocimetry (from ‐1 d to 1.5 d) and pressure recordings on the blades were made, but no far wake data are available. The measurements are evaluated and compared to simulations in the (currently ongoing) MexNext project, IEA Annex 29 [ http://www.mexnext.org ]. Second, the experiments from Chamorra and Porté‐Agel with hot‐wire measurements in an atmospheric wind tunnel provide high‐resolution spatial and temporal measurements that can be used to validate LES models. Thirdly, there are the experiments on wake meandering in the ABL wind tunnel at the University of Orléans, with wind turbines modeled as porous disks. Besides wind tunnel data, full‐scale field tests are a valuable source of experimental data; they have been cataloged in a database by IEA, Annex XVIII. They are the only way to obtain data on full‐scale wake development and wake‐wake interaction. A common problem of field measurements is that the inflow is often not fully known because the number of measurement stations (giving point data) is limited; recent advances in the LIDAR technique can significantly increase the knowledge of the entire flow field. Well‐known single‐wake experiments are the Tjæreborg turbine, the Nibe turbine and the Sexbierum measurements. Well‐known sites for multiple‐wake experiments are Horns Rev, Nysted, Middlegrunden, Vindeby and the ECN Wind Turbine Test Site Wieringermeer. Valuable comparative studies between different wake models and experiments at these sites have been obtained in the ENDOW project and, more recently, the UPWIND project. Some results of the latter project are shown in Figure , where experimental data and numerical simulations of the Horns Rev wind farm are compared. The numerical simulations have been made with engineering models (WAsP from Risø, WindFarmer from Garrad Hassan) and RANS models, namely a parabolic code (ECN's WAKEFARM) and an elliptic code (NTUA). LES of the same wind farm were already shown in Figure . First, it should be noted that the standard deviation in experimental data is relatively large, as indicated by the error bars. Taken this into account, the agreement between models and experiment is good, typically within 15% of the experiments, and it is remarkable that the engineering models are at least as accurate as the RANS and LES models, while having much lower computational costs. Comparison of models and experiments in predicting the power output of Horns Rev for different wind directions. At 270°, the wind is aligned with the turbines. Reproduced from Barthelmie et al. In the recent Bolund experiment, linearized models, RANS models and LES models were compared in predicting the flow over a hill, without the presence of wind turbines. Linearized models were found to be unfavorable because of the steepness of the hill. RANS models performed best when compared with experimental data, with similar trends between a large range of different codes and turbulence models. A much larger spread was observed in the results obtained by LES models, and it was concluded that LES has not matured enough for widespread application to atmospheric flow simulation. 6 CONCLUSIONS A growing number of researchers are using CFD to study wind turbine wake aerodynamics. CFD has become a mature tool for predicting a wide range of flows; however, one important ongoing challenge is the accurate representation of turbulence. Since in wake aerodynamics turbulence is a dominating factor, affecting both blade performance and wake behavior, most efforts are now directed at improving turbulence modeling. In the near wake, full 3D Navier‐Stokes solvers are used with overset grid approaches to resolve the boundary layers on the rotating blades. Such computations are able to predict separation and transition to reasonable accuracy, although the investigation of the effect of different turbulence models is an ongoing work. This direct approach is not yet suitable for optimizing turbine blades or optimizing farm layouts. The current contribution of the direct modeling approach is to increase the understanding of the flow physics on the blades and to improve simpler models such as the BEM method. To reduce computational requirements for wake simulations, the presence of the blades is in general taken into account with the generalized actuator disk approach, in which the rotor is represented by forces. A careful numerical treatment is necessary to correctly handle these forces. Three methods are currently in use, being the actuator disk, actuator line and actuator surface methods. The more expensive actuator line and surface techniques use more detailed aerodynamic blade characteristics than the actuator disk, and as a result, they are more accurate. The actuator disk and line techniques are maturing, and current research is on the actuator surface technique, but because of the lower computational effort, the actuator disk remains the most widely used method for multiple‐wake simulations. The use of the actuator approach seems justified because in the far wake, differences with the direct modeling approach are shown to be small. However, a detailed study of the effect of the blade geometry on the wake and a comparison with actuator line simulations are still lacking. For modeling turbulence in the wake, RANS will most likely prevail to be the engineer's choice, even though many eddy viscosity‐based models like k −ϵ proved to be too diffusive. Ongoing research is now in different adaptations and corrections to reduce diffusive effects in the near wake. Increased interest is in the application of LES, a method that naturally takes into account the unsteadiness and the anisotropy of the turbulence in the atmosphere, which is not properly done in most RANS methods. In single‐wake computations, LES shows indeed an improved agreement with experimental data compared with RANS, but the increase in the required computational effort, which is about two orders of magnitude, still prohibits its widespread use for farm applications, so that its use is currently confined to gain physical insight and improve simpler methods. Because of the high computational costs associated with fully resolved LES, mainly under‐resolved simulations are available nowadays. Thorough mesh refinement and subgrid‐scale model studies are so far still open problems. The use of accurate wind tunnel data in subgrid‐scale model studies is indispensable; wake measurements on large turbines in atmospheric tunnels are a necessity. An impression of the magnitude of subgrid modeling errors and discretization errors is important since it gives the opportunity to (partially) quantify the uncertainty in CFD results. Indeed, in a time where CFD codes are becoming readily available for wind turbine wake computations and increased computer power enables larger problems and more refined modeling, a possible next step is the quantification of uncertainties in the computations. These uncertainties originate not only from discretization and turbulence modeling errors but also from the description of the inflow, terrain geometry, rotor geometry, etc. A quantification of uncertainties would make the comparison with experimental data more fair and will give a guideline in which areas the CFD of wind turbine wakes has to be improved. In this respect, the reduction of the uncertainty in the prescription of atmospheric inflow conditions remains a great challenge. More research on the effect of stratification on power production is to be expected. A coupling with mesoscale atmospheric models can shed more light on appropriate boundary conditions and, at the same time, can be used to investigate the effect of wind farms on the local meteorology. ACKNOWLEDGEMENTS Colleagues A.J. Brand, P.J. Eecen, J.G. Schepers and H. Snel from ECN are acknowledged for their help in supplying literature and correcting the manuscript. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Wind Energy Wiley

Review of computational fluid dynamics for wind turbine wake aerodynamics

Wind Energy , Volume 14 (7) – Oct 1, 2011

Loading next page...
 
/lp/wiley/review-of-computational-fluid-dynamics-for-wind-turbine-wake-0N6TkCsYEG

References (132)

Publisher
Wiley
Copyright
Copyright © 2011 John Wiley & Sons, Ltd.
ISSN
1095-4244
eISSN
1099-1824
DOI
10.1002/we.458
Publisher site
See Article on Publisher Site

Abstract

1 INTRODUCTION During the last decades, wind turbines have been installed in large wind farms. The grouping of turbines in farms introduces two major issues: reduced power production, because of wake velocity deficits, and increased dynamic loads on the blades, because of higher turbulence levels. Depending on the layout and the wind conditions of a wind farm, the power loss of a downstream turbine can easily reach 40% in full‐wake conditions. When averaged over different wind directions, losses of approximately 8% are observed for onshore farms and 12% for offshore farms (see e.g. Barthelmie et al. ). When studying power losses and blade loading, wind turbine wakes are typically divided into near and far wakes. The near wake is the region from the turbine to approximately one or two rotor diameters downstream, where the turbine geometry directly affects the flow, leading to the presence of distinct tip vortices. Tip and root vortices lead to sharp gradients in the velocity and to peaks in the turbulence intensity. For very high tip speed ratios, the tip vortices form a continuous vorticity sheet: a shear layer. The turbine extracts momentum and energy from the flow, causing a pressure jump and consequently an axial pressure gradient, an expansion of the wake and a decrease of the axial velocity. In the far wake, the actual rotor shape is only felt indirectly, by means of the reduced axial velocity and the increased turbulence intensity. Turbulence is the dominating physical process in the far wake, and three sources can be identified: atmospheric turbulence (from surface roughness and thermal effects), mechanical turbulence (from the blades and the tower) and wake turbulence (from tip vortex breakdown). Turbulence acts as an efficient mixer, leading to the recovery of the velocity deficit and a decrease in the overall turbulence intensity. Far downstream, the velocity deficit becomes approximately Gaussian, axisymmetric and self‐similar. Wake meandering, the large‐scale movement of the entire wake, might further reduce the velocity deficit, although it can considerably increase fatigue and extreme loads on a downwind turbine. It is believed to be driven by the large‐scale turbulent structures in the atmosphere. The distinction between near and far wakes is also apparent when classifying the existing numerical models for wind turbine wake aerodynamics (see Table ). The first and simplest approach is an analytical method that exploits the self‐similar nature of the far wake to obtain expressions for the velocity deficit and the turbulence intensity. The second, Blade Element Momentum (BEM) theory, uses a global momentum balance together with two‐dimensional (2D) blade elements to calculate aerodynamic blade characteristics. The vortex‐lattice and vortex‐particle methods assume inviscid, incompressible flow and describe it with vorticity concentrated in sheets or particles. Panel methods similarly describe an inviscid flow field, but the blade geometry is taken into account more accurately, and viscous effects can be included with a boundary‐layer code; the wake follows as in vortex‐wake methods. These four methods have been extensively discussed in previous reviews, such as Vermeer et al. , Crespo et al. , Snel and Hansen et al. The last two methods, the generalized actuator disk method and the direct method, are relatively new and are commonly called computational fluid dynamics (CFD) methods. In this survey, we will review these methods, discuss their ability to predict wind turbine wakes and give an outlook into possible future developments. Classification of models. Method Blade model Wake model Kinematic Thrust coefficient Self‐similar solutions BEM Actuator disk + blade element Quasi one‐dimensional momentum theory Vortex lattice, vortex particle Lifting line/surface + blade element Free/fixed vorticity sheet, particles Panels Surface mesh Free/fixed vorticity sheet Generalized actuator Actuator disk/line/surface Volume mesh, Euler/RANS/LES Direct Volume mesh Volume mesh, Euler/RANS/LES RANS, Reynolds‐averaged Navier‐Stokes; LES, large eddy simulation. The paper is organized as follows. First, we will consider the Navier–Stokes equations and discuss their use to predict turbulent flows ( Section 2 ). Section 3 then discusses rotor modeling, Section 4 discusses wake modeling and Section 5 deals with the question how to verify and validate wind turbine wake CFD codes. 2 GOVERNING EQUATIONS 2.1 The incompressible Navier–Stokes equations It is reasonable to assume that the flow field in wind turbine wakes is incompressible, since the velocities upstream and downstream of a turbine placed in the atmosphere are typically in the range of 5–25 m s −1 . Only when calculating the aerodynamics at blade tips that compressibility effects may be important. Since in most calculations of wind turbine wakes the rotor is not modeled directly (which will be discussed later), the incompressible Navier–Stokes equations are a suitable model to describe the aerodynamics of wind turbine wakes: ∇ ⋅ u = 0 , ∂ u ∂ t + ( u ⋅ ∇ ) u = − 1 ρ ∇ p + ν ∇ 2 u , supplemented with initial and boundary conditions, which will be discussed in Section 2.3 . In the case of a non‐neutral atmosphere, the Boussinesq approximation is typically employed to account for buoyancy effects, and an extra equation for the temperature has to be solved. The effect of the rotation of the Earth, given by the Coriolis term, is typically neglected in many wake studies but can have an effect when computations involve large wind turbines and wind farms (see e.g. Porté‐Agel et al. ). Although this set of equations provides a complete model for the description of turbulent flows, it is not easily solved. The difficulty associated with turbulent flows is the presence of the non‐linear convective term, which creates a wide range of time and length scales. For example, in the atmospheric boundary layer (ABL), the largest turbulent scales are of the order of 1 km, whereas the smallest scales are of the order of 1 mm. Inside the blade boundary layers, the scales are even smaller. The range of scales depends on the Reynolds number ( Re), the dimensionless parameter that indicates the ratio of convective forces to viscous forces in the flow. Large values of the Reynolds number, encountered in blade and wake calculations, lead to a large range of scales, making computer simulations extremely expensive. Resolving all scales in the flow, so‐called direct numerical simulation (DNS), is therefore not feasible. Turbulence models need to be constructed, modeling the effect of the unresolved small scales based on the behavior of the large scales. However, even with the cost reduction provided by a turbulence model, one cannot resolve both the boundary layers on the turbine blades and the turbulent structures in the wake. This necessitates a simplified representation of the wind turbine in case of wake calculations and a simplified representation of the wake in case of blade calculations. 2.2 Turbulence modeling A large number of turbulence models have been constructed in the last decennia (see e.g. Wilcox, Sagaut and Geurts ). This section will discuss the two most important methodologies in turbulence modeling for wind turbine wakes, RANS and LES, their applicability and their limitations. 2.2.1 RANS Reynolds‐averaged Navier–Stokes methods aim for a statistical description of the flow. Flow quantities such as velocity and pressure are split in an average and a fluctuation, the so‐called Reynolds decomposition: u ( x , t ) = u ¯ ( x ) + u ′ ( x , t ) The averaging procedure, ensemble averaging, is such that u ¯ ¯ ( x ) = u ¯ ( x ) and u ′ ¯ ( x , t ) = 0 . The Reynolds decomposition [equation ] is substituted into the Navier–Stokes equations, which are then averaged, resulting in ∂ u ¯ ∂ t + u ¯ ⋅ ∇ u ¯ = − 1 ρ ∇ p ¯ + ν ∇ 2 u ¯ − ∇ ⋅ u ′ u ′ ¯ The term u ′ u ′ ¯ is called the Reynolds stress tensor, which appears as a consequence of the non‐linearity of the convective term, and represents the averaged momentum transfer because of turbulent fluctuations. The Reynolds stresses can be interpreted as turbulent diffusive forces. In wind turbine wakes, they are much larger than the molecular diffusive forces v ∇ 2 ū , except near solid boundaries. In order to close the system of equations, a model is needed to express the Reynolds stresses in terms of mean flow quantities. A widely adopted approach of modeling the Reynolds stresses exploits the Boussinesq hypothesis (not to be confused with the Boussinesq approximation mentioned earlier). Based on an analogy with laminar flow, it states that the Reynolds stress tensor can be related to the mean velocity gradients via a turbulent ‘eddy’ viscosity ν T : u ′ u ′ ¯ = v T ∇ u ¯ + ∇ u ¯ T so that the RANS equation becomes ∂ u ¯ ∂ t + u ¯ ⋅ ∇ u ¯ = − 1 ρ 1 ∇ p ¯ + ∇ v + v T ∇ u ¯ + ∇ u ¯ T This approach of modeling the effect of turbulence as an added viscosity is widely used for turbulent flow simulations. It is very useful as an engineering method, because the computational time is only weakly dependent on the Reynolds number. However, the validity of the Boussinesq hypothesis is limited. In contrast to ν,ν T is not a property of the fluid but rather a property of the type of flow in question. Since eddies are fundamentally different from molecules, there is no sound physical basis for equation , and DNS calculations have indeed not shown a clear correlation between u ′ u ′ ¯ and ∇ ū . The Boussinesq hypothesis is therefore inadequate in many situations, for example, for flows with sudden changes in mean strain rate (e.g. the shear layer of the wake), anisotropic flows (e.g. the atmosphere) and three‐dimensional (3D) flows. Many different methods have been suggested to calculate ν T , typically called zero‐equation (algebraic closure, such as mixing length), one‐equation and two‐equation models (see e.g. Wilcox and the references therein). The k ‐ϵ modelis an example of a two‐equation model often encountered in wind energy wake applications; the k ‐ω model [with shear stress transport (SST) limiter] is more convenient near blade surfaces. In the k ‐ϵ model, two additional partial differential equations are introduced, one for the turbulent kinetic energy k and one for the turbulent diffusion ϵ.They contain a number of constants that have been determined by applying the model to some very general flow situations (e.g. isotropic turbulence decay or flow over a flat plate). A fundamentally different approach is the Reynolds stress model (RSM), also called differential second‐moment closure model, which does not rely on the Boussinesq hypothesis. In the RSM, all the components of the Reynolds stress tensor are modeled, which makes it suitable for anisotropic flows. However, it leads to six additional partial differential equations (PDEs), making the approach expensive. Moreover, these PDEs contain terms that have to be modeled again, and often, closure relations resembling the Boussinesq hypothesis are still employed. Lastly, the disappearance of the (stabilizing) eddy viscosity term can lead to numerical problems. 2.2.2 LES In recent years, LES is receiving more attention in the wind energy wake community, because of its ability to handle unsteady, anisotropic turbulent flows dominated by large‐scale structures and turbulent mixing. This is a significant advantage over RANS methods, but the drawback is that the computational requirements for LES are much higher than for RANS. In LES, the large eddies of the flow are calculated, whereas the eddies smaller than the grid are modeled with a subgrid‐scale model. This is based on the assumption that the smallest eddies in the flow have a more or less universal character that does not depend on the flow geometry. Mathematically, this scale separation is carried out by spatially filtering the velocity field, splitting it in a resolved (also called large scale, simulated or filtered) velocity and an unresolved (small scale) part. In general, this filtering operation is defined as a convolution integral: ũ ( x , t ) = ∫ u ( ξ , t ) G ( x − ξ , Δ ) dξ where G ( x − ξ , Δ ) is the convolution kernel, depending on the filter width Δ . The subgrid velocity is then defined as the difference between the flow velocity and the filtered velocity: ũ ( x , t ) = ( x − t ) − ũ ( x , t ) This decomposition resembles that of Reynolds averaging but with the difference that in general u ˜ ˜ ≠ u ˜ and u ′ ˜ ≠ 0 . Applying the filtering operation to the Navier‐Stokes equations (and assuming certain properties of the filter) leads to the following equation: ∂ u ˜ ∂ t + u ˜ ⋅ ∇ u ˜ = − 1 ρ ∇ p ˜ + v ∇ 2 u ˜ − ∇ ⋅ u u ˜ − u ˜ u ˜ As in the RANS equation a new term appears, the subgrid‐scale (SGS) stresses. These stresses represent the effect of the small (unresolved) scales on the large scales. A widely used model to calculate these stresses is the Smagorinsky model, which again employs the Boussinesq hypothesis: τ SGS = u u ˜ − u ˜ u ˜ = − v SGS ∇ u ˜ + ∇ u ˜ T Possible ways to calculate the subgrid‐scale eddy viscosity ν SGS are to use an analogy of the mixing‐length formulation and to use one‐equation or two‐equation models involving kinetic energy and turbulent dissipation. Because of the use of the Boussinesq hypothesis, similar limitations as in RANS are encountered. A large number of other subgrid‐scale models have therefore been proposed, for example, dynamic models, regularization models and variational multiscale models (see e.g. Sagaut and Geurts ). In contrast to RANS, where the computational cost is only weakly dependent on Re, the computational cost of LES scales roughly with Re 2 . Near solid boundaries, where boundary layers are present, LES is extremely expensive because it requires refinement in three directions, whereas RANS only requires refinement in the direction normal to the wall. A possibility is to employ a hybrid approach: RANS to resolve the attached boundary layers and LES outside the wall region, so‐called detached eddy simulation (DES). Since equations and both have a similar form, this switch between RANS and LES can be made by changing the eddy viscosity based on a wall‐distance function. 2.3 Boundary conditions One of the ongoing challenges in CFD simulations of wind turbine wakes, especially when comparing with experimental data, is the prescription of inflow conditions that mimic all relevant characteristics of the atmosphere, such as the sheared velocity profile, the anisotropy of the turbulence and the instationary nature of the inflow. In early CFD simulations, uniform, laminar inflow profiles were used, but LES simulations showed later that both the presence of the shear inflow profile (the ABL) as well as the turbulence in the incoming flow have a pronounced effect on the flow field behind the rotor. For RANS simulations, Monin‐Obukhov similarity theory (see e.g. Panofsky and Dutton ) can be used to prescribe profiles for the velocity components and turbulence quantities, which are independent of time. For LES, unsteady inflow data are necessary, and basically, two different types of methods exist: synthesized inlet methods and precursor simulation methods. The synthesis technique was employed in the models of Mann, Veers and Winkelaar. They are based on the construction of spectral tensors, which model the frequency content of the wind field, and have the advantage that certain parameters like turbulent length or time scales can be directly specified. Mann's model combines the rapid distortion theory (based on a linearization of the Navier‐Stokes equations) and the von Kármán spectral tensor to generate a 3D incompressible atmospheric turbulence field, which is homogeneous, stationary, Gaussian and anisotropic but does not include the effect of the ground boundary. Troldborg used this model to create a turbulent inflow profile. This profile is introduced close to the rotor (one rotor radius upstream) to prevent the decay of turbulent fluctuations before they reach the rotor. In order to ensure a divergence‐free velocity field, the turbulent fluctuations are introduced via body forces in the momentum equation, following Mikkelsen et al. In the second category, a separate precursor simulation is made as an inflow model for a successor simulation. The advantages of this method are that the ground boundary is taken into account and that the resulting turbulent velocity field is a solution to the non‐linear Navier‐Stokes equations. However, it cannot be easily manipulated to obtain certain desired turbulence characteristics, and it is computationally more expensive. In the work of Bechmann, Meyers et al. and Stovall et al. , such a method was used. In Stovall et al. , roughness blocks were introduced to generate the turbulence. When a desired turbulence level is achieved, the blocks are removed, and the turbines are introduced in the domain. The model of Bechmann can handle the anisotropy of the atmosphere, except near the surface, where it switches to a RANS solver. As was stated earlier by Vermeer, Bechmann observed again that ‘imposing realistic and even experimentally obtained inflow conditions was found extremely difficult’. At other boundaries of the computational domain, the prescription of boundary conditions is not straightforward either. The most physical model for the presence of the ground is a body‐fitted mesh with a no‐slip boundary condition, but such a condition does not allow the prescription of a ground roughness and would require very fine grids near the surface. Therefore, a common approach in both RANS and LES models is to use wall functions to prescribe wall friction. Porté‐Agel et al. prescribed the surface shear stress and heat flux by adopting a time‐dependent variant of Monin‐Obukhov similarity theory. An alternative to no‐slip conditions are slip conditions, so that a laminar shear profile can be imposed without depending on the building‐up of a boundary layer associated with the Reynolds number of the atmosphere. A turbulent inflow is difficult to prescribe, because the turbulent viscosity of the ABL is much (10 6 times) larger than the molecular viscosity, making the effective Reynolds number of the blades much too low. At the upper boundary, symmetry conditions are often prescribed, but care should be taken when the height of the ABL is of the same order as the turbine height, something which is not unusual for modern‐size turbines, especially at night time. In streamwise direction, periodicity conditions can be used to simulate infinite wind farms or to generate a turbulence field with the precursor method discussed before, so that an ABL can form without the need for a very large domain. Periodicity conditions are also convenient from a numerical point of view because they allow the application of fast Poisson solvers like the fast Fourier transform. 3 ROTOR MODELING To solve the RANS equation or LES equation in the near and far wakes of a wind turbine, a representation of the blades is necessary. Basically, two approaches exist: the generalized actuator disk approach, in which the blades are represented by a body force ( Section 3.1 ), and the direct approach, in which the presence of the blades is taken into account by discretizing the actual blades on a computational mesh ( Section 3.2 ). 3.1 Generalized actuator disk modeling In many near and far wake calculations, the rotor is represented by an actuator disk or an actuator line. Such a representation circumvents the explicit calculation of the blade boundary layers, reducing computational cost and easing mesh generation. The actuator disk exerts a force on the flow, acting as a momentum sink. This force is explicitly added to the momentum equation : ∫ Ω ∂ u ∂ t d Ω + ∫ ∂ Ω u u ⋅ n d S = − ∫ ∂ Ω 1 ρ n d S + ∫ ∂ Ω v ∇ u ⋅ n d S + ∫ A ∩ Ω f d A which are written in weak form, because the force leads to a discontinuity in pressure. Apart from this momentum sink, one should also introduce sources of turbulence corresponding to the mechanical turbulence generated by the blades. Currently, three different approaches for prescribing the force term f exist: the actuator disk, the actuator line and the actuator surface models (see Figure ). Illustration of the actuator disk (AD), line (AL) and surface (AS) concept. 3.1.1 Actuator disk In case of a uniformly loaded actuator disk, f acts on the rotor‐disk surface A and is usually expressed in terms of the thrust coefficient C T only: ρ f = 1 2 ρ V ref 2 C T The determination of the reference velocity V ref in order to calculate C T is not obvious. For a turbine facing the undisturbed flow, V ref is evidently V ∞ , but for a turbine in the wake of an upstream turbine or in complex terrain, this is not the case. Prospathopoulos et al. proposed an iterative procedure to obtain the reference velocity and the thrust coefficient for downstream turbines modeled as actuator disks; starting with a certain V ref , one determines the thrust coefficient, from which the axial induction a follows, and then a new reference velocity based on the local flow field is computed: V ref = V local /(1 − a ). This procedure is repeated until convergence is achieved. Meyers and Meneveau and Calaf et al. used a similar approach by taking the local velocity and axial induction factor to determine the reference wind speed, but not in an iterative manner. The local velocity was obtained by disk averaging and time filtering (a LES model was used). For non‐uniformly loaded disks, the force is depending on radial position but constant over an annulus (see Figure ). Sectional lift and drag coefficients ( c l and c d ) are then used to find the local forces on the blades, like in BEM: ρ f 2 D = 1 2 ρ V rel 2 c c 1 e L + c d e D where e L and e D are unit vectors in the direction of lift and drag; c l and c d are functions of the Reynolds number and the angle of attack α. The force integral in equation then becomes a line integral; the disk is recovered as a time average of line forces. The relative velocity V rel at the radial positions is found by interpolating the velocity field in the surrounding computational cells. This is different from BEM, where V rel is found from an iterative procedure that employs a global momentum balance. Another difference with BEM is the application of a tip‐loss correction. The assumption of an infinite number of blades is corrected in BEM by locally changing the induced velocity. In Navier‐Stokes computations, this is not necessary because the flow field will notice the presence of the disk so that the induction changes automatically. However, the use of 2D airfoil data still requires a correction to obtain the right flow angle and flow velocity. A related problem is the determination of the local α to find c l and c d . Shen et al. developed a technique with which α can be determined based on information slightly upstream of the rotor. Rajagopalan et al. were one of the first to use the actuator‐type approach in a CFD code, for the calculation of vertical‐axis turbines. Time‐averaged forces were prescribed, and a finite‐difference laminar flow solver was used to solve the steady Navier‐Stokes equations. Ammara et al. and Masson et al. followed the time‐averaging approach of Rajagopalan in a control‐volume finite‐element setting. They included a second grid around the turbine in order to evaluate the surface force integral [see equation ]accurately and linearized the force term in an iterative procedure to a steady‐state solution. In Masson et al. , the lift and drag coefficients were obtained with a dynamic‐stall model. In this work, the tower was, like the rotor, also modeled as a porous surface. The forces on this surface were obtained from experimental data on the drag of a cylinder and were enforced by imposing a pressure discontinuity instead of adding the surface force in the equations. The presence of oscillations as a result of the use of collocated methods was mentioned (i.e. storing pressure and velocity variables at the same location), which was resolved by storing two different pressure values for the points located on the disk surface. Unsteady computations with the actuator disk approach were made by Sørensen et al. by using cylindrical coordinates in a rotor‐fixed reference frame. In Sørensen and Myken, a finite‐difference method was employed to solve the unsteady Euler equations in a vorticity‐stream function formulation (for advantages and disadvantages of such a formulation, see Vermeer et al. . A constant rotor loading was specified, but because of numerical difficulties at the disk edge, it was replaced by an elliptic distribution with equivalent total loading. In Sørensen and Kock, viscous terms were included, but only to stabilize the solution. The force on the non‐uniformly loaded actuator disk was obtained from tabulated airfoil data, and a non‐linear filter was applied to suppress oscillations in the vicinity of the disk. In Sørensen et al. , the vorticity‐velocity formulation was used to study different wake states, including windmill brake, turbulent wake, vortex ring and hover state. This formulation has the advantage of avoiding pressure‐velocity coupling issues, although a regularization kernel is still necessary to smoothly distribute the loading from the disk to the surrounding mesh points: f mesh = f disk * η ∈ The regularization function η ϵ is, for example, a Gaussian. The choice for ϵ, indicating the amount of smearing, is often a trade‐off between stability and accuracy. With this approach, the wiggles present in Sørensen and Kock disappeared. Madsen investigated both uniformly and non‐uniformly loaded actuator disks and the effect of turbulent mixing to show the validity of the BEM theory. It was found that BEM, with the application of a tip correction, gives a good correlation with the CFD results. In later work, Madsen et al. proposed a correction to BEM based on the comparison with actuator disk simulations, which showed that BEM overestimates the induction at the inboard part of the rotor and underestimates it at the outboard part. Apart from the uniform or non‐uniform axial loading described above, one can also introduce tangential forces on the disk surface to account for rotational effects. Meyers and Meneveau applied this in a LES context and showed that the effect of the tangential forces on the wake and extracted power appears to be negligible in the case of moderate power coefficient and high tip speed ratio. However, Porté‐Agel et al. and Wu and Porté‐Agel showed that the inclusion of rotation and non‐uniform loading leads to significant improvement in the prediction of the mean velocity and the turbulence intensity with respect to the uniformly loaded disk (see Figure ). This is especially apparent in the center of the near wake, where the uniformly loaded disk leads to an underestimation of the wake deficit and the turbulence intensity. Further downstream, the effect of rotation and non‐uniform loading disappears. Comparison of actuator disk without (dashed line) and with rotation (solid line), actuator line (solid black dots) and experiments (open red dots). Average horizontal velocity as function of height at different downstream positions. Reproduced from Porté‐Agel et al. Another approach to describe an actuator disk is the actuator shape model by Réthoré. The common cells of two different grids, one for the computational domain and one for the actuator (one dimension lower), were determined, and based on the intersecting polygons, forcing was applied to a cell. A comparison between this actuator model and a full‐rotor computation shows that modeling the wake by using forces is a good approximation for the mean flow quantities at distances larger than a rotor diameter from the wind turbine. It was observed that 10 cells per rotor diameter are sufficient; a similar number is typically found in other studies as well. However, the forces fail to represent the mechanical turbulence generated at the blade location. This turbulence can therefore be added at the disk location, independently of the actuator force, but its effect on the far wake in comparison with atmospheric and wake turbulence is small. As was mentioned, in actuator disk simulations, the boundary layers are not explicitly simulated, but their effect is taken into account via the lift and drag coefficients. Correctly simulating the blade Reynolds number is then less important, reducing the required computer resources considerably. It was shown in Sørensen et al. that the velocity field does not change noticeably when the Reynolds number is larger than 1000 (see also Snel, Ivanell et al. and Sørensen and Shen ). This corresponds roughly to results obtained by research on the nature of the interface of turbulent wakes and the outer flow, which changes around Re ∼ 10 4 . Whale et al. also found that the behavior of the wake is possibly rather insensitive to the blade Reynolds number. 3.1.2 Actuator line As an extension of the non‐uniformly loaded actuator disk approach, Sørensen and Shen introduced the actuator line approach (see Figure ). The line forces are not averaged over the disk but depend on time. Whereas in the actuator disk model vorticity is shed into the wake as a continuous sheet, in the actuator line model distinct tip vortices can be calculated. As for the non‐uniformly loaded actuator disk, the actuator line method requires knowledge of the lift and the drag on the blades. Corrections for Coriolis, centrifugal and tip effects are necessary when 2D airfoil data are used. In order to transfer the rotating line forces to the stationary mesh, a regularization kernel similar to equation is used. Mikkelsen investigated the actuator line method in detail and implemented it in EllipSys3D, a finite volume code for the solution of the incompressible Navier‐Stokes equations in pressure‐velocity formulation in general curvilinear coordinates. Howard and Pereira used point forces to represent the blades but did not take into account the distance between the force location and the cell centers, leading to high‐frequency noise in the power output. They modeled the tower as a square cylinder and found that it generates a wake that partially destroys the blade tip vortices. They recommend for future work to model the tower by point forces as well. 3.1.3 Actuator surface Shen et al. extended the actuator line method to an actuator surface method and applied it to vertical‐axis wind turbines. Whereas in the actuator line model the blade is represented by a line, in the actuator surface model it is represented by a planar surface (see Figure ). This requires more accurate airfoil data; instead of c l and c d , knowledge of the pressure and skin‐friction distribution on the airfoil surface is needed: f 2D AS ξ = f 2D F dist ξ where ξ is directed along the chord and F dist (ξ) is determined by fitting empirical functions to chordwise pressure distributions. These were obtained in Shen et al. with Xfoil, a highly accurate tool to compute pressure‐friction and skin‐friction profiles on airfoils. A comparison with 2D RANS calculations on a body‐fitted mesh around an airfoil shows that the pressure field can be accurately computed approximately one chord away from the airfoil. Furthermore, it was shown that the use of a dynamic‐stall model in vertical‐axis wind turbine calculations significantly improves the agreement with experiments. In Shen et al. , the method was applied to 3D turbine calculations and compared with the actuator line technique. Some improvements were seen in the representation of tip vortices and the flow behavior near the airfoil surface. In the actuator surface approach of Dobrev et al. , the pressure distribution was represented by a piecewise linear function that is directly calculated from the lift and drag coefficients. As in Shen et al. , the flow field was found to be more realistic compared to the actuator line approach, but the method was not yet able to model the wake of an airfoil since the shear forces on the blade were not accounted for. Leclerc, Watters and Masson used a slightly different approach with their actuator surface concept. Inviscid aerodynamic theory was used to relate vorticity distributions to both pressure and velocity discontinuities across a porous surface. This bears resemblance with vortex methods and lifting‐line theory, which describe an inviscid flow field by concentrated vortex sheets or lines. Given the lift coefficient, the circulation was computed, and subsequently, the velocity discontinuity over the surface was obtained by assuming a parabolic or cubic distribution. Viscous drag was not taken into account, but in three dimensions, the actuator surface could still extract energy from the flow because of induced drag. In two dimensions, an actuator surface could not extract energy, and a rotor disk was then modeled by prescribing the vorticity distribution of the slipstream surface, instead of the momentum loss through the disk area. The use of surface forces instead of volume forces was found to be the reason that the solution did not exhibit spurious oscillations. 3.1.4 Conclusions To conclude, a hierarchy of actuator models has emerged, from ‘simple’ actuator disks to more accurate actuator line and surface models. Gaining accuracy is accompanied by higher computational cost and the need for more detailed airfoil data: from C T (uniform actuator disk) to c l and c d (non‐uniform actuator disk, actuator line) to C p and C f (actuator surface). It would seem that the unsteady nature of actuator line and surface methods makes them most suitable for LES simulations and that the steady nature of actuator disk methods limits their application to RANS simulations. However, the computational costs limit the use of the actuator line technique to single‐wake investigations, and most LES simulations of wind farms are being performed with actuator disks. 3.2 Direct modeling The complete or direct modeling of the rotor by constructing a body‐fitted grid is physically the most sound method to compute the flow around a turbine. The focus is here on work that has contributed to the knowledge of wake aerodynamics. Compared to the generalized actuator disk approach, the blade is represented ‘exactly’, instead of a disk/line/surface approximation. However, accurately simulating the boundary layer on the blades, including possible transition, separation and stall, is computationally very expensive. Additionally, compressibility effects at blade tips can require the solution of the compressible Navier‐Stokes equations, whereas the wake remains essentially incompressible. Furthermore, the generation of a high‐quality moving mesh is not trivial. Mesh generation is therefore commonly done with so‐called ‘overset’ or ‘chimera’ grids: different overlapping grids, often structured, that communicate with each other. An example of such a grid is shown in Figure . Layout of overset grids around tower, blade and far field. Reproduced from Zahle et al. The first simulations (for wind turbine applications) with direct modeling were done by Sørensen and Hansen, employing a rotating reference frame and the k −ω SST model. The rotor power is predicted well for wind speeds below 10 m s −1 , but at higher speeds, the power is underpredicted. The strongly separated flow on the blade is not correctly captured at these speeds, which the authors attribute to insufficient mesh resolution and limitations of the turbulence model. Duque et al. used the compressible thin‐layer Navier‐Stokes equations with the Baldwin‐Lomax (algebraic) turbulence model. A comparison with the NREL phase II rotor showed good agreement of the pressure distributions, but the rotor‐tower interaction was not predicted well. In later work, the full equations were solved with the Baldwin‐Barth (one‐equation) turbulence model and compared with the NREL phase VI rotor, showing good agreement with experiments, including stalled flow and cross flow. Different transition models were investigated by Xu and Sankar with a hybrid Navier‐Stokes potential model. In this approach, the viscous compressible flow equations are solved only in a small region surrounding the rotor, whereas the rest of the flow field is modeled using an inviscid free‐wake method. The position of the transition line on the blades was shown to clearly depend on the transition model. Benjanirat and Sankar used an overset grid model to study the effect of different turbulence models. All the models predicted the out‐of‐plane forces and associated bending moments well, but difficulties were found in predicting the in‐plane forces and therefore in predicting power generation. This can probably be explained by the fact that lift prediction is less sensitive to turbulence modeling than drag prediction. Sørensen (with the k −ω SST model) and Johansen (DES) performed simulations for the NREL phase VI rotor with a rotor‐fixed reference frame. No evidence was found that the DES computations improved the RANS results in predicting blade characteristics. Madsen et al. compared direct modeling with a generalized actuator approach and concluded that the local flow angle is generally better predicted by the direct model. In computations by Johansen and Sørensen, the full 3D solution was used to extract airfoil characteristics. A method was developed to find the local angle of attack, based on earlier work of Hansen et al. The resulting c l −α and c d −α curves can be used as an input for BEM calculations and can improve actuator line or surface computations. Shen et al. devised a method to extract the angle of attack for more general flow conditions. For further improvement, BEM calculations in shear were performed, and it was found that although the power variation over a revolution was small, the variation of individual blade forces was much higher. The effect of the rotor on the velocity profile was investigated, and it was shown that disturbances in the velocity profile were still noticeable at five rotor diameters (5 D ) upstream. In later work, the same computations were repeated with an overset grid that indicated that the induction was not felt anymore more than 2 d upstream, showing that the rotor‐fixed reference frame approach is less appropriate when dealing with sheared inflow. Zahle and Sørensen investigated the influence of the tip vortices on the velocity at the rotor plane and the influence of the tower in downwind and upwind configurations, using the k −ω SST model and a fully turbulent flow. A reduction in thrust and torque of 1 − 2% because of the tower shadow was reported, as well as signs of dynamic‐stall behavior on the blades when exiting the shadow. Their results showed the breakdown of the root vortices after 3 d and the tip vortex breakdown after 4–6 d. Furthermore, it turned out that fine grids are needed to track the tip vortices and that the rotor thrust and torque are hardly affected by the far wake resolution. Investigations on transition modeling were done by Sørensen using the k −ω model and the Langtry‐Menter transition model. It was found that transitional computations lead to better agreement with experiments than fully turbulent simulations. A similar comparison with the same turbulence and transition model was made by Laursen et al. This study also showed that the application of the transition model leads to a more realistic performance of the blade: increased lift and lowered drag. Near‐wake calculations without transition model were made by Bechmann and Sørensen on the MEXICO rotor (see Section 5 ). The axial velocity in the near wake corresponds well with measurements, but in the far wake (from 2.5 d), coarsening of the mesh leads to an unphysical velocity increase. The radial velocity at 1 d downstream compares reasonably well with measurements. The blade pressure distributions are used to extract 3D airfoil characteristics (i.e. sectional c l and c d when 3D flow effects are present) that can, as mentioned before, be used for different wake models. Calculations by Wußow et al. were aimed at wake characteristics and resolved the entire wake but used a coarse resolution on the blade to limit computational efforts. An exception to the structured overset grid approach is the work of Sezer‐Uzol et al. who used a moving unstructured grid. This circumvents the interpolation problem inherent to the (structured) overset grid approach and allows for grid coarsening or refinement, but the continuous re‐meshing can be expensive. 4 WAKE MODELING By having discussed the possible ways to model the rotor in the previous section, this section discusses the different ways of modeling the wake. As mentioned before, the simulation of turbulence is of major importance in CFD of wind turbine wakes, and therefore, the work done in this area will be grouped according to the two turbulence techniques encountered in wake modeling, being RANS and LES. 4.1 RANS Although the application of the averaging procedure to the Navier‐Stokes equations, resulting in the RANS equations, and the use of the generalized actuator approach lead to a significant reduction in computational effort, solving equation is still a formidable job for wake calculations. One of the main reasons for this is the divergence‐free constraint in equation , which requires that the pressure is calculated implicitly via an elliptic equation, the pressure Poisson equation. Two simplifications to this approach, parabolization and linearization, and the original elliptic models are discussed in this section. 4.1.1 Parabolic models The first wake calculations based on RANS used axial‐symmetric and parabolic forms of the Navier‐Stokes equations. Since most parabolic models were already reviewed quite extensively in Vermeer et al. and Crespo et al. , only the main developments will be mentioned here. The parabolic model, also called boundary‐layer wake model, was obtained by neglecting both the diffusion and the pressure gradient in streamwise direction. This allowed for a fast solution by space marching, but the (pressure driven) expansion of the wake cannot be predicted properly. Ainslie used a parabolic eddy viscosity model with axial symmetry, steady flow and neglected the pressure gradients in the outer flow. Since the parabolicity assumption is not valid in the near wake, the calculations started at 2 d downstream of the rotor, where the wake is assumed to be fully developed. At that location, a Gaussian velocity profile was used as a starting condition for the far wake, and the solution was computed by advancing in the wind direction. Improvements have been made by using a variable‐length near wake and by refining the initial condition for the far wake computation by linking it to a BEM method. Reasonable results were obtained when compared to wind tunnel experiments, but because of the assumption of axial symmetry, the downward movement of the wake centerline and the upwind movement of the location of maximum turbulence intensity cannot be predicted. A parabolic model that does not assume axial symmetry is UPMWAKE, by Crespo et al. The k −ϵ turbulence model was employed, with the constants adapted for neutral ABL flow. The presence of the turbine was taken into account by imposing perturbation values of velocity, temperature, kinetic energy and dissipation at the turbine location. The UPMWAKE model was the foundation for the WAKEFARM model of the Energy research Centre of the Netherlands (ECN). This model was subsequently improved by including the axial pressure gradient via a boundary‐layer analogy and coupling it to the near wake via a vortex‐wake method, so that the expansion of the wake is taken into account. A version of UPMWAKE that accounts for the anisotropic nature of the turbulence in the atmosphere and the wake, UPMANIWAKE, was made by Gómez‐Elvira et al. An explicit algebraic model for the components of the turbulent stress tensor was used. The near wake was assumed to be fully turbulent, and the analysis was carried out for a neutral atmosphere. Improved results with respect to isotropic models were obtained, and analytical expressions for the added turbulence intensity in the wake as function of thrust coefficient were derived. They suggested that future improvements can be made by using LES. The UPMWAKE model was extended to UPMPARK to calculate parks with many turbines. In a parabolic setting, this is a sequence of single‐wake computations, where the outflow of one wake forms the inflow condition of the next. 4.1.2 Elliptic models Crespo and Hernández extended the parabolic code to a fully elliptic version to resolve discrepancies with experimental results that are mainly found in the near‐wake region, but no significant improvements were observed in single‐wake computations. Rajagopalan et al. used an elliptic approach to study interaction of vertical‐axis wind turbines in a farm with a laminar solver, from which it is concluded that the mutual influence of turbines is significant. Ammara et al. used an algebraic closure (Prandtl' mixing length) with values from Panofsky and Dutton. In single‐wake computations, a good prediction of the velocity deficit is found. In farm simulations, the existence of a ‘Venturi’ effect is demonstrated, positively influencing the power output of a farm. Different researchers found that the standard k −ϵ and k −ω models result in too diffusive wakes: the velocity deficits are too small, and the turbulence intensity does not show the distinct peaks observed in experiments. The reasons for the failure of these standard turbulence models in wind turbine wakes were explained by Réthoré and are basically caused by the limited validity of the Boussinesq hypothesis, as mentioned before. To improve correspondence with experimental data, several adaptations of the k −ϵ model have been suggested. These adaptations all strive to reduce (directly or indirectly) the eddy viscosity and hence the diffusion, in the near wake. Firstly, El Kasmi and Masson added an extra term to the transport equation for the turbulent energy dissipation (ϵ) in a region around the rotor, based on the work of Chen and Kim. This method leads to the introduction of two additional parameters: a model constant and the size of the region where the model is applied. Compared with experimental data, significant improvements over the original k −ϵ model were observed. Cabez‐n et al. and Rados et al. have also applied this approach and found acceptable agreement with other experimental data. Réthoré noticed however that increasing the dissipation proportionally to the production of turbulence contradicts the LES results. A second approach is the realizability model. In this model, the eddy viscosity is reduced to enforce that Reynolds stresses respect the so‐called realizability conditions (see e.g. Schumann ). According to Réthoré, these conditions are not satisfied in the near wake by the eddy viscosity‐based k −ϵ model, because of the large strain rate at the edge of the wake, where the Boussinesq hypothesis is inadequate. Both Réthoré and Cabezón used a realizable model based on the work of Shih. Compared with the standard k −ϵ model, better results are obtained with the realizable model, although the prediction of both wake deficit and wake spreading and turbulence intensity remains unsatisfactory. More adaptations are under investigation. Prospathopoulos et al. adapted the constants of the k −ϵ and k −ω model for atmospheric flows and investigated the effect of complex terrain. Rados et al. used Freedman's model to change the constants of the ϵ equation in case of stable atmospheric conditions. Prospathopoulos et al. obtained good results with Durbin's correction, in which the turbulent length scale is bounded. Another possibility comes from an analogy with forest and urban canopy simulations, where obstacles are also modeled as body forces. Réthoré adapted the k and ϵ equations to take into account the extraction of turbulent kinetic energy by the actuator and found significant improvement with respect to the standard k −ϵ model. The RSM approach, not relying on the Boussinesq hypothesis, was employed by Cabezón et al. and gave more accurate results than both the El Kasmi‐Masson model and the realizability model. However, in predicting the velocity deficit in the near wake, the parabolic UPMPARK code outperformed the elliptic models because in the parabolic code, the streamwise diffusive effects are totally neglected. In the far wake, the velocity deficit is predicted similarly to the elliptic models, but the turbulence intensity is overestimated and does not show the distinct peaks. In case of stratification, the production of turbulence because of buoyancy has to be taken into account in turbulence closure models. Most researchers use the constants of the k −ϵ model proposed by Crespo et al. for neutral atmospheric flows. For non‐neutral stratification, Alinot and Masson changed the closure coefficient of the k −ϵ model related to buoyancy production to match the basic flow described by Monin–Obukhov similarity theory; no wind turbines are present. With respect to Crespo and Hernández, this gives improved results for stable atmospheric conditions, especially when considering profiles of turbulent kinetic energy and dissipation. Sanz et al. found that CFD modeling with Monin–Obukhov theory is not valid for stable conditions in offshore situations. 4.1.3 Linearized models Ott et al. used the k −ω model employing the canopy approach mentioned before, in combination with the fast linearized model of Corbett et al. The linearization of the Navier‐Stokes equations, combined with a spectral method and lookup tables (containing a.o. velocity profiles), allows for a very fast solution of entire wind farms. In such a spectral method, the governing equations are transformed to Fourier space, which has the advantage that no computational grid and therefore no artificial diffusion is present, although restrictions on the boundary conditions apply (typically periodicity in lateral direction). The results for the velocity field are most accurate in the far wake, where the assumption of small perturbations (linearization) is justified. However, compared to a full‐model prediction, the turbulent kinetic energy shows serious mismatches in both the near and far wakes. 4.2 LES As mentioned before, LES has the advantage over most RANS models in that it is better able to predict the unsteady, anisotropic turbulent atmosphere. Jimenez et al. used a dynamic subgrid‐scale model with the rotor represented by a uniformly loaded actuator disk. A good comparison with experiments was found. The same model was used in a paper by Jimenez et al. to study the spectral coherence in the wake and in another paper by Jimenez et al. to study the wake deflection because of yaw; in both cases, reasonable agreement with experiments and an analytical model was observed. Wußow et al. used a Smagorinsky‐Lilly subgrid‐scale model with a wind field generated with a von Kármán spectrum to study loading of wind turbines in a farm. It was shown that when calculating the mean velocity field, the averaging time must be long enough to obtain a good comparison with 10 min averaged experimental data. Better meshing strategies led to improvements of velocity deficit and turbulence intensity. Troldborg et al. used a vorticity‐based mixed‐scale subgrid‐scale model with the actuator line technique to study the influence of shear and inflow turbulence on wake behavior. The ABL was imposed with the force‐field technique of Mikkelsen, and atmospheric turbulence was generated with Mann's method (see Section 2.3 ). It was found that shear and ground effects cause an asymmetric development of the wake with larger expansions upward and sideward than downward. Inflow turbulence, when compared with laminar inflow, destabilized the vortex system in the wake, resulting in an improved recovery of the velocity deficit but also an increased turbulence intensity. Furthermore, it was found that the wake is also unstable and becomes turbulent for uniform, laminar inflow, especially at low tip speed ratios. In the work of Ivanell, the same LES model was used for farm simulations. Twenty non‐uniformly loaded actuator disks, in combination with periodic boundary conditions, were used to simulate the 80 turbines that comprise the Horns Rev wind farm (Figure (a)). The power production of downstream turbines agreed reasonably well with experimental data (Figure (b)). LES was also used to study the stability of tip vortices. By inserting harmonic perturbations at blade tips, an exponential relationship between the perturbation frequency (related to turbulence intensity) and the length of the stable part of the wake was observed. LES of the flow through 20 actuator disks, representing the Horns Rev wind farm. Reproduced from Ivanell. Meyers and Meneveau performed LES with a (non‐dynamic) Smagorinsky model to study infinite arrays of wind turbines in staggered and non‐staggered arrangement, and it was observed that the staggered arrangement had a higher power production. Porté‐Agel et al. and Wu and Porté‐Agel et al. employed a Lagrangian scale‐dependent SGS model, a parameter‐free model that performs better in ABL simulations than the traditional Smagorinsky and standard dynamic models. Similar to what Réthoré showed for C μ in the k −ϵ model, Wu and Porté‐Agel showed that the Smagorinsky coefficient C S is not constant in the wake but increases in the center of the wake and decreases near the ground and in the shear layer. With an actuator‐type approximation for the turbine, this model provides a very good agreement with atmospheric wind tunnel data. Furthermore, the effect of a wind turbine on a stable ABL was investigated, demonstrating that the momentum and buoyancy flux at the surface are reduced and as a result possibly influence the local meteorology. A comparison between RANS (standard k −ϵ) and LES (method of Bechmann ) was made by Réthoré in order to explain the aforementioned discrepancies between the k −ϵ model and the experimental data. The LES results are clearly superior to RANS when comparing both mean velocity profiles and turbulence quantities, but the computational time increases from hours (RANS) to days (LES) for a single‐wake case. In the study of Stovall et al. , the standard k −ϵ model and LES with a one‐equation SGS model were compared. Again, it was observed that the LES results are closer to experimental data in a single‐wake calculation and that the wake recovery because of turbulent mixing with the outer flow is much better captured by LES in a multiple‐wake situation. The difference in computational time was a factor 60, when the same grid for both RANS and LES was used. However, a well‐resolved LES should require a much finer grid than RANS, making it more expensive, as will be explained in the following. An estimate of the necessary grid resolution can be obtained based on the Reynolds number. For a sufficiently resolved LES, one needs cell sizes far enough in the inertial subrange, i.e. of the order of the Taylor microscale. The fact that the Taylor microscale scales with Re −1/2 and the Reynolds number based on the diameter of a modern‐size wind turbine is O 10 8 leads to the requirement of cell sizes around 1 cm 3 . This is somewhat conservative; for atmospheric flows (based on an integral length scale of 1 km), the Taylor microscale was estimated at 10 cm. Most of the LES computations to date use cell sizes of approximately 1–10 m, which might be too coarse to resolve enough scales. Computations require typically 10 7 –10 8 grid points and run on supercomputers for several days or weeks, even for single‐wake situations. With the continuous advance in computer power, finer meshes will be possible, but the associated increase in data analysis will leave this approach unattractive for many engineering purposes. 4.3 Numerical issues Most RANS codes use second‐order accurate finite volume schemes on structured meshes, with upwind discretization of the convective terms and central discretization of the diffusive terms, leading to stable and robust schemes. Implicit methods are normally used to find a steady‐state solution. In LES, the temporal and spatial discretization of the convective terms should be done more carefully. The use of upwind schemes for the spatial discretization can influence the energy cascade from large to small scales because of the introduction of numerical dissipation. Several authors have indeed mentioned premature turbulence decay and related it to possible numerical dissipation (see e.g. Bechmann, Stovall et al. and Troldborg et al. . This forms the rationale for the application of central or spectral schemes (see e.g. Bechmann and Meyers et al. ). Apart from the spatial discretization, the temporal discretization can also introduce numerical dissipation. In Bechmann, this limits the time step, partly removing the advantages of an implicit method. High‐order low‐dissipative explicit schemes, such as the standard four‐stage Runge‐Kutta method, used in Calaf et al. can then be an attractive alternative. Lately, high‐order (typically fourth order) schemes have also been employed for the spatial discretization, reducing interference between the subgrid‐scale model and the discretization. A blend with third‐order upwind schemes is sometimes made to ensure numerical stability. However, it should be realized that the formal order of accuracy of a discretization scheme is only obtained in the limit of very high spatial resolution, something typically not encountered in LES of wind turbine wakes. Mesh refinement studies are often computationally too expensive, and one has to rely on the energy spectrum to check if the inertial subrange is captured well. It should be noted that partial cancelation of subgrid modeling and numerical errors may occur, leading to the counterintuitive conclusion that high‐order accurate schemes, improved subgrid models or finer meshes, may lead to worse results. 5 VERIFICATION AND VALIDATION The process of testing a CFD model should consist of two steps: verification (‘solving the equations right’) and validation (‘solving the right equations’). The first step, verification, is necessary to evaluate the errors that are made by approximating the continuous equations by a discrete set. This gives an indication of the mesh size necessary to reach a required accuracy. Highly accurate numerical simulations (benchmark tests) or analytical solutions may serve as verification data. An example of such tests is the actuator disk simulation performed by Réthoré: the flow through lightly loaded (2D and 3D ) and heavily loaded actuator disks. The second step, validation, requires a comparison with experimental data to investigate the validity of the model and to assess modeling errors. Generally, the turbulence model, the blade approximation and the inflow conditions are the main sources. Wind tunnel tests on rotating wind turbine models yield very suitable experimental data because they take place in a constant and controllable environment. However, the wind tunnel conditions are typically different from atmospheric conditions (scaling and blocking effects, turbulence intensity). One can therefore only conclude how well the CFD model simulates the wind tunnel conditions, not how well it simulates the ‘real’ flow through turbines standing in the field. An elaborate overview of the near and far wake wind tunnel experiments is given in Vermeer. Since that review, three important wind tunnel tests have been performed. First is the tests in the German‐Dutch wind tunnel DNW in 2006, on a 4.5 m diameter rotor, the so‐called MEXICO project. Near‐wake measurements with particle image velocimetry (from ‐1 d to 1.5 d) and pressure recordings on the blades were made, but no far wake data are available. The measurements are evaluated and compared to simulations in the (currently ongoing) MexNext project, IEA Annex 29 [ http://www.mexnext.org ]. Second, the experiments from Chamorra and Porté‐Agel with hot‐wire measurements in an atmospheric wind tunnel provide high‐resolution spatial and temporal measurements that can be used to validate LES models. Thirdly, there are the experiments on wake meandering in the ABL wind tunnel at the University of Orléans, with wind turbines modeled as porous disks. Besides wind tunnel data, full‐scale field tests are a valuable source of experimental data; they have been cataloged in a database by IEA, Annex XVIII. They are the only way to obtain data on full‐scale wake development and wake‐wake interaction. A common problem of field measurements is that the inflow is often not fully known because the number of measurement stations (giving point data) is limited; recent advances in the LIDAR technique can significantly increase the knowledge of the entire flow field. Well‐known single‐wake experiments are the Tjæreborg turbine, the Nibe turbine and the Sexbierum measurements. Well‐known sites for multiple‐wake experiments are Horns Rev, Nysted, Middlegrunden, Vindeby and the ECN Wind Turbine Test Site Wieringermeer. Valuable comparative studies between different wake models and experiments at these sites have been obtained in the ENDOW project and, more recently, the UPWIND project. Some results of the latter project are shown in Figure , where experimental data and numerical simulations of the Horns Rev wind farm are compared. The numerical simulations have been made with engineering models (WAsP from Risø, WindFarmer from Garrad Hassan) and RANS models, namely a parabolic code (ECN's WAKEFARM) and an elliptic code (NTUA). LES of the same wind farm were already shown in Figure . First, it should be noted that the standard deviation in experimental data is relatively large, as indicated by the error bars. Taken this into account, the agreement between models and experiment is good, typically within 15% of the experiments, and it is remarkable that the engineering models are at least as accurate as the RANS and LES models, while having much lower computational costs. Comparison of models and experiments in predicting the power output of Horns Rev for different wind directions. At 270°, the wind is aligned with the turbines. Reproduced from Barthelmie et al. In the recent Bolund experiment, linearized models, RANS models and LES models were compared in predicting the flow over a hill, without the presence of wind turbines. Linearized models were found to be unfavorable because of the steepness of the hill. RANS models performed best when compared with experimental data, with similar trends between a large range of different codes and turbulence models. A much larger spread was observed in the results obtained by LES models, and it was concluded that LES has not matured enough for widespread application to atmospheric flow simulation. 6 CONCLUSIONS A growing number of researchers are using CFD to study wind turbine wake aerodynamics. CFD has become a mature tool for predicting a wide range of flows; however, one important ongoing challenge is the accurate representation of turbulence. Since in wake aerodynamics turbulence is a dominating factor, affecting both blade performance and wake behavior, most efforts are now directed at improving turbulence modeling. In the near wake, full 3D Navier‐Stokes solvers are used with overset grid approaches to resolve the boundary layers on the rotating blades. Such computations are able to predict separation and transition to reasonable accuracy, although the investigation of the effect of different turbulence models is an ongoing work. This direct approach is not yet suitable for optimizing turbine blades or optimizing farm layouts. The current contribution of the direct modeling approach is to increase the understanding of the flow physics on the blades and to improve simpler models such as the BEM method. To reduce computational requirements for wake simulations, the presence of the blades is in general taken into account with the generalized actuator disk approach, in which the rotor is represented by forces. A careful numerical treatment is necessary to correctly handle these forces. Three methods are currently in use, being the actuator disk, actuator line and actuator surface methods. The more expensive actuator line and surface techniques use more detailed aerodynamic blade characteristics than the actuator disk, and as a result, they are more accurate. The actuator disk and line techniques are maturing, and current research is on the actuator surface technique, but because of the lower computational effort, the actuator disk remains the most widely used method for multiple‐wake simulations. The use of the actuator approach seems justified because in the far wake, differences with the direct modeling approach are shown to be small. However, a detailed study of the effect of the blade geometry on the wake and a comparison with actuator line simulations are still lacking. For modeling turbulence in the wake, RANS will most likely prevail to be the engineer's choice, even though many eddy viscosity‐based models like k −ϵ proved to be too diffusive. Ongoing research is now in different adaptations and corrections to reduce diffusive effects in the near wake. Increased interest is in the application of LES, a method that naturally takes into account the unsteadiness and the anisotropy of the turbulence in the atmosphere, which is not properly done in most RANS methods. In single‐wake computations, LES shows indeed an improved agreement with experimental data compared with RANS, but the increase in the required computational effort, which is about two orders of magnitude, still prohibits its widespread use for farm applications, so that its use is currently confined to gain physical insight and improve simpler methods. Because of the high computational costs associated with fully resolved LES, mainly under‐resolved simulations are available nowadays. Thorough mesh refinement and subgrid‐scale model studies are so far still open problems. The use of accurate wind tunnel data in subgrid‐scale model studies is indispensable; wake measurements on large turbines in atmospheric tunnels are a necessity. An impression of the magnitude of subgrid modeling errors and discretization errors is important since it gives the opportunity to (partially) quantify the uncertainty in CFD results. Indeed, in a time where CFD codes are becoming readily available for wind turbine wake computations and increased computer power enables larger problems and more refined modeling, a possible next step is the quantification of uncertainties in the computations. These uncertainties originate not only from discretization and turbulence modeling errors but also from the description of the inflow, terrain geometry, rotor geometry, etc. A quantification of uncertainties would make the comparison with experimental data more fair and will give a guideline in which areas the CFD of wind turbine wakes has to be improved. In this respect, the reduction of the uncertainty in the prescription of atmospheric inflow conditions remains a great challenge. More research on the effect of stratification on power production is to be expected. A coupling with mesoscale atmospheric models can shed more light on appropriate boundary conditions and, at the same time, can be used to investigate the effect of wind farms on the local meteorology. ACKNOWLEDGEMENTS Colleagues A.J. Brand, P.J. Eecen, J.G. Schepers and H. Snel from ECN are acknowledged for their help in supplying literature and correcting the manuscript.

Journal

Wind EnergyWiley

Published: Oct 1, 2011

There are no references for this article.