An entrainment model for fully-developed wind farms: effects of atmospheric stability and an ideal limit for wind farm performance

TL;DR

This paper develops a model for large wind farm performance considering atmospheric stability, providing a new upper performance limit and validating predictions with field data and simulations.

Contribution

It introduces a novel entrainment-based model for fully-developed wind farms that accounts for atmospheric stability effects and establishes an ideal performance upper bound.

Findings

Model predictions agree well with field data and simulations.

Atmospheric stability significantly influences wind farm power output.

An upper performance limit is an order of magnitude higher than current turbines.

Abstract

While a theoretical limit has long been established for the performance of a single turbine, no corresponding upper bound exists for the power output from a large wind farm, making it difficult to evaluate the available potential for further performance gains. Here we build a model describing the essential features of a large array of turbines with arbitrary design and layout, by considering a fully-developed wind farm whose upper edge is bounded by a self-similar boundary layer. The exchanges between the wind farm, the overlaying boundary layer, and the outer flow are parameterized by means of the classical entrainment hypothesis. We obtain a concise expression for the wind farm's power density (corresponding to power output per unit planform area), as a function of three coefficients, which represent the array thrust and the turbulent exchanges at each of the two interfaces. We assess the performance of our model by comparing the predicted power density to field data, laboratory measurements and large-eddy simulations for the fully-developed regions of wind farms, finding good agreement. We extend our model to include the effect of atmospheric stability on power output, by using an established parameterization relating entrainment coefficients to local Froude numbers. Our predictions for power variation with atmospheric stability are in broad agreement with field measurements. To the best of our knowledge, this constitutes the first quantitative comparison between an atmospheric-stability-dependent model and field data. Finally, we consider an ideal limit for array operation, whereby turbines are designed to maximize momentum exchange with the overlying boundary layer. This enables us to obtain an upper bound for the performance of large wind farms, which we determine to be an order of magnitude larger than the output of contemporary turbine arrays.

Tables

Table 1. Table 1: Wind farm data used in the preparation of Fig.  6 . Values of P / P 1 𝑃 subscript 𝑃 1 P/P_{1} for Horns Rev and Nysted correspond to case ‘ER’, from Fig. 5 of Ref.  [ 4 ] . Lillgrund data for P / P 1 𝑃 subscript 𝑃 1 P/P_{1} is adapted from Fig. 30 of Ref.  [ 34 ] . Power data for the last row is used in each case. † † \dagger assumes C t ≃ 0.75 similar-to-or-equals subscript 𝐶 𝑡 0.75 C_{t}\simeq 0.75 , ‡ ‡ \ddagger assumes C p ≃ 0.4 similar-to-or-equals subscript 𝐶 𝑝 0.4 C_{p}\simeq 0.4 . § § \S uses U o / U ∞ subscript 𝑈 𝑜 subscript 𝑈 U_{o}/U_{\infty} from Ref.  [ 40 ] , where the precursor method is discussed.

Source	$P/P_1$	$C_p$	$s_x$	$s_y$	$h_f/H$	$\frac{U_{o,\infty }}{U_{\infty }}$	$\frac{U_o}{U_{o,\infty }}$	$c_{ft}^{\prime }$	$c_{fp}\times {10}^3$	Ref.	Case	Symbol
Horns Rev	0.63	0.44	7	7	N/A	1.11	1	0.0249	$3.24$	[4]	“Aligned”	$\bullet$
Nysted	0.61	0.45	10.3	5.8	N/A	1.11	1	0.0233	$2.63$	[4]	“Aligned”	$\bullet$
Lillgrund	0.30	0.42	4.3	3.3	N/A	1.12	1	0.0863	$5.03$	[34]	“Aligned”	$\bullet$
Experiment	0.558	0.328	6	3	0.15	1.11	1.10	0.0582†	4.34	[35]	“Uniform”	-
Experiment	0.595	0.316	6	3	0.15	1.11	1.10	0.0582†	4.46	[35]	“Row-by-row”	+
Experiment	0.664	0.305	6	3	0.15	1.11	1.10	0.0582†	4.81	[35]	“Col.-by-col.”	*
Experiment	0.623	0.305	6	3	0.15	1.11	1.10	0.0582†	4.51	[35]	“Checkerboard”	$\times$
Full-farm LES	0.550	0.4‡	7.85	3.49	0.075	1.11§	1.07	0.0382	3.82	[36]	“Aligned”	$▷$
Full-farm LES	0.569	0.4‡	7.85	3.49	0.075	1.11§	1.07	0.0382	3.95	[36]	“Staggered”	$\diamond$
Full-farm LES	0.619	0.4‡	7.85	5.23	0.075	1.11§	1.06	0.0255	2.91	[36]	“Aligned”	$▷$
Full-farm LES	0.651	0.4‡	7.85	5.23	0.075	1.11§	1.06	0.0255	3.07	[36]	“Staggered”	$\diamond$
Full-farm LES	0.549	0.445	7	7	0.105	1.13	1.07	0.0263	2.20	[37]	“$\pm 1^{\circ }$”	$△$
Full-farm LES	0.655	0.4‡	7.85	5.24	0.075	1.11§	1.05	0.0255	3.13	[38]	“A3, aligned”	$\square$
Full-farm LES	0.580	0.4‡	5.24	5.24	0.075	1.11§	1.06	0.0381	4.07	[38]	“C4, aligned”	$\square$
Full-farm LES	0.462	0.4‡	5.24	3.49	0.075	1.11§	1.07	0.0573	4.74	[38]	“D3, aligned”	$\square$
Full-farm LES	0.660	0.4‡	7.85	5.24	0.075	1.11§	1.05	0.0255	3.16	[39]	“A3, staggered”	$▽$
Full-farm LES	0.563	0.4‡	5.24	5.24	0.075	1.11§	1.06	0.0382	3.96	[39]	“C4, staggered”	$▽$
Full-farm LES	0.542	0.4‡	3.49	7.85	0.075	1.11§	1.06	0.0382	3.81	[39]	“F3, staggered”	$▽$
Full-farm LES	0.453	0.4‡	3.49	5.24	0.075	1.11§	1.07	0.0573	4.66	[39]	“G3, staggered”	$▽$
Full-farm LES	0.783	0.4‡	7.85	7.85	0.075	1.11§	1.05	0.0170	2.55	[39]	“H3, staggered”	$▽$

Table 2. Table 2: Main sources of uncertainty, and their impact on c f ​ p subscript 𝑐 𝑓 𝑝 c_{fp} , for field measurements. The variation in P / P 1 𝑃 subscript 𝑃 1 P/P_{1} is from the published data for these turbine sites [ 4 , 34 ] . For the wind exponent α 𝛼 \alpha and height z o / h f subscript 𝑧 𝑜 subscript ℎ 𝑓 z_{o}/h_{f} for estimating U o subscript 𝑈 𝑜 U_{o} , an uncertainty of 20 % percent 20 20\% relative to reference values is assumed. Overall uncertainties in c f ​ p subscript 𝑐 𝑓 𝑝 c_{fp} are estimated by the root-sum-square method [ 41 ] .

Site	$[\frac{\delta \left(P/P_1\right)}{P/P_1}$	$\frac{\delta c_{fp}}{c_{fp}}]$	$[\frac{\delta \alpha }{\alpha }$	$\frac{\delta c_{fp}}{c_{fp}}]$	$[\frac{\delta \left(z_o/h_f\right)}{z_o/h_f}$	$\frac{\delta c_{fp}}{c_{fp}}]$	Overall $\frac{\delta c_{fp}}{c_{fp}}$	Ref.
Horns Rev	[0.29	0.29]	[0.2	0.03]	[0.2	0.48]	0.56	[4]
Nysted	[0.19	0.19]	[0.2	0.03]	[0.2	0.48]	0.51	[4]
Lillgrund	[0.085	0.085]	[0.2	0.03]	[0.2	0.48]	0.49	[34]

