Analytical model of fully developed wind farms in conventionally neutral atmospheric boundary layers

Abstract

The wind energy industry relies on computationally efficient engineering-type models to design wind farms. Typically these models do not account for the effect of atmospheric stratification in either the boundary layer or the free atmosphere. This study proposes a new analytical model for fully developed wind-turbine arrays in conventionally neutral atmospheric boundary layers frequently encountered in nature. The model captures the effect of the free-atmosphere stratification, Coriolis force, wind farm layout and turbine operating condition on the wind farm performance. The model is developed based on the physical insight derived from large-eddy simulations. We demonstrate that the geostrophic drag law (GDL) for flow over flat terrain can be extended to flow over fully developed wind farm arrays. The presence of a vast wind farm significantly increases the wind farm friction velocity compared with flow over flat terrain, which is modelled by updated coefficients in the GDL. The developed model reliably captures the vertical wind speed profile inside the wind farm. Furthermore, the power production trends observed in simulations are reliably reproduced. The wind farm performance, normalized by the geostrophic wind speed, decreases as the free-atmosphere thermal stability increases or the Coriolis force decreases. In addition, we find that the optimal wind farm performance is obtained at a lower thrust coefficient than the Betz limit, which indicates that optimal operating conditions for turbines in a wind farm are different than for a single turbine.

1. Introduction

Extended wind farms lead to the formation of a fully developed wind-turbine array boundary layer in which the power extraction from the wind farm originates from vertical energy transport from the flow aloft (Stevens & Meneveau 2017; Porté-Agel, Bastankhah & Shamsoddin 2020). As wind-turbine arrays are covering increasingly larger areas, the fully developed wind farm flow has been extensively studied using wind-tunnel experiments (Chamorro, Arndt & Sotiropoulos 2011; Bossuyt, Meneveau & Meyers 2017; Segalini & Chericoni 2021), large-eddy simulations (LES, Calaf, Meneveau & Meyers 2010; Lu & Porté-Agel 2011; Goit & Meyers 2015) and analytical models (Frandsen 1992; Calaf et al. 2010; Abkar & Porté-Agel 2013). These studies have deepened our fundamental understanding of the turbulent energy transport in extended wind farms (Meneveau 2019), which is relevant to accurately representing the impact of large wind farms in climate model simulations (Keith et al. 2004; Wang & Prinn 2010; Li et al. 2018) and models that capture regional meteorology (Baidya-Roy, Pacala & Walko 2004). Note that two main types of wind farm parametrizations are used in weather and climate models (Porté-Agel et al. 2020): In the first type, the wind farm is considered as an increased surface roughness length; and in the second type, the wind farm is considered as an elevated sink of momentum. While the former was used early both in climate and mesoscale models, the latter is becoming more popular in mesoscale models.

The most popular analytical method to describe fully developed wind-turbine arrays is the top-down approach, which models the impact of a wind farm on the atmospheric flow above by determining the effective surface roughness height that corresponds to the wind farm. It has been inspired by the impact of plant canopies on atmospheric flows and has been considered to model wind farm dynamics since the pioneering studies of Templin (1974) and Newman (1977). In particular, Frandsen (1992) obtained an explicit formulation for effective surface roughness by assuming that the vertical profile of the horizontally averaged wind speed can be split into two logarithmic layers, one below and one above the turbine hub height. The model relates the corresponding friction velocities by balancing the momentum fluxes and assuming that the wind speed is continuous at hub height. The existence of the two logarithmic layers, with a wake layer in between, has been confirmed by the wind farm simulation of Calaf et al. (2010). They also introduced the concept of wake eddy viscosity to represent the wake layer and showed that it gives improved predictions compared with the model of Frandsen (1992). In the models of Frandsen (1992) and Calaf et al. (2010), the horizontally averaged wind speed at hub height is considered to be representative of the velocity at the turbine locations. However, the appropriateness of this assumption depends on the wind farm layout. To account for the local variations in the wind speed at hub height, Yang, Kang & Sotiropoulos (2012) introduced the inhomogeneity factor $\beta$, which is the ratio of the velocity at the turbine location over the horizontally averaged velocity. Using one-dimensional mass and momentum conservation, they determined $\beta$ by approximating the mean streamwise velocity distributions at hub height. Stevens, Gayme & Meneveau (2015, 2016b) introduced an ‘effective wake area coverage’ $w_f$ in their generalized coupled wake boundary-layer model, which is equivalent to setting $\beta ^{2} = {1/w_{f}}$. More recently, Zhang et al. (2021) used a revised model of Jensen (1983) to capture the effect of the turbine positioning on the velocity distribution in the wind farm and the model of Calaf et al. (2010) to simulate the effect of the vertical kinetic energy flux on the wind farm performance. These two models were coupled such that $\beta$ can be solved iteratively. Note that these models are all explicitly based on the logarithmic law. In contrast, Nishino & Dunstan (2020) derived a two-scale momentum theory directly from the law of momentum conservation without using the logarithmic law, so that the theory may account for the effect of large-scale motions of the atmospheric boundary layer (ABL) in a time-dependent (rather than statistical) manner.

