1. Introduction
The global wind energy production is set to rise in the next few decades. To achieve this, wind farm clusters are expected to be built which are an order of magnitude larger than existing ones (Maas & Raasch Reference Maas and Raasch2022). To optimise their design, it is important to accurately predict the total farm power as well as aerodynamic loads on each turbine. However, it is very difficult to model the aerodynamics of wind farms because of the multi-scale nature of the wind farm flows (Porté-Agel, Bastankhah & Shamsoddin Reference Porté-Agel, Bastankhah and Shamsoddin2020). In 2019, Ørsted, one of the largest offshore wind farm developers, announced that its wind farms were producing less power than expected (Ørsted 2019). The underprediction of farm power was attributed to two effects: wake and farm blockage effects.
Behind every turbine, there is a turbulent wake which has a reduced wind speed. When turbine wakes interact with other turbines within a farm, this can cause substantial power losses. This is known as the ‘wake effect’ and has been measured to reduce the power of turbines in an existing wind farm by up to 40 % in the worst case (Barthelmie et al. Reference Barthelmie, Pryor, Frandsen, Hansen, Schepers, Rados, Schlez, Neubert, Jensen and Neckelmann2010). This effect has been extensively investigated in the literature using large-eddy simulations (LES) (Porté-Agel, Wu & Chen Reference Porté-Agel, Wu and Chen2013; Wu & Porté-Agel Reference Wu and Porté-Agel2015; Stevens, Gayme & Meneveau Reference Stevens, Gayme and Meneveau2016a). The ‘farm blockage’ is a recently observed effect in large wind farms (Bleeg et al. Reference Bleeg, Purcell, Ruisi and Traiger2018). The flow resistance caused by a wind farm reduces the wind speed upstream as well inside the farm. Hence, the total power produced by the wind farm is reduced compared to the ideal situation where the upstream wind speed is not affected by the farm (Nishino & Dunstan Reference Nishino and Dunstan2020).
Wind farm aerodynamics have been traditionally modelled using ‘wake’ models, which predict the velocity behind a turbine (e.g. Jensen Reference Jensen1983; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2014). To account for interactions between multiple turbines, the wake velocity deficits are superposed (e.g. Katic, Hojstrup & Jensen Reference Katic, Hojstrup and Jensen1986; Zong & Porté-Agel Reference Zong and Porté-Agel2020). However, these models do not account for the response of the atmospheric boundary layer (ABL) to the farm. As such, they tend to perform poorly for large wind farms (Stevens et al. Reference Stevens, Gayme and Meneveau2016a). A different approach to modelling wind farms is to use ‘top-down’ models (e.g. Frandsen Reference Frandsen1992; Frandsen et al. Reference Frandsen, Barthelmie, Pryor, Rathmann, Larsen, Højstrup and Thøgersen2006; Calaf, Meneveau & Meyers Reference Calaf, Meneveau and Meyers2010). They model the response of an idealised ABL (which follows a logarithmic law) to an infinitely large wind farm, but cannot take into account the details of turbine–wake interactions explicitly. Hence, these wind farm models cannot correctly capture the two-way interactions of turbine-scale and farm-scale flow effects which determine the performance of large wind farms.
To capture this two-way interaction, it has been proposed to couple wake and top-down models; for example, by adjusting parameters in both models to match the hub-height-averaged velocities (Stevens, Gayme & Meneveau Reference Stevens, Gayme and Meneveau2016b; Starke et al. Reference Starke, Meneveau, King and Gayme2021). Examples of two-way coupling can also be found in existing engineering software, e.g. the ‘Deep Array Wake Model’ (Brower & Robinson Reference Brower and Robinson2012). However, these models involve the coupling of low-order flow models and are therefore limited by the assumptions made by the constituent models, e.g. wake superposition or a log-law wind profile. To account for the effects of more realistic flow physics, it would be beneficial to use an approach based on more fundamental laws of fluid mechanics.
The optimal design of a large wind farm under realistic atmospheric conditions remains a challenge as it requires consideration of the complex atmospheric response to the farm. It is often too expensive to run a large number of simulations which resolve both individual turbines and the atmospheric response to the farm. In addition, to find an optimal design for the long-term performance of a wind farm, such as the annual energy production (AEP), the range of timescales we would need to consider is too wide. As such, Nishino & Dunstan (Reference Nishino and Dunstan2020) proposed the ‘two-scale momentum theory’ to split the multi-scale flow problem into ‘internal’ turbine/array-scale and ‘external’ farm/atmospheric-scale problems. West & Lele (Reference West and Lele2020) performed LES of infinitely large wind farms and the results showed a good agreement with the two-scale momentum theory. However, their LES study was limited to ‘fully aligned’ turbine layouts, and their discussion was also limited to the special case where the momentum supplied by the atmosphere to the wind farm site was fixed. In reality, the strength of atmospheric response to the wind farm resistance depends on mesoscale weather patterns (Patel, Dunstan & Nishino Reference Patel, Dunstan and Nishino2021) as well as atmospheric stability and gravity waves (Allaerts & Meyers Reference Allaerts and Meyers2017, Reference Allaerts and Meyers2018, Reference Allaerts and Meyers2019).
The aim of the present study is to better understand the fluid mechanics processes which determine the power production of large wind farms, using a combined theoretical and computational approach. First, we will perform a large suite of LES of infinitely large wind farms with different turbine layouts and wind directions. We then combine the results of LES with the two-scale momentum theory to investigate and explain expected performance of large finite-sized wind farms with a realistic range of atmospheric response strengths. Using this approach allows the combined effects of turbine-scale and farm-scale flow characteristics on wind farm power to be determined.
In § 2, we summarise the definitions of key wind farm parameters in the two-scale momentum theory (Nishino & Dunstan Reference Nishino and Dunstan2020). Section 3 details the methodology of the LES and wind turbine implementation. In § 4, we present the results including validation of the LES code. These results are discussed in § 5 and concluding remarks are given in § 6.
2. Theory
2.1. Two-scale momentum theory
Figure 1 shows a pair of control volumes for a given farm site. We first consider the momentum balance of the control volume without the turbines (figure 1a). For this scenario, the following equation can be derived for a ‘short-time-averaged’ flow:
\begin{equation} \frac{\partial [\rho U_0]}{\partial t}=X_0-C_0-\left[\frac{\partial{p_0}}{\partial{x_{F0}}}\right]-\frac{\langle\tau_{w0}\rangle S_F}{V_{CV}}, \end{equation}
where 
$U$ is the velocity in the hub-height wind direction (i.e. streamwise direction) 
$x_F$, 
$X$ represents the net streamwise momentum injection through the top and side boundaries of the control volume (due to advection and Reynolds stress), 
$C$ is the streamwise component of the Coriolis force averaged over the control volume, 
$\partial p / \partial {x_F}$ is the pressure gradient in the direction 
$x_F$, 
$\tau _w$ is the bottom shear stress, 
$S_F$ is the wind farm area and 
$V_{CV}$ is the volume of the control volume. The subscript 
$0$ refers to values without the turbines present, 
$[\square ]$ refers to control-volume-averaged values and 
$\langle \square \rangle$ refers to farm-area-averaged values.
Control volume for an entire wind farm site: (a) without turbines and (b) with turbines.
By considering the control volume with the turbines present (figure 1b), the following equation can be derived:
\begin{equation} \frac{\partial [\rho U]}{\partial t}=X-C-\left[\frac{\partial{p}}{\partial{x_F}}\right]-\frac{\langle\tau_w\rangle S_F+\sum_{i=1}^n T_i}{V_{CV}}, \end{equation}
where 
$T_i$ is the thrust of turbine 
$i$ in the farm and 
$n$ is the total number of turbines in the farm. Equations (2.1) and (2.2) can be combined to obtain the non-dimensional farm momentum (NDFM) equation (Nishino & Dunstan Reference Nishino and Dunstan2020):
\begin{equation} C_T^* \frac{\lambda}{C_{f0}} \beta^2 + \beta^\gamma = M, \end{equation}
where 
$\beta$ is the farm wind-speed reduction factor defined as 
$\beta \equiv U_F/U_{F0}$ (with 
$U_F$ defined as the average wind speed in the nominal farm-layer of height 
$H_F$, and 
$U_{F0}$ is the farm-layer-averaged speed without the presence of the turbines); 
$\lambda$ is the array density defined as 
$\lambda \equiv nA/S_F$ (where 
$A$ is the rotor swept area); 
$C_T^*$ is the (farm-averaged) ‘local’ or ‘internal’ turbine thrust coefficient defined as 
$C_T^*\equiv \sum _{i=1}^{n}T_i/\frac {1}{2}\rho U_F^2nA$; 
$C_{f0}$ is the natural friction coefficient of the surface defined as 
$C_{f0}\equiv \langle \tau _{w0} \rangle /\frac {1}{2}\rho U_{F0}^2$; 
$\gamma$ is the bottom friction exponent defined as 
$\gamma \equiv \log _\beta (\langle \tau _w \rangle / \langle \tau _{w0} \rangle )$; and 
$M$ is the momentum availability factor given by
\begin{equation} M = \frac{X-C-\left[\dfrac{\partial p}{\partial x_F}\right] - \dfrac{\partial [\rho U]}{\partial t}}{X_0-C_0-\left[\dfrac{\partial p_0}{\partial x_{F0}}\right] - \dfrac{\partial [\rho U_0]}{\partial t}}. \end{equation}Note that the first term on the left-hand side of (2.3) can be extended to include the impact of support structure drag (see Ma, Nishino & Antoniadis Reference Ma, Nishino and Antoniadis2019), but its impact is usually small and is therefore neglected in this study for simplicity.
The height of the farm-layer, 
$H_F$, is used to define the reference velocities 
$U_F$ and 
$U_{F0}$. The value of 
$H_F$ is typically between those of 
$2H_{hub}$ and 
$3H_{hub}$ (where 
$H_{hub}$ is the turbine hub-height); the NDFM (2.3) is valid as long as the same 
$H_F$ value is used for both ‘internal’ and ‘external’ problems. In this study, we use a fixed definition of 
$H_F=2.5H_{hub}$. This is discussed further in Appendix A.
Equation (2.3) helps the analysis of large wind farm aerodynamics. This is because 
$C_T^*$ and 
$\gamma$ on the left-hand side are expected to depend primarily on the turbine/array-scale flow physics or ‘internal’ conditions, for example, the turbine layout, operating conditions and local wind conditions, whereas 
$M$ on the right-hand side is expected to depend largely on ‘external’ conditions. Following Nishino & Dunstan (Reference Nishino and Dunstan2020), in this study, we assume that the ‘internal’ problem (to be modelled using LES in § 3 to calculate 
$C_T^*$ and 
$\gamma$) can be modelled without explicitly considering the effects of ‘external’ conditions such as wind farm size and location, and the response strength of the atmosphere.
The ‘external’ problem is to determine the parameter 
$M$, which represents how much the amount of momentum available to the farm site differs from its ‘natural’ value. This problem is largely independent of the small-scale flow features and can be modelled using a numerical weather prediction (NWP) model with a wind farm parameterisation, i.e. without resolving individual turbines. Patel et al. (Reference Patel, Dunstan and Nishino2021) used such an NWP model to demonstrate that, for most cases, 
$M$ varied almost linearly with 
$\beta$ (for a realistic range of 
$\beta$ between 1 and 0.8). Therefore, 
$M$ can be approximated by
\begin{equation} M=1+\zeta (1-\beta), \end{equation}
where 
$\zeta$ is called the ‘momentum response’ factor or ‘wind extractability’ factor. Patel et al. (Reference Patel, Dunstan and Nishino2021) found 
$\zeta$ to vary between 5 and 25 for a typical offshore wind farm site. Note that 
$\zeta =0$ corresponds to the case where momentum available to the farm site is assumed to be fixed, i.e. 
$M=1$. An infinitely large value of 
$\zeta$ corresponds to the case where there is no wind speed reduction in the farm-layer, i.e. 
$\beta =1$. Preliminary results of an extended study from Patel et al. (Reference Patel, Dunstan and Nishino2021) show that 
$\zeta$ changes according to atmospheric conditions and decreases exponentially with increasing farm size (see Appendix B). Although the details of how 
$\zeta$ changes with weather conditions are still unclear and need to be clarified in future studies, in the present study, we take 
$0<\zeta <25$ as a typical range of the wind extractability for large offshore wind farms.
Using 
$\beta$ obtained from (2.3) for a set of 
$C_T^*$, 
$\gamma$ and 
$\zeta$, the following equation can be used to calculate the power coefficient of the turbines within the farm:
\begin{equation} C_p = \beta^3 C_p^*, \end{equation}
where 
$C_p$ is the (farm-averaged) turbine power coefficient defined as 
$C_p\equiv \sum _{i=1}^{n}P_i/\frac {1}{2}\rho U_{F0}^3nA$ (
$P_i$ is power of turbine 
$i$ in the farm) and 
$C_p^*$ is the (farm-averaged) ‘local’ or ‘internal’ turbine power coefficient defined as 
$C_p^*\equiv \sum _{i=1}^{n}P_i/\frac {1}{2}\rho U_F^3nA$.
2.2. Analytical model of ideal wind farm performance
In this study, we consider arrays of actuator discs (or aerodynamically ideal turbines operating below the rated wind speed). For an actuator disc, 
$C_p^*=\alpha C_T^*$, where 
$\alpha$ is the turbine-scale wind speed reduction factor defined as 
$\alpha \equiv U_T/U_F$ (
$U_T$ is the streamwise velocity averaged over the rotor swept area). We can estimate 
$\alpha$ using the expression 
$\alpha =\sqrt {C_T^*/C_T'}$, where 
$C_T'\equiv T/\frac {1}{2}\rho U_T^2 A$ is a turbine resistance coefficient describing the turbine operating conditions (noting that this is strictly valid only for infinitely large regular arrays of turbines where the farm-averaged turbine thrust is identical to the thrust of each individual turbine). The theoretical 
$C_p$ of an actuator disc is therefore given by
\begin{equation} C_p = \beta^3\alpha C_T^* = \beta^3 {C_T^*}^{3/2}{C_T'}^{-{1}/{2}}, \end{equation}
where 
$C_T^*$ may be predicted using a simple analytical model (Nishino Reference Nishino2016) given by
\begin{equation} C_T^* = 4\alpha (1 - \alpha) = \frac{16C_T'}{(4+C_T')^2} \end{equation}
using the expression 
$C_T'=C_T^*/\alpha ^2$ to express 
$C_T^*$ as a function of 
$C_T'$. The model predicts 
$C_T^*$ as a function of turbine-scale wind-speed reduction by using an analogy to the classical actuator disc theory. This simple analytical model will be compared with LES results later in § 4. Using the analytical model of 
$C_T^*$, (2.8), and the linear approximation of 
$M$, (2.5), (2.3) and (2.7) can be solved to give a theoretical prediction of 
$C_p$, which we will call 
$C_{p,Nishino}$. Note that West & Lele (Reference West and Lele2020) also introduced this, but only for the special case with 
$\zeta =0$. As shown by Nishino & Dunstan (Reference Nishino and Dunstan2020), 
$C_{p,Nishino}$ is sensitive to 
$\zeta$ but much less sensitive to 
$\gamma$. If we assume that 
$\gamma =2.0$, then we can obtain an analytical expression for 
$C_{p,Nishino}$, i.e.
\begin{equation} C_{p,Nishino} = \frac{64C_T'}{(4+C_T')^3} \left[ \dfrac{-\zeta + \sqrt{\zeta^2 + 4\left( \dfrac{16C_T'}{(4+C_T')^2}\dfrac{\lambda}{C_{f0}}+1\right)(1+\zeta)}}{2\left( \dfrac{16C_T'}{(4+C_T')^2}\dfrac{\lambda}{C_{f0}}+1\right)} \right]^3. \end{equation} It is worth noting that the power coefficient of an isolated turbine, 
$C_{p,Betz}$, is given by
\begin{equation} C_{p,Betz} = \frac{64C_T'}{(4+C_T')^3}, \end{equation}
which takes the well-known maximum value of 
$16/27$ at 
$C_T'=2$. This equation can be obtained by substituting (2.8) into (2.7) with 
$\beta =1$ (i.e. assuming flow mechanisms as described by the classical actuator disc theory and no farm-scale wind speed reduction). Note that this is the same as solving (2.9) for two special cases: (1) with 
$\lambda /C_{f0}=0$; and (2) with an infinitely large value of 
$\zeta$.
3. LES modelling
3.1. Governing equations of the flow
We performed LES of flow over periodic turbine arrays using the MetOffice/NERC Cloud (MONC) Model (Brown et al. Reference Brown, Weiland, Hill and Shipway2018). The flow is driven by an imposed pressure gradient and is neutrally stratified. The flow is governed by the incompressible Navier–Stokes equations, i.e.
$$\begin{gather} \frac{\partial{u_i}}{\partial{x_i}}=0, \end{gather}$$
$$\begin{gather}\frac{\partial{u_i}}{\partial{t}}+u_j\frac{\partial{u_i}}{\partial{x_j}} ={-} \frac{\partial}{\partial{x_i}}\left(\frac{p'}{\rho}\right) + \frac{1}{\rho}\frac{\partial{\tau_{ij}}}{\partial{x_j}} + f_i - \frac{1}{\rho}\frac{\partial{p_{\infty}}}{\partial{x_i}}, \end{gather}$$
where 
$u_i$ is the resolved velocity in the 
$i$ direction, 
$p'$ is the pressure perturbation from the reference state, 
$\rho$ is the reference density, 
$\tau _{ij}$ is the subgrid stress term, 
$f_i$ is the force added to model the wind turbines and 
$\partial {p_{\infty }}/\partial {x_i}$ is the imposed pressure gradient.
The subgrid stress model is a standard Smagorinsky model (Smagorinsky Reference Smagorinsky1963) given by 
$\tau _{ij}=\rho \nu S_{ij}$, where 
$\nu$ is the subgrid-scale eddy viscosity and 
$S_{ij}$ is the rate of strain tensor. The eddy viscosity is given by a mixing length model 
$\nu =l^2S$, where 
$l$ is the mixing length scale and 
$S$ is the modulus of the rate of strain tensor 
$S=\|S_{ij}\|/\sqrt {2}$. Near the bottom boundary, the mixing length scale 
$l$ is damped using the function 
$1/l^2 = 1/(l_0)^2 + 1/[\kappa (z+z_0)]^2$ described by Brown, Derbyshire & Mason (Reference Brown, Derbyshire and Mason1994), where 
$l_0$ is the basic mixing length scale. Here, 
$l_0$ is given by 
$l_0=c_s \varDelta$ with a coefficient of 
$c_s=0.23$ and a grid spacing given by 
$\varDelta =\max (\Delta x,\Delta y)$ (Brown et al. Reference Brown, Derbyshire and Mason1994).
All velocity components are set to zero at the bottom boundary. The shear stress at the surface is parameterised by specifying 
$\nu$ following the classical Monin–Obukhov similarity theory. The horizontal boundary conditions are periodic for all prognostic quantities. The top boundary has a zero vertical velocity boundary condition and a damping layer for the top 200 m of the domain (which was not necessary in the present study for neutrally stratified flows, but still included for future studies to explore the effect of atmospheric stability).
3.2. Actuator disc implementation
We model individual turbines as actuator discs following the methodology used by the KULeuven code described by Calaf et al. (Reference Calaf, Meneveau and Meyers2010). The approach uses a Gaussian convolution filter to apply the turbine force from the rotor plane onto the LES grid. This allows the position and orientation of turbines to be changed easily. The thrust force exerted by a single turbine is given by
\begin{equation} F ={-}\frac{1}{2} \rho C_T' \widehat{U_T}^2 \frac{\rm \pi}{4} D^2, \end{equation}
where 
$\widehat {U_T}$ is the time-filtered disc-averaged velocity and 
$D$ is the turbine diameter. This turbine thrust force is spatially distributed using a normalised indicator function 
$\mathcal {R}(\boldsymbol {x})$, defined as
\begin{equation} \mathcal{R}(\boldsymbol{x}) = \frac{4}{{\rm \pi} D^2} \iint{G(\boldsymbol{x}-\boldsymbol{x'})}\,{\rm d}^2\boldsymbol{x'}, \end{equation}
where 
$G(\boldsymbol {x})$ is a filtering kernel. This integral is calculated over the surface of the disc. We divide the disc area into 10 segments in the radial and angular directions, respectively. This was sufficient for 
$\mathcal {R}(\boldsymbol {x})$ to be independent of the number of segments. MONC uses a staggered grid where the 
$u$ and 
$v$ velocities are evaluated at different points. As such, two different indicator functions, 
$\mathcal {R}_x(\boldsymbol {x})$ and 
$\mathcal {R}_y(\boldsymbol {x})$, are calculated for the 
$x$ and 
$y$ directions. We use the same filtering kernel as described by Shapiro, Gayme & Meneveau (Reference Shapiro, Gayme and Meneveau2019),
\begin{equation} G(\boldsymbol{x}) = \left(\frac{6}{{\rm \pi} \delta^2}\right)^{3/2} \exp\left(-\frac{6\|\boldsymbol{x}\|}{\delta^2}\right), \end{equation}
where 
$\delta$ is the filter width, which, following the approach of Shapiro et al. (Reference Shapiro, Gayme and Meneveau2019), is given by 
$\delta = 1.5\sqrt {\Delta x^2 + \Delta y^2 + \Delta z^2}$.
The force per unit density at a given grid point 
$\boldsymbol {x}$ in the direction 
$i$ is given by
\begin{equation} f_i = \frac{1}{2}C_T' \widehat{U_T}^2\frac{\rm \pi}{4} D^2 \mathcal{R}_i(\boldsymbol{x}). \end{equation} The disc-averaged turbine velocity 
$U_T$ is calculated using the indicator function 
$\mathcal {R}(\boldsymbol {x})$ as a weighting function,
\begin{equation} U_T = \iiint{u(\boldsymbol{x})\mathcal{R}_x(\boldsymbol{x})\cos{\theta}\,{\rm d}^3\boldsymbol{x}} + \iiint{v(\boldsymbol{x})\mathcal{R}_y(\boldsymbol{x})\sin{\theta}\,{\rm d}^3\boldsymbol{x}}, \end{equation}
where 
$\theta$ is the wind direction relative to the 
$x$ direction. We use a constant value for 
$\theta$ which is the direction of the pressure gradient forcing. Note that 
$u$ refers to the velocity in the 
$x$ direction whereas 
$U$ describes velocities in the wind direction.
The spatially averaged velocity 
$U_T$ is then temporally averaged using a one-sided exponential time filter with a time window of 10 minutes to calculate 
$\widehat {U_T}$. To calculate 
$C_T^*$ from the LES, we use the following relationship:
\begin{equation} C_T^*= \frac{C_T'}{nU_F^2} \sum_{i=1}^n(\widehat{U_T}^2)_i, \end{equation}
noting that turbine velocity 
$U_T$ is time filtered before being squared and then averaged over all 
$n$ discs. Here, 
$U_F$ is calculated by integrating the streamwise velocity (
$u\cos \theta +v\sin \theta$) between the surface and 
$2.5H_{hub}$ across the entire domain. Unlike 
$\widehat {U_T}$, no time filter is used to calculate 
$U_F$. The 
$C_T^*$ varies with time during the LES so is time-averaged over a long period to give a single value of 
$\overline {C_T^*}$.
To calculate the (farm-averaged) turbine power coefficient from the LES, we use the expression,
\begin{equation} C_p= \frac{C_T'}{nU_{F0}^3} \sum_{i=1}^n(\widehat{U_T}^3)_i, \end{equation}
where 
$U_{F0}$ is the farm-layer-averaged velocity in an LES without turbines. The 
$C_p$ varies with time so 
$\overline {C_p}$ is calculated by time averaging over a long period.
4. Results
4.1. LES code validation
We first validate our LES framework with the new actuator disc implementation by comparing with the benchmark cases reported by Calaf et al. (Reference Calaf, Meneveau and Meyers2010). We then investigate the sensitivity of our results to horizontal resolution, domain size and pressure solver.
For the validation cases summarised in table 1, we use a surface roughness length of 
$z_0=0.1$ m and a pressure gradient of 
$(1/\rho )\,{\rm d}p_{\infty }/{{\rm d}\kern0.7pt x} = 1\times 10^{-3}$ m s
$^{-2}$. The turbines all have a hub height of 100 m and a diameter of 100 m. We use a turbine resistance of 
$C_T'=1.33$ and the same turbine spacing as for Case A1 of Calaf et al. (Reference Calaf, Meneveau and Meyers2010) (
$S_x=7.85D$ and 
$S_y=5.23D$). The surface roughness length, pressure gradient, turbine design and resistance are chosen to match the values used by Calaf et al. (Reference Calaf, Meneveau and Meyers2010). Validation cases V-1, V-2 and V-3 use a fast Fourier transform (FFT) pressure solver, whereas V-4 uses an iterative pressure solver. Cases V-1, V-3 and V-4 have a domain size 
$L_x \times L_y \times L_z$ of 
$3.14 \times 3.14 \times 1$ km (with 24 turbines) and case V-2 has a domain size of 
$6.28 \times 6.28 \times 1$ km (with 96 turbines). All validation cases were run for 100 000 seconds and flow data averaged between 
$t= 30\,000$ and 100 000 s. The convergences of two flow statistics for case V-1 are shown in figure 2. Figure 3 shows the instantaneous streamwise velocity plotted on a cross-streamwise plane 
$2.5D$ behind a row of turbines in the validation case with a double horizontal resolution.
Time series of 10-minute-time-averaged (a) internal turbine thrust coefficient 
$C_T^*$ and (b) friction velocity 
$u_*$ for validation case V-1. Note that the time-averaging here is different from the time-filtering used to calculated 
$\widehat {U_T}$.
Instantaneous streamwise velocity behind a row of turbines in validation case V-3.
Figure 4(a) shows the time and horizontally averaged streamwise velocity for the four validation cases. The velocity profiles agree well with the results of Calaf et al. (Reference Calaf, Meneveau and Meyers2010). Figure 4(b) shows the profiles of total shear stress and the results reported by Calaf et al. (Reference Calaf, Meneveau and Meyers2010). The shear stress profiles match well except for the region near the bottom surface. Note that the shear stress in our LES does not approach zero at the bottom surface as it includes both modelled and resolved components; the latter of which approaches zero as in the results of Calaf et al. (Reference Calaf, Meneveau and Meyers2010). The velocity profiles are insensitive to the horizontal domain size and the pressure solver.
(a) Mean velocity profiles and (b) mean shear stress profiles for all validation cases, compared with Case A1 of Calaf et al. (Reference Calaf, Meneveau and Meyers2010).
In wind farm LES using the actuator disc method, the disc-averaged velocity is usually overpredicted (Shapiro et al. Reference Shapiro, Gayme and Meneveau2019). This is because at coarse resolutions, the vorticity shed from the disc edge is not fully captured. This can be seen in our results by considering cases V-1 and V-3 in figure 4(a). In the coarse grid case V-1, the disc-averaged velocity is overpredicted so the turbine thrust applied is greater. This results in a slightly lower 
$\langle \bar {u}\rangle$ throughout the entire domain. This effect is also seen in the higher uncorrected 
$\overline {C_T^*}$ value in case V-1 compared to V-3 (table 1), which is due to the overprediction of the disc-averaged velocity.
To correct the overprediction of disc velocity in coarse LES, the following correction factor was proposed by Shapiro et al. (Reference Shapiro, Gayme and Meneveau2019):
\begin{equation} N = \left( 1+ \frac{C_T'}{2} \frac{1}{\sqrt{3{\rm \pi}}}\frac{\delta}{D}\right)^{{-}1}. \end{equation}
We apply this correction factor to the disc-averaged velocity by multiplying our uncorrected 
$\overline {C_T^*}$ values by 
$N^2$. The correction is applied after the simulation and not during. After correction, the horizontal resolution used in cases V-1 and V-3 only had a small impact on 
$\overline {C_T^*}$, suggesting that this correction factor can be successfully applied to a periodic array of actuator discs. For the value of 
$C_T'$ used here, the analytical model of 
$C_T^*$ in (2.8) gives a value of 0.75. All the validation cases in table 1 have a lower 
$\overline {C_T^*}$ than this because of wake interactions among turbines.
We also consider the effect of resolution on the wake velocity deficit. Figure 5 shows the average wake profiles for each of our validation cases in table 1. The wake velocity profiles are normalised by the farm-averaged velocity 
$\overline {U_F}$ for each validation case. This shows the far wake velocity deficit does not vary with the domain size, horizontal resolution and pressure solver. Comparing the wake profiles for V-1 and V-3 at 2
$D$ downstream (figure 5a) shows a small difference in the velocity deficit at the centre. This is because V-3 uses a different filter size for the projection of the turbine area (see § 3.2). When comparing the wake profiles 4
$D$ and 6
$D$ downstream, this difference is negligible. This shows that the far wake velocity profile is insensitive to the filter size used for the turbine area projection.
Summary of validation cases.
Normalised wake velocity deficit for each validation case (a) 2
$D$, (b) 4
$D$ and (c) 6
$D$ downstream.
To validate the capability of the code to simulate different wind directions, we also performed a simulation with a wind direction of 45
$^{\circ }$. The turbine layout in this simulation corresponds to Case K of Calaf et al. (Reference Calaf, Meneveau and Meyers2010) and is shown in figure 6(b). We used a resolution of 
$\Delta x/D \times \Delta y/D \times \Delta z/D$ of 
$0.227\times 0.227\times 0.0787$ and a domain size of 
$L_x \times L_y \times L_z$ of 
$3.61 \times 3.61 \times 1$ km (with 32 turbines). The horizontally averaged streamwise velocity is shown in figure 6(a). There is an excellent agreement with the results of Calaf et al. (Reference Calaf, Meneveau and Meyers2010) for the same turbine layout, demonstrating that our new actuator disc implementation for various wind directions (see § 3.2) is valid.
Validation case with 
$45^\circ$ wind direction: (a) mean profile of the horizontally average streamwise velocity and (b) turbine layout.
4.2. LES results
The LES cases for this study have the same setup as the validation case V-4 (see § 4.1) except for the domain size, the surface roughness length and the streamwise pressure gradient. To model offshore wind farms, we use a surface roughness length of 
$1 \times 10^{-4}$ m and a pressure gradient of 
$(1/\rho )\,{\rm d}p_{\infty }/{\rm d}\kern0.7pt x_F = 8.38 \times 10^{-5}$ m s
$^{-2}$ in the wind direction 
$x_F$, which results in 
$U_{F0}=10.103$ m s
$^{-1}$ and 
$C_{f0} = 1.607 \times 10^{-3}$ for a fixed nominal farm-layer height of 
$H_F=250$ m (both obtained from LES with no turbines). All cases were run for 100 000 s and flow data averaged between 
$t= 60\,000$ and 100 000 s. Note that we adopted a different spin-up and averaging period (compared to § 4.1) because of the different pressure gradient forcing.
We performed a suite of 50 simulations with different turbine layouts which are described by the parameters 
$S_x$, the turbine spacing in the 
$x$ direction, 
$S_y$, the turbine spacing in the 
$y$ direction and 
$\theta$, the wind direction relative to the 
$x$ direction (see figure 7a). The turbine operating conditions are the same for all simulations and are given by 
$C_T'=1.33$. We consider a realistic range of turbine layouts and wind directions: 
$S_x \in [5 D, 10 D]$, 
$S_y \in [5 D, 10 D]$ and 
$\theta \in [0^\circ, 45^\circ ]$. We only consider regular arrays and so, by symmetry, we only need to consider wind directions up to 
$45^\circ$. We adopt the minimum possible horizontal domain size (
$L_x$ and 
$L_y$) within the range between 3.14 and 6.28 km (depending on 
$S_x$ and 
$S_y$) as the validation results presented in § 4.1 suggest that the results would be insensitive to the domain size within this range.
Design of numerical experiments: (a) input parameters and (b) maximin design of 50 wind farm layouts.
We use a space filling maximin design (Johnson, Moore & Ylvisaker Reference Johnson, Moore and Ylvisaker1990; Santner, Williams & Notz Reference Santner, Williams and Notz2018) to select different turbine layouts in the parameter space (
$S_x$, 
$S_y$, 
$\theta$). The maximin algorithm iteratively selects a point which maximises the minimum distance to other points and to the boundaries of the parameter space. Figure 7(b) shows the 50 different turbine layouts selected in the parameter space.
Figure 8 shows the time-averaged flow fields from 4 of 50 cases. Figure 8(a) is for a case where the wind direction is almost perfectly aligned with a relatively small streamwise spacing between turbines of 5.76
$D$. This case gives a low 
$\overline {C_T^*}$ value of 0.585 due to strong wake effects. High speed regions between rows of turbines are formed because of the large cross-streamwise turbine spacings and aligned wind direction. Figure 8(b) is for a case with a high turbine density and the wind direction almost aligned along the diagonal. This arrangement is similar to a staggered layout. The streamwise spacings between turbines is larger than for figure 8(a) so the 
$\overline {C_T^*}$ has a higher value of 0.669 because of the increased wake recovery between turbines.
Time-averaged streamwise velocity at the turbine hub height for: (a) 
$S_x=5.76D$, 
$S_y=8.51D$, 
$\theta =1.32^\circ$; (b) 
$S_x=5.27D$, 
$S_y=5.07D$, 
$\theta =43.8^\circ$; (c) 
$S_x=7.59D$, 
$S_y=5.47D$, 
$\theta =16.7^\circ$; (d) 
$S_x=6.00D$, 
$S_y=6.21D$, 
$\theta =37.6^\circ$.
The flow field for a case with an intermediate wind angle is shown in figure 8(c). The turbine wakes are mostly misaligned with downstream turbines which minimises wake effects and gives a high 
$\overline {C_T^*}$ value of 0.752. This result agrees qualitatively with Stevens, Gayme & Meneveau (Reference Stevens, Gayme and Meneveau2014) in which it was found that the maximum farm power was produced by an intermediate wind direction. The results also give further evidence that the analytical model of 
$C_T^*$ proposed by Nishino (Reference Nishino2016) can be used to predict an upper bound to wind farm performance as it gives 
$C_T^*=0.75$ in this case. Figure 8(d) shows the streamwise velocity for a partially waked turbine layout. The partial wake effects cause the 
$\overline {C_T^*}$ to be reduced slightly to 0.713.
Figure 8 also shows the effect of the turbine layout on the farm-averaged wind speed 
$U_F$. The farm shown in figure 8(a) has a low array density and so has a high farm-averaged wind speed of 
$\beta =0.347$ (shown by the brighter colour). Figure 8(b) shows a farm with a high array density which resulted in a low farm-averaged speed of 
$\beta =0.248$. Figures 8(c) and 8(d) have similar intermediate array densities and so had intermediate farm-averaged wind speeds of 
$\beta =0.289$ and 
$\beta =0.284$.
Figure 9 shows that 
$\overline {C_T^*}$ was not a strong function of 
$S_x$ or 
$S_y$. The 
$\overline {C_T^*}$ was found to be a much stronger function of 
$\theta$ (see figure 10). The lowest 
$\overline {C_T^*}$ values were for small values of 
$\theta$ because of the high degree of turbine–wake interactions. When 
$\theta$ was very small 
$\overline {C_T^*}$ was also sensitive to the value of 
$S_x$ (see figure 10a). As 
$\theta$ increases, 
$\overline {C_T^*}$ increases rapidly until the maximum value around 15
$^\circ$. As 
$\theta$ increases further, 
$\overline {C_T^*}$ slowly decreases. This is because turbines start to become aligned along the diagonal (similar to the layout shown in figure 8b). When 
$\theta$ is greater than 15
$^\circ$, the minimum 
$\overline {C_T^*}$ value tends to be observed when 
$S_x\tan (\theta )/S_y$ is close to 1 (figure 10b). This corresponds to layouts where turbines are aligned along the diagonal (similar to the layout shown in figure 8b).
Scatter plots of 
$\overline {C_T^*}$ against (a) turbine spacing in the 
$x$ direction 
$S_x/D$ and (b) turbine spacing in the 
$y$ direction 
$S_y/D$ with the colour given by the wind direction 
$\theta$.
Scatter plots of 
$\overline {C_T^*}$ against wind direction with the colour given by (a) turbine spacing in the 
$x$ direction 
$S_x/D$ and (b) 
$S_x\tan (\theta )/S_y$.
Figure 11 shows that the range of 
$\bar {\gamma }$ from the wind farm LES is small (varying between 1.7 and 1.8). Nishino (Reference Nishino2016) suggested that 
$\gamma$ would be slightly less than 2 because the presence of turbines would increase the turbulence intensity in the ABL. These results, along with the findings of Dunstan, Murai & Nishino (Reference Dunstan, Murai and Nishino2018), provide evidence for this. Figure 11 shows that there is a slight variation of 
$\bar {\gamma }$ with wind direction and effective array density.
Scatter plot of 
$\bar {\gamma }$ against the wind direction 
$\theta$ with the colour given by the effective array density 
$\lambda /C_{f0}$.
4.3. Prediction of wind farm performance
Now we compare the (farm-averaged) turbine power coefficient 
$\overline {C_p}$ from the wind farm LES with the analytical model derived from the two-scale momentum theory, 
$C_{p,Nishino}$, (2.9).
To make a fair comparison of 
$C_p$ between the theory and the LES, we need to consider the fact that the coarse LES resolution caused the turbine thrust to be overpredicted. For 
$\overline {C_T^*}$, this was corrected for in post-processing using the correction factor 
$N$, (4.1). However, for 
$\overline {C_p}$, there are two simultaneous factors that need to be corrected. The first is that the disc velocity (relative to the farm-layer velocity) has been overpredicted, which increases the turbine power. The second is that the farm-layer velocity has been underpredicted (as shown earlier in figure 4a), which reduces the turbine power. The first effect can be corrected by using the correction factor 
$N$ (i.e. multiplying the raw 
$\overline {C_p}$ from the LES by 
$N^3$), but the second effect cannot be corrected in this manner.
To adjust for the second effect, we estimate the farm wind-speed reduction that would be obtained if a sufficiently fine resolution was used for the LES, 
$\beta _{fine,LES}$. This should be higher than the value obtained from the coarse resolution LES, 
$\beta _{coarse,LES}$. To calculate 
$\beta _{fine,LES}$, we assume
\begin{equation} \frac{\beta_{fine,LES}}{\beta_{coarse,LES}} = \frac{\beta_{fine,theory}}{\beta_{coarse,theory}}, \end{equation}
where 
$\beta _{fine,theory}$ and 
$\beta _{coarse,theory}$ are the predictions from the two-scale momentum theory. Here, 
$\beta _{fine,theory}$ is calculated by solving (2.3) (with 
$M=1$) using 
$\overline {C_T^*}$ from the LES and the corresponding array density 
$\lambda /C_{f0}$. To calculate 
$\beta _{coarse,theory}$, we do the same but with the uncorrected turbine thrust (i.e. 
$\overline {C_T^*}/N^2$). Since (2.3) is derived from the law of momentum conservation, the only assumption we are making in (4.2) is that the flow is independent of Reynolds number (which is a reasonable assumption as the change in Reynolds number between the fine and the coarse resolution cases is only of the order of 10 %). We can therefore estimate the turbine power from a fine resolution LES, using the expression
\begin{align} \frac{\overline{C_p}_{fine,LES}}{\overline{C_p}_{coarse,LES}} &= \left( \frac{{U_T}_{fine,LES}}{{U_T}_{coarse,LES}} \right)^3 = \left( \frac{{U_T}_{fine,LES}/{U_F}_{fine,LES}}{{U_T}_{coarse,LES}/{U_F}_{coarse,LES}} \right)^3\nonumber\\ &\quad \times \left( \frac{{U_F}_{fine,LES}}{{U_F}_{coarse,LES}}\right)^3 =N^3 \left(\frac{\beta_{fine,LES}}{\beta_{coarse,LES}} \right)^3, \end{align}
where 
$\overline {C_p}_{coarse,LES}$ is from the coarse resolution wind farm LES, (3.9).
The green symbols marked by 
$\zeta$=0 in figure 12 show the corrected 
$\overline {C_p}$ for each of the 50 wind farm LES runs. These are plotted against the effective array density 
$\lambda /C_{f0}$. The blue line shows the 
$C_p$ prediction using the theoretical model (2.9). The theory matches remarkably well with the LES results, despite that the theory does not account for the turbine layout parameters, namely 
$S_x$, 
$S_y$ and 
$\theta$. The reason for this excellent agreement for 
$\zeta =0$ will be discussed later in § 5.
Average turbine power coefficient 
$\overline {C_p}$ from the 50 wind farm LES adjusted for different wind extractability factors (
$\zeta$). The 
$C_p$ predicted by the two-scale momentum theory ((2.9) with 
$C_T'=1.33$) is shown by the blue lines for comparison.
Next, we use the results from LES of infinitely large wind farms to estimate the average power coefficients for large but finite-sized wind farms, following the concept of the two-scale momentum theory. We combine the infinite wind farm LES results with the simple linear model of the momentum availability factor 
$M$, (2.5). This does not capture the finite-size effects observed near the edge of the farm, but reveals general trends of turbine layout effects with different atmospheric responses. We assume that the farms are still sufficiently large that the flow over the farm is mostly fully developed (or more specifically, the dependency of 
$C_T^*$ on 
$S_x$, 
$S_y$ and 
$\theta$ is still approximately the same as that in the corresponding infinitely large farm). ‘Sufficiently large’ depends on atmospheric conditions (Wu & Porté-Agel Reference Wu and Porté-Agel2017), but it is likely to be of the order of 10 km.
First, we consider the balance of the pressure gradient forcing (PGF) with surface stress and turbine thrust in our LES,
\begin{equation} \rho (S_FC_f^*+nAC_T^*)U_F^2/2 = S_F L_z \Delta p/\Delta x, \end{equation}
where 
$C_f^*$ is the ‘local’ or ‘internal’ friction coefficient 
$C_{f}^*\equiv \langle \tau _{w} \rangle /\frac {1}{2}\rho U_{F}^2$ and 
$\Delta p/\Delta x$ is the PGF applied to the LES of an infinitely large wind farm.
If we assume 
$C_f^*$ and 
$C_T^*$ are independent of the PGF (as the Reynolds number of the flow is very high), then (4.4) shows that 
$U_F^2$ should be proportional to 
$\Delta p / \Delta x$; hence,
\begin{equation} \frac{U_{F,finite}^2}{U_F^2} = \frac{\Delta p_{finite}}{\Delta x} \left/\, \frac{\Delta p}{\Delta x},\right. \end{equation}
where 
$U_F$ is the farm-layer-averaged wind speed in the LES and 
$\Delta p / \Delta x$ is the corresponding PGF, and 
$U_{F,finite}$ is the new farm-layer-averaged wind speed that would result from an unknown PGF for a finite-sized wind farm, 
$\Delta p_{finite} / \Delta x$.
As the flow is quasi-steady and horizontally periodic with no Coriolis force in this study (meaning that we ignore the contributions of all terms except the PGF in (2.4)), the right-hand side of (4.5) is equivalent to the momentum availability factor 
$M$ in (2.4). This will not be true for real wind farms under real atmospheric conditions. However, this is valid for our simplified analysis where the flow across the farm is mostly fully developed and is driven purely by a PGF. As noted earlier, it was found by Patel et al. (Reference Patel, Dunstan and Nishino2021) that 
$M$ can be well approximated by 
$M=1+\zeta (1-\beta )$ and 
$\zeta$ was typically between 5 and 25 for an offshore wind farm site. Substituting this expression for 
$M$ into (4.5) gives
\begin{equation} \frac{U_{F,finite}^2}{U_F^2} = 1 + \zeta \left(1-\frac{U_{F,finite}}{U_{F0}}\right). \end{equation}
Since 
$U_F$ and 
$U_{F0}$ are both known from LES results, 
$U_{F,finite}$ can be calculated analytically for a given 
$\zeta$. Finally, we assume (again based on the Reynolds number independency) that the disc-averaged velocity 
$U_T$ scales with 
$U_F$ and as such
\begin{equation} C_{p,finite}=C_p \left(\frac{U_{F,finite}}{U_F}\right)^3, \end{equation}
where 
$C_{p,finite}$ is the average turbine power coefficient of the large finite-sized wind farm.
Figure 12 also shows the average power coefficients 
$\overline {C_p}$ estimated for three different 
$\zeta$ values. This shows that the atmospheric response can significantly impact the farm power. The 
$\overline {C_p}$ for 
$\zeta =25$ are roughly an order of magnitude higher than 
$\overline {C_p}$ for the same layout with 
$\zeta =0$. Note that these results are with a constant turbine operating condition of 
$C_T'=1.33$ which may not be optimal for a given 
$\zeta$. Figure 4 of Nishino & Dunstan (Reference Nishino and Dunstan2020) shows how the theoretically optimal turbine power coefficient, 
$C_{p,max}$, varies with the strength of atmospheric response.
Of particular interest in figure 12 is that, for a given 
$\lambda /C_{f0}$, the variation of 
$\overline {C_p}$ (due to different turbine layouts) increases with 
$\zeta$. To better understand this trend, the 
$\overline {C_p}$ values for the 50 turbine layouts are presented together with their 
$\overline {C_T^*}$ values in figure 13 for six different 
$\zeta$ values separately. Figure 13(a) shows the results for infinitely large farms with 
$\zeta$=0. This shows that the average turbine power coefficient is a strong function of the effective array density and insensitive to 
$\overline {C_T^*}$ meaning that turbine–wake interactions have little effect on the farm power. However, for finite-sized farms (figure 13b–f), the layouts with high wake interactions and a low 
$\overline {C_T^*}$ value tend to produce less power. These are shown by the darker plots which fall well below the theoretical prediction (blue line). Interestingly, some of the layouts seem to produce slightly higher power than predicted by the two-scale momentum theory. These layouts have a 
$\overline {C_T^*}$ slightly higher than 0.75 (the value for an isolated turbine) and this seems to be due to locally accelerated flow caused by the local blockage effect (Nishino & Draper Reference Nishino and Draper2015; Ouro & Nishino Reference Ouro and Nishino2021). These results suggest that both the array density and turbine–wake interactions are important for the performance for large finite-sized wind farms.
Turbine power coefficient 
$\overline {C_p}$ for (a) infinite wind farms (
$\zeta =0$) and (b)–( f) finite wind farms with the colour giving the 
$\overline {C_T^*}$ recorded for each layout. The 
$C_p$ predicted by the two-scale momentum theory is shown by the blue line.
Figure 13 also suggests that finite-sized wind farms are less sensitive to the effective array density than infinitely large farms. Figure 13(a) shows a roughly inverse relationship between 
$\overline {C_p}$ and 
$\lambda /C_{f0}$, whereas figure 13(b–f) show a more linear decrease of 
$\overline {C_p}$. Overall, the analytical model (2.9) predicts the variation of 
$\overline {C_p}$ with the effective array density well for all of the atmospheric responses. The theoretical model does not predict the effect of turbine–wake interactions because the model of 
$C_T^*$ (2.8) is a function of turbine operating conditions only. These results support the argument that the two-scale momentum theory can be used to provide an approximate upper limit on wind farm performance.
5. Discussion
The analytical model of 
$C_T^*$ (2.8) has been shown in this study to provide an approximate upper bound to wind farm performance. The model is based on the classical actuator disc theory with an upstream velocity of 
$U_F$. Equation (2.8) would provide accurate predictions of 
$C_T^*$ if the following two conditions are met: (1) the wind speed upstream of each turbine in the farm is 
$U_F$; and (2) the mechanism of the flow around each turbine is the same as that around an isolated actuator disc. In reality, the wind speed upstream of each turbine is often lower than 
$U_F$ due to wake effects, reducing 
$C_T^*$ (and thus 
$C_p$). Conversely, local blockage effects (Nishino & Draper Reference Nishino and Draper2015) may increase 
$C_T^*$ because: (1) it creates locally accelerated flows which may allow the upstream velocity of most turbines to be higher than 
$U_F$ (see e.g. Ouro & Nishino Reference Ouro and Nishino2021); and (2) it changes the mechanism of the flow around each turbine, allowing for a higher 
$U_T$ (for a given 
$C_T'$) than that predicted by the classical actuator disc theory. However, such a positive effect of local blockage can be exploited only when the layout is carefully optimised for a specific wind direction. As such, (2.8) can be used to predict the upper bound of farm performance or the performance of an ideal wind farm without the negative effect of turbine–wake interactions.
The results from this study also show that the 
$\overline {C_p}$ of infinitely large wind farms depends mainly on the effective array density (see figure 13a), i.e. the farm power is insensitive to the turbine-scale flow interactions. A closer look at the results suggest that the limited impact of turbine-scale interactions is due to the fact that the turbine drag is typically much greater than the surface drag. The ratio of total turbine drag to surface drag is given by
\begin{equation} \frac{\sum_{i=1}^n T_i}{\langle \tau_w \rangle S_F}=\frac{\frac{1}{2}\rho C_T^* n A U_F^2}{\frac{1}{2}\rho C_f^* S_F U_F^2} = \frac{\lambda C_T^*}{C_f^*}. \end{equation}
We measured the time-averaged value of 
$\lambda C_T^*/C_f^*$ for all 50 turbine layouts simulated. The mean value was 5.22, the minimum value was 
$2.93$ and the maximum was 
$8.59$. Therefore, the turbine drag was typically five times greater than the surface drag. Consider an offshore farm where 
$\lambda C_T^*/C_f^*=5$, 
$\lambda /C_{f0}=10$ (both typical values) and 
$\zeta =0$. This means the turbine drag is five times greater than the surface drag and the momentum supplied by the atmosphere does not change in response to the farm. The composition of the total drag (normalised by 
$\langle \tau _{w0} \rangle S_F$) for this scenario is shown by the bar 1 in figure 14(a). This state 1 is an equilibrium state, i.e. the normalised drag is balanced by 
$M$ as in (2.3). If there is a sudden small change in wind direction which increases the degree of wake interactions, then this will decrease 
$C_T^*$ (or the ratio of 
$T$ to 
$U_F^2$) for the farm. In this example, the turbine drag is now only four times greater than the surface drag (bar 2 in figure 14a). This corresponds to 
$C_T^*$ decreasing from 
$0.75$ to 
$0.6$ (which is close to the largest difference observed in the 50 wind farm LES). This state 2 is a non-equilibrium state, i.e. the normalised drag is not balanced by 
$M$ (due to the sudden small change of wind direction). However, if the momentum supplied by the atmosphere is unchanged, the wind speed in the farm 
$U_F$ will eventually increase to compensate for the reduced turbine drag. The surface and turbine drag both scale with 
$U_F^2$ so they are both expected to increase at the same rate. Bar 3 in figure 14(a) shows the new equilibrium state with a new composition of turbine and surface drag after the increase in 
$U_F$. Comparing states 1 and 3, the total turbine drag has been reduced only slightly because of the constant amount of momentum supplied by the atmosphere to the farm site. Therefore, for aerodynamically ideal turbines (or turbines operating below the rated wind speed), the average turbine power coefficient 
$C_p$ is also insensitive to turbine-scale flow interactions when the momentum supplied by the atmosphere is constant (
$\zeta =0$).
Examples of offshore wind farm drag composition at 3 different states: (1) equilibrium state before wind direction change; (2) non-equilibrium state immediately after a small change of wind direction; and (3) new equilibrium state after atmospheric response, for (a) 
$\zeta =0$ and (b) 
$\zeta =15$. (c) Schematic of the three states on the 
$M$ versus 
$(1-\beta )$ plot.
Figure 14(b) explains why turbine-scale flow interactions become more important as 
$\zeta$ increases. Now, the momentum supplied by the atmosphere changes with the farm wind-speed reduction factor 
$\beta$ according to 
$M=1+\zeta (1-\beta )$. The ratios between turbine and surface drag for states 1 and 2 are exactly the same as for the 
$\zeta =0$ case in figure 14(a). However, now when the wind speed increases in response to the reduced turbine drag, the momentum supplied by the atmosphere to the farm site changes. As the wind speed 
$U_F$ increases, 
$\beta$ increases so the momentum supplied by the atmosphere decreases (see figure 14c). Therefore, the sum of the turbine and surface drag for state 3 decreases compared to state 1. This explains why there is a much greater reduction in turbine drag (and thus power) for the non-zero 
$\zeta$ case than for 
$\zeta =0$. As 
$\zeta$ increases, there will be a larger decrease in the total drag in response to an increased wind speed. Therefore, the power losses due to turbine-scale flow interactions also increase as 
$\zeta$ increases.
To better understand the factors which determine the power output of wind farms, we propose three power loss factors. First, the turbine-scale loss factor 
$\varPi _T$ which is defined by
\begin{equation} \varPi_T \equiv 1 - \frac{C_p}{C_{p,Nishino}}, \end{equation}
where 
$C_{p,Nishino}$ is given by (2.9). Here, 
$\varPi _T$ represents the power losses due to turbine-scale (or internal) flow interactions only, separate from the losses due to the farm-scale atmospheric response. For real turbines, 
$\varPi _T$ would also include turbine design losses (i.e. power losses due to a non-ideal rotor design for a given 
$C_T'$). Figure 15(a) shows that the losses caused by turbine interactions are small for infinitely large farms, typically less than 5 %. For some turbine layouts, 
$\varPi _T$ is negative, meaning that 
$C_p$ exceeds 
$C_{p,Nishino}$. These layouts have 
$\overline {C_T^*}$ values greater than 0.75 which is the value given by (2.8). As discussed earlier, this is likely to be due to local blockage effects increasing the turbine incident velocity above the farm-layer velocity 
$U_F$. Figure 15(b–f) show that the turbine-scale losses are greater in finite-sized farms with the same turbine layout. Under different atmospheric conditions, the same turbine layout can give different turbine-scale losses. Across a realistic range of wind extractability factors, the turbine-scale losses from the same layout can vary significantly. As an example, the losses from one layout vary from 13 % to 22 % as 
$\zeta$ changes from 5 to 25. Figure 15 suggests the maximum losses due to turbine-scale flow interactions is likely to be approximately 20 % for large offshore wind farms.
Turbine-scale loss factor 
$\varPi _T$ for (a) infinite wind farms (
$\zeta =0$) and (b)–( f) finite wind farms with the colour giving the 
$\overline {C_T^*}$ recorded for each layout.
The farm-scale loss factor, 
$\varPi _F$, represents the power loss due to the atmospheric response to the whole farm, and is defined by
\begin{equation} \varPi_F \equiv 1 - \frac{C_{p,Nishino}}{C_{p,Betz}}, \end{equation}
where 
$C_{p,Betz}$ is given by (2.10). Note that in this study, 
$C_T'=1.33$ and hence 
$C_{p,Betz}=0.563$ (
$C_T'=2$ would give the optimal performance for an isolated turbine of 
$C_{p,Betz}=16/27$). Here, 
$\varPi _F$ represents the power loss accompanied by the reduction of the farm-average wind speed. Figure 16 shows that the farm-scale losses are typically more than twice the turbine-scale losses (i.e. 
$\varPi _T/\varPi _F$ is generally less than 0.5). This suggests that for large offshore wind farms, the atmospheric response to the array density is more important than the turbine-scale interactions (i.e. wake effects). Similarly, the total power loss factor, 
$\varPi$, can also be defined as
\begin{equation} \varPi \equiv 1 - \frac{C_p}{C_{p,Betz}} = 1 - (1-\varPi_T)(1-\varPi_F). \end{equation}Ratio of turbine-scale to farm-scale loss factors 
$\varPi _T/\varPi _F$ for (a) infinite wind farms (
$\zeta =0$) and (b)–( f) finite wind farms with the colour giving the 
$\overline {C_T^*}$ recorded for each layout.
It is important to note that the turbine-scale loss factor discussed above is smaller than what is typically referred to as ‘wake losses’ (see figure 17). ‘Wake losses’ are traditionally evaluated by comparing the farm power with the power produced by the first row of turbines (
$C_{p,1}$ in figure 17). However, this includes not only the losses due to wake interactions between turbines, but also the atmospheric response to the array density. In contrast, 
$\varPi _T$ represents the power losses solely due to the interactions between turbines. For the same reason, the farm-scale loss factor, 
$\varPi _F$, is larger than what is often referred to as ‘farm blockage losses’. These two different classifications of power losses (one using 
$C_{p,Nishino}$ and the other using 
$C_{p,1}$ as a point of reference) are both useful in different ways. The latter classification is straightforward when the value of 
$C_{p,1}$ is known; however, sometimes 
$C_{p,1}$ cannot be defined unambiguously (e.g. when there is no regular ‘front row’ which is perpendicular to the wind direction). The merit of using 
$C_{p,Nishino}$ is that this can be predicted analytically using (2.9).
Comparison of turbine-scale loss (TSL) and farm-scale loss (FSL) with what is known as wake loss (WL) and farm blockage loss (FBL). Here, 
$C_{p,1}$ is the power coefficient recorded by the first row of turbines in a farm.
It should be noted that figure 17 is for actuator discs (i.e. ideal turbines). Real turbines will experience additional power losses due to practical (non-ideal) turbine design.
As demonstrated earlier in figures 13, 15 and 16, the strength of the atmospheric response alters the wind farm performance. The importance of array density and turbine layout vary with the large-scale atmospheric conditions. The results shown in Appendix B suggest that the value of 
$\zeta$ tends to decrease as the wind farm size increases (at least within the range between 10 and 30 km; note that the values of 
$\zeta$ reported by Patel et al. (Reference Patel, Dunstan and Nishino2021) were for a fixed wind farm size of 20 km). These results seem to suggest that larger wind farms will be less sensitive to turbine-scale wake effects and more sensitive to farm-scale losses. This could have significant implications for the design of future wind farms.
It is also worth noting that similar trends of large wind farm aerodynamics have already been reported in the literature, e.g. in the wind farm LES performed by Wu & Porté-Agel (Reference Wu and Porté-Agel2017) in which they also used a relatively low surface roughness length of 0.05 m. They performed LES of finite and infinite wind farms under a weak and strong free-atmosphere stratification, using a rotational actuator disc model (representing Vestas V-80 2MW turbine) with an aligned and a staggered turbine layout with the same array density. The varying free-atmosphere stratifications change the strength of atmospheric response, which affects the wind farm blockage. Figure 11 of Wu & Porté-Agel (Reference Wu and Porté-Agel2017) shows the power production of different turbine rows for the different wind farms. For both atmospheric conditions, the layout of the infinite wind farm did not change the farm power. This agrees with our finding that the power output of infinite wind farms is insensitive to turbine interactions. The stronger free-atmosphere stratification induced a larger pressure gradient across the farm which implies a larger 
$\zeta$ value. The power output of the finite farm under the strong free-atmosphere stratification was more sensitive to the turbine layout. Under a weak stratification, the layout had a smaller impact on farm power. This is further confirmation that the large-scale atmospheric response changes the importance of turbine-scale interactions. Other LES studies have also generally found that the power output of finite wind farms (e.g. Archer, Mirzaeisefat & Lee Reference Archer, Mirzaeisefat and Lee2013; Porté-Agel et al. Reference Porté-Agel, Wu and Chen2013; Stevens et al. Reference Stevens, Gayme and Meneveau2014) is more sensitive to turbine layout than infinite farms (e.g. Yang, Kang & Sotiropoulos Reference Yang, Kang and Sotiropoulos2012; Abkar & Porté-Agel Reference Abkar and Porté-Agel2013; Yang & Sotiropoulos Reference Yang and Sotiropoulos2014). However, these previous studies used surface roughness lengths typical of rough onshore locations rather than offshore, meaning that the ratio of total turbine drag to surface drag (5.1) was smaller than the present study.
The impact of turbine layout on wind farm performance is typically investigated by reporting the normalised turbine power, where the turbine power is normalised by either the power of a turbine in the first row or that of a standalone turbine. The normalised power is often used as a measure of turbine–wake interactions within the wind farm. However, this could be misleading for a large wind farm where the power is also reduced due to the ABL response to the farm resistance. The turbine-scale loss factor, 
$\varPi _T$, gives the power losses due to wake interactions separated from farm-scale effects. Our results suggest that the power losses in large farms with the same layout can change with the strength of the large-scale atmospheric response. Hence, the relative performance of different turbine layouts can also vary with the strength of the atmospheric response. Layouts which produce more power with a specific large-scale atmospheric response may not perform as well under different responses.
A limitation of this study is that the wind farm LES consider only neutrally stratified atmospheric boundary layers. One uncertainty is how different stratifications would affect wake recovery within large wind farms, which could affect the 
$\overline {C_T^*}$ value for a given turbine layout. However, since the two-scale momentum theory has been derived from the principle of momentum conservation, which is valid for all atmospheric conditions, we expect that the trends found in this study (on how the importance of turbine-scale interactions changes with the strength of the atmospheric response) would be observed generally. The study also used a fixed turbine operating condition given by 
$C_T'=1.33$. This is typical for the current wind turbines operating below their rated wind speed (at which the farm blockage effects are most significant). Therefore, the results of this study give the trends for large farms if the operating conditions are similar to those currently used. As discussed earlier, the trends found in this study have also been observed in LES of finite wind farms with different atmospheric conditions and a different turbine model, supporting the generality of our findings.
6. Conclusions
In this study, we performed a combined theoretical and computational analysis of large wind farms, using the two-scale momentum theory (Nishino & Dunstan Reference Nishino and Dunstan2020) and new LES of infinitely large wind farms with different turbine layouts and wind directions. To consider a range of wind directions, we used a new implementation of the actuator disc model which is not aligned with the structured grid. First, we validated the LES against the results published by Calaf et al. (Reference Calaf, Meneveau and Meyers2010) and found an excellent agreement. We also confirmed that when the correction factor proposed by Shapiro et al. (Reference Shapiro, Gayme and Meneveau2019) is applied, the average internal turbine thrust coefficient 
$\overline {C_T^*}$ (representing the turbine thrust relative to the square of the average wind speed across the wind farm layer) is insensitive to the grid resolution. We then performed 50 wind farm LES with different turbine spacings and wind directions. The 
$\overline {C_T^*}$ was found to be a strong function of wind direction with a weaker dependence on turbine spacings.
This study adopted an approach to adjust the average power coefficient 
$\overline {C_p}$ of an infinitely large wind farm to estimate the power of large but finite-sized farms with the same layout, following the concept of the two-scale momentum theory. A full validation of this approach is currently infeasible as it would require a large set of finite-size wind farm LES, which will be the subject of future studies. However, we have provided detailed theoretical arguments on why we can estimate the power of large finite-size farms from the results of infinitely large farm LES. For infinitely large farms, 
$\overline {C_p}$ was found to be a strong function of the array density and agree remarkably well with the analytical model (2.9). The power produced by large finite-sized farms was found to depend on both the array density and turbine layout (including the effect of wind direction), although the analytical model still provides an approximate upper limit of the power (for a given array density). These results confirm that to model large but finite-sized wind farms, it is important to consider both the effect of array density and the turbine layout.
The impact of array density and turbine layout on farm power varies with the strength of atmospheric response or the ‘wind extractability’. The analytical model seems to predict very well the impact of array density for a given farm-scale atmospheric response. As such, we propose a new classification of power losses in large wind farms by introducing new metrics, namely a turbine-scale loss factor 
$\varPi _T$, a farm-scale loss factor 
$\varPi _F$ and a total loss factor 
$\varPi$. The turbine-scale loss factor describes the power losses due to turbine-scale wake interactions only (which is smaller than what is commonly called ‘wake losses’). Importantly, the turbine-scale loss factor varies with the strength of the large-scale atmospheric response. As an example, 
$\varPi _T$ varied from 0.13 to 0.22 for a given layout across a realistic range of the extractability factor (from 
$\zeta =5$ to 
$25$ in this study). Although further studies are required to better quantify the extractability factor, the results obtained to date suggest that as wind farms get larger, the proportion of power losses due to turbine–wake interactions will decrease. This is because (1) having a better turbine layout (with less wake interactions and higher 
$C_T^*$) does not substantially increase the momentum available to the farm site when 
$\zeta$ is small (figure 14) and (2) 
$\zeta$ seems to decrease as the wind farm size increases. The farm-scale loss factor describes the power loss due to the atmospheric response to the array density (which is larger than what is often called ‘farm blockage losses’). The 
$\varPi _F$ was typically more than twice as large as 
$\varPi _T$ for the 50 offshore wind farms considered in this study with a realistic range of atmospheric responses. This suggests that farm-scale flow effects have a greater impact than turbine-scale flow effects on the performance of large offshore wind farms.
Acknowledgements
A.K. acknowledges the NERC-Oxford Doctoral Training Partnership in Environmental Research (NE/S007474/1) for funding and training. We acknowledge use of the Monsoon2 system, a collaborative facility supplied under the Joint Weather and Climate Research Programme, a strategic partnership between the Met Office and the Natural Environment Research Council. We also thank Dr François-Xavier Briol for helpful discussions and assistance with the design of numerical experiments (figure 7b) and Joseph Ashton for assistance with the analysis of NWP data shown in Appendix B.
Funding
This work was supported by the Natural Environmental Research Council (NERC) award NE/S007474/1.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available at https://github.com/AndrewKirby2/wind-farm-modelling. This includes the wind farm LES results and the code used to calculate the data for figures 13, 15 and 16.
Author contributions
T.N. derived the theory, A.K. and T.D.D. performed the simulations, A.K. and T.N. analysed the data. A.K. wrote the paper with corrections from T.N. and T.D.D.
Appendix A. Sensitivity of results to nominal farm-layer height 
$H_F$
The height of the farm-layer, 
$H_F$, can be defined in several ways. This study uses the definition 
$H_F=2.5H_{hub}$, where 
$H_{hub}$ is the turbine hub height. Differences in the value of 
$H_F$ will change the values of 
$U_F$, 
$U_{F0}$, 
$\overline {C_T^*}$ and 
$C_{f0}$ of the LES results presented in § 4. However, as noted by Nishino & Dunstan (Reference Nishino and Dunstan2020), the conservation of momentum arguments used to derive (2.3) and to eventually predict 
$C_p$ from (2.6) are valid irrespective of the 
$H_F$ value, provided that the same 
$H_F$ value is used in both ‘internal’ and ‘external’ sub-problems. The question here is how much the results of the ‘internal’ sub-problem (i.e. LES results in this study) and ‘external’ sub-problem (i.e. the value of 
$\zeta$ obtained from NWP simulations) change with the choice of 
$H_F$. Although the ‘external’ sub-problem is outside the scope of this paper, Patel et al. (Reference Patel, Dunstan and Nishino2021) reported that the sensitivity of their NWP results to the choice of 
$H_F$ was minor. In the following, we present how the sensitivity of the LES results to the choice of 
$H_F$ eventually affects the value of 
$C_p$ obtained (assuming that 
$\zeta$ is not affected by 
$H_F$).
To investigate the sensitivity of the results to 
$H_F$, we repeat the analysis with 
$\zeta =25$ using 
$H_F=3H_{hub}$. Figure 18 compares the 
$\overline {C_p}$ values obtained for 
$H_F=2.5H_{hub}$ and 
$3H_{hub}$. Changing the value of 
$H_F$ has no effect on the theoretical predictions of 
$C_p$ for a given 
$\zeta$. Increasing 
$H_F$ decreases the LES-based results of 
$\overline {C_p}$ for each farm because the reference velocity 
$U_{F0}$ slightly increases. The increase in 
$U_{F0}$ is the same for all wind farm cases. Therefore, increasing 
$H_F$ to 
$3H_{hub}$ reduces the 
$\overline {C_p}$ of all wind farm cases by the same factor of 
$1.043$ in this study. The absolute values of 
$\overline {C_p}$ are affected by the value of 
$H_F$, but the trends reported in this study remain unchanged. The 
$\overline {C_T^*}$ values obtained from the LES are also reduced when the farm-layer height is increased. This is because the reference velocity 
$U_F$ slightly increases. Nevertheless, the two-scale momentum theory still captures the trend well and can provide an upper bound estimate to wind farm performance.
LES and theoretical results for 
$\overline {C_p}$ with 
$\zeta =25$ for (a) 
$H_F=250$ m and (b) 
$H_F=300$ m.
Nishino & Dunstan (Reference Nishino and Dunstan2020) proposed a definition for 
$H_F$ based on the undisturbed velocity profile 
$\overline{U_0(z)}$. Here, 
$H_F$ is defined as the farm-layer height with which the farm-layer average of the undisturbed velocity is equal to that of the rotor average, i.e.
\begin{equation} \frac{\displaystyle\int_0^{H_F}\overline{U_0}\,{\rm d}z}{H_F}=\frac{\displaystyle\int\overline{U_0}\,{\rm d}A}{A}. \end{equation}
Using this definition, the exact value of 
$H_F$ will vary with turbine design and (undisturbed) ABL profile. Figure 19 shows the value of 
$H_F$ according to (A1) for a wide range of velocity profiles and turbine designs. We assume that the velocity profile is given by the logarithmic law for a wide range of surface roughness lengths that includes onshore and offshore locations. For the turbine design, we assume that the distance between the surface and bottom of the rotor is fixed at 50 m. The range of 
$D/H_{hub}$ used corresponds to turbine diameters from 100 m up to 300 m. Across a wide range of turbine design and velocity profiles, the value of 
$H_F$ obtained from (A1) only varies between 
$2.5H_{hub}$ and 
$3H_{hub}$.
Variation of normalised wind farm-layer height 
$H_F/H_{hub}$ with surface roughness length 
$z_0$ and turbine design parameter 
$D/H_{hub}$ for a log law wind profile.
The average turbine power coefficient 
$\overline {C_p}$ is changed only slightly when 
$H_F$ is increased from 
$2.5H_{hub}$ to 
$3H_{hub}$. The trends reported in this study are unaffected by a change in 
$H_F$. The 
$H_F$ is relatively insensitive to velocity profiles and turbine design and is expected to be in the range 
$2.5H_{hub}$ to 
$3H_{hub}$. Therefore, the results of the present study are insensitive to the exact value of 
$H_F$.
Appendix B. Effect of wind farm size on the wind extractability factor $\zeta$
To explore the effect of wind farm size on the wind extractability factor 
$\zeta$, we performed additional twin NWP simulations using the same methodology as Patel et al. (Reference Patel, Dunstan and Nishino2021). Within the simulations, the diameter of a hypothetical circular wind farm (located off the east coast of Scotland) was varied from 10 to 30 km. Example results for a 24-hour period with a relatively constant wind speed across the farm site (corresponding to Case B of figure 9 of Patel et al. Reference Patel, Dunstan and Nishino2021) are shown in figure 20. It can be seen that the value of 
$\zeta$ tends to decrease with increasing the wind farm size for all the time. Further simulations are required in future studies to derive a correction function for the effect of wind farm size on 
$\zeta$, but our preliminary results suggest that an exponential relationship may exist between the wind farm size and 
$\zeta$. Such an exponential correction would satisfy the expected asymptotic behaviour of 
$\zeta$ (i.e. 
$\zeta$ approaches to +
$\infty$ and 0 as the farm size approaches to 0 and +
$\infty$, respectively). Figure 20 also shows that 
$\zeta$ varies significantly with time due to changing atmospheric conditions. These preliminary results show that 
$\zeta$ depends both on wind farm size and atmospheric conditions.
Variation of 
$\zeta$ throughout a 24 h period (corresponding to Case B in figure 9 of Patel et al. Reference Patel, Dunstan and Nishino2021) for different wind farm diameters.
Open access