Table 3. Table 3: . Dependence of power output on stability conditions. Values of P / P 1 𝑃 subscript 𝑃 1 P/P_{1} and corresponding intervals of L / h f 𝐿 subscript ℎ 𝑓 L/h_{f} are calculated from Ref.  [ 14 ] . For Horns Rev, h f = 110 subscript ℎ 𝑓 110 h_{f}=110  m.

Conditions	$P/P_1$	$L/h_f$ interval	$U_o/U_{\infty }$	$c_{ft}^{\prime }$	$c_{fp}\times {10}^3$
Very stable	0.30	$[0.091,0.45]$	1.57	0.0249	$1.01$
Stable	0.61	$[0.45,1.8]$	1.47	0.0249	$1.22$
Neutral/Unstable	0.59	$[1.8,\infty ]$	1.11	0.0249	$3.24$

Equations

$${\overline{w}}_E=-E|U_b-U_o|,$$

$$(U_b-U_o)h_b$$

$$(U_b-U_o{)}^2{h}_b$$

$$\frac{\partial \overline{u}}{\partial x}+\frac{\partial \overline{w}}{\partial z}$$

$$\frac{\partial {\overline{u}}^2}{\partial x}+\frac{\partial \overline{u}\overline{w}}{\partial z}$$

$${\int }_0^H\frac{\partial \overline{u}}{\partial x}dz+{\overline{w}|}_{z=z_o}-{\overline{w}|}_{z=0}$$

$$\frac{d}{dx}\left(h_b{U}_b\right)=-\left({\overline{w}|}_{z=z_o}-U_o\frac{d{h}_b}{dx}\right).$$

$$\frac{d}{dx}\left(h_b{U}_b\right)=E|U_o-U_b|.$$

$${\int }_0^{z_o}\frac{\partial {\overline{u}}^2}{\partial x}dz$$

$$\frac{d}{dx}[h_b(U_b-U_o{)}^2]+2U_o\frac{d}{dx}[h_b(U_b-U_o)]=-U_o{\overline{w}|}_{z=z_o}-\nu {\frac{\partial u}{\partial z}|}_{z=0}.$$

$$\frac{d}{dx}[h_b{U}_b^2]=-U_o\left({\overline{w}|}_{z=z_o}-U_o\frac{d{h}_b}{dx}\right)-\nu {\frac{\partial u}{\partial z}|}_{z=0}.$$

$$\frac{d}{dx}[h_b{U}_b^2]=E{U}_o|U_o-U_b|-\nu {\frac{\partial u}{\partial z}|}_{z=0}.$$

$$\frac{\partial ⟨\overline{u}⟩}{\partial x}+\frac{\partial ⟨\overline{w}⟩}{\partial z}$$

$$\frac{\partial ⟨\overline{u}{⟩}^2}{\partial x}+\frac{\partial ⟨\overline{u}⟩⟨\overline{w}⟩}{\partial z}$$

$$⟨u⟩(x,y,z)=\frac{1}{{\lambda }_x{\lambda }_y}{\int }_{x-{\lambda }_x/2}^{x+{\lambda }_x/2}{\int }_{y-{\lambda }_y/2}^{y+{\lambda }_y/2}u(x^{\prime },y^{\prime },z)d{x}^{\prime }d{y}^{\prime },$$

$$0=\frac{\partial }{\partial z}\left(-⟨\overline{u^{\prime }{w}^{\prime }}⟩-⟨{\overline{u}}^{\prime \prime }{\overline{w}}^{\prime \prime }⟩+\nu \frac{\partial ⟨\overline{u}⟩}{\partial z}\right)+f.$$

$${\left(⟨\overline{u^{\prime }{w}^{\prime }}⟩+⟨{\overline{u}}^{\prime \prime }{\overline{w}}^{\prime \prime }⟩\right)}_{z=h_f}={\int }_0^{h_f}fdz-\nu {\frac{\partial ⟨\overline{u}⟩}{\partial z}|}_{z=0}.$$

$$h_f{U}_f={\int }_0^{h_f}⟨\overline{u}⟩dz.$$

$${\int }_0^{h_f}fdz=-\frac{c_{ft}^{\prime }}{2}{U}_f^2$$

$$c_{ft}^{\prime }=\frac{T}{\frac{1}{2}\rho U_f^2{A}_{plan}},$$

$$-{\left(⟨\overline{u^{\prime }{w}^{\prime }}⟩+⟨{\overline{u}}^{\prime \prime }{\overline{w}}^{\prime \prime }⟩\right)}_{z=h_f}=C_M(U_b-U_f{)}^2$$

$$\frac{c_{ft}^{\prime }}{2}{U}_f^2=C_M{\left(U_b-U_f\right)}^2.$$

$$\frac{d}{dx}\left(h_b{U}_b\right)$$

$$\frac{d}{dx}(h_b{U}_b^2)$$

$$C_M(U_b-U_f{)}^2=E(U_o-U_b{)}^2.$$

$$U_b=\frac{(C_M/E{)}^{\frac{1}{2}}{U}_f+U_0}{1+(C_M/E{)}^{\frac{1}{2}}}.$$

$$U_f$$

$$U_b$$

$$\frac{d{h}_b}{dx}$$

$$c_{fp}=\frac{P_{farm}/A_{farm}}{\frac{1}{2}\rho U_o^3}=c_{ft}^{\prime }{\left(\frac{U_f}{U_o}\right)}^3$$

$$c_{fp}=\frac{c_{ft}^{\prime }}{{\left[\left(C_M^{-\frac{1}{2}}+E^{-\frac{1}{2}}\right){\left(\frac{c_{ft}^{\prime }}{2}\right)}^{\frac{1}{2}}+1\right]}^3}.$$

$$U_f$$

$$U_b$$

$$c_{fp}=\frac{c_{ft}^{\prime }}{{\left[\left(C_M^{-\frac{1}{2}}+E^{-\frac{1}{2}}\right){\left(\frac{c_{ft}^{\prime }+c_d^{\prime }}{2}\right)}^{\frac{1}{2}}+1\right]}^3}.$$

Full text

An entrainment model for fully-developed wind farms: effects of atmospheric stability and an ideal limit for wind farm performance

Paolo Luzzatto-Fegiza111 fegiz@engineering.ucsb.edu and Colm-cille P. Caulfieldb,c

a Department of Mechanical Engineering, University of California, Santa Barbara, CA 93106, USA

b BP Institute for Multiphase Flow, Madingley Rise, Madingley Road, Cambridge, CB3 0EZ, UK

c Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Rd, CB3 0WA, UK

Abstract

While a theoretical limit has long been established for the performance of a single turbine, no corresponding upper bound exists for the power output from a large wind farm, making it difficult to evaluate the available potential for further performance gains. Recent work involving vertical-axis turbines has achieved large increases in power density relative to traditional wind farms (Dabiri, J.O., J. Renew. Sust. Energy 3, 043104 (2011)), thereby adding motivation to the search for an upper bound. Here we build a model describing the essential features of a large array of turbines with arbitrary design and layout, by considering a fully-developed wind farm whose upper edge is bounded by a self-similar boundary layer. The exchanges between the wind farm, the overlaying boundary layer, and the outer flow are parameterized by means of the classical entrainment hypothesis. We obtain a concise expression for the wind farm’s power density (corresponding to power output per unit planform area), as a function of three coefficients, which represent the array thrust and the turbulent exchanges at each of the two interfaces. Before seeking an upper bound on farm performance, we assess the performance of our simple model by comparing the predicted power density to field data, laboratory measurements and large-eddy simulations for the fully-developed regions of wind farms, finding good agreement. Furthermore, we extend our model to include the effect of atmospheric stability on power output, by using a parameterization (which had been previously developed in the context of geophysical fluid dynamics) relating entrainment coefficients to local Froude numbers. Our predictions for power variation with atmospheric stability are in broad agreement with field measurements. To the best of our knowledge, this constitutes the first quantitative comparison between an atmospheric-stability-dependent model and field data. Finally, we consider an ideal limit for array operation, whereby turbines are designed to maximize momentum exchange with the overlying boundary layer. This enables us to obtain an upper bound for the performance of large wind farms, which we determine to be an order of magnitude larger than the output of contemporary turbine arrays.