The top-down models are derived from truly neutral ABLs driven by a mean pressure gradient. However, ABLs are seldom truly neutral but often capped by a stably stratified free atmosphere, which is known as conventionally neutral ABLs (Zilitinkevich & Esau 2002). The mean pressure gradient is balanced by the Coriolis force associated with the geostrophic wind when the flow is in geostrophic equilibrium. This leads to the geostrophic drag law (GDL) for flat terrain, which was originally derived by Rossby & Montgomery (1935) from a turbulence closure model and later by Blackadar & Tennekes (1968) from a more general asymptotic similarity. For extended wind farms, the equivalent wind farm surface roughness height is commonly used in conjuncture with the GDL to model the effect of the geostrophic wind forcing (Frandsen et al. 2006; Calaf et al. 2010; Meyers & Meneveau 2012),

(1.1)\begin{equation} \left( \frac{\kappa G}{u_{*,2}} \right)^{2} = \left[ \ln \left( \frac{u_{*,2}}{|f| z_{0,2}} \right) - A \right]^{2} + B^{2}, \end{equation}

where $\kappa =0.4$ is the von Kármán constant, $G$ is the geostrophic wind speed, $u_{*,2}$ is the wind farm friction velocity, $z_{0,2}$ is the equivalent wind farm surface roughness, $f=2\varOmega \sin \phi$ is the Coriolis parameter, $\varOmega$ is the angular speed of the Earth and $\phi$ is the latitude. Typical values for the GDL coefficients $A$ and $B$ reported from meteorological observations are $A=1.3\unicode{x2013}1.9$ and $B=4.4\unicode{x2013}4.9$ (Hess & Garratt 2002). It has been shown that, for conventionally neutral ABLs, the coefficients $A$ and $B$ only depend on the Zilitinkevich number $Zi = N/f$ (Esau 2004; Kadantsev, Mortikov & Zilitinkevich 2021; Liu, Gadde & Stevens 2021a,b; Liu & Stevens 2022). Here $N=\sqrt {g \varGamma /\theta _0}$ is the free-atmosphere Brunt–Väisälä frequency, $g$ is the gravity acceleration, $\varGamma$ is the free-atmosphere lapse rate and $\theta _0$ is the reference potential temperature. However, it is not a priori clear whether the similarity properties of the GDL apply to atmospheric flow above extended wind farms and, thus, whether a functional dependence that only depends on the Zilitinkevich number can describe the GDL coefficients $A$ and $B$.

To account for the effect of free-atmosphere stratification, Abkar & Porté-Agel (2013) modified the model of Frandsen (1992) by adding a term $a_u N z$ to the logarithmic laws below and above hub height, where $a_u\approx 0.3$ is an empirical constant and $z$ is the altitude. The model system was closed by assuming the wind speed profile to be valid up to the boundary-layer height, where the wind speed equals the geostrophic wind speed. In agreement with their LES results, the model predicts that the power extracted by each turbine decreases with increasing stability of the free atmosphere. However, this model can not distinguish different wind farm layouts nor accurately capture the wind speed profile as the additional term $a_u N z$ will increase the wind speed at the hub height while the turbines impose an additional drag at this height. Also Allaerts & Meyers (2015) proposed a simple analytical model to obtain the wind farm power output for the fully developed regime. With the assumption that the wind farm friction velocity is much larger than the surface friction velocity, the authors derived an approximate relation between the wind farm power output $P$ and the geostrophic drag coefficient $u_{*,2}/G$. The latter was determined by further assuming that the GDL of (1.1) was still valid above extended wind farms. Their model system is closed by estimating the equivalent wind farm surface roughness height $z_{0,2}$ from previous models (e.g. Frandsen 1992; Calaf et al. 2010) and determining explicit expressions for the GDL coefficients $A$ and $B$ with some specific eddy viscosity profiles (e.g. Csanady 1974; Nieuwstadt 1983). The resultant model is thus a function of the ABL height and the Rossby number. However, as remarked previously, the validity of the GDL above extended wind farms is not a priori clear and needs validation. In addition, as pointed out by the authors, the model of Allaerts & Meyers (2015) is not valid when the wind farm thrust coefficient becomes very small, as in this case the wind farm friction velocity is no longer much larger than the surface friction velocity.

In the present work, we use the LES of wind farms to establish that insight from the GDL applies to fully developed wind-turbine array boundary layers in conventionally neutral ABLs. After providing parametrization of the GDL coefficients $A$ and $B$ for extended wind farms, we combine this updated GDL in conjuncture with an adjusted version of the model of Frandsen (1992) to model the wind farm production in conventionally neutral ABLs. Using this model, we investigate the effects of free-atmosphere lapse rate and turbine thrust coefficient on the wind farm power output and compare the results with LES data. The good agreement between the model prediction and LES data confirms the validity of the present model.

