The Effect of Free-Atmosphere Stratification on Boundary-Layer Flow and Power Output from Very Large Wind Farms

Abstract

Large-eddy simulation is used to study the influence of free-atmosphere stratification on the structure of atmospheric boundary-layer flow inside and above very large wind farms, as well as the power extracted by the wind turbines. In the simulations, tuning-free Lagrangian scale-dependent dynamic models are used to model the subgrid-scale turbulent fluxes, while the turbine-induced forces are parameterized with an actuator-disk model. It is shown that for a given surface cover (with and without turbines) thermal stratification of the free atmosphere reduces the entrainment from the flow above compared with the unstratified case, leading to lower boundary-layer depth. Due to the fact that in very large wind farms vertical energy transport associated with turbulence is the only source of kinetic energy, lower entrainment leads to lower power production by the wind turbines. In particular, for the wind-turbine arrangements considered in the present work, the power output from the wind farms is reduced by about 35% when the potential temperature lapse rate in the free atmosphere increases from 1 to 10 K/km (within the range of values typically observed in the atmosphere). Moreover, it is shown that the presence of the turbines has significant effect on the growth of the boundary layer. Inspired by the obtained results, a simple one-dimensional model is developed to account for the effect of free-atmosphere stability on the mean flow and the power output from very large wind farms.

1. Introduction

Within a very large and sufficiently dense wind farm, the atmospheric boundary-layer (ABL) flow asymptotes to the fully-developed regime. In this regime, the performance of the wind turbines is not affected by the farm entrance region, and the power extraction from the wind farm is mainly due to the vertical energy transport from the flow above [1,2,3]. In other words, the kinetic energy must be entrained from the free atmosphere to balance that extracted by the farm, leading to an increased boundary-layer depth [4,5].

In the absence of wind turbines, entrainment and boundary-layer growth are controlled by different factors such as Earth’s rotation, surface momentum and buoyancy fluxes, and static stability of the free atmosphere immediately above the ABL. The shear at the surface and the positive surface buoyancy fluxes facilitate the growth of the boundary layer, while the Earth’s rotation, the negative surface buoyancy fluxes, and the static stability of the free atmosphere limit the boundary-layer growth [6]. Since the power extracted by the turbines in very large wind farms is directly linked to the entrainment of kinetic energy from the external environment, it is of great importance to investigate how the above mentioned parameters (buoyancy fluxes, free-atmosphere stability, and earth’s rotation) can affect the performance of wind farms.

Most of the previous numerical studies on the interaction between atmospheric boundary layer and wind farms, using computational fluid dynamics (CFD) [2,5] and one-dimensional models [1,3], have focused on the purely neutral ABL. In that case, the effect of surface buoyancy fluxes and free-atmosphere stratification can be ignored. However, it has been shown that the purely neutral ABLs are rarely observed in the atmosphere, and most of the real ABLs that are classified as neutral are almost always conventionally neutral [7]. Conventionally neutral ABLs are defined as neutrally-stratified ABLs capped by the stably-stratified free atmosphere. The main effect of the free-atmosphere stratification is to reduce the entrainment from the free atmosphere, leading to lower boundary-layer depth. The smaller boundary-layer depth limits the size of the largest eddies, which have a relatively important contribution to the turbulent kinetic energy and fluxes. As a result, the free-atmosphere stability limits the turbulent transport away from the surface compared with the unstratified case [8,9,10,11].

In order to isolate the effects of free-atmosphere stratification on the turbulent flow through wind farms and the turbines performance, in this study we focus on the interaction between conventionally neutral ABLs and very large wind farms. In this regard, a suit of large-eddy simulations (LES) of fully-developed wind-farm ABL flow is performed including the effect of earth’s rotation and free-atmosphere stability. Different values of the potential temperature lapse rate ($\Gamma =\frac{\partial \mathsf{\theta }}{\partial \text{z}}$) in the free-atmosphere within the range of 1 to 10 K/km, typically observed in the atmosphere [12], are considered. In the simulations, tuning-free Lagrangian scale-dependent dynamic models [13] are used to model the subgrid-scale fluxes, while the turbine-induced forces are parameterized with an actuator disk model [14]. It should be highlighted that considering the Coriolis forces in the governing equations allows investigating how the wind direction changes inside and above the farm due to the presence of the turbines, which is not possible in the unidirectional boundary layer flow resulting from an imposed pressure gradient. Moreover, the important effect of the wind farm on the boundary-layer growth can be explicitly resolved and studied. In Section 2, the LES framework used in this work is described. The influence of stably-stratified free atmosphere on the growth of the boundary layer and the surface-layer scaling in the absence of wind turbines is examined in Section 3. Afterwards, the results obtained from the LES of ABLs over very large wind farms are presented. In Section 4, our simulation results are also employed to develop a simple one-dimensional model to account for the effect of static stability of the free atmosphere on the power extracted by the turbines, the effective roughness experienced by the ABL and also the boundary-layer depth. Finally, a summary and conclusions are provided in Section 5.

2. Large-Eddy Simulation Framework

2.1. LES Governing Equations

LES solves the filtered continuity equation, the filtered Navier-Stokes equations (written here in rotational form and using the Boussinesq approximation), and the filtered transport equation for potential temperature:

$$\begin{gathered} \frac{\partial {\tilde{\text{u}}}_{\text{i}}}{\partial {\text{x}}_{\text{i}}}=0 \\ \frac{\partial {\tilde{\text{u}}}_{\text{i}}}{\partial \text{t}}+{\tilde{\text{u}}}_{\text{j}}\left(\frac{\partial {\tilde{\text{u}}}_{\text{i}}}{\partial {\text{x}}_{\text{j}}}-\frac{\partial {\tilde{\text{u}}}_{\text{j}}}{\partial {\text{x}}_{\text{i}}}\right)=-\frac{\partial {\tilde{\text{p}}}^{\text{*}}}{\partial {\text{x}}_{\text{i}}}-\frac{\partial {\mathsf{\tau }}_{\text{ij}}^{\text{d}}}{\partial {\text{x}}_{\text{j}}}+{\mathsf{\delta }}_{\text{i}3}\text{g}\frac{\tilde{\mathsf{\theta }}-\langle \tilde{\mathsf{\theta }}\rangle }{{\mathsf{\theta }}_0}+{\text{f}}_{\mathrm{c}}{\mathsf{\epsilon }}_{\text{ij}3}{\tilde{\text{u}}}_{\text{j}}+{ℱ}_{\text{i}} \\ \frac{\partial \tilde{\mathsf{\theta }}}{\partial \text{t}}+{\tilde{\text{u}}}_{\text{j}}\frac{\partial \tilde{\mathsf{\theta }}}{\partial {\text{x}}_{\text{j}}}=-\frac{\partial {\text{q}}_{\text{j}}}{\partial {\text{x}}_{\text{j}}} \end{gathered}$$

(1)