1 Introduction

While it has been argued that wind energy presents a large untapped potential [1], and that it could realistically mitigate climate change [2], a persistent shortcoming of wind farms is that their power density (defined as power output per unit land area) is low, reaching typically only a few W/m2 in a large farm [3, 4]. For this reason, the production of substantial amounts of electricity from wind requires the construction of extensive wind farms, which can comprise hundreds of units (an example is given by the Horns Rev facility, shown in Fig. 1A). The flow physics of these large turbine arrays are fundamentally different from those of a single turbine operating in isolation; in a large wind farm, most of the fluid kinetic energy enters not from the front, but rather through turbulent motions at the top of the array [5]. Remarkably, while a theoretical limit has long been established for the performance of a single turbine [6, 7], no corresponding theory appears to exist for a general, large-scale energy extraction array. Recently, Ref. [8] obtained numerical results showing that power increases of the order of 10% may be possible, by applying optimal control techniques to traditional turbine designs. Work with novel layouts of vertical-axis turbines has demonstrated large performance gains relative to conventional turbine arrays [3]. Somewhat surprisingly, the power density that has been measured is even larger than the value that would be achieved if a traditional wind farm could be operated without interference losses (i.e. without changes in turbine spacing). These results add urgency to the search for a theoretical upper bound.

Previous work on theoretical modeling of turbine arrays has focused on two main methodologies. The most commonly used wind-farm theories are based on superposing empirical models for wakes of individual turbines [5, 10]. As such, these do not consider effects arising from nonlinear wake interactions, which can be significant in affecting momentum fluxes in large arrays (as is the case, for example, in vegetation canopies [11]). In sufficiently large wind farms, the wakes eventually merge laterally, and the flow is often approximated as a canonical turbulent boundary layer. The wind farm is treated by resorting to the classic concept of effective roughness [5, 12], whereby the wind farm drag is imposed through a local increase in the effective roughness height. This provides a concise approximation that can be readily implemented in atmospheric models, provided the new roughness height is known. However, this approach is constrained by the need for semi-empirical formulae relating the effective roughness height (which is different from the actual obstacle size) to the physical turbine design and layout [12]. These limitations make current wind farm models unsuitable for examining arrays of energy-extraction devices of arbitrary design and layout.

A further issue that has received attention involves the thermal coupling between the atmosphere and the array. While much work has focused on the effect of large wind farms on surface temperatures [13], recent observations have shown that power output is increased when the atmosphere transitions from a stable condition to a neutral or unstable state (for a given wind velocity [14, 15]). Modeling of these atmospheric stability effects has focused on extending classical effective-roughness models (by making use of Monin-Obukhov similarity theory [16, 17]), and on estimating corresponding corrections for models describing individual wakes [17]. However, to the best of our knowledge, there are no investigations of atmospheric stability effects that have provided a quantitative comparison between theory and field data.

In this paper, we develop an entrainment-based model for power output from a large wind farm, and employ this theory to address the outstanding questions raised above. In Sec. 2, we introduce our entrainment-based wind farm model, and compare its predictions to field measurements from full-scale wind farms. In Sec. 3, we extend the model to encompass the effect of atmospheric stability, and test its predictions against available full-scale data. In Sec. 4 we consider the limit of an ‘ideal’ wind farm. The results are discussed briefly in Sec. 5. Finally, conclusions are presented in Sec. 6.

2 An entrainment model for large wind farms

In this section, we employ the classical entrainment hypothesis (generally attributed to G. I. Taylor, see e.g. Refs. [18], [19]) to develop a model for flow in a fully-developed canopy. Since the entrainment hypothesis has not previously been used in the context of canopy flows or wind farm aerodynamics (to the best of our knowledge), the next subsection briefly introduces the key elements of this turbulence closure.

2.1 The entrainment hypothesis, and its implications for boundary layers

The fundamental assumption underlying the entrainment hypothesis concerns the rate at which the interface between turbulent and nonturbulent fluid advances into the nonturbulent region. While the assumption is usually discussed in the context of plumes or gravity currents, here we consider a turbulent boundary layer, with outer velocity $U_{o}$. The (Reynolds-averaged) velocity $\overline{w}_{E}$ at which low-turbulence fluid crosses into the turbulent interface, measured in a frame of reference moving with the interface, is assumed to be [18, 19]:

[TABLE]

where the minus sign corresponds to the fact that turbulent fluid is below the low-turbulence fluid, in this example. Here, $U_{b}$ is a characteristic velocity for the interior of the boundary layer (defined below), an overline denotes Reynolds averaging, and $E$ is the entrainment coefficient, which is an empirical, nondimensional parameter. Much effort has been devoted to measuring $E$ in laboratory experiments, as well as in the ocean and atmosphere, as a function of Reynolds number and ambient stratification [20]. In this parameterization, $U_{b}$ is defined by the relations [19]

[TABLE]

where $\bar{u}$ is the Reynolds-averaged $x$-component of velocity, $z$ is the upward-pointing direction, and the upper bound of integration $H$ is chosen so that $\overline{u}(z=z_{o})=U_{o}$. One should choose $z_{o}$ to be large enough to be outside the wind farm’s boundary layer, but not so large that the resulting $U_{o}$ is no longer representative of the flow velocity driving the wind farm. Note that, in general, the value of $U_{o}$ will depend on atmospheric conditions. To obtain a corresponding set of depth-integrated equations, which can incorporate (1), we start from the Reynolds-averaged mass and momentum equations, assuming the flow is steady and statistically independent of the transverse direction $y$, and retaining only the leading Reynolds and viscous stresses:

[TABLE]

where the Reynolds decomposition for each quantity is denoted by $u=\overline{u}+u^{\prime}$, such that $\overline{u}$ is the ensemble-average (often implemented as a time-average) and $u^{\prime}$ denotes a fluctuation from the ensemble, satisfying $\overline{u^{\prime}}=0$. In (5), $\nu$ is the kinematic viscosity and we consider a zero pressure gradient flow. We integrate (4) from $z=0$ to $z=z_{o}$ to obtain

[TABLE]

where we have used the no-through-flow condition at $z=0$, and have assumed ${\partial U_{o}}/{\partial x}=0$ and $\left.\overline{u^{\prime}w^{\prime}}\right|_{z={z_{o}}}=0$. We rearrange (6) to obtain

[TABLE]

Note that, by the chain rule of differentiation, $U_{o}\frac{dh_{b}}{dx}=\frac{dh_{b}}{dt}$ is the vertical velocity of the turbulent interface. Therefore the term in parentheses on the right-hand-side of (7) is the vertical fluid velocity in a frame of reference moving with the turbulent interface, that is, $\overline{w}_{E}$. Using (1), we obtain the integral form of the continuity equation (4), under the entrainment hypothesis:

[TABLE]

Similarly, we integrate the momentum equation (5) as follows. Note that the first term on the left-hand-side of (5) becomes

[TABLE]

such that the momentum equation is

[TABLE]

Expanding the terms in the left-hand side and rearranging we have

[TABLE]

Using the definition of the entrainment velocity, we finally obtain the integral version of the momentum equation under the entrainment hypothesis:

[TABLE]

When discussing stresses in turbulent boundary layers, it is common to consider the idealized case of a flow that is very slowly evolving in the $x$-direction [21], such that the total stress (from turbulent fluctuations and viscosity) is constant. In such a case, (12) implies $E\,U_{o}|U_{o}-U_{b}|\simeq\nu\left.\frac{\partial u}{\partial z}\right|_{z=0}=\tau$, where $\tau$ is the wall stress. Since this is equal to the magnitude of the Reynolds stress in the boundary layer, we have shown that the first term on the right-hand-side of (12) corresponds to the Reynolds stress.

2.2 A two-interface entrainment model for fully-developed flow in wind farms

We distinguish three horizontal regions involving qualitatively different dynamics, as shown in Fig. 2. Starting from the bottom, these layers comprise the wind farm, a boundary layer and an outer region far above the farm. This implies that fluxes are governed by processes at two interfaces, which we label as the ‘farm’ and ‘outer’ interfaces. We will assume that the properties of the upper interface are primarily determined by atmospheric conditions, whereas those of the lower interface are directly affected by the wind farm design, as described below.

To model the flow through the farm, we consider the double-averaged canopy equations, in a zero pressure gradient, where we retain only the leading-order stresses [22, 23]

[TABLE]

These equations can be obtained, for example, by starting with the ensemble-averaged equations, and further decomposing $\overline{u}=\langle\overline{u}\rangle+\overline{u}^{\prime\prime}$, where the angle brackets denote an horizontal average over one canopy wavelength in $x$ and $y$, that is

[TABLE]

and $\overline{u}^{\prime\prime}$ is defined such that $\langle\overline{u}^{\prime\prime}\rangle=0$. Substituting $\overline{u}=\langle\overline{u}\rangle+\overline{u}^{\prime\prime}$ into (4, 5), and averaging over the canopy wavelength, one obtains (13, 14). The resulting dispersive stress $\langle{\overline{u}^{\prime\prime}\overline{w}^{\prime\prime}}\rangle$ accounts for exchanges of momentum due to flow features smaller than the spacing between canopy elements. In the momentum equation (14), $f$ is a body force corresponding to the drag exerted by the canopy, per unit mass.

We consider the fully-developed regime, such that the velocity in the farm is independent of downstream distance $x$. This assumption is consistent with power measurements in large wind farms [4, 14], as also shown in Fig. 1(b), and corresponds, inside the farm, to $\partial\langle\overline{u}\rangle/\partial x=0$. By continuity, this implies $\partial\langle\overline{w}\rangle/\partial z=0$. Since $\langle\overline{w}\rangle=0$ at the wall, we must have $\langle\overline{w}\rangle=0$ throughout the wind farm. Then the momentum equation (14) is

[TABLE]

Integrating this from $z=0$ to $z=h_{f}$, where $h_{f}$ is height of the wind farm:

[TABLE]

We now wish to express the stress and force terms in terms of depth-averaged velocities. To this end, let us define the canopy velocity as

[TABLE]

Assuming $U_{f}>0$, the vertically-integrated wind farm drag is written as (this step is similar to existing theories for vegetation canopies, and to other wind farm models [24, 11, 5])

[TABLE]

The coefficient $c_{ft}^{\prime}$ accounts for the turbines’ individual thrust and layout, and is defined as [22, 12]

[TABLE]

where $\rho$ is the air density, $T$ is the thrust from each turbine, the planform area associated with one turbine is $A_{plan}=s_{x}s_{y}D^{2}$, and $s_{x},s_{y}$ are the array spacings (normalized by the turbine diameter $D$) in the two horizontal directions. The prime in $c_{ft}^{\prime}$ emphasizes the fact that a local fluid velocity (rather than freestream) is used (consistently with notation, for example, in Ref. [25]).

In accordance with studies of vegetation canopies [26, 27], we assume that the combination of Reynolds and dispersive stresses can be parameterized as

[TABLE]

where $C_{M}$ is a momentum-exchange coefficient, which is expected to vary with the geometry of the array, the design of the turbines, Reynolds number, and ambient stratification (as encoded by an appropriate Froude number). A brief discussion of this parameterization, including possible appropriate values for $C_{M}$ (in a wind farm and in other canopy flows) is provided below in Sec. 2.3. We will also develop a simple heuristic model of momentum exchange, for the specific case in which the farm interface consists of a classical mixing layer.

Neglecting viscous stresses, the integral form of the momentum equation inside the wind farm is therefore

[TABLE]

Following a derivation almost identical to the one described in the previous section, the flow in the overlaying boundary layer is described by the following integral forms of the mass and momentum equations (assuming $U_{o}>U_{b}>U_{f}$)

[TABLE]

To summarize, in this fully-developed regime, the horizontal momentum that is drained by the wind farm (and ultimately converted into electrical power) is replenished by vertical turbulent exchange processes with the boundary layer above the turbine array, which therefore take place at the ‘farm interface’, as sketched in Fig. 2.

Since $U_{f}$ does not change with downstream distance, (22) implies that $U_{b}$ must also be independent of $x$. As the boundary layer grows, outer fluid is entrained at the upper edge of the boundary layer, and the vertical extent $h_{b}$ of the boundary layer is expected to increase. To get an explicit expression for $U_{f}$ as a function of $(c_{ft}^{\prime},E,C_{M})$, we first take $U_{b}$ out of the derivatives on the left-hand sides. Then multiplying (23) by $U_{b}$ and subtracting (24) to eliminate $d{h_{b}}/dx$, we obtain

[TABLE]

Since we assumed that $U_{o}>U_{b}>U_{f}$, taking the square root and rearranging we get a linear relation between $U_{b}$ and $U_{f}$

[TABLE]

Similarly, taking the square root of (22) yields $U_{b}=U_{f}\{1+[c_{ft}^{\prime}/(2C_{M})]^{1/2}\}$. Subtracting this from (26) to eliminate $U_{b}$, and rearranging for $U_{f}$, we obtain

[TABLE]

Note that, since $U_{o}$ and $U_{b}$ are both independent of $x$, (29) implies that the boundary layer height increases linearly. As a first check of our theory, we note that this linear growth, for $h_{b}$, is consistent with experimental results for vegetation canopies (see, for example, Fig. 6 of Ref. [28]).

The power density (corresponding to the power extracted per unit land area) is given by $(c_{ft}^{\prime}/2)\rho U_{f}^{3}$, to a first approximation. We define a nondimensional power density coefficient, for the whole farm, by normalizing the power density (say, $P_{farm}/A_{farm}$) using the energy flux in the outer region, such that

[TABLE]

and, using (27), we obtain the following concise expression:

[TABLE]

Finally, we note that while the theory presented so far neglects the effect of bottom friction (as may be reasonable to do for offshore wind farms [29]), it is not difficult to introduce ground dissipation effects, by adding a force, within the canopy, equal to $(c_{d}^{\prime}/2)\rho U_{f}^{2}$, where $c_{d}^{\prime}$ is a bottom drag coefficient. Then it is sufficient to replace $c_{ft}^{\prime}$ with $c_{ft}^{\prime}+c_{d}^{\prime}$ in equations (27) and (28) for $U_{f}$ and $U_{b}$, and so

[TABLE]

whereas $c_{fp}$ remains equal to $c_{ft}^{\prime}(U_{f}/U_{o})^{3}$ (where $U_{f}$ is of course now given by (32)), leading to

[TABLE]

To assess the accuracy of our expression for $c_{fp}$, we need to first estimate $E$, $C_{M}$, $c_{ft}^{\prime}$ and $c_{d}^{\prime}$ for existing wind farms. For the entrainment coefficient $E$, we note that a wide range of measurements in turbulent flows have found that, at very large Reynolds numbers, the entrainment coefficient is of the order $E\sim 0.1$ (see for example the large dataset compiled in Ref. [20] for ocean overflows, which is reproduced in Sec. 3). The calculations needed to estimate $C_{M}$ and to evaluate $(c_{ft}^{\prime},c_{fp})$ from field data are described in the next two subsections.

2.3 Estimating $C_{M}$ for an interface comprising a classical mixing layer

As noted earlier, parameterizations of the form $C_{M}(U_{b}-U_{f})^{2}$, for the vertical momentum flux, have been widely employed in the context of modeling canopy flows [27]. One way to argue that this is a relevant scaling, specifically in a turbulent shear flow, is to use the classic closure ${-\overline{u^{\prime}w^{\prime}}\propto|\partial\overline{u}/\partial z|(\partial\overline{u}/\partial z)}$ and note that the velocity scale for the gradient is $(U_{b}-U_{f})$, such that

[TABLE]

To estimate $C_{M}$, we need a more specific physical model for the flow dynamics at the farm interface (which separates the farm and the boundary layer). Here we present a simple model, which may also be applicable to the study of other canopy flows. Recent field measurements suggest that the farm interface comprises vortical structures which closely resemble a mixing layer [30]. We therefore assume a classical mixing layer structure, with thickness much smaller than $h_{b}$ or $h_{f}$, such that the mixing layer dynamics might be considered separately from the overall wind farm flow, and subsumed in a particular choice of $C_{M}$. We expect that, after each turbine, the mixing layer grows downstream, entraining fluid from the farm and from the boundary layer, and mixing together the entrained fluid from these two regions, as sketched in Fig. 3(a). The characteristic mixing layer velocity is taken as $U_{mix}=(U_{b}+U_{f})/2$, in accordance with classical mixing-layer models [21]. As the mixing layer grows, its thickness $\delta_{mix}$ and momentum flow $M_{mix}$ increase with distance downstream, as the layer entrains fluid from each side (see Fig. 3):

[TABLE]

Note that the mixing layer overlaps with both the farm and boundary layer regions. We assume that, to a first approximation, the mixing layer protrudes by the same distance into the two regions, and, when encountering the following turbine, is partitioned equally between the two domains, as sketched in Fig. 3(b), thereby preventing the mixing layer from growing continuously along the length of the farm. The net momentum flux into the farm can therefore be expressed by the loss to the mixing layer, together with a gain corresponding to half the overall momentum flux into the mixing layer, i.e.

[TABLE]

This provides a simple heuristic argument suggesting that the momentum flux from the boundary layer to the farm is proportional to $(U_{b}-U_{f})^{2}$. Incidentally, the functional dependence that we obtained from the entrainment hypothesis can also be shown to be consistent with the well-known experimental result [31], for mixing layers at high Reynolds numbers, that the mixing layer thickness (say, $\delta_{mix}$) grows in proportion to $(U_{b}-U_{f})/(U_{b}+U_{f})$. To see this, since $dU_{mix}/dx=0$, (36) implies

[TABLE]

such that expanding $U_{mix}=\frac{1}{2}(U_{b}+U_{f})$ and rearranging gives

[TABLE]

consistently with the classic experimental results of Ref. [31]. Heuristically, one may expect $C_{M}$ to be independent of $c_{ft}^{\prime}$, to a first approximation, since the immediate effect of increasing $c_{ft}^{\prime}$ (without otherwise changing the geometry, or drastically changing the Reynolds number) is to reduce $U_{f}$, and therefore increase the velocity difference across the interface. However, for mixing layers without stratification, experimental evidence indicates that $E$ is largely independent of the velocity difference [31].

Furthermore, comparing (37) to the momentum flux parameterization $C_{M}\rho(U_{b}-U_{f})^{2}$ used in Sec. 2, we can deduce that $C_{M}=E/4$ in contemporary wind farms, provided the farm interface resembles a mixing layer. Since $E\sim 0.1$, we note that the value of $C_{M}$ from our estimate is also consistent with empirical values for vegetation canopies [27]. Note, however, that canopies comprising obstacles with more complex geometries and layouts may have different values of $C_{M}$. Indeed, we believe that a key goal in designing new wind farms should be that of maximizing $C_{M}$, as discussed later in Sec. 5.

2.4 Using field measurements to calculate $c_{ft}^{\prime}$ and $c_{fp}$ for large wind farms

Following classical wind turbine theory, the thrust from one turbine is:

[TABLE]

where $C_{t}$ is the turbine’s thrust coefficient and $A_{rotor}=\pi D^{2}/4$. Note that the conventional definition of $C_{t}$ makes use of the freestream velocity upstream of the turbine (at rotor height), which we denote as $U_{\infty}$. As noted in Ref. [12], to relate this to the velocity at the rotor, one can use classical actuator disk theory:

[TABLE]

Consistently with wind farm models based on the effective-roughness concept [12, 5], we assume that $U_{rotor}\simeq U_{f}$. We can combine (40) and (41) with (20) to obtain

[TABLE]

which enables the estimation of $c_{ft}^{\prime}$ from field measurements.

Now consider $c_{fp}$, which we defined in Sec. 2 as

[TABLE]

where $P$ is the power output from one turbine. In studies of turbine arrays, $P$ is typically reported as a fraction of $P_{1}$, which is the power output of the first turbine in the row. This makes it convenient to write $P=(P/P_{1})P_{1}$. Using classical turbine theory, $P_{1}$ can be written:

[TABLE]

where $C_{p}$ is the power coefficient. Therefore, $c_{fp}$ is given by:

[TABLE]

Unfortunately, $U_{o}$ (the velocity far above the wind farm) is not easily measured in the field. Here we use the classical power-law approximation $U(z)=U_{\infty}(z/h_{hub})^{\alpha}$. To decrease sensitivity to the choice of $z$, we take $U_{o}$ as the average between $z=h_{f}$ and $z=2h_{f}$, such that $U_{o}/U_{\infty}=(2^{\alpha+1}-1)/(\alpha+1)$.

We first calculate $c_{fp}$ from the field measurements of Refs. [4], [32]. Since these consist of long-time averages, and the corresponding locations exhibit predominantly neutral conditions [4, 33], we neglect the effect of ambient stratification at this stage, and set $\alpha\simeq 0.12$ (see Ref. [33]). The key parameters used in calculating $c_{ft}^{\prime}$ and $c_{fp}$ are summarized in Table 1; the resulting values are displayed by the symbols in Fig. 6.

Note that there is appreciable uncertainty in estimating $U_{o}/U_{\infty}$. Table 2 provides a breakdown of the main sources of uncertainty in the calculation required to obtain $c_{fp}$. For the uncertainty $\delta(P/P_{1})$ we use available published data, whereas variations of $20\%$ are assumed for the wind exponent and upper bound of integration of the wind power-law to obtain $U_{o}$. The overall uncertainty is calculated by the root-sum-square method, and is shown by the error bars in Fig. 6, which combine the uncertainty from the original measurement together with the uncertainty that is introduced through the calculation of $c_{fp}$.

[FIGURE:]

2.5 Estimating a blockage correction for experiments and full-farm LES

For wind tunnel measurements, even in an empty test section, boundary layers effectively reduce the cross sectional area available to the flow, causing the outer flow velocity to gradually increase with downstream distance. This effect is of course stronger for arrays of obstacles, since their drag causes the flow to slow down locally, thereby forcing a substantial amount of fluid out of the array and into the outer flow, such that $U_{o}$ is increased. Consistently with classic work on wings and bluff-bodies [42], the resulting blockage is greater than the value that would be obtained by just considering the contraction associated with the cross-sectional area of the objects. This is true in experiments [35] as well as simulations that employ a no-through-flow boundary condition at the top of the domain [38, 36, 39, 37]. The outer flow velocity $U_{o}$ is typically not reported, but the inflow profile is described. This must be related to the outer flow velocity $U_{o}$ above the wind farm, at a location where the farm flow is fully developed. While $U_{o}$ increases with downstream distance, we assume that $U_{o}$ varies slowly, by comparison to the distance over which turbulent fluxes equilibrate, so that the steady model from Sec. 2 still applies.

In the spirit of classic blockage models [42], and consistently with the model in Sec. 2, we use a layer-wise approach, illustrated in Fig. 4. Continuity implies

[TABLE]

where $H$ is the height of the test section. Equation (46) is written as

[TABLE]

such that, dividing both sides by $h_{f}$ and solving for the blockage correction factor $U_{o}/U_{o,\infty}$

[TABLE]

To evaluate (48), we use (32-33) for $U_{f}/U_{o},U_{b}/U_{o}$, and approximate $h_{b}/h_{f}\simeq(dh_{b}/dx)(x/h_{f})$, where $dh_{b}/dx$ is found from (29) and $x/h_{f}$ depends on the dataset. The blockage correction is reported in Table 1. Note that while a value of $E$ must be assumed, the blockage correction has only very weak sensitivity to the entrainment coefficient, since changing $E$, say, from $0.1$ to $0.16$ yields a variation of less than 1% in (48).

2.6 Estimating the bottom drag coefficient $c_{d}^{\prime}$

From (40,18), the stress due to the combined effect of the wind farm and of bottom drag is

[TABLE]

Consider the flow ahead of the wind farm, where $c^{\prime}_{ft}=0$. For $z<h_{f}$, the log-law is a reasonable approximation [35], such that

[TABLE]

where $\kappa$ is the von Kármán constant (approximately equal to 0.4) and $z_{0}$ is the equivalent roughness height for the bottom surface, which is assumed known. Substituting (50) into (49), integrating, and solving for $c_{d}^{\prime}$ yields

[TABLE]

Values of $c^{\prime}_{d}$ do not vary significantly in experiments and simulations. For example, [35] report $z_{0}/h_{f}\simeq 5.56\times 10^{-4}$, implying $c_{d}^{\prime}\ \simeq 0.0076$, whereas [40] have $z_{0}/h_{f}\simeq 9.77\times 10^{-4}$, yielding $c_{d}^{\prime}\ \simeq 0.0091$. For simplicity, in our theory, we take $c_{d}^{\prime}=0.008$.

2.7 Comparing the two-interface model with measurements

In order to perform a comparison between our theory and measurements, we still need to select a value for the entrainment coefficient $E$. For large-scale turbulent flows (such as gravity currents, jets and plumes) the entrainment coefficient is typically of the order $E\sim 0.1$ (e.g. Ref. [20]). There appear to be no precise published values of $E$ for turbulent boundary layers. To determine a value of $E$, we use published data for the growth of a boundary layer over wind turbine arrays. Since boundary layer thickness is not measured in full-scale wind farms (to the best of our knowledge), we compare our solution to the wind tunnel results of Ref. [43] for aligned and staggered wind turbines. In order to show an example at larger $c^{\prime}_{ft}+c^{\prime}_{d}$ also, we include the boundary layer for a flume experiment with a model canopy comprising vertical cylinders, from Ref. [27]. To express the boundary layer growth in terms of a standard quantity for boundary layers, consider the definition of the displacement thickness $\delta^{*}$ [44], and use equation (2)

[TABLE]

Differentiate by $x$, using the fact that $dU_{b}/dx=0$ in fully-developed flow

[TABLE]

where $U_{b}/U_{o}$ and $dh_{b}/dx$ are found from (33,29). Figure 5 shows that the two-interface theory with $E=0.16$ and $C_{M}=E/4=0.04$ is in agreement with experimental data.

We now use equation (34) and the values $E=0.16$ and $C_{M}=E/4=0.04$ to obtain an approximate value of $c_{fp}$, as a function of $c_{ft}^{\prime}$, for existing wind farms, and compare it to the data summarized in Table 1. The result is shown by the dashed line in Fig. 6($a$), where we take $c_{d}^{\prime}=0.008$, following the analysis in Sec. 2.6. The effect of varying $E$ and $C_{M}$ by $\pm 20\%$ is shown by the shaded region. In spite of the simplicity of our model, it is in agreement with available data. The plot of $c_{fp}$ versus $c_{ft}^{\prime}$ has a few notable features. Firstly, $c_{fp}$ is a very small number (of the order of $0.1\%$), expressing the fact that a ground-based wind farm can extract a small fraction of the power flux overhead. Secondly, as shown in Fig. 6($b$), at larger $c_{ft}^{\prime}$ the $c_{fp}$ curve reaches a maximum, which corresponds to optimizing over the possible turbine spacings and thrust settings. The value of $c_{ft}^{\prime}$ that maximizes power output is found by differentiating (34) and equating it to zero, yielding

[TABLE]

where $\zeta=(C_{M}^{-\frac{1}{2}}+E^{-\frac{1}{2}})^{-1}$. Here, the optimal ${c_{ft}^{\prime}}\simeq 0.179$, which corresponds to $c_{fp}\simeq 5.0\times 10^{-3}$. Assuming $C_{t}=8/9$ (corresponding to the Betz limit) and $s=s_{x}=s_{y}$, this would imply a turbine spacing of 3.0$D$, which is somewhat smaller than the value at the Lillgrund site (which has $s_{x}=4.3$, $s_{y}=3.3$), as illustrated in Fig. 6($b$). However, since the optimum in Fig. 6($b$) is very broad, there is little reason to seek such small turbine spacings. For example, the theory predicts that increasing $s$ by 50% from the optimal value will result in a power decrease of only about 6%.

Note that, according to (20), $c_{ft}^{\prime}$ is inversely proportional to the square of the turbine spacing, such that, for a given land extent, an increase in $c_{ft}^{\prime}$ corresponds to an increase in the number of installed turbines. This implies that increasing the number of turbines beyond a certain value can actually be detrimental to power output (with given land area), as the lower values of $U_{f}$ that result do not make up for the large values of $c_{ft}^{\prime}$.

3 Modeling the effect of atmospheric stability on wind farm performance

3.1 Linking $E$ and $C_{M}$ to the Obukhov length

Here we extend our theory to include the effect of atmospheric stability, by considering the influence of ambient conditions on the parameters $E$ and $C_{M}$. This approach accounts for the effect of stratification on the turbulent motions responsible for mass and momentum exchanges. We employ a recently introduced parameterization [20], which describes $E$ as an empirical function of the Reynolds number $Re$ and Froude number $Fr$, based on a large compilation of laboratory experiments and field observations involving ocean overflows; here we refer to this fit for $E$ using the expression $E_{ca10}(Re,Fr)$. For completeness, the dataset on which the parameterization is based is shown in Fig. 7, and the fit equations are (using the same notation as Ref. [20])

[TABLE]

where $E_{min}=4\times 10^{-5}$, $E_{max}=1$, $A=3.4\times 10^{-3}$, $B=243.52$, $\alpha=7.18$, $\beta=0.5$, $Fr_{0}=0.51$. In Ref. [20], the Reynolds and Froude numbers are based on the gravity current height and on the velocity difference across the turbulent interface, and are defined as

[TABLE]

where $g$ is the gravitational acceleration, $h$ is the vertical length-scale of the flow (corresponding to the vertical extent of the gravity current in ocean measurements), $\Delta U$ and $\Delta\rho$ are the velocity and density differences across the interface, and $\rho_{0}$ is a reference density.

From Fig. 7, it is apparent that for $Re$ larger than about $10^{4}$, $E_{ca10}$ can exceed 0.16 for large $Fr$. However, while there are atmospheric measurements that support $E$ as high as 0.2 [45], published data do not seem to support significantly larger values of $E$, as also noted in Ref. [46]. For this reason, we propose a modified version of the fit, which saturates at a value $E_{sat}$. While we could simply alter $E_{max}$ in (55), this has the unwanted effect of altering $E_{ca10}$ for all $Re$. At fixed $Re$, we therefore simply evaluate (55) as a function of $Fr$. Based on the value of the resulting $E_{ca10}$, our $E$ is given by

[TABLE]

where we set E_{cut}=0.8E_{sat}{\color[rgb]{0,0,0}=0.128} and $Fr_{cut}$ is the value of $Fr$ at which $E_{ca10}=E_{cut}$, which is found by interpolation. For $Re=10^{8}$, $Fr_{cut}\simeq 1.95$. This provides a smooth algebraic transition between the parameterization of Ref. [20] and $E_{sat}$, as shown by the red continuous and dashed lines in Fig. 7.

In our flow of interest, the characteristic scale $h$ is the farm height $h_{f}$, and we use the classic atmospheric approximation $\Delta\rho/\rho_{o}=-\Delta\theta/\theta_{o}$, where $\theta$ is the potential temperature. For each interface, we consider the corresponding differences in potential temperature and velocity to calculate the relevant $Fr$. Therefore we write, for each interface:

[TABLE]

where the subscripts to $Fr$ specify the interface considered (‘outer’ or ‘farm’).

We assume that a constant specific heat flux $q$ is applied at the bottom of the farm. By analogy with the momentum flux parametrization, we assume that the exchange of heat at the farm interface is proportional to $C_{\theta}(U_{b}-U_{f})(\theta_{b}-\theta_{f})$, where $C_{\theta}$ is an exchange coefficient analogous to $C_{M}$. Using the same subscripts as before to denote each region, the equations governing potential temperature are

[TABLE]

In the above, $c_{p}$ is the specific heat of air, which has been assumed constant. Consistently with our fully-developed assumption, we take the temperature in each region to be constant with downstream distance (such that the left-hand side of (60) is $U_{b}\theta_{b}dh_{b}/dx$, whereas the left-hand side of (61) vanishes). Then (60) and (61) are two linear equations for $\theta_{f}$, $\theta_{b}$, with constant coefficients. Using (29) for $dh_{b}/dx$, these expressions can be manipulated to yield the nondimensional temperature differences relative to the outer flow:

[TABLE]

We now need to relate the nondimensional heat flux $q/(c_{p}U_{o}\theta_{o})$ to quantities that are reported through field measurements. For example, as discussed later in the next subsection, the wind farm measurements of Ref. [14] are sorted using estimates of the Obukhov length $L$ [47]

[TABLE]

where $\tau$ is the shear stress in the boundary layer, and $\kappa\simeq 0.4$ is the von Kármán constant. Physically, $L$ (given by (64) above) expresses the height above which damping from the stratification exceeds the production of turbulent kinetic energy. Since turbulent production peaks close to the ground, a smaller $L$ corresponds to an increasingly stable atmosphere [21, 47]. Since here $\tau=\frac{1}{2}\rho(c_{ft}^{\prime}+c_{d}^{\prime})U_{f}^{2}$, we obtain

[TABLE]

With given values of $c_{ft}^{\prime},c_{d}^{\prime}$ and $L/h_{f}$ (as well as $gh_{f}/U_{o}^{2}$), the resulting set of equations (27, 28, 59$a$,$b$, 62, 63, 65) for $U_{f}$, $U_{b}$, $\theta_{f}$, $\theta_{b}$, $E$, $C_{M}$ and $q$ can be solved by iteration. We find that a simple relaxation method is sufficient to obtain convergence. For simplicity, we set $C_{\theta}=C_{M}$, which essentially corresponds to the assumption (common in turbulent flows [47]) that the turbulent Prandtl number is close to unity.

3.2 Power dependence on atmospheric stability in theory and field measurements

Measurements of atmospheric stability at wind turbine sites are currently quite limited. We are aware of two such published studies, namely Refs. [14, 48]. The former provides data for the well-known Horns Rev site, whereas the latter acquired data at a wind farm in the U.S. midwest (the specific site, layout, and turbine model do not appear to have been described). No data on ambient stability is available for the Nysted and Lillgrund sites, to the best of our knowledge. To avoid complications associated with neighboring topography, here we focus on the measurements [14] from the Horns Rev site. These used estimates of the Obukhov length $L$ to categorize atmospheric conditions. The field measurements of Ref. [14] found that as the atmosphere transitioned from stable to neutral, power output increased [14]. However, as the atmosphere shifted from neutral to unstable, power output did not increase further; for this reason, results were reported by binning together measurements involving nearly stable, neutral and unstable conditions [14].

Recall from Sec. 2.4 that extracting $c_{fp}$ from field data requires estimating $U_{o}$ above the wind farm. However, the value of $\alpha$ in the power-law approximation for $U_{o}$ depends on the Obukhov length $L$. Here we use the measurements of Ref. [24], for flow over relatively smooth terrain (equivalent roughness $z_{0}\simeq 0.01\,$m), which report $\alpha$ as a function of atmospheric stability category. To match these values of $\alpha$ to the corresponding $L$, for Horns Rev, we use the field measurements of Ref. [14] (which associate $L$ intervals to stability conditions), thereby obtaining the following approximate $(\alpha,L)$ pairings: (0.53, 50 m), (0.34, 200 m), (0.12, 500 m). An empirical function $\alpha(L)$ is then defined by using a piecewise cubic Hermite interpolant between these values, with constant value outside of the interval $L=[50\,\mathrm{m},500\,\mathrm{m}]$. The associated data are reported in Table 3, whereas the results are displayed in Fig. 8, as discussed below.

Our prediction for $c_{fp}$ versus $L/h_{f}$ is shown by the dashed line in Fig. 8, whereas field data [14] are shown by the shaded bars. The gray shaded region shows the sensitivity of the results to a change of $\pm 20\%$ in $E_{par}$, with given Froude numbers. The theory captures the trend in $c_{fp}$, and is in good quantitative agreement with records for stable and neutral/unstable atmospheres. Our model seems to underpredict $c_{fp}$ for very stable conditions. This might be due to the parameterization used for $E$, which assumes that the entrainment tends to a small value as the atmosphere becomes strongly stratified, as is the case in flows without large obstacles that can maintain fluid mixing. However, in a wind farm, the presence of turbines would drive the flow towards a larger $E$. In addition, the turbulent Prandtl number might be deviating from unity, as the flow becomes more strongly stratified (this could make $C_{\theta}$ depart significantly from $C_{M}$). Interestingly, the literature on vegetation canopies indicates that traditional Monin-Obukhov theory also underestimates momentum fluxes for plants in a stable stratification, although there does not seem an accepted explanation for this underestimate; see for example Ref. [49] and references therein. However, for wind farms, since very stable conditions were relatively rare in the field measurements [14] (corresponding to only a few percent of observations), this model limitation might not have a very significant practical impact.

4 An ideal limit for wind farm performance

We now turn to the problem of finding an upper bound for the relative power density, as the design of a wind farm is varied. The maximum of $c_{fp}$ over $c_{ft}^{\prime}$ was found earlier as (54). Since we focus on idealized performance, assume $c_{d}^{\prime}\ll c_{ft}^{\prime}$, such that

[TABLE]

and (31) yields

[TABLE]

The parameter $E$ controls the transfer of momentum at the interface between the outer flow and the boundary layer, and is unlikely to be affected by the design of the wind farm. However, we propose that it may be possible to redesign a wind farm to increase $C_{M}$, by promoting exchanges between the farm and the boundary layer. To establish an upper bound on $c_{fp}$, we take the limit as $C_{M}$ becomes much larger than $E$. In this idealized situation, the boundary layer mixes perfectly with the farm flow, and the performance is limited only by the dynamics at the interface between the boundary layer and the outer flow. The limiting value of $c_{fp}$ is then:

[TABLE]

If $E\simeq 0.16$, the maximum $c_{fp}$ is approximately $0.047$, which is an order of magnitude larger than the output of current wind farms, which is around $c_{fp}\simeq 0.005$ (as can be verified, for example, by inspecting values in Fig. 8).

In practice, however, reaching this limit requires $C_{M}^{1/2}\gg E^{1/2}$, implying $C_{M}$ would need to be of order ten, which is several hundred times current values (which we estimated around 0.04 in Sec. 2.3). Nevertheless, it is instructive to consider the change in $c_{fp}$ that would be accomplished by increasing $C_{M}$ through an order of magnitude, as plotted in Fig. 9. The regime accessed by current wind farms corresponds to $C_{M}\simeq 0.04$, $c_{fp}\simeq 0.005$. Our theory predicts, for example, that increasing $C_{M}$ to 0.4 would yield a $c_{fp}$ of approximately 0.018, corresponding to over three times its current value. This indicates that even moderate increases in the momentum exchange coefficient $C_{M}$ could lead to practically significant gains in turbine performance.

5 Discussion

Our model emphasizes turbulent and dispersive transport at two planar interfaces as the mechanism enabling energy extraction by the array. This appears to be a reasonable assumption in the fully developed regime; however, it would not hold near the front of a turbine array, where the turbine wakes are separated by regions of relatively clean air. This issue would likely result in strong disagreements between predicted and measured power outputs (which should scale with the average of the cube of $u$, rather than the cube of the average, which we used in defining $c_{fp}$).

A point worth noting is that the wind farm boundary layer is relatively thin; for example, for the Lillgrund wind farm (which has the fastest-growing boundary layer, on account of its larger $c_{ft}^{\prime}$), $dh_{b}/dx$ from our model is only 0.042. If we consider flow along the longest direction of the Lillgrund farm, we can estimate, at the end of the farm, $h_{b}\sim h_{f}$, implying that the top of the boundary layer reaches a height of the order $\sim 2h_{f}$. Consistently with this estimate, a recent simulation of the complete Lillgrund wind farm [33] has obtained excellent agreement with field measurements while using a computational domain having an overall vertical height of only $2.5D\simeq 2.1h_{f}$.

According to Table 1, removing wind shadowing effects between turbines (such that all turbines would output as much power as the first one) in a wind farm with dense spacing such as Lillgrund (sketched in Fig. 10$a$) would yield a power output that is roughly 3.3 times that of existing arrays (in the fully developed regime). This is in contrast with the ten-fold improvement identified by the ideal limit in Sec. 4. This ideal limit proposes that the array be redesigned to enhance the downward transport of momentum, relative to a traditional boundary layer (which constrains the operation of even a single turbine). As noted in Sec. 4, approaching the ideal wind farm limit requires a $C_{M}$ of around 10, which may not be achievable; it would be interesting to investigate bounds for the value of $C_{M}$, as this would help establishing a sharper ideal wind farm limit. Nevertheless, it is interesting to note that, if one could hypothetically achieve wake recovery over a length of approximately $4D$, the order-of-magnitude improvement in power density could be obtained by turbine layouts that are only slightly denser than for the Lillgrund array, as sketched in Fig. 10($b,c$).

We should briefly also comment here on recent works [3, 50] involving bio-inspired layouts of counter-rotating, vertical-axis turbines (sketched in Fig. 10$c$), which have reported increases in power density of one order of magnitude, relative to the fully-developed regime of Horns Rev and other large wind farms. Our model suggests that such an order-of-magnitude increase in power density might be difficult to achieve in the fully developed regime. Concerning this issue, it would be interesting to examine in detail the contribution of developing-flow effects to the large power output described in Refs. [3, 50]. It remains of great practical interest to develop simple and accurate theoretical models of developing flows in wind farms and other canopies. It would also be valuable to construct a general framework linking the kinematics of an array of obstacles (such as the rotating-cylinder configurations recently examined in Refs. [51, 52]) to a predicted value of the momentum exchange coefficient $C_{M}$, and therefore to achievable power output.

6 Conclusions

In this paper, we introduce an entrainment-based model for the power output from the fully-developed region of large wind farms. The model comprises a fully-developed canopy whose upper edge is bounded by a self-similar boundary layer, and can represent arrays of turbines with arbitrary design and layout. We obtain a concise closed-form expression for the power density, which is in agreement with available field data, experiments, and large-eddy-simulations. Our theory ultimately frames the performance of a large wind farm in terms of two parameters associated with turbulent and dispersive transport, namely $E$ (the entrainment coefficient at the edge of the boundary layer) and $C_{M}$ (the momentum exchange coefficient at the wind farm-boundary layer interface). A valuable feature of this approach is that it enables the use of existing parameterizations from the study of geophysical flows, and from canonical studies of turbulent flows.

By introducing equations governing the potential temperature in the farm and in the boundary layer, we extend our theory to describe the effect of atmospheric stability on power output, finding agreement with the majority of available field measurements (with the exception of very stable conditions, where the theory underpredicts the power output). The theory’s predictions for very stable conditions might be improved by new parameterizations that account for the effect of strong bottom roughness in stratified flows. We build on our model and define an ideal performance limit, for a large wind farm, by considering a turbine array that has been optimized to maximise energy exchanges at the interface between the array and the boundary layer (such that $C_{M}$ is relatively large). In this idealised situation, the power output is limited by the rate at which energy is entrained at the upper edge of the boundary layer. The limiting value of power density (normalized by the power flux far above the wind farm) is $8E/27$, where $E$ is the entrainment coefficient (typically of order 0.1 in high-$Re$, unstratified flows). This simple result constitutes a reference value against which one can measure the performance of existing and proposed wind farm designs, and corresponds to an order of magnitude increase relative to current turbine arrays.

In practice, reaching this ideal $8E/27$ limit requires an increase in $C_{M}$ by a factor of 400. Nevertheless, even moderate increases in $C_{M}$ could have a significant impact; for example, our theory predicts that increasing $C_{M}$ by 20% (above its current value) would yield a corresponding 13% increase in power output. Our work points to a new system-wide approach for increasing wind farm performance, focusing on maximizing the momentum exchange coefficient $C_{M}$ for the whole farm, rather than on optimizing individual turbine performance. We suggest here that this may be achieved by redesigning wind farms to control the large-scale turbulent structures governing momentum fluxes. This is in direct contrast with current practice, which relies on optimizing the layout of turbines, each of which is designed in isolation. While it must be emphasized that the link between power density and electricity cost can be complex [25], we are hopeful that this approach might enable marked decreases in the unit price of wind-generated electricity.

Acknowledgments

P. L.-F. gratefully acknowledges support from a Junior Research Fellowship from Churchill College, Cambridge. We are also grateful to Claudia Cenedese for sharing with us the data used to construct the entrainment parameterization, to Heidi Nepf for sharing data on canopy measurements, and to Antonio Segalini for helpful comments on a draft of this manuscript.