The organization of the paper is as follows. In § 2 the simulation method is explained. In § 3 the GDL for extended wind farms is parameterized based on the LES data. In § 4 we construct an analytical solution of the fully developed wind farm in conventionally neutral ABLs. In §§ 5 and 6 we investigate the effects of free-atmosphere lapse rate and turbine thrust coefficient, respectively. Finally, § 7 concludes the paper.

2. Large-eddy simulations of wind farms

We use an LES to simulate the conventionally neutral ABL flow over an infinite wind farm. In LES all turbulent structures larger than the filter scale are resolved, and the contributions of all smaller scales are parameterized using subgrid scale models. We integrate the spatially filtered governing equations (Albertson 1996; Liu & Stevens 2021a)

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} = 0, \end{gather}$$

(2.2)$$\begin{gather}\partial_t \tilde{\boldsymbol{u}} + \tilde{\boldsymbol{\omega}} \times \tilde{\boldsymbol{u}} = \boldsymbol{f}_{wt} + f \boldsymbol{e}_z \times (\boldsymbol{G} - \tilde{\boldsymbol{u}}) + g ( \tilde \theta / \langle \tilde \theta \rangle - 1) \boldsymbol{e}_z - \boldsymbol{\nabla} \tilde{p} / \rho - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}, \end{gather}$$

(2.3)$$\begin{gather}\partial_t \tilde{\theta} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\theta} ={-} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{q}}, \end{gather}$$

where $\tilde {{\cdot }}$ denotes the spatial filtering, $\langle \cdot \rangle$ represents the horizontal averaging, $\boldsymbol{e}_z$ is the unit vector in the vertical direction, $\tilde{\boldsymbol {u}}$ is the velocity, $\tilde {\boldsymbol {\omega } }= \boldsymbol {\nabla } \times \tilde {\boldsymbol {u}}$ is the vorticity, $\tilde p$ is the modified pressure departure from equilibrium, $\rho$ is the fluid density, $\boldsymbol {f}_{wt}$ is the force exerted by the wind turbines, $\tilde \theta$ is the potential temperature, $g$ is the gravity acceleration, $\boldsymbol {G}$ is the geostrophic wind velocity, $\boldsymbol {\tau }$ is the deviatoric part of the subgrid scale shear stress $\widetilde {\boldsymbol {u} \boldsymbol {u}} - \tilde {\boldsymbol {u}} \, \tilde {\boldsymbol {u}}$ and $\boldsymbol {q} = \widetilde {\boldsymbol {u} \theta } - \tilde {\boldsymbol {u}} \tilde \theta$ is the subgrid scale heat flux. Viscous terms in the governing equations are neglected as the Reynolds number is very high. The subgrid scale shear stress and heat flux are parameterized by the anisotropic minimum dissipation model (Abkar & Moin 2017; Gadde, Stieren & Stevens 2021).

The computational grid is uniform in the horizontal and vertical directions and staggered in the wall-normal direction. The first vertical velocity grid plane is located at the ground, while the first horizontal velocities grid plane is located at half a vertical grid distance above the ground. We use a pseudo-spectral discretization, and, thus, periodic boundary conditions, in the horizontal directions and a second-order finite difference method in the vertical direction. At the top boundary, we enforce a zero vertical velocity and zero shear stress. At the ground, we employ the classic wall model based on the logarithmic law of the wall (Moeng 1984; Bou-Zeid, Meneveau & Parlange 2005). Time integration is performed using a second-order Adams–Bashforth method, and the projection method is used to ensure that the velocity field is divergence free. A proportional-integral-derivative controller is used such that the wind at hub height is always directed in the same way to ensure the wind farm layout is not changed due to the change of flow conditions in all simulations (Sescu & Meneveau 2014; Liu & Stevens 2021a). In addition, a Rayleigh damping layer is applied in the top 25 % of the domain to reduce the effects of wind farm induced gravity waves (Klemp & Lilly 1978; Gadde & Stevens 2021; Liu & Stevens 2021a).

The wind turbines are modelled using an actuator disk approach (Calaf et al. 2010) in which the undisturbed upstream velocity $U_\infty$ is used to calculate the turbine force $F_{wt}$,

(2.4)\begin{equation} F_{wt} ={-} \frac{1}{2}\rho C_T U_\infty^{2} \frac{\rm \pi}{4} D^{2}, \end{equation}

where $C_T$ is the thrust coefficient based on $U_\infty$ and $D$ is the turbine diameter. However, when a turbine operates in the wake of upstream turbines $U_\infty$ is not readily available. Therefore, $F_{wt}$ is actually calculated as (Calaf et al. 2010)

(2.5)\begin{equation} F_{wt} ={-} \frac{1}{2}\rho C_T' U_d^{2} \frac{\rm \pi}{4} D^{2}, \end{equation}

where $U_d$ is the disk-averaged velocity and $C_T'$ is the thrust coefficient based on $U_d$. From the momentum theory (e.g. Burton et al. 2001) it follows that $C_T=4a(1-a)$ and $C_T'=4a/(1-a)$, where $a$ is the axial induction factor. To avoid numerical instability, the turbine force is commonly filtered via a Gaussian smoothing kernel and this can result in a small overprediction of the wind-turbine production (Munters & Meyers 2017). It has been shown by Shapiro, Gayme & Meneveau (2019) that this overprediction can be prevented by correcting the disk-averaged velocity $U_d$ in (2.5) by $U_d^{corr} = M U_d$, where $1/M=1+(C_T' \varDelta ) / (2 \sqrt {3{\rm \pi} } D)$ is the inverse correction factor and $\varDelta$ is the filter width. Thus, this correction method is always adopted in our simulations.

2.1. Conventionally neutral ABL simulations

The simulations are performed using an LES framework that has been validated extensively for flow in conventionally neutral ABLs (Liu et al. 2021a,b) and wind farms (Stevens, Martínez-Tossas & Meneveau 2018; Liu & Stevens 2020, 2021b). Two sets of simulations are performed for aligned and staggered wind farm layouts. In the first set, we keep the thrust coefficient $C_T'=4/3$ constant to study the effect of the free-atmosphere lapse rate $\varGamma$ and latitude $\phi$ (table 1). In the second set, we keep the free-atmosphere lapse rate $\varGamma =4$ K km$^{-1}$ constant to study the effect of the turbine thrust coefficient $C_T'$ and latitude $\phi$ (table 2). The computational domain size is $6\ {\rm km} \times 3\ {\rm km} \times 3\ {\rm km}$ in the streamwise, spanwise and vertical directions, and the corresponding grid resolution is $256\times 256\times 384$. Here streamwise and spanwise are relative to the mean wind direction at the hub height. We consider wind turbines with a rotor diameter $D=100$ m and hub height $z_h=100$ m. The turbine spacing normalized by the turbine diameter in the streamwise and spanwise directions is $s_x=s_y=6$, and, hence, the number of turbines is $50$ in each simulation. The initial potential temperature profile is $\theta (z) = \theta _0 + \varGamma z$ (Abkar & Porté-Agel 2013), where $\theta _0 = 300$ K is the potential temperature at the surface, $z$ is the altitude and $\varGamma = 0.2\sim 10$ K km$^{-1}$. The initial velocity field is $\boldsymbol {u} = G \boldsymbol {e}_x$, where $G=12$ m s$^{-1}$ is the geostrophic wind speed and $\boldsymbol {e}_x$ is the unit vector in the streamwise direction. To spin-up turbulence, small random perturbations are added to the initial fields of $\boldsymbol {u}$ and $\theta$ near the surface ($z \le 100$ m). The surface roughness length is $z_{0,1}= 10^{-4}$ m, which is representative for offshore conditions, and the latitude is $\phi =30^{\circ }, 50^{\circ }$ and $80^{\circ }$. The statistics are averaged for the dimensionless time $ft\in [10, 10+2{\rm \pi} ]$, where the flow has reached its quasi-equilibrium state. In this state the mean potential temperature is nearly constant up to the capping inversion height (Abkar & Porté-Agel 2013). In the inversion layer there is a strong potential temperature gradient, which recovers to the free atmospheric lapse rate above in the free atmosphere. In total, 66 high-resolution simulations are performed and statistics are summarized in tables 1 and 2.

Table 1. Summary of the first set of simulations with thrust coefficient $C_T'=4/3$.

Table 2. Summary of the second set of simulations with lapse rate $\varGamma =4$ K km$^{-1}$.

2.2. Boundary-layer height

Figure 1 shows the normalized boundary-layer height $u_{*,2}^{2}/(fh)^{2}$ as a function of the $Zi$ number, which reveals clearly that the boundary-layer height decreases as the stability of the free atmosphere increases (see also table 1). Here $u_{*,2}=(u_{*,1}^{2}+0.5 c_{ft}' U_d^{2})^{1/2}$ is the wind farm friction velocity, with $u_{*,1}$ being the surface friction velocity and $c_{ft}'={\rm \pi} C_T'/(4 s_x s_y)$ the wind farm thrust coefficient, and $h$ is the boundary-layer height, which is determined as the height where the total shear stress $\tau =0.05 u_{*,2}^{2}$ (Abkar & Porté-Agel 2013). The expression that describes well the LES data is $u_{*,2}^{2}/(fh)^{2}=Zi/C_N^{2}$, with $C_N=1.61$ being an empirical constant. This is very similar to the situation over flat terrain (Zilitinkevich, Esau & Baklanov 2007; Zilitinkevich et al. 2012), indicating that the similarity as defined in the GDL applies to fully developed wind-turbine array boundary layers. However, we note that the boundary-layer height above a wind farm is significantly higher than over flat terrain as the friction velocity $u_{*,2}$ in the former is much higher than in the latter. In addition, we note that the inversion height $z_i$, which is determined as the height where the total heat flux reaches its minimum value (Abkar & Moin 2017), can also be described well by $u_{*,2}^{2}/(fz_i )^{2}=Zi/C_N'^{2}$, where the empirical constant $C_N'=1.45< C_N=1.61$. Thus, the inversion height $z_i$ is generally shallower than the ABL height $h$, and their ratio is approximately $z_i/h \approx 0.9$.

Figure 1. The $Zi$-dependence of the inverse square of the dimensionless boundary-layer height $u_{*,2}^{2}/(fh)^{2}$. Filled symbols: LES data of table 1 for aligned (circles) and staggered (diamonds) wind farms; open symbols: the corresponding LES data of table 2; solid line: model prediction of $u_{*,2}^{2}/(fh)^{2}=Zi/C_N^{2}$ with $C_N=1.61$.

Figure 2 shows the velocity distribution at hub height for aligned and staggered wind farms normalized by the geostrophic wind speed. It is demonstrated that the wind velocity in this plane is non-uniform due to the wind-turbine wakes. High-velocity wind speed channels are formed between the turbines for the aligned wind farm configuration. At the same time, for the staggered layout, the wakes have a relatively long distance to recover. As a result, the velocity at the turbine location $U_d$ compared with the horizontally averaged velocity at hub height $U_h$ is relatively high for the staggered wind farms and relatively low for the aligned wind farms. Because the power output $P$ is proportional to $U_d^{3} = (1-a)^{3} \beta ^{3} U_h^{3}$ rather than $U_h^{3}$, it is important to model $\beta$ to account for the local velocity variations compared with the horizontal mean.

Figure 2. Temporally averaged streamwise velocity at hub height normalized by the geostrophic wind speed for (a,c,e) aligned and (b,d, f) staggered wind farms. The free-atmosphere lapse rate is $\varGamma =4$ K km$^{-1}$ and the latitude is (a,b) $\phi =30^{\circ }$, (c,d) $\phi =50^{\circ }$ and (e, f) $\phi =80^{\circ }$. The short lines denote the location of turbines.

3. Geostrophic drag law for extended wind farms

Figure 3 shows the $Zi$-dependence of the GDL coefficients

(3.1a,b)\begin{equation} A = \ln{\left(\frac{u_{*,2}}{f z_{0,2}}\right)} - \frac{\kappa G}{u_{*,2}} \cos \alpha_0, \quad B = \frac{\kappa G}{u_{*,2}} \sin \alpha_0, \end{equation}

where $z_{0,2}$ is determined by $z_{0,2}=2z_h \exp (-\kappa U_{2h}/u_{*,2})$, with $U_{2h}$ being the horizontally averaged wind speed at $z=2z_h$ (Calaf et al. 2010; Zhang et al. 2021), and $\alpha _0$ is the total wind angle change across the entire boundary layer (see table 1). It is demonstrated that $A$ and $B$ are well parametrized by a function that depends only on $Zi$, similar to what is observed in the classical GDL for flow over flat terrain (Zilitinkevich & Esau 2005; Liu et al. 2021a). This confirms that the similarity arguments apply for flow over extended wind farm arrays, even though their physics are different. Because the friction velocity $u_{*,2}$ is significantly increased compared with the flat terrain case, the value of $B$ is lower, and the value of $A$ is higher for the wind farm case than for the traditional flat terrain case. Here we parameterize $A$ and $B$ for fully developed wind-turbine arrays in conventionally neutral ABLs using a least-square fitting to the data as follows:

(3.2a,b)\begin{equation} A=1.54+0.18 \ln {Zi}, \quad B=1.74+0.011 Zi. \end{equation}

The corresponding results are shown in figure 3, which confirm that this model indeed captures the LES data well.

Figure 3. The $Zi$-dependence of the (a) $A$ and (b) $B$ coefficients in the GDL. Filled symbols: LES data of table 1 for aligned (circles) and staggered (diamonds) wind farms; open symbols: the corresponding LES data of table 2; solid line: the least-square fitting using (3.2a,b); dashed line: coefficients over flat terrain (Liu et al. 2021a).

It is worth pointing out that a complete description of $A$ and $B$ in extended wind farms involves many parameters, such as the hub height $z_h$, rotor diameter $D$, turbine spacing $(s_x,s_y)$, surface roughness height $z_{0,1}$, geostrophic wind speed $G$, turbine thrust coefficient, wind farm layout, latitude and free-atmosphere lapse rate. For simplicity, we fixed $(z_h,D,s_x,s_y,z_{0,1},G)$ in our simulations. As shown in figure 3, the data collapse well for different wind farm layouts and turbine thrust coefficients, and only depend on the Zilitinkevich number. However, these results cannot guarantee the values of $A$ and $B$ are unchanged when $(z_h,D,s_x,s_y,z_{0,1},G)$ are varied, which should be kept in mind for future studies.

4. A new top-down model for conventionally neutral ABLs

Similar to the models of Frandsen (1992) and Abkar & Porté-Agel (2013), we assume that the vertical profiles of the horizontally averaged wind speed $U$ is characterized by two logarithmic velocity regions, i.e. one below hub height and one above. Figure 4 shows a sketch of the corresponding wind speed profile, which is modelled as

(4.1) \begin{equation} U = \left\{\begin{array}{@{}ll} \displaystyle\dfrac{u_{*, 1}}{\kappa} \ln \left( \dfrac{z}{z_{0, 1}} \right) - a_u \beta^{2} u_{*,1} \dfrac{z}{z_h}, & z \le z_h,\\ \displaystyle\dfrac{u_{*, 2}}{\kappa} \ln \left( \dfrac{z}{z_{0, 2}} \right), & z \ge z_h. \end{array}\right. \end{equation}

Here $a_u=a_u(C_T')$ is an empirical constant that can be determined from LES data by setting $z=z_h$ in the first formula of (4.1). For the cases considered here, we find that $a_u \approx 4.3$ for $C_T'=4/3$; see figures 5(a) and 8(a). The factor $\beta =U_d/[(1-a)U_h]$, which accounts for the flow inhomogeneity at hub height, is incorporated into the first line of (4.1). To determine $\beta$ analytically, Zhang et al. (2021) introduced a new coupled model that combines a revised model of Jensen (1983) and the top-down model of Calaf et al. (2010). Here we adopted the model of Zhang et al. (2021) and added ghost turbines (Lissaman 1979; Ge et al. 2021) to further account for the ground effects. We find that $\beta =0.973$ for an aligned wind farm and $\beta =1.102$ for a staggered wind farm. Figure 5(b) shows that these values agree within $4\,\%$ of the values determined from LES. There may be various reasons for the observed small differences. First of all, the model of Zhang et al. (2021) neglects the effects of the Coriolis force and thermal stratification. Furthermore, the performance of the model may depend on the parameters such as $s_x,s_y,z_{0,1}$, the wind farm layout, the employed wake superposition method and the boundary-layer height parametrization.

Figure 4. A sketch showing the different parameters used in the present model.

Figure 5. The $Zi$-dependence of (a) the empirical constant $a_u$ and (b) the factor $\beta$ for aligned (circles) and staggered (diamonds) wind farms. Solid line: (a) $a_u=4.3$ and (b) analytical predictions for aligned ($\beta =0.973$) and staggered ($\beta =1.102$) wind farms, see details in text.

The wind speed profile should be continuous at hub height, and from (4.1) we thus obtain

(4.2)\begin{equation} U_h = \frac{u_{*, 1}}{\kappa} \ln \left( \frac{z_h}{z_{0, 1}} \right) - a_u \beta^{2} u_{*,1} =\frac{u_{*, 2}}{\kappa} \ln \left( \frac{z_h}{z_{0, 2}} \right). \end{equation}

Furthermore, the surface friction velocity $u_{*,1}$, the wind farm friction velocity $u_{*,2}$ and the hub height wind speed $U_h$ can be related by the horizontally averaged momentum balance above and below the turbines. When neglecting the Coriolis effect in the turbine wake region, the relation reads

(4.3a,b)\begin{equation} u_{*, 2}^{2} = u_{*, 1}^{2} + \frac{1}{2} c_{ft} \beta^{2} U_h^{2}, \quad c_{ft} = \frac{{\rm \pi} C_T}{4 s_x s_y}. \end{equation}

In addition, $u_{*,2}$ and $z_{0,2}$ are related by the GDL (i.e. (1.1)) with $A$ and $B$ being parameterized by (3.2a,b). This means that the system of equations ((1.1), (4.2) and (4.3a,b)) can be solved for the three unknowns $(u_{*, 1}, u_{*, 2}, z_{0, 2})$.

Figure 6 compared the vertical profiles of the horizontally averaged velocity obtained from LES with the above corresponding profiles obtained from the model (i.e. (4.1)). Very close to the wall, the LES data agrees well with the logarithmic law that originates from the wall. However, around hub height the influence of the wake effects, which are modelled by the additional term in (4.1), becomes pronounced. The present model also predicts the wind speed profile above hub height very well. The only exception is the region around the inversion layer, where the effects of the Coriolis force and thermal stratification result in the formation of a low-level jet.

Figure 6. The normalized horizontally averaged wind speed $U/G$ profile for aligned (a,c,e) and staggered (b,d, f) wind farms. The free-atmosphere lapse rate $\varGamma =4$ K km$^{-1}$ and the latitude (a,b) $\phi =30^{\circ }$, (c,d) $\phi =50^{\circ }$ and (e, f) $\phi =80^{\circ }$. Solid line: prediction of the present wind farm model (4.1); dashed-dotted line: the logarithmic law for the flat terrain very close to the ground. The vertical dashed lines denote the location of turbine rotor area.

5. Effect of free-atmosphere lapse rate

Figure 7(a,b) compares the $\varGamma$-dependence of the normalized hub height wind speed $U_h/G$ obtained from LES with the model prediction. The LES data shows that the normalized wind speed at hub height increases as the free-atmosphere lapse rate $\varGamma$ decreases or the latitude $\phi$ increases (see also table 1). In addition, it is revealed that the wind farm layout has a significant effect as $U_h$ for the aligned layout is about $6\sim 8\,\%$ higher than for the staggered layout. The reason is that staggered wind farms extract more energy from the flow and this increased drag results in a lower normalized velocity at hub height. Panels (c,d) show the corresponding normalized total wind farm power output per surface area $10^{3} P/G^{3}$, where $P= 0.5 c_{ft}' U_d^{3}$ is the wind farm power output per surface area, $c_{ft}'={\rm \pi} C_T'/(4 s_x s_y)$ is the normalized wind farm thrust coefficient and $U_d=(1-a)\beta U_h$ is the velocity at the turbine location. These panels show that the predicted power output agrees well with the LES data. As expected, the power output is slightly higher in the staggered wind farms than in the aligned ones.

Figure 7. The $\varGamma$-dependence of the (a,b) normalized wind speed at hub height $U_h/G$ and the corresponding (c,d) normalized wind farm power output per surface area for (a,c) aligned and (b,d) staggered wind farm. Lines: present model; filled symbols: LES data.

It should be pointed out that the differences in wind farm performance between aligned and staggered layouts are strongly dependent on the spanwise spacing (among other parameters). Thus, in the fully developed regime, a staggered layout does not always lead to a higher power production. For example, the present study focuses on a uniform turbine spacing of $6D$ in the streamwise and spanwise directions and finds a higher power output for a staggered layout than for an aligned layout. However, previous LES and wind-tunnel experiments (Stevens, Gayme & Meneveau 2014, 2016a; Wu & Porté-Agel 2017; Bossuyt, Meneveau & Meyers 2018) have shown that, for certain conditions with smaller spanwise spacings (${\le }5D$), staggered and aligned layouts can lead to the same asymptotic power output in the fully developed regime.

A closer inspection of figure 7 reveals that the predictions for the staggered layout are generally better than the aligned layout. This is reasonable as the flow in the staggered wind farm is closer to the horizontally homogeneous condition than the aligned layout (Stevens et al. 2016a), see figure 2, which is one of the underlying assumptions of the model. Furthermore, it should be appreciated that wind farm designs of extended wind farms would be closer to a staggered layout as it optimizes wind farm performance. The figure also shows that the power decreases as the free atmospheric lapse rate increases. The reason is that kinetic energy entrainment, which is the only source of kinetic energy that balances the power extracted by the turbines, decreases as the stability in the free atmosphere increases. As a result, less power can be extracted from the flow by the wind turbines (Abkar & Porté-Agel 2013, 2014). In addition, one notes that the predictions for the hub height velocity tend to be closer to the LES values than the corresponding power production estimates. The reason is that the power production scales as $\beta ^{3} U_h^{3}$, making the power more difficult to predict than just the velocity $U_h$.

6. Effect of turbine thrust coefficient

To examine the validity of the proposed theoretical framework, we consider the second set of simulations in which the free-atmosphere lapse rate $\varGamma =4$ K km$^{-1}$ is fixed, while the latitude $\phi$, the wind farm layout and the thrust coefficient $C_T'$ are varied. In particular, the thrust coefficient $C_T'$ changes between $C_T' = 0.1$ and $C_T' = 2$, i.e. the Betz limit; see table 2. Figure 8(a) shows that the empirical constant $a_u$ in (4.1) increases from zero at $C_T'=0$ and asymptotically approaches $a_u \approx 4.3$ for $C_T'=4/3$. The lower limit is consistent with the physical assumptions as it indicates that the wake correction term in (4.1) disappears when wind turbines are turned off. The fact that $a_u$ asymptotes to a constant value indicates that the deviation from the logarithmic law below the turbine hub height caused by wake effects saturates for a sufficiently large thrust coefficient, say $C_T' \ge 4/3$, as the energy entrainment from above is limited. For the analytical predictions, the value of $a_u$ is parametrized as $a_u=4.3\tanh (2C_T')$, which is in good agreement with the LES data as shown in figure 8(a). Moreover, figure 8(b) shows the comparison of the analytical predictions and LES measurements of the flow inhomogeneity factor $\beta$ for different thrust coefficients, of which the overall agreement is good for both the aligned and staggered wind farms. It is also presented that the value of $\beta$ is nearly independent of the latitude, indicating that our assumption of neglecting the Coriolis force in the wake region is reasonable.

Figure 8. The $C_T'$-dependence of (a) the empirical constant $a_u$ and (b) $\beta$ for aligned (circles) and staggered (diamonds) wind farms. The solid line in (a) $a_u=4.3\tanh (2C_T')$ and (b) $\beta$ is obtained from model calculations; see § 4 for details.

Figure 9(a,b) shows that the normalized wind speed at hub height $U_h/G$ reduces monotonically as the thrust coefficient $C_T'$ increases. As staggered wind farms extract more kinetic energy from the ABL flow than aligned wind farms under the same conditions, $U_h/G$ is smaller for the staggered wind farms. Furthermore, $U_h/G$ is larger at higher latitudes. This phenomenon originates from the increasing Coriolis force in the free atmosphere that drives the conventionally neutral ABL flow. With increasing latitude $\phi$, the Coriolis parameter $f=2\varOmega \sin \phi$ increases, allowing for an increased normalized wind speed at hub height $U_h/G$. It is worthwhile to point out that increasing $C_T'$ also increases the drag coefficient $u_{*,2}/G$ and the equivalent wind farm surface roughness $z_{0,2}$ perceived by the mean flow aloft, a trend supported well by the results presented in table 2.

Figure 9. The $C_T'$-dependence of the (a,b) normalized wind speed at hub height $U_h/G$ and the corresponding (c,d) normalized power output per surface area $10^{3}P/G^{3}$ for (a,c) aligned and (b,d) staggered wind farms. Lines: present model; filled symbols: LES data.

Figure 9(c,d) shows the corresponding normalized wind farm power outputs per surface area. When considering a single turbine in idealized conditions, the optimal operating condition of the turbine corresponds to $C_T'=2$, corresponding to the Betz limit (Burton et al. 2001). However, this figure shows that the normalized production of fully developed wind farms reaches its maximum value at $C_T'\approx 4/3$. A similar phenomenon has also been observed for the fully developed wind farm in truly neutral ABLs by Goit & Meyers (2015). The reason is that the power production $P/G^{3}$ is proportional to the thrust coefficient $C_T'$ and the cubed hub height velocity $(U_h/G)^{3}$. As the thrust coefficient $C_T'$ increases, the hub height velocity $U_h/G$ decreases due to the wake effects (see table 2). For a low thrust coefficient, the benefit of increasing $C_T'$ is more significant than the decrease of $U_h/G$, leading to increased power output. However, the benefit of increasing $C_T'$ is overwhelmed by the reducing hub height velocity for a high thrust coefficient.

7. Conclusions

Using LES of fully developed wind farms in conventionally neutral ABLs, we demonstrate that the normalized boundary-layer height and the GDL coefficients $A$ and $B$ can be parameterized as a function of just the Zilitinkevich number $Zi$. This confirms that the similarity arguments from the GDL apply for flow over extended wind farms, even though their physics are different. The coefficients $A$ and $B$ are changed with respect to flow over flat terrain as the wind farm friction velocity is higher than over flat terrain. The canopy and wind farm flow similarity suggests that these findings may be relevant for canopy flow.

Inspired by the classical work of Frandsen (1992), we develop a new top-down model that predicts the performance of fully developed wind farm arrays as a function of the geostrophic wind forcing in conventionally neutral ABLs. The model captures the effect of the free-atmosphere lapse rate $\varGamma$, the latitude $\phi$, the wind farm layout and the thrust coefficient $C_T$ on the wind farm performance. We demonstrate that the highest wind farm power output is obtained for the thrust coefficient well below the Betz limit. This means that the optimal wind farm performance is obtained when the wind farm is operated as a complete system and not when each turbine optimizes its performance. Furthermore, the model predicts that the normalized power production per surface area decreases with the free-atmosphere lapse rate, which is consistent with the finding of Abkar & Porté-Agel (2013). This phenomenon originates from the decreasing kinetic energy entrainment with increasing stability in the free atmosphere. In contrast, the power production per surface area normalized by the cubed geostrophic wind speed increases with the latitude, which is consistent with the finding of Antonini & Caldeira (2021) and results from the increasing Coriolis force in the free atmosphere with increasing latitude.

Future research will focus on developing a similar model to include the effect of thermal stratification in the boundary layer. To do this, two main adjustments may be needed. The first one is to include a thermal stratification correction to the initially assumed logarithmic profile (Sescu & Meneveau 2015). The second one is to extend the GDL approach towards thermally stratified flows (Tong & Ding 2020; Kadantsev et al. 2021).