where the tilde represents a spatial filtering at scale $\widetilde{∆}$, $\text{t}$ is time; ${\tilde{\text{u}}}_{\text{i}}$ is instantaneous resolved velocity in the $\text{i}$-direction (with $\text{i}=1,2,3$ corresponding to the streamwise (x), spanwise (y) and vertical (z) direction, respectively); $\tilde{\mathsf{\theta }}$ denotes the resolved potential temperature; ${\mathsf{\theta }}_0$ is the reference temperature, the angle brackets represent a horizontal average; $\text{g}$ refers to the gravitational acceleration; ${\text{f}}_{\mathrm{c}}$ is the Coriolis parameter; ${\mathsf{\delta }}_{\text{ij}}$ is the Kronecker delta; ${\mathsf{\epsilon }}_{\text{ijk}}$ denotes the alternative unit tensor; ${\tilde{\text{p}}}^{\text{*}}=\tilde{\text{p}}/\mathsf{\rho }+\frac{1}{2}{\tilde{\text{u}}}_{\text{i}}{\tilde{\text{u}}}_{\text{i}}+\frac{1}{3}{\mathsf{\tau }}_{\text{kk}}$ is the modified kinematic pressure; ${ℱ}_{\text{i}}$ is a forcing term (e.g., wind-turbine induced forces); ${\text{q}}_{\text{j}}=\widetilde{{\text{u}}_{\text{j}}\mathsf{\theta }}-{\tilde{\text{u}}}_{\text{j}}\tilde{\mathsf{\theta }}$ denotes the SGS heat flux; ${\mathsf{\tau }}_{\text{ij}}=\widetilde{{\text{u}}_{\text{i}}{\text{u}}_{\text{j}}}-{\tilde{\text{u}}}_{\text{i}}{\tilde{\text{u}}}_{\text{j}}$ represents the kinematic SGS stress; and ${\mathsf{\tau }}_{\text{ij}}^{\text{d}}$ is its deviatoric part. Note that ${\mathsf{\tau }}_{\text{ij}}^{\text{d}}$ and ${\text{q}}_{\mathrm{j}}$ are unknown and need to be parameterized as a function of filtered (resolved) fields.

In LES, all turbulent structures larger than the filter scale ($\widetilde{∆}$) are resolved and the contribution of the unresolved eddies or small-scale structures on the resolved field is parameterized using a subgrid-scale (SGS) model. A common parameterization strategy in LES consists of computing the deviatoric part of the SGS stress with an eddy-viscosity model [15]:

$${\mathsf{\tau }}_{\text{ij}}^{\text{d}}={\mathsf{\tau }}_{\text{ij}}-\frac{1}{3}{\mathsf{\tau }}_{\text{kk}}{\mathsf{\delta }}_{\text{ij}}=-2{\widetilde{∆}}^2{\text{C}}_{\mathrm{S}}^2\left|\tilde{\text{S}}\right|{\tilde{\text{S}}}_{\text{ij}}$$

(2)

and the SGS heat flux with an eddy-diffusivity model:

$${\mathrm{q}}_{\text{j}}=-{\widetilde{∆}}^2{\mathrm{C}}_{\text{S}}^2{\text{Pr}}_{\text{sgs}}^{-1}\left|\tilde{\text{S}}\right|\frac{\partial \tilde{\mathsf{\theta }}}{\partial {\text{x}}_{\text{j}}}$$

(3)

where ${\tilde{\text{S}}}_{\text{ij}}=\left(\partial {\tilde{\text{u}}}_{\text{i}}/\partial {\text{x}}_{\text{j}}+\partial {\tilde{\text{u}}}_{\text{j}}/\partial {\text{x}}_{\text{i}}\right)/2$ is the resolved strain rate tensor and $\tilde{\text{S}}=\sqrt{2{\tilde{\text{S}}}_{\text{ij}}{\tilde{\text{S}}}_{\text{ij}}}$ is the strain rate magnitude. ${\text{C}}_{\mathrm{s}}$ represents the Smagorinsky coefficient and ${\text{C}}_{\mathrm{S}}^2{\text{Pr}}_{\text{sgs}}^{-1}$ is the lumped coefficient, where ${\text{Pr}}_{\text{sgs}}$ is the SGS Prandtl number. Here, we employ the scale-dependent Lagrangian dynamic models [13] to compute the local optimized value of the model coefficients without any ad hoc tuning. In contrast with the traditional dynamic models [16,17], the scale-dependent dynamic models compute dynamically not only the value of the model coefficients in the eddy-viscosity and eddy-diffusivity models, but also the dependence of these coefficients with scale. More details on the formulation of scale-dependent dynamic models for the SGS stress and the SGS scalar fluxes can be found in Porté-Agel et al. [18], Porté-Agel [19] and Stoll and Porté-Agel [13]. The scale-dependent Lagrangian dynamic model has been successfully applied to neutral homogeneous and heterogeneous ABLs [13,20], to the homogeneous stable boundary layer (SBL) [21], to boundary-layer flow over topography [22,23] and boundary-layer flow through the wind farm [4,24].

2.2. Wind-Turbine Parameterization

To parameterize the turbine-induced forces, the actuator-disk model with rotation [14,24] is used. Through this model, the lift and drag forces acting on the turbines are parameterized using the blade element momentum theory, and are distributed over the rotor area. Unlike the standard actuator-disk model, which assumes the loads are distributed uniformly over the rotor disk and acting only in the axial direction, the actuator-disk model with rotation includes the effect of turbine-induced flow rotation as well as the non-uniform force distribution. Figure 1 shows a cross-sectional element of radius r in the $\left(\mathsf{\theta },\text{x}\right)$ plane, where $\text{x}$ is the axial direction. Different forces, velocities and angles are shown in this figure. ${\text{V}}_{\mathrm{a}}={\text{V}}_{\mathrm{a}}(\text{r},\mathsf{\theta })$ and ${\text{V}}_{\mathsf{\theta }}={\text{V}}_{\mathsf{\theta }}(\text{r},\mathsf{\theta })$ are axial and tangential velocities of the incident flow at the blades, respectively, in the inertial frame of reference. The local velocity relative to the rotating blade is defined as ${\text{V}}_{\text{rel}}=({\text{V}}_{\mathrm{a}},{\text{V}}_{\mathsf{\theta }}-\mathsf{\Omega }\text{r})$, where $\mathsf{\Omega }$ is the turbine angular velocity. The angle of attack is defined as $\mathsf{\alpha }=\mathsf{\phi }-\mathsf{\gamma }$, where $\mathsf{\phi }={\tan }^{-1}(\frac{{\text{V}}_{\text{a}}}{\mathsf{\Omega }\text{r}-{\text{V}}_{\mathsf{\theta }}})$ is the angle between ${\text{V}}_{\text{rel}}$ and the rotor plane and $\mathsf{\gamma }$ is the local pitch angle. The resulting force is given by:

$${\text{f}}_{\text{disk}}=\frac{\text{dF}}{\text{dA}}=\frac{1}{2}{\mathsf{\rho }\text{V}}_{\text{rel}}^2\frac{\text{Bc}}{2\mathsf{\pi }\text{r}}\left({\text{C}}_{\text{L}}{\text{e}}_{\text{L}}+{\text{C}}_{\text{D}}{\text{e}}_{\text{D}}\right)$$

(4)

where an annular area of differential size is $\text{dA}=2\mathsf{\pi }\text{rdr}$; $\text{B}$ is the number of blades; ${\text{C}}_{\mathrm{L}}={\text{C}}_{\mathrm{L}}\left(\mathsf{\alpha },\text{Re}\right)$ and ${\text{C}}_{\mathrm{D}}={\text{C}}_{\mathrm{D}}\left(\mathsf{\alpha },\text{Re}\right)$ are lift and drag coefficient, respectively; $\text{c}$ is the chord length; and ${\text{e}}_{\mathrm{L}}$ and ${\text{e}}_{\mathrm{D}}$ denote the unit vectors in the directions of the lift and the drag, respectively. Through the simulation, the tangential (${\text{V}}_{\mathsf{\theta }}$) and the axial (${\text{V}}_{\mathrm{a}}$) velocities at the rotor plane are known. Hence, the local velocity relative to the rotating blade (${\text{V}}_{\text{rel}}$) and the angle of attack ($\mathsf{\alpha }=\mathsf{\phi }-\mathsf{\gamma }$) can be computed. Then, the resulting force is obtained using Equation (4). Porté-Agel et al. [25] provided details of the SWT-2.3-93 wind turbine in a simulation of an operational wind farm.

Figure 1. A cross-sectional aerofoil element.

Figure 1. A cross-sectional aerofoil element.

2.3. Numerical Setup

The LES code used in this study is a modified version of the one described by Albertson and Parlange [26], Porté-Agel et al. [18], Stoll and Porté-Agel [13] and Wu and Porté-Agel [14]. In the simulations, tuning-free Lagrangian scale-dependent dynamic models are used to model the subgrid-scale turbulent fluxes, and the turbine-induced forces are parameterized using the above-mentioned actuator-disk model. Since the Reynolds number of the ABL is very high, no near-ground viscous processes are resolved, and the viscous term is neglected in the momentum equation. The vertical derivatives are approximated with second-order central differences and the horizontal directions are discretized pseudo-spectrally. Periodic boundary conditions are applied horizontally; as a result, a very large (infinite-size) wind farm is simulated. Full dealiasing of the nonlinear terms is obtained by padding and truncation according to the 3/2 rule [27]. The time advancement is carried out using a second-order-accurate Adams–Bashforth scheme [28]. The finite-difference scheme in the vertical direction requires specification of boundary conditions at the top and bottom of the domain. The upper boundary condition consists of a stress-free condition. A Rayleigh damping region is also used following the GABLS case description [4,20] to limit gravity-wave reflection from the top of the domain. At the bottom surface, the standard formulation based on local application of Monin-Obukhov similarity theory [29] is used. Although this theory was developed for mean quantities, it is a common practice in LES of atmospheric flows to compute the instantaneous filtered surface momentum flux [30,31] as follows:

$${{\mathsf{\tau }}_{\text{i}3}|}_{\text{w}}=-{\text{u}}_{*}^2\frac{{\tilde{\text{u}}}_{\text{i}}}{{\tilde{\text{u}}}_{\text{r}}}=-{\left(\frac{{\tilde{\text{u}}}_{\text{r}}\mathsf{\kappa }}{\text{ln}\left(\frac{\text{z}}{{\text{z}}_{\mathrm{o}}}\right)}\right)}^2\frac{{\tilde{\text{u}}}_{\text{i}}}{{\tilde{\text{u}}}_{\text{r}}}(\text{i}=1,2)$$

where ${{\mathsf{\tau }}_{\text{i}3}|}_{\text{w}}$ is the instantaneous local wall stress; ${\text{u}}_{\mathrm{*}}$ is the friction velocity; ${\text{z}}_{\mathrm{o}}$ is the roughness length; $\mathsf{\kappa }$ is the von Kármán constant; and ${\tilde{\text{u}}}_{\text{r}}$ is the local filtered horizontal velocity at the first level ($\text{z}=\Delta \text{z}/2$). In this study, we focus on the conventionally neutral ABL and, therefore, the surface heat flux is set to zero.

The boundary layer is driven by an imposed uniform geostrophic wind (G) of 10 m/s; the Coriolis parameter is set to ${\text{f}}_{\mathrm{c}}=1\times {10}^{-4}$ rad/s. The initial mean vertical profiles are specified as linearly increasing potential temperature with a prescribed depth-constant temperature gradient $\Gamma$, and a constant wind velocity (equal to geostrophic wind magnitude) at all levels. The initial flow is laminar with imposed very small perturbations at the first 100 m from the surface. It should be mentioned the angle of the geostrophic wind is set so that, when the flow reaches the quasi-steady state, the mean velocity direction is aligned with the turbine axes at the hub-height level. The code is run for a long-enough time to guarantee that quasi-steady conditions are reached. The domain is divided into ${\text{N}}_{\mathrm{x}}$, ${\text{N}}_{\mathrm{y}}$, and ${\text{N}}_{\mathrm{z}}$ uniformly spaced grid points in streamwise, spanwise and wall-normal directions, respectively. The grid planes are staggered in the vertical direction, with the first vertical velocity plane at a distance ${∆}_{\text{z}}={\text{L}}_{\mathrm{z}}/\left({\text{N}}_{\text{z}}-1\right)$ from the surface, and the first horizontal velocity plane ${∆}_{\text{z}}/2$ from the surface. The filter width is computed using the common formulation $\widetilde{\mathrm{∆}}={\left({\mathrm{∆}}_{\text{x}}{\mathrm{∆}}_{\text{y}}{\mathrm{∆}}_{\text{z}}\right)}^{1/3}$, where ${\mathrm{∆}}_{\text{x}}={\mathrm{L}}_{\text{x}}/{\mathrm{N}}_{\text{x}}$ and ${\mathrm{∆}}_{\text{y}}={\text{L}}_{\mathrm{y}}/{\mathrm{N}}_{\text{y}}$. The Siemens SWT-2.3-93 wind turbines, with rotor diameter ($\text{D}$) of 93 m and a hub-height (${\text{z}}_{\mathrm{h}}$) of 80 m, are “immersed” in the flow. In order to investigate the layout effects, the framework is also applied to study several cases of aligned and staggered wind farms with different streamwise (${\mathrm{s}}_{\text{x}}$) and spanwise (${\mathrm{s}}_{\text{y}}$) spacings. The key parameters of the various LES cases are summarized in Table 1.

Table 1. Key parameters of the various LES cases.

Case	zo (m)	Abbreviation	Γ (K/km)	Number of wind turbines (Ntx × Nty)	sx × sy	Turbine arrangement	Lx × Ly × Lz (m3)	Nx × Ny × Nz
A1	0.1	−	1	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
A2	0.01	−	1	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
A3	0.1	−	10	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
A4	0.01	−	10	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
B1	0.1	s5 × 5 − Γ1	1	8 × 7	5 × 5	Staggered	3720 × 3255 × 1687	126 × 175 × 128
B2	0.1	a5 × 5 − Γ1	1	8 × 7	5 × 5	Aligned	3720 × 3255 × 1687	126 × 175 × 128
C1	0.1	s5 × 5 − Γ10	10	8 × 7	5 × 5	Staggered	3720 × 3255 × 1045	126 × 175 × 80
C2	0.1	a5 × 5 − Γ10	10	8 × 7	5 × 5	Aligned	3720 × 3255 × 1045	126 × 175 × 80
D1	0.1	s7 × 7 − Γ1	1	6 × 5	7 × 7	Staggered	3906 × 3255 × 1687	126 × 175 × 128
D2	0.1	a7 × 7 − Γ1	1	6 × 5	7 × 7	Aligned	3906 × 3255 × 1687	126 × 175 × 128
E1	0.1	s7 × 7 − Γ10	10	6 × 5	7 × 7	Staggered	3906 × 3255 × 1045	126 × 175 × 80
E2	0.1	a7 × 7 − Γ10	10	6 × 5	7 × 7	Aligned	3906 × 3255 × 1045	126 × 175 × 80

Table 1. Key parameters of the various LES cases.

Case	zo (m)	Abbreviation	Γ (K/km)	Number of wind turbines (Ntx × Nty)	sx × sy	Turbine arrangement	Lx × Ly × Lz (m3)	Nx × Ny × Nz
A1	0.1	−	1	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
A2	0.01	−	1	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
A3	0.1	−	10	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
A4	0.01	−	10	No farm	−	−	3906 × 3255 × 1045	126 × 175 × 80
B1	0.1	s5 × 5 − Γ1	1	8 × 7	5 × 5	Staggered	3720 × 3255 × 1687	126 × 175 × 128
B2	0.1	a5 × 5 − Γ1	1	8 × 7	5 × 5	Aligned	3720 × 3255 × 1687	126 × 175 × 128
C1	0.1	s5 × 5 − Γ10	10	8 × 7	5 × 5	Staggered	3720 × 3255 × 1045	126 × 175 × 80
C2	0.1	a5 × 5 − Γ10	10	8 × 7	5 × 5	Aligned	3720 × 3255 × 1045	126 × 175 × 80
D1	0.1	s7 × 7 − Γ1	1	6 × 5	7 × 7	Staggered	3906 × 3255 × 1687	126 × 175 × 128
D2	0.1	a7 × 7 − Γ1	1	6 × 5	7 × 7	Aligned	3906 × 3255 × 1687	126 × 175 × 128
E1	0.1	s7 × 7 − Γ10	10	6 × 5	7 × 7	Staggered	3906 × 3255 × 1045	126 × 175 × 80
E2	0.1	a7 × 7 − Γ10	10	6 × 5	7 × 7	Aligned	3906 × 3255 × 1045	126 × 175 × 80

3. LES Results

3.1. No-Farm Case

Figure 2a,b and 2d show the vertical profiles of the horizontally-averaged velocity magnitude (in linear and semi-log scales) and potential temperature, respectively, for two different values of $\Gamma$ and two different land-surface roughnesses ${\text{z}}_{\mathrm{o}}$. As expected, for a given surface cover, the stronger stratification (i.e., larger value for $\Gamma$) leads to lower entrainment from the free atmosphere and, consequently, shallower boundary-layer depth. It is also evident that, for a given $\Gamma$, increasing the surface roughness leads to an increase in the ABL depth. It is due to the fact that, by increasing the surface roughness, the shear production near the surface increases, and as a result, turbulent transport away from the surface is enhanced. In addition, as displayed in Figure 2c, for a given ${\text{z}}_{\mathrm{o}}$, the wind direction (arctangent of the ratio of the ageostrophic to the geostrophic velocity components) inside the ABL increases when the free-atmosphere stability increases.

Figure 2. Vertical profiles of horizontally-averaged velocity magnitude in M (a) linear scale; and (b) semi-log scale; (c) wind direction; and (d) horizontally-averaged potential temperature (Θ) inside the ABL for two different values of Γ and ${\text{z}}_{\mathrm{o}}$. The horizontal dotted lines show the top-tip and bottom-tip heights.

Figure 2. Vertical profiles of horizontally-averaged velocity magnitude in M (a) linear scale; and (b) semi-log scale; (c) wind direction; and (d) horizontally-averaged potential temperature (Θ) inside the ABL for two different values of Γ and ${\text{z}}_{\mathrm{o}}$. The horizontal dotted lines show the top-tip and bottom-tip heights.

The vertical profiles of the total shear stress are presented in Figure 3. In the absence of turbines, for a given ${\text{z}}_{\mathrm{o}}$, increasing the free-atmosphere stability reduces the surface momentum flux. This reduction of the surface shear stress is mainly due to the fact that the smaller boundary-layer depth limits the size of the largest eddies, which have an important contribution to the turbulent kinetic energy and fluxes [9]. It is also clear that, for a given value of $\Gamma$, the surface shear stress increases by increasing the surface aerodynamic roughness, leading to larger boundary-layer depth.

Figure 3. Vertical profile of total shear stress $\left(\sqrt{{\mathsf{\tau }}_{\text{xz}}^2+{\mathsf{\tau }}_{\text{yz}}^2}\right)$ for two different values of $\Gamma$ and ${\text{z}}_{\mathrm{o}}$.

Figure 3. Vertical profile of total shear stress $\left(\sqrt{{\mathsf{\tau }}_{\text{xz}}^2+{\mathsf{\tau }}_{\text{yz}}^2}\right)$ for two different values of $\Gamma$ and ${\text{z}}_{\mathrm{o}}$.

3.1.1. The Effect of Free-Atmosphere Stability on the ABL Depth

Zilitinkevich et al. [7] showed that static stability of the free atmosphere is a key parameter that needs to be taken into account to estimate the height of the boundary layer. They showed that the stably-stratified free atmosphere can reduce the height of the boundary layer more than 10 times compared with the purely neutral (Ekman) condition. For the purely neural condition, the height of the boundary layer is estimated using the classical Rossby and Montgomery [32] formula as follows:

$${\mathsf{\delta }}_{\text{bl}}={\text{C}}_{\mathrm{R}}\frac{{\text{u}}_{\text{*}}}{\left|\text{f}\right|}$$

(5)

where ${\text{u}}_{\mathrm{*}}$ is the friction velocity. As mentioned before, purely neutral ABLs are rarely observed over the land and never over the ocean, where the ABL is always conventionally neutral [33]. For the conventionally neutral ABL, Zilitinkevich and Esau [8] generalized Rossby and Montgomery’s expression for the height of the boundary layer by including Brunt-Vӓisӓlӓ frequency ($\text{N}=\sqrt{\frac{\text{g}}{{\mathsf{\theta }}_0}\frac{\partial \mathsf{\theta }}{\partial \text{z}} }$) in the free atmosphere as an additional parameter as follows:

$${\mathsf{\delta }}_{\text{bl}}={\text{C}}_{\mathrm{R}}{\left(1+{\text{C}}_{\text{N}}\frac{\text{N}}{\left|\text{f}\right|}\right)}^{-\frac{1}{2}}\frac{{\text{u}}_{\text{*}}}{\text{f}}={\text{C}}_{\mathrm{R}}^{\prime }\frac{{\text{u}}_{\text{*}}}{\text{f}}$$

(6)

where ${\text{C}}_{\mathrm{R}}$ and ${\text{C}}_{\mathrm{N}}$ are empirical dimensionless constants. The values of ${\text{C}}_{\mathrm{R}}$ and ${\text{C}}_{\mathrm{N}}$ depend upon the definition of the boundary-layer height. One possible definition for the height of the boundary layer is the height at which the total momentum flux reaches 5% of the surface value (i.e., ${\left|\mathsf{\tau }\right|}_{\text{z }={\mathsf{\delta }}_{\text{bl}}}=0.05\times {\text{u}}_{\mathrm{*}}^2$). Based upon this definition, Zilitenkevich and Esau [33] found ${\text{C}}_{\mathrm{R}}\approx 0.5$ and ${\text{C}}_{\mathrm{N}}~0.05-0.2$ from the atmospheric data. Our LES results, displayed in Figure 3, give ${\text{C}}_{\mathrm{R}}\approx 0.5$ and ${\text{C}}_{\mathrm{N}}\approx 0.11$, which are in good agreement with the above estimations. Figure 4 shows the comparison of ${\mathsf{\delta }}_{\text{bl}}$ computed using Equation (6) with the one directly obtained from LES. From that figure, it is clear that Equation (6) provides a reasonable estimation for the height of the boundary layer. It is worth mentioning that the value of Brunt-Väisälä frequency is $5.8\times {10}^{-3}{\text{ s}}^{-1}$ for $\Gamma =1\text{ K}/\text{km}$ and $1.8\times {10}^{-2}{\text{ s}}^{-1}$ for $\Gamma =10\text{ K}/\text{km}$.

Figure 4. Correlation between the ABL heights: ${\mathsf{\delta }}_{\text{bl}}(\text{theory})$ calculated from Equation (6) (with ${\text{C}}_{\mathrm{R}}=0.5$ and ${\text{C}}_{\mathrm{N}}=0.11$) and ${\mathsf{\delta }}_{\text{bl}}(\text{LES})$ directly obtained from LES.

Figure 4. Correlation between the ABL heights: ${\mathsf{\delta }}_{\text{bl}}(\text{theory})$ calculated from Equation (6) (with ${\text{C}}_{\mathrm{R}}=0.5$ and ${\text{C}}_{\mathrm{N}}=0.11$) and ${\mathsf{\delta }}_{\text{bl}}(\text{LES})$ directly obtained from LES.

Another possible definition for the ABL height (denoted here as ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}$) is the lowest elevation at which the horizontally-averaged velocity equals the geostrophic value. In that case, the value of ${\text{C}}_{\mathrm{R}}\approx 0.16$ has been reported for the purely neutral ABL [2]. In conventionally neutral conditions, the LES results shown in Figure 2a give ${\text{C}}_{\mathrm{R}}\approx 0.16$ and ${\text{C}}_{\mathrm{N}}\approx 0.02$. Figure 5 illustrates the comparison of ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}$ using Equation (6) (with ${\text{C}}_{\mathrm{R}}\approx 0.16$ and ${\text{C}}_{\mathrm{N}}\approx 0.02$) with the one directly obtained from LES. It is also evident that Equation (6) provides a good estimation of the boundary-layer height ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}$. As will be explained in Section 4, estimating the ABL height based upon the velocity profile is preferable as it allows to develop a simple analytical model to predict the interaction of the ABL and the wind farms while accounting for the effect of free-atmosphere stability.

Figure 5. Correlation between ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}(\text{theory})$ calculated from Equation (6) (with ${\text{C}}_{\mathrm{R}}=0.16$ and ${\text{C}}_{\mathrm{N}}=0.02$) and ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}(\text{LES})$ directly obtained from LES.

Figure 5. Correlation between ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}(\text{theory})$ calculated from Equation (6) (with ${\text{C}}_{\mathrm{R}}=0.16$ and ${\text{C}}_{\mathrm{N}}=0.02$) and ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}(\text{LES})$ directly obtained from LES.

3.1.2. The Effect of Free-Atmosphere Stability on the Surface-Layer Parameterization

There are different methods to take into account the non-local effect of the free-atmosphere stratification on the surface-layer parameterization [34,35,36]. For instance, Blackadar [34] proposed a model based upon redefinition of the mixing length in the log-law. Zilitenkevich et al. [37] proposed another method which does not need to change the mixing length. Their model is based on adding an additional non-local term to the log-law in order to take into account the effect of free-atmosphere stability on the surface-layer scaling. Esau [9] performed LES to evaluate the performance of different surface-layer parameterizations to predict the surface drag coefficient, and showed that Zilitenkevich’s method was the most accurate parameterization in the conventionally neutral condition. That model is as follows:

$$\frac{\text{dU}}{\text{dz}}=\frac{{\text{u}}_{\text{*}}}{\mathsf{\kappa }\text{z}}+{\text{a}}_{\mathrm{u}}\text{N}\to \text{U}\left(\text{z}\right)=\frac{{\text{u}}_{\mathrm{*}}}{\mathsf{\kappa }}\ln \left(\frac{\text{z}}{{\text{z}}_{\text{o}}}\right)+{\text{a}}_{\mathrm{u}}\text{Nz}$$

(7)

where ${\text{a}}_{\mathrm{u}}$ is an empirical constant. We found that ${\text{a}}_{\mathrm{u}}=0.3$ provides the best fit of this model to the LES velocity profiles. The same value of ${\text{a}}_{\mathrm{u}}$ has also been reported by Esau [9]. It should be mentioned that this method asymptotes to the standard log-law when $\text{N}\to 0$ for purely neutral condition.

3.2. Wind-Farm Case

Figure 6 presents the vertical profiles of velocity magnitude, potential temperature and wind direction in a very large wind farm for two different wind-turbine spacings (${\text{s}}_{\mathrm{x}}={\text{s}}_{\mathrm{y}}=5$ and ${\text{s}}_{\mathrm{x}}={\text{s}}_{\mathrm{y}}=7$) and two different values for $\Gamma$. It is evident that the presence of the turbines has significant effect on the growth of the boundary layer. It is also observed from Figure 6c that the ratio of the ageostrophic to geostrophic velocity component increases due to the presence of the wind turbines. In other words, the wind direction changes inside the wind farm due to the presence of the turbines. Recently, a similar trend was reported by Johnstone and Coleman [5]. They performed numerical simulation of the neutral turbulent Ekman layer over an infinite wind farm, and showed that the Ekman spiral becomes more pronounced when the wind turbines are present. It is important to note that the Reynolds number used in the mentioned work was extremely small in comparison to that of the real atmospheric boundary layer.

In a very large wind farm, as observed in Figure 7 and also reported in previous studies [2,5,24], the surface shear stress decreases due to extraction of momentum by the turbines compared with the no-farm case. Besides, the total shear stress has a peak at the turbine-top level, where the strong wind shear occurs, and has a higher value compared with the surface shear stress in the absence of turbines.

It should be mentioned that the previous studies were limited to the purely neutral condition, thus ignoring the influence of free-atmospheric stratification. As presented in Figure 7, applying the plane-averaged momentum balance between the top-tip and the bottom-tip heights shows that the extracted momentum by the farm increases when the static stability of the free atmosphere decreases. It can be concluded that more kinetic energy can be entrained from the flow above when the potential temperature lapse rate decreases, leading to higher power production by the farm.

Figure 6. Vertical profiles of horizontally-averaged velocity magnitude in (a) linear scale; and (b) semi-log scale; (c) wind direction; and (d) horizontally-averaged potential temperature (Θ) through very large wind farms for two different wind-turbine spacings, and two different values of Γ.

Figure 6. Vertical profiles of horizontally-averaged velocity magnitude in (a) linear scale; and (b) semi-log scale; (c) wind direction; and (d) horizontally-averaged potential temperature (Θ) through very large wind farms for two different wind-turbine spacings, and two different values of Γ.

Figure 7. Vertical profile of total shear stress $\left(\sqrt{{\mathsf{\tau }}_{\text{xz}}^2+{\mathsf{\tau }}_{\text{yz}}^2}\right)$ for two different wind-turbine spacings, and two different values of $\Gamma$ through very large wind farms.

Figure 7. Vertical profile of total shear stress $\left(\sqrt{{\mathsf{\tau }}_{\text{xz}}^2+{\mathsf{\tau }}_{\text{yz}}^2}\right)$ for two different wind-turbine spacings, and two different values of $\Gamma$ through very large wind farms.

It is interesting to mention that the above-mentioned formula used to estimate the height of the ABL in the absence of wind turbines can also be used to predict the ABL depth for a very large wind farm. Figure 8 shows the comparison of ${\mathsf{\delta }}_{\text{bl}}$ and ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}$ using Equation (6) written for ${\text{u}}_{*}={\text{u}}_{\text{*}2}$ (i.e., friction velocity at the turbine-top level) with the data obtained from LES. It should be noted that in the case of the wind farm, the top-tip height needs to be added to the value obtained from Equation (6) (i.e., ${\mathsf{\delta }}_{\text{bl}}={\text{C}}_{\mathrm{R}}^{\prime }\frac{{\text{u}}_{\text{*}2}}{\text{f}}+{\text{z}}_{\mathrm{h}}+\frac{\text{D}}{2}$). This is because the magnitude of the momentum flux (shear stress), and thus the production of TKE, has a peak at the turbine-top level, and it decreases almost linearly above that height.

Figure 8. Correlation between ${\mathsf{\delta }}_{\text{bl}}$ and ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}$ theory and their counterparts directly obtained from LES.

Figure 8. Correlation between ${\mathsf{\delta }}_{\text{bl}}$ and ${\mathsf{\delta }}_{\text{bl}}^{\text{*}}$ theory and their counterparts directly obtained from LES.

3.3. Layout Effect

In this study, we also investigate the effect of wind farm layout (aligned vs. staggered) on the flow and the power extracted by the farm. Figure 9 shows vertical profiles of the horizontally-averaged velocity and the total shear stress. Note that, above the wind-turbine region, the profiles are almost identical for different layouts. The main effect of configuration is found below the turbine region, where the staggered wind farms extract more momentum compared with the aligned ones. Using the plane-averaged momentum balance between the top-tip and the bottom-tip heights, one can show that the momentum and, consequently, the power extracted by the farm is higher when the layout changes from aligned to staggered.

In Table 2 we summarize the time-averaged power output from the wind turbines for the different cases. As stated before and also presented in this table, the power extracted by the turbines is higher for the staggered wind farms compared with the aligned ones. In particular, for the wind farms considered here, the power output from the staggered wind farms is about 10% higher than the aligned ones. It should be noted that the power output is defined as $\text{P}={{\displaystyle \int }}_0^{2\mathsf{\pi }}{{\displaystyle \int }}_{{\text{r}}_{\text{n}}}^{\text{D}/2}{\text{r}\mathsf{\Omega }\mathsf{\delta }\text{F}}_{\mathsf{\theta }}\left(\mathsf{\theta },\text{r}\right)\text{rd}\mathsf{\theta }\text{dr}$, where ${\mathsf{\delta }\text{F}}_{\mathsf{\theta }}$ is the local tangential force per unit area of the disk and ${\text{r}}_{\mathrm{n}}$ is the nacelle radius.

Figure 9. Vertical profiles of (a) horizontally-averaged velocity magnitude M; and (b) total shear stress for two different values of and two different layouts.

Figure 9. Vertical profiles of (a) horizontally-averaged velocity magnitude M; and (b) total shear stress for two different values of and two different layouts.

Table 2. The time-averaged power output for different cases.

Case	Power/Turbine (MW)
s5 × 5 − Γ1	0.3069
a5 × 5 − Γ1	0.2835
s5 × 5 − Γ10	0.1995
a5 × 5 − Γ10	0.1841
s7 × 7 − Γ1	0.4303
a7 × 7 − Γ1	0.3811
s7 × 7 − Γ10	0.2993
a7 × 7 − Γ10	0.2690

Table 2. The time-averaged power output for different cases.

Case	Power/Turbine (MW)
s5 × 5 − Γ1	0.3069
a5 × 5 − Γ1	0.2835
s5 × 5 − Γ10	0.1995
a5 × 5 − Γ10	0.1841
s7 × 7 − Γ1	0.4303
a7 × 7 − Γ1	0.3811
s7 × 7 − Γ10	0.2993
a7 × 7 − Γ10	0.2690

Figure 10 displays the contour plots of mean and instantaneous streamwise velocity (at the hub-height level) for the aligned and staggered layout. As shown in this figure, in aligned farms, the presence of high-speed “channels” between the turbine rows and low-speed regions along the turbine axes is clear. In contrast, in the staggered farm, the flow is more uniform in the spanwise direction and no high-speed channels are clearly visible. Similar flow pattern has been shown by Meyer and Meneveau [38], Markfort et al. [39] and Wu and Porté-Agel [24]. In addition, in the staggered farms, the turbine wakes have a longer distance to recover before the next downwind turbine, which results in a higher velocity at the turbines and, consequently, larger power output compared with the aligned counterparts. For this reason, although the power extracted by the staggered farms is higher compared with the aligned ones, the horizontally-averaged velocity is lower for the staggered cases (see Figure 9a).

Figure 10. Contours of mean (a) and instantaneous (b) streamwise velocity (m/s) at the hub-height level for the different cases. Only a section of the domain is shown.

Figure 10. Contours of mean (a) and instantaneous (b) streamwise velocity (m/s) at the hub-height level for the different cases. Only a section of the domain is shown.

4. One-Dimensional (1D) Modeling of Very Large Wind Farms in Conventionally Neutral Conditions

One of the first one-dimensional models to study the interactions of very large wind farms with the turbulent atmospheric boundary layer was proposed by Frandsen [1]. He formulated a model for the surface roughness induced by very large wind farms. Calaf et al. [2] improved Frandsen’s model by including the effects of turbine wake mixing. Yang et al. [40] refined Calaf et al.’s surface-roughness model by taking into account the effects of streamwise and spanwise turbine spacings in infinite aligned wind farms. Meneveau [3] and Meyers and Meneveau [41] combined Calaf et al.’s surface-roughness model with the outer-layer defect law proposed by Tennekes and Lumley [42] to predict the velocity and, consequently, the power in very large wind farms in purely neutral conditions. It is important to mention that in all the above-mentioned models the effect of free-atmosphere stratification is ignored. In this section, we follow the methodology of Frandsen [1] to derive a simple one-dimensional model to predict the mean flow and the power output from very large wind farms under conventionally neutral conditions. Following Frandsen’s approach, the horizontally-averaged velocity is divided into two equilibrium layers, one below and the other above the hub height (Figure 11). They are characterized by the local friction velocities ${\text{u}}_{\text{*}1}$ and ${\text{u}}_{\text{*}2}$, respectively. For the layer below the hub height, one can write:

$$\frac{{\text{dU}}_1}{\text{dz}}=\frac{{\text{u}}_{\text{*}1}}{\mathsf{\kappa }\text{z}}+{\text{a}}_{\text{u}}\text{N}\to {\text{U}}_1(\text{z})=\frac{{\text{u}}_{\text{*}1}}{\mathsf{\kappa }}\ln \left(\frac{\text{z}}{{\text{z}}_{\text{o}1}}\right)+{\text{a}}_{\text{u}}\text{Nz }$$

(8)

and for the layer above the hub height:

$$\frac{{\text{dU}}_2}{\text{dz}}=\frac{{\text{u}}_{\text{*}2}}{\mathsf{\kappa }\text{z}}+{\text{a}}_{\text{u}}\text{N}\to {\text{U}}_2(\text{z})=\frac{{\text{u}}_{\text{*}2}}{\mathsf{\kappa }}\ln \left(\frac{\text{z}}{{\text{z}}_{\text{o}2}}\right)+{\text{a}}_{\text{u}}\text{Nz }$$

(9)

where $\mathsf{\kappa }=0.4$ is von Karman constant; and ${\text{z}}_{\text{o}1}$ and ${\text{z}}_{\text{o}2}$ are the aerodynamic surface roughness and wind-farm-induced roughness, respectively. Assuming that the horizontally-averaged velocity is continuous at the hub height level, leads to:

$$\frac{{\text{u}}_{\text{*}1}}{\mathsf{\kappa }}\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}1}}\right)+{\text{a}}_{\text{u}}{\text{Nz}}_{\text{h}}=\frac{{\text{u}}_{\text{*}2}}{\mathsf{\kappa }}\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}2}}\right)+{\text{a}}_{\text{u}}{\text{Nz}}_{\text{h}}⟹{\text{u}}_{\text{*}1}={\text{u}}_{\text{*}2}\frac{\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}2}}\right)}{\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}1}}\right)}$$

(10)

Using the plane-averaged momentum balance above and below the turbines, one can write:

$${\text{u}}_{\text{*}2}^2={\text{u}}_{\text{*}1}^2+\frac{1}{2}{\mathrm{C}}_{\text{T}}\left(\frac{\mathsf{\pi }}{4}{\text{D}}^2\right){\text{u}}_{\text{h}}^2\left(\frac{1}{{\text{s}}_{\text{x}}{\text{s}}_{\text{y}}{\text{D}}^2}\right)$$

(11)

where ${\text{s}}_{\text{x}}\text{D}$ and ${\text{s}}_{\text{y}}\text{D}$ are streamwise and spanwise spacing of the wind turbines, respectively; ${\mathrm{C}}_{\text{T}}$ is the turbine’s trust coefficient; and ${\text{u}}_{\text{h}}$ is horizontally-averaged velocity at the hub-height level. In a very large wind farm, it is common to define the trust coefficient of the farm according to ${\text{c}}_{\text{ft}}=\frac{{\mathsf{\pi }\text{C}}_{\text{T}}}{4{\text{s}}_{\mathrm{x}}{\text{s}}_{\mathrm{y}}}$. Thus, the equation above can be rewritten as follows:

$${\text{u}}_{\text{*}2}^2={\text{u}}_{\text{*}1}^2+\frac{1}{2}{\text{c}}_{\text{ft}}{\text{u}}_{\text{h}}^2$$

(12)

Substituting Equation (10) in Equation (12) gives:

$${\text{u}}_{\text{*}2}^2={\left({\text{u}}_{\text{*}2}\frac{\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}2}}\right)}{\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}1}}\right)}\right)}^2+\frac{1}{2}{\text{c}}_{\text{ft}}{\left(\frac{{\text{u}}_{\text{*}2}}{\mathsf{\kappa }}\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}2}}\right)+{\text{a}}_{\text{u}}{\text{Nz}}_{\text{h}}\right)}^2$$

(13)

For the purely neutral condition ($\text{N}=0)$, Equation (13) can be solved for ${\text{z}}_{\text{o}2}$ and one can get Frandsen’s formula for the effective roughness of a very large wind farm as follows:

$${\text{z}}_{\text{o}2,\text{Frandsen}}={\text{z}}_{\text{h}}\text{exp}(-\frac{\mathsf{\kappa }}{\sqrt{\frac{1}{2}{\text{c}}_{\text{ft}}+{\left[\frac{\mathsf{\kappa }}{\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}1}}\right)}\right]}^2}}).$$

(14)

Unlike the purely neutral condition, in the conventionally neutral case the problem is not closed due to the presence of ${\text{a}}_{\text{u}}{\text{Nz}}_{\text{h}}$ term in Equation (13). Hence, it is not possible to derive an explicit equation for wind farm effective roughness. We can rewrite Equation (12) by using Equation (8) as follows:

$${\text{u}}_{\text{*}2}^2={\left[\frac{\mathsf{\kappa }}{\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}1}}\right)}\left({\text{u}}_{\text{h}}-{\text{a}}_{\text{u}}{\text{Nz}}_{\text{h}}\right)\right]}^2+\frac{1}{2}{\text{c}}_{\text{ft}}{\text{u}}_{\text{h}}^2$$

(15)

To close the problem, we can rewrite Equation (9) by integrating this equation from the hub height to the height where the velocity equals the geostrophic value (i.e., $\text{z}={\mathsf{\delta }}_{\text{bl}}^{\text{*}}={\text{C}}_{\text{R}}^{\prime }\frac{{\text{u}}_{\text{*}2}}{\left|\text{f}\right|}+{\text{z}}_{\text{h}}+\frac{\text{D}}{2}$) as follows:

$$\text{G}-{\text{u}}_{\text{h}}=\frac{{\text{u}}_{\text{*}2}}{\mathsf{\kappa }}\ln \left(\frac{{\mathsf{\delta }}_{\text{bl}}^{\text{*}}}{{\text{z}}_{\text{h}}}\right)+{\text{a}}_{\text{u}}\text{N}\left({\mathsf{\delta }}_{\text{bl}}^{\text{*}}-{\text{z}}_{\mathrm{h}}\right)$$

(16)

By substituting Equation (6) in Equation (16) we have:

$$\text{G}-{\text{u}}_{\text{h}}=\frac{{\text{u}}_{\text{*}2}}{\mathsf{\kappa }}\ln \left(\frac{{\text{C}}_{\mathrm{R}}^{\prime }\frac{{\text{u}}_{\text{*}2}}{\left|\text{f}\right|}+{\text{z}}_{\text{h}}+\frac{\text{D}}{2}}{{\text{z}}_{\text{h}}}\right)+{\text{a}}_{\text{u}}\text{N}\left({\text{C}}_{\mathrm{R}}^{\prime }\frac{{\text{u}}_{\text{*}2}}{\left|\text{f}\right|}+\frac{\text{D}}{2}\right)$$

(17)

Finally, combining Equations (15) and (17), we can write the following two equations:

$${\text{u}}_{\text{*}2}={\left[{\left(\frac{\mathsf{\kappa }}{\ln \left(\frac{{\text{z}}_{\text{h}}}{{\text{z}}_{\text{o}1}}\right)}\left({\text{u}}_{\text{h}}-{\text{a}}_{\text{u}}{\text{Nz}}_{\text{h}}\right)\right)}^2+\frac{1}{2}{\text{c}}_{\text{ft}}{\text{u}}_{\text{h}}^2\right]}^{\frac{1}{2}}$$

(18)

and:

$${\text{u}}_{\text{h}}=\text{G}-\frac{{\text{u}}_{\text{*}2}}{\mathsf{\kappa }}\ln \left(\frac{{\text{C}}_{\mathrm{R}}^{\prime }\frac{{\text{u}}_{\text{*}2}}{\left|\text{f}\right|}+{\text{z}}_{\text{h}}+\frac{\text{D}}{2}}{{\text{z}}_{\text{h}}}\right)-{\text{a}}_{\text{u}}\text{N}\left({\text{C}}_{\mathrm{R}}^{\prime }\frac{{\text{u}}_{\text{*}2}}{\left|\text{f}\right|}+\frac{\text{D}}{2}\right)$$

(19)

By solving Equations (18) and (19) numerically, we can compute the values for ${\text{u}}_{\text{*}2}$ and ${\text{u}}_{\text{h}}$. Then, by using Equations (8) and (9), one can compute ${\text{u}}_{\text{*}1}$ and ${\text{z}}_{\text{o}2}$. It is important to note that when $\text{N}\to 0$ the obtained value for ${\text{z}}_{\text{o}2}$ is identically equal to the value obtained from Frandsen’s formula [Equation (14)].

Figure 11. A schematic of wind profile inside a very large wind farm.

Figure 11. A schematic of wind profile inside a very large wind farm.

Next, the proposed model is used to evaluate the effect of free-atmosphere stability on the extracted power by the turbines, the boundary-layer growth and the effective roughness of the wind farm for various wind turbine spacings $\text{s}$ (where $\text{s}=\sqrt{{\text{s}}_{\mathrm{x}}{\text{s}}_{\mathrm{y}}}$) and underlying surface roughnesses ${\text{z}}_{\text{o}}$. In order to evaluate the performance of the proposed model, we use the same values for geostrophic velocity ($\text{G}$), Coriolis parameter ($\text{f}$), rotor diameter ($\text{D})$, and hub height (${\text{z}}_{\text{h}}$) as those used in the numerical simulations described in Section 3. To estimate the thrust coefficient (${\mathrm{C}}_{\text{T}}$) of the turbines, first we calculate another thrust coefficient (denoted as ${\mathrm{C}}_{\text{T}}^{\prime }$) based on the average velocity at the rotor (${\text{u}}_{\text{d}}$):

$${\mathrm{C}}_{\text{T}}^{\prime }=\frac{\text{T}}{\frac{1}{2}{\mathsf{\rho }\text{u}}_{\text{d}}^2{\text{A}}_{\mathrm{d}}}$$

(19)

where T is the total trust over the rotor area obtained from the drag and lift forces computed in the actuator disk model; and ${\text{A}}_{\text{d}}$ is the area of the disk. Then, using the one-dimensional momentum theory, it is possible to relate ${\mathrm{C}}_{\text{T}}^{\prime }$ to the axial induction factor ($a$) as follows:

$${\mathrm{C}}_{\text{T}}^{\prime }=\frac{4a}{1-a}$$

(20)

We can compute the axial induction factor ($a$) using the equation above. Then, the trust coefficient (${\mathrm{C}}_{\text{T}}$) used in the 1D model is computed using the one-dimensional momentum theory as follows:

$${\mathrm{C}}_{\text{T}}=4a\left(1-a\right)$$

(21)

The values of $a$, ${\mathrm{C}}_{\text{T}}^{\prime }$ (directly obtained from LES) and ${\mathrm{C}}_{\text{T}}$ [using Equation (21)] are presented in Table 3.

Table 3. The values of $a$, ${\mathrm{C}}_{\text{T}}^{\prime }$, and ${\mathrm{C}}_{\text{T}}$ for the different cases.

Case	a	$C_T^{\prime }$	$C_T$
$\text{s}5\times 5-\Gamma 1$	$0.197$	$0.98$	$0.63$
$\text{a}5\times 5-\Gamma 1$	$0.199$	$1.00$	$0.64$
$\text{s}5\times 5-\Gamma 10$	$0.203$	$1.02$	$0.64$
$\text{a}5\times 5-\Gamma 10$	$0.207$	$1.04$	$0.65$
$\text{s}7\times 7-\Gamma 1$	$0.195$	$0.97$	$0.63$
$\text{a}7\times 7-\Gamma 1$	$0.195$	$0.97$	$0.63$
$\text{s}7\times 7-\Gamma 10$	$0.192$	$0.95$	$0.62$
$\text{a}7\times 7-\Gamma 10$	$0.195$	$0.97$	$0.63$

Table 3. The values of $a$, ${\mathrm{C}}_{\text{T}}^{\prime }$, and ${\mathrm{C}}_{\text{T}}$ for the different cases.

Case	a	$C_T^{\prime }$	$C_T$
$\text{s}5\times 5-\Gamma 1$	$0.197$	$0.98$	$0.63$
$\text{a}5\times 5-\Gamma 1$	$0.199$	$1.00$	$0.64$
$\text{s}5\times 5-\Gamma 10$	$0.203$	$1.02$	$0.64$
$\text{a}5\times 5-\Gamma 10$	$0.207$	$1.04$	$0.65$
$\text{s}7\times 7-\Gamma 1$	$0.195$	$0.97$	$0.63$
$\text{a}7\times 7-\Gamma 1$	$0.195$	$0.97$	$0.63$
$\text{s}7\times 7-\Gamma 10$	$0.192$	$0.95$	$0.62$
$\text{a}7\times 7-\Gamma 10$	$0.195$	$0.97$	$0.63$

Figure 12 shows the effect of free-atmosphere stability on the power extracted by the turbines, for different wind turbine spacings. The power output directly obtained from LES is also shown for comparison. It should be mentioned that, in the 1D model, the power output is defined as $\text{P}=\frac{1}{2}{\mathsf{\rho }\text{C}}_{\text{p}}{\text{u}}_{\text{h}}^3{\text{A}}_{\text{d}}$, where ${\text{C}}_{\text{p}}$ is the power coefficient which is defined as ${\text{C}}_{\text{p}}=4a{\left(1-a\right)}^2$. The results show that the power extracted by each turbine decreases by increasing the stability of the free atmosphere. It is also observed that the magnitude of the differential change ($\text{dP}/\text{d}\Gamma$) is relatively larger for low values of $\Gamma$ and it becomes smaller for higher potential temperature lapse rates.

As mentioned before, in a fully developed wind farm, the power extraction is connected to the kinetic energy entrainment from the flow above. This effect leads to an increase in the height of the boundary layer. Figure 13 shows the effect of free-atmosphere stability on the height of the boundary layer using Equation (6). It is clear that, by increasing the stability of the free atmosphere, the height of the boundary layer decreases. It means that less energy can be entrained from the flow above and, consequently, less power can be extracted. It is also clear that the wind farms with the highest density of wind turbines have the strongest impact on the ABL growth. Consistent with the above result, the rate of change (${\text{d}\mathsf{\delta }}_{\mathrm{bl}}/\text{d}\Gamma$) is higher for lower values of $\Gamma$ and it decreases when the potential temperature lapse rate increases.

Figure 12. Effect of free-atmosphere stability on the power extracted by the turbines.

Figure 12. Effect of free-atmosphere stability on the power extracted by the turbines.

Figure 13. Influence of free-atmosphere stability on the height of the boundary layer.

Figure 13. Influence of free-atmosphere stability on the height of the boundary layer.

As explained in Section 3, in the absence of turbines, the height of the boundary layer is a function of surface shear stress (${\text{u}}_{\text{*}}$) and the external stability parameter (${\mathsf{\mu }}_{\text{N}}=\frac{\text{N}}{\text{f}}$). Inside a very large wind farm, the maximum shear stress happens at the turbine-top level and has a dominant effect on the height of the boundary layer. Figure 14 shows the effect of the free-atmosphere stability on the friction velocity at the top level of the farm. It is observed that ${\text{u}}_{\text{*}2}$ decreases with increasing the stability of the free atmosphere.

The influence of free-atmosphere stability on the effective roughness (${\text{z}}_{\text{o}2}$) is plotted in Figure 15. As shown in this figure, free-atmosphere stratification increases the wind farm-induced surface roughness. The impact of the free-atmosphere stability on the effective roughness is higher when the density of the wind farm increases.

Figure 14. Influence of free-atmosphere stability on the friction velocity at the turbine-top level.

Figure 14. Influence of free-atmosphere stability on the friction velocity at the turbine-top level.

Figure 15. Influence of free-atmosphere stability on wind farm effective roughness.

Figure 15. Influence of free-atmosphere stability on wind farm effective roughness.

The effect of the underlying surface roughness on the power output from the turbines is provided in Figure 16 for values of $\Gamma$ from 0 to 20 K/km. The power extraction from the wind farm increases with decreasing underlying surface roughness. A similar trend has been reported recently by Meneveau [3] for a purely neutral ABL. As shown in this figure, when the stability of the free atmosphere increases, a smaller amount of energy can be extracted by the turbines.

Figure 16. Effect of underlying surface roughness and free-atmosphere stability on the power output for $\text{s}=7.$

Figure 16. Effect of underlying surface roughness and free-atmosphere stability on the power output for $\text{s}=7.$

5. Conclusions

The present work focuses on the influence of thermal stratification of the free atmosphere on the structure of ABL flow, as well as the power extracted by the wind turbines. A suite of large-eddy simulations of fully-developed wind-farm ABL flow has been performed including the effect of earth’s rotation and free-atmosphere stability. In the simulations, tuning-free Lagrangian scale-dependent dynamic models [13] are used to model the subgrid-scale turbulent fluxes, while the turbine-induced forces are parameterized with the actuator disk model with rotation [14].

The large-eddy simulations performed in this study demonstrate that the presence of the turbines has a significant effect on the boundary-layer depth. In very large wind farms, the vertical energy transport associated with turbulence is the only source of kinetic energy. As a result, the kinetic energy must be entrained from the external environment, leading to an increase in the ABL height. On the other hand, the results presented here indicate that the thermal stratification of the free atmosphere reduces the growth of the boundary-layer height, leading to lower kinetic energy entrainment from the flow above and, consequently, lower power production by the turbines. In particular, for the wind-turbine arrangements considered in the present study, the power output from the wind farms is reduced by about 35% when the potential temperature lapse rate in the free atmosphere increases from 1 to 10 K/km.

Inspired by the LES results, we have developed a simple one-dimensional model to include the effect of free-atmosphere stability on the prediction of the mean flow and the power output from the wind farm. The proposed model is a modified version of the model initially developed by Frandsen [1]. It should be noted that Frandsen’s model does not account for thermal stratification in the free-atmosphere and, thus, it is only valid for the purely neutral condition. In contrast, the proposed new model is able to capture the effect of thermal stratification of the free atmosphere. The proposed model was used to systematically investigate the effect of free-atmosphere stability, aerodynamic surface roughness and turbine density on the boundary-layer depth and the power extracted by the turbines.

Future research will focus on developing a similar model to include the effect of surface buoyancy fluxes in stable and convective boundary layers on the power production in large wind farms. Furthermore, the LES framework will also be used to investigate the effect of free-atmosphere stability on all the terms in the budgets of mean and turbulent kinetic energy inside and above infinite as well as finite-size wind farms.