On the spread and decay of wind turbine wakes in ambient turbulence

The decay of the downstream wake of a wind turbine plays an important role in the performance of wind farms. The spread and decay of a wake depend both on wake meandering (advection of the wake as a whole) and wake diffusion (widening of the wake within its meandering frame of reference). Both of these effects depend strongly on the intensity of the ambient turbulence relative to the velocity deficit in the wake, and on the integral length scale of the turbulence relative to the wake width. Recent theory, which we review here, shows how intense large-scale turbulence can lead to a rapid x−2 decay in the time-averaged centreline velocity deficit, as compared to a x−1 decay for smaller scale turbulence, where x is distance downstream. We emphasise in this paper that common wind farm models do not predict this rapid decay. We present new experimental measurements of the velocity deficit downstream of a porous disc in relatively large-scale ambient turbulence which corroborate predictions of a x−2 decay, and we show theoretically that the commonly used k-ϵ model does not capture this effect. We further show that a commercial CFD package, configured to match our experiments and employing the k-ϵ model, fails to predict such rapid decay. We conclude that steady simulations of wind turbine wake dynamics are insufficient for informing wind farm layout optimisation.


Introduction
Wind farm developers must trade off the cost of spacing turbines further apart, with the associated benefits of increased energy capture due to increased wake recovery. A sound understanding of the spread and decay of a wake downstream of a turbine is therefore critical. The rate at which the far wake decays depends on the ambient turbulence, but most current wind farm models do not predict this interaction correctly. The 'far' wake begins a few diameters downstream of the turbine and we define it here as the region after which the details of the initial flow disturbance have become insignificant, i.e. tip vortices have disintegrated and the flow profile has smoothened out.
Broadly speaking there are two methods of analysing the aerodynamics of wind farms: 'theoretical' models, also known as explicit or kinematic models; and 'numerical' models, also known as implicit or field models [1]. In theoretical models, the decay of the time-averaged centreline velocity deficit, u m , with distance downstream, x, is often modelled as where d is the rotor diameter, U ∞ the upstream velocity and C and n are dimensionless constants [2]. Theoretical models use various superposition theories for multiple wakes and are convenient for engineers because they are simple and quick to use (for reviews see [1,2,3,4]). For high Reynolds number (turbulent) axisymmetric wakes in an irrotational external flow the exponent is n = 2/3 and in a laminar wake or a wake that grows by diffusion due to small-scale ambient turbulence n = 1 (observations find 0.75 < n < 1.25) [2]. However, in large-scale turbulence the exponent is n = 2 [5] corresponding to much more rapid mixing; we refer to this as 'ballistic' mixing and in this paper we discuss when it is expected and its relevance to wind farms. By contrast numerical models, which are more physically comprehensive but computationally demanding, solve equations describing the fluid motion throughout the region of a wind farm. The wide range of geometrical length scales of the wind farm and turbine blades, and range of scales within the flow, means that some reduced modelling approach is required to represent the turbines -for example actuator discs -and the fluid dynamics.
Turbulence models separate the mean flow from the fluctuating flow and use an eddy viscosity model to account for their dynamic interaction. A large number of numerical simulations of wind farm aerodynamics have used the k-turbulence model in some form (see [1,2,6] for reviews) as part of a Reynolds Averaged Navier-Stokes (RANS) model. In more recent years some large eddy simulations (LES) have been used, in which the larger scales are computed explicitly [7,8,9,10]. LES models capture more of the physics of the flow but are computationally expensive.
To predict the power output of and loads on a turbine in a wind farm the probability density function (PDF) of the wind speed at each turbine is required, together with further characteristics of the turbulence. The PDF is necessary because power output and loading do not depend linearly on velocity, hence it is incorrect to compute a time-averaged velocity field and use this to calculate a time-averaged power production.
The aim of this paper is to highlight a limitation of common wind farm models, namely that they do not capture the fundamentally important physical effect of ballistic wake decay. In §2 we review theoretical models of wake spreading in turbulent flow and we show theoretically that a k-turbulence closure model is incapable of predicting ballistic spreading -a weakness not previously highlighted in reviews such as [6]. In §3 we present experimental measurements from a water flume, with large-scale ambient turbulence, which show rapid decay of the wake of a porous disc. In §3 we also present a simulation of this experiment using a commercial software package that employs the k-closure model, showing that this model does not capture the rapid decay. In §4 we give a general disussion and conclusion.

Theory
To calculate the velocity deficit within a wind farm we consider each wind turbine to be a momentum sink; a schematic of the problem considered is shown in Figure 1. The momentum sink is characterised by a mean drag force,F D , which is related to the momentum deficit flux in the mean flow byF whereC D andū are the time-averaged drag coefficient and streamwise flow speed respectively, A is the area swept by the rotor and ρ is the fluid density. The coordinate system x = (x, y, z) is Cartesian, with x being the streamwise component. This is the classic formula proposed by Betz [11] and holds only sufficiently far downstream that the streamwise pressure gradient is negligible and u m U ∞ , where u m (x) = U ∞ − u(x, 0, 0). Regardless of the mechanism of wake spreading, this integral constraint applies.

Wake dynamics in the presence of ambient turbulence
The turbulence within a wind farm is composed of the pre-existing turbulence in the atmospheric boundary layer and additional turbulence caused by the presence of the turbines. Everywhere the turbulence is strong and significantly influences the evolution of the wakes. Small-scale turbulence acts to diffuse (widen) the wake in its moving frame of reference, whereas very large-scale ambient turbulence acts to advect the wake as a whole -a phenomenon known as wake meandering [12,13]. This is illustrated in Figure 1, where we take wake meandering to mean transverse movement of the wake centreline and wake diffusion to mean widening of the wake relative to its centreline. A useful analogy is that to single-and twoparticle passive dispersion; the mean displacement of a single particle with respect to a fixed reference point (cf. R c ) is influenced primarily by the large energy-containing length scales of the turbulence, whereas the increase in the mean separation of two particles (cf. R ∆ ) is mainly influenced by the scales of the order of the separation (Batchelor [14,15]).
There is no clear-cut borderline between the length scales that contribute to meandering and those that contribute to diffusion. The turbulence is comprised of a continuous distribution of scales and only those in the very large or very small limits relative to the wake width will contribute solely to meandering or diffusion respectively. Nevertheless, evidence is emerging [5,16,17,18] to suggest that, whereas small-scale deformation and diffusion depend on dynamic interactions within the wake, overall advection (meandering) of flow structures due to large-scale turbulence is a process of nearly passive transport -i.e. it can be modelled by a consideration of only the kinematics of the large-scale ambient turbulence.
The phenomenon of wake meandering has been known for a long time [19,12], and was more recently treated by Larsen et al. [13] who were the first to propose a model of meandering based on the conjecture that the advection of the wake deficit by large-scale air movements is passive.

Model for the evolution of weak wakes
We recently presented a theoretical model to explain how an axisymmetric [5] and a planar [16] wake -behind a cylinder and a sphere respectively -grow in the presence of strong ambient turbulence. We take 'strong turbulence' to mean I t > u m /U ∞ , where I t = u * /U ∞ is the turbulence intensity and u * = v * = w * are the RMS of the velocity components. In the limiting cases of the ambient turbulence having a very large or very small integral length scale relative to the wake width, we determined analytical expressions for the time-averaged wake growth and decay of the velocity deficit, and we summarise those results here. Since the analysis is focused on establishing time-averaged quantities, the relevant wake width is R w , which is defined later in this section; R c and R ∆ do not enter explicitly. Consider a momentum sink introduced into a fluid flow with velocity field u = (U ∞ , 0, 0) + (u , v , w ), where the two parts are the uniform mean flow and a homogeneous and isotropic turbulent velocity field respectively. The integral length scale is where an overline denotes the time average. If the ambient turbulence is strong, the wake will effectively be advected passively, which can be understood by treating the wake as a continuous release of passive fluid particles possessing velocity deficit. Following Taylor [20], the equation governing the spreading is where Y and Z are the Lagrangian coordinates of the 'wake particles'. We consider the limits of small times and large times relative to the Lagrangian integral timescale, where analytical solutions to (4) can be established. Using Taylor's hypothesis [21], the time-dependence can be recast into dependence on distance downstream, letting dx = U ∞ dt. These assumptions yield, after substituting R 2 w = Y 2 + Z 2 , Combining (5 -6) with (2) and assuming a self-similar Gaussian velocity profile, the timeaveraged centreline velocity deficit decays as The existence of the first (ballistic) regime has been confirmed by experiments and numerical simulations [22,23,24,25]. The hypothesis that strong, large-scale turbulence is capable of passively advecting flow structures in general is also supported by other studies -for example [18,26,27] on wind turbine wakes and [17] on wing-tip vortices. In order to predict the PDF of the velocity deficit within the wake we require both information on wake meandering and wake diffusion; the knowledge of how R w and the time-averaged velocity deficit evolve downstream is not enough. However, the limit of R w L is effectively that for which the fluctuations in R c are negligible and the growth is dominated by the increase in R ∆ . In the limit of R w L the downstream increase in R w is primarily due to fluctuations in R c , although R ∆ may also increase significantly with downstream distance. An analysis of a more general case is beyond the scope of the present paper -the authors intend to address this in future work.

The k-model
The k-model is one of the most widely used eddy-viscosity models [6]; it is largely based on dimensional analysis and is used to approximate the turbulent viscosity, ν t , from the turbulent kinetic energy, k, and the mean dissipation rate, , by the relationship ν t = C µ k 2 / . The transport equations for k and are ∂k/∂t +ū · ∇k = ∇ · ((ν + ν t /σ k )∇k) + G − ; where G = 2ν t S ij S ij , S ij = 1 2 ∇u + (∇u) T being the strain rate tensor. Typical values for the coefficients are C µ = 0.09, C 1 = 1.44, C 2 = 1.92, σ k = 1 and σ = 1.3 [28].
In the absence of a rigid body, the turbulent kinetic energy decays due to dissipation, where d /dk = C 2 /k, Substituting into (9), which reduces to U ∞ dk/dx = − , and integrating gives (11) The decay of the RMS velocity with x (from 11) is consistent with experimental observations for grid turbulence. For weak wakes, u m u * , the contribution of wake generated turbulence to ν t is small, in which case the above relationships hold approximately. The momentum equation then reduces to where r = y 2 + z 2 is the radial (transverse) coordinate and the turbulent viscosity ν t = ν t 0 (k/k 0 ) 2−C 2 depends weakly on x (since C 2 ≈ 2) and is much larger than ν. The growth of the wake width is where λ > 0 is a constant. Assuming a Gaussian profile, the velocity deficit calculated from (2) is which yields a centreline deficit decay of Ultimately, u m ∼ x −1 , which is dramatically different to (7); the k-model is not able to account for the rapid time-averaged wake spreading in the presence of intense large-scale ambient turbulence that results from meandering. This is of course to be expected since it treats the influence of the external turbulence on the wake as extra diffusion.  Table 1. Summary of the grid independence study for the numerical simulation. The table compares the calculated centreline velocity deficit with calculations using 4 × 10 6 elements.

Comparison with experimental and numerical results
Here we compare theoretical predictions and arguments with results from a laboratory study and a numerical study performed using commercial CFD software that implements the k-model.

Experimental study
The experiments were conducted in a recirculating water channel, 1.2 m wide and 0.7 m deep. Ambient turbulence was generated by the pumps and swept along the channel; the turbulence intensity was I t ≈ 0.16. A porous disc of 0.2 m diameter was held in the flow at centre of the channel by a thin strut. To match the drag coefficient to that of a turbine (C D = 8/9), and based on [29], a disc of 58% porosity was used (4mm circular holes spaced by 5mm in a hexagonal pattern). The mean freestream velocity was U ∞ = 0.17 m s −1 , hence the Reynolds number based on the diameter of the disc was 3.4×10 4 . Point measurements of the three velocity components were taken at a number of positions downstream of the disc using an Acoustic Doppler Velocimeter (ADV). Samples at 25 Hz were taken for 360 s, ensuring statistical convergence. The measured velocity was an average of data sampled in a cylindrical volume, 3.1 mm long and 6 mm in diameter, with the cylinder axis parallel to the mean flow. The upstream integral length scale was L 0 ≈ 0.2 m. Details of apparatus and diagnostic tools are in [16].

Numerical study using the k-model
The commercial software package CFX was used to simulate a flow configuration similar to that of the experiment using an unstructured grid. The model domain was 2.5 m long, and extended 0.5 m upstream of the porous disc, where the inlet profile was defined using experimental data; a power-law shear profile was fitted through the measured points, using a 1/7th exponent. The k-model was used to calculate the turbulent viscosity as described in §2. The porous disc was represented by a thin pressure loss, defined as ∆p = 0.5ρU 2 d κ, where U d is the fluid speed at the disc and κ = 2. The computational model requires both k and to be specified on the inlet surface, based on the turbulent intensity, I t , and integral lengthscale, L 0 , defined by Half of the laboratory tank was simulated, cut in a vertical plane parallel to the mean flow down the centre and symmetry was imposed on that plane. The free surface was modelled using a free slip wall condition; for sub-critical Froude numbers this is a good approximation [30] and it is computationally efficient. A qualitative assessment during the experiments confirmed negligible surface elevation changes. The grid density was increased around the actuator disc (i.e. where gradients were highest) and Table 1 summarises a grid independence study. Figure 2 shows the velocity deficit for the experimental measurements and numerical calculations. The experimental measurements clearly show a high rate of decay with an exponent  Figure 2. Decay of the mean velocity deficit with distance downstream from an actuator disc in an environment with intense large-scale ambient turbulence.

Results
close to 2 in the far wake, while the numerical calculations show an eventual exponent closer to 1. For the theoretical prediction of (7), the initial wake width is approximated using the experimental measurements. The experiments show a slightly slower decay than that predicted by the theory, which is to be expected since the integral length scale of the ambient turbulence was comparable to the disc diameter. The results clearly show that the k-model does not account for the effects of large, intense ambient turbulence.

Discussion and conclusion
We have addressed the influence of ambient turbulence on the spreading and decay of wind turbine wakes. When the length scale of the turbulence is large the wake will meander, which affects the rate of spreading of the time-averaged wake and subsequently the rate of decay of the time-average velocity deficit. The decay rate changes from diffusive (x −1 ) to ballistic (x −2 ), the latter of which has previously been shown experimentally [23,16], numerically [22,24,25], and with theoretical arguments [5,16] and can exist within significant regions of a wind turbine's wake when R w < L. Theoretical models of wind farms do not account for ballistic spreading. We also note that the calculation of power production also requires information about the unsteadyness of the flow (e.g. velocity PDF) in addition to the time averaged velocity field, since the relationship between power and velocity is not linear. When wake meandering is significant this becomes an important consideration but is not a part of current theoretical models.
Time-averaged numerical models are also incapable of modelling ballistic spreading. We have shown theoretically that the common k-formulation can only predict a diffusive spreading. We also tested this theory by performing experiments with large, intense turbulence and a porous disc, showing rapid wake decay, and simulated the experiment with a commercial CFD code that uses the k-model. This test confirmed that the k-model does not capture the rapid decay of the velocity deficit. More comprehensive numerical calculations using large eddy simulations (LES) should in principle capture the more rapid wake decay associated with large-scale turbulence. This was shown for the wake of a sphere by [24]. Studies of the wake of wind turbines or actuator discs have also been made using LES methods [7,8,9,10] showing significant wake meandering caused by large-scale turbulence. This implies ballistic decay rates but we are not aware of an LES study that has specifically identified this; such a study would be valuable.
A number of authors have made experimental measurements in the wakes of wind turbines, for example [18,26,27], showing that significant wake meandering is common. Estimates of the rate of decay of the average velocity deficit using this data would be of much interest. Ideally a theoretical model would include ballistic and diffusive spreading within a broader model for wind farms such as reviewed by [2,4]. This will require a new understanding of the interaction of wake meandering and diffusion, and wake interaction with the atmospheric boundary layer, and this should be the focus of future efforts applied to theoretical modeling of wind turbine wake decay.