A fast-running physics-based wake model for a semi-infinite wind farm

Abstract

This paper presents a new generation of fast-running physics-based models to predict the wake of a semi-infinite wind farm, extending infinitely in the lateral direction but with finite size in the streamwise direction. The assumption of a semi-infinite wind farm enables concurrent solving of the laterally averaged momentum equations in both streamwise and spanwise directions. The developed model captures important physical phenomena such as vertical top-down transport of energy into the farm, variable wake recovery rate due to the farm-generated turbulence and also wake deflection due to turbine yaw misalignment and Coriolis force. Of special note is the model's capability to predict and shed light on the counteracting effect of Coriolis force causing wake deflections in both positive and negative directions. Moreover, the impact of wind farm layout configuration on the flow distribution is modelled through a parameter called the local deficit coefficient. Model predictions were validated against large-eddy simulations extending up to 45 km downstream of wind farms. Detailed analyses were performed to study the impacts of various factors such as incoming turbulence, wind farm size, inter-turbine spacing and wind farm layout on the farm wake.

1. Introduction

Offshore wind is projected to experience rapid expansion in the coming decades, emerging as a significant global renewable energy source (Veers et al. 2019). To achieve this goal, many new offshore wind farms are anticipated to be erected in specific and promising geographical areas, particularly in regions like the North Sea, where strong and consistent winds are present. Consequently, the interaction among neighbouring offshore wind farms, as their wakes affect each other, has become an essential and pressing subject of research. Recent satellite images and field measurements have revealed that wakes of wind farms can last for many kilometres (Christiansen & Hasager 2005; Nygaard & Christian Newcombe 2018; Ahsbahs et al. 2020). Significant power degradation and fatigue loads can thus occur for a wind farm subject to wakes of adjacent wind farms (Stevens & Meneveau 2017). Beyond technical complexities, interactions between adjacent wind farms may lead to legal and financial disputes between operators of neighbouring facilities. As a result, accurate and reliable modelling of wind farm wake effects becomes of great importance for optimising future wind farms in increasingly competitive offshore environments.

High-fidelity numerical simulations such as large-eddy simulation (LES) are powerful tools for modelling complex turbulent wake flows, offering detailed insights into flow dynamics and wake interactions (Porté-Agel, Bastankhah & Shamsoddin (2020), and references therein). However, simulating a cluster of wind farms in congested areas such as the North Sea with LES is computationally intensive and time-consuming, making it impractical for real-time or large-scale studies. To address this challenge and enable more efficient simulations, there is a clear demand for fast-running engineering wake models striking a balance between accuracy and computational cost. Major advantages of these models are their ease of use and low computational costs, allowing for quicker assessments of various scenarios and aiding in optimisation of wind farm layouts and real-time control. Below, we attempt to classify the engineering wake models developed in the literature.

The typical method for modelling airflow distribution within wind farms involves predicting the wake generated by each individual turbine. A superposition method is then applied to consider the combined impact of these wake effects. These individual wake models range mainly from top-hat models (Jensen 1983; Katić, Højstrup & Jensen 1986) to Gaussian-type models (Bastankhah & Porté-Agel 2014). The Jensen top-hat model (also known as the Park model) has been extended in recent works to account for variable wake recovery rate due to turbine-generated turbulence (Nygaard et al. 2020). Over time, Gaussian wake models have also been refined and extended in several studies to more accurately describe the near-wake region (e.g. Keane et al. 2016; Shapiro et al. 2019; Blondel & Cathelain 2020; Schreiber, Balbaa & Bottasso 2020), to better capture wake expansion and its asymmetric shape (e.g. Abkar & Porté-Agel 2015; Xie & Archer 2015; Pedersen et al. 2022; Vahidi & Porté-Agel 2022) or to capture effects of yaw angle (e.g. Bastankhah & Porté-Agel 2016; King et al. 2020; Bastankhah et al. 2022; Bay et al. 2023) and wind veer (Abkar, Sørensen & Porté-Agel 2018; Mohammadi et al. 2022; Narasimhan, Gayme & Meneveau 2022). Moreover, a variety of wake superposition methods exist, aiming to model cumulative wake effects in wind farms (e.g. Lissaman 1979; Voutsinas, Rados & Zervos 1990; Niayifar & Porté-Agel 2016; Zong & Porté-Agel 2020; Bastankhah et al. 2021; Lanzilao & Meyers 2022). Some of these methods are solely empirical in nature, while others have a foundation in flow physics. See Bastankhah et al. (2021) for a detailed discussion of different wake superposition methods. This simple approach has proven to be very useful in providing detailed information on the flow field within small-sized wind farms and has been extensively used in wind farm layout optimisation and real-time flow control; see the review of Meyers et al. (2022) and references therein. However, this modelling approach cannot properly describe the interaction of wind farms with the atmospheric boundary layer (ABL) which involves scales that are comparable to the size of the entire wind farm or the ABL thickness. Most notably, these models fall short of capturing the crucial vertical transport of kinetic energy from higher-altitude layers of the atmosphere into the wind farm/wind farm wake (Stevens & Meneveau 2017). This becomes especially problematic as the size of wind farms grows, or if we seek information about the wake of an entire wind farm several kilometres downstream.

Capturing large-scale wind farm physics may be more readily achieved using infinite wind farm models. In this approach, the wind farm is assumed to be infinitely large in both lateral and streamwise directions, and the whole wind farm is modelled as an area with an increased aerodynamic surface roughness. Unlike single-wake modelling, this approach is able to capture the vertical transport of energy caused by turbulent fluxes, which is in balance with the energy extracted by wind turbines in infinite wind farms. The interested reader is referred to the seminal works of Frandsen (1992) and Calaf, Parlange & Meneveau (2011) and other subsequent studies (e.g. Frandsen et al. 2006; Meneveau 2012; Meyers & Meneveau 2012; Yang, Kang & Sotiropoulos 2012; Abkar & Porté-Agel 2013; Stevens, Gayme & Meneveau 2016) for more information. Despite the great advantage of these models in capturing the farm–atmosphere interaction, the concept of an infinite wind farm can be only regarded as an asymptotic case that resembles what very large wind farms may tend to approach. More importantly, these models fail to offer any insight into the wake of the wind farm due to their core assumption that the wind farm extends infinitely in the streamwise direction.

The other group of existing models, which we classify within the broad category of multi-scale models, strive to leverage the benefits of both large-scale farm and small-scale single-turbine modelling. Within this category, different approaches have been adopted to model wind farm flows. Stevens, Gayme & Meneveau (2015) and subsequent works (e.g. Shapiro et al. 2019; Starke et al. 2021) coupled the infinite-farm approach with the single-turbine approach by matching the predicted mean velocity at the turbine hub height. Other studies coupled the wind farm scale with the turbine scale through a parameter called the farm induction factor (e.g. Nishino & Dunstan 2020; Kirby, Nishino & Dunstan 2022) with more recent works modelling blockage effects as well (e.g. Legris et al. 2022). In another type of multi-scale model, the exchange of energy between the layer consisting of wind turbines and the overlaying boundary layer was parametrised using the classical entrainment theory (e.g. Luzzatto-Fegiz & Colm-cille 2018; Bempedelis, Laizet & Deskos 2023). Other multi-scale models characterised the farm–atmosphere interaction and farm-scale blockage effects caused by meso-scale phenomena such as gravity waves (Allaerts & Meyers 2019; Stipa et al. 2023). The coupling between the different scales in these models usually involves an iterative process or the numerical solution of governing equations. Moreover, the focus of the majority of the models discussed above is farm power production or the flow field within the wind farm, and less attention has been paid to the wake of the entire farm.

In this work, we propose a new category of wake models by considering a semi-infinite wind farm, i.e. a wind farm that extends infinitely in the lateral direction but has a finite size in the streamwise direction. The infinite lateral extent of the wind farm allows us to perform lateral averaging, which significantly simplifies the flow's governing equations and leads to a closed-form explicit solution without the need for using an iterative approach. The finite length of the wind farm also makes it possible to systematically model the wake of the entire farm. A schematic of the semi-infinite farm modelling in comparison with single-turbine modelling (i.e. finite approach) as well as infinite-farm modelling is shown in figure 1. A particular focus of this work is given to the prediction of the deflection of the farm wake. Predicting the magnitude of the farm wake deficit is important but not sufficient. The wake deflection also needs to be quantified to determine whether the wake of a wind farm may impinge on a downstream farm. In general, the wake deflection is mainly caused by (i) meso-scale phenomena such as Coriolis force (and its by-product wind veer) and (ii) yaw misalignment. The latter has recently received a great deal of attention because of its importance in wake steering strategies (e.g. Fleming et al. 2017; Bastankhah & Porté-Agel 2019; Howland, Lele & Dabiri 2019; Campagnolo et al. 2020). The former, however, is mostly overlooked in prior modelling works. While the deflection of a single-turbine wake due to Coriolis force is expected to be negligible (Mohammadi et al. 2022), several studies mainly based on numerical simulations have underpinned the importance of the farm wake deflection caused by Coriolis force (e.g. van der Laan et al. 2015; Abkar, Sharifi & Porté-Agel 2016; Allaerts & Meyers 2017; Eriksson et al. 2019; Gadde & Stevens 2019). Interestingly, there has not been a universal agreement in the literature with regards to the direction of the wake deflection caused by the Coriolis force. van der Laan & Sørensen (2017) argued that this is due to the fact that the Coriolis force has two effects on the wake deflection. The direct effect turns the wake in the anticlockwise direction (seen from top) in the northern hemisphere, while the indirect effect (through wind veer) rotates the wake in the clockwise direction. Gadde & Stevens (2019) explained this phenomenon based on the direction of the vertical turbulent fluxes in the entrance region in comparison with those in the wake. One of our objectives with this new modelling framework is to capture the conflicting influence of the Coriolis force on the farm wake. This is achieved by concurrently solving momentum equations in both streamwise and spanwise directions. The outcome is a simple one-dimensional model that predicts both laterally averaged streamwise and spanwise velocities at the turbine hub height within and downwind of the wind farm.

Figure 1. Different approaches used to model wind farm flows. (a) Modelling each individual wind turbine wake and then using superposition techniques to account for cumulative wake effects. (b) Modelling a wind farm that is extended to infinity in both streamwise and spanwise directions as an added aerodynamic surface roughness. (c) Modelling a wind farm that is extended to infinity in the lateral direction but has a finite length in the streamwise direction.

The rest of the paper is organised as follows. Section 2 develops the laterally averaged Reynolds-averaged Navier–Stokes equations. Section 3 describes the high-fidelity numerical simulations used in this study to validate the developed model. Section 4 discusses the budget analysis that is conducted to identify dominant terms in the momentum equations. The farm wake model is then developed in § 5. Results are discussed in § 6, and finally a summary is provided in § 7.

2. Streamwise and spanwise laterally averaged Reynolds-averaged Navier–Stokes equations

We start by writing the steady-state Reynolds-averaged Navier–Stokes equation for high-Reynolds-number flows (i.e. negligible friction forces) using Einstein notation. For simplicity, we non-dimensionalise all variables and equations throughout this paper using a selection of scales based on the incoming flow and turbine characteristics. All spatial dimensions are normalised by the turbine rotor diameter $D$. All velocities $u_i$ are normalised by the incoming velocity $\mathcal {U}_h$ at the turbine hub height $z_h$. Static pressure $p$ is normalised by $\rho \mathcal {U}_h^2$, where $\rho$ is the air density. The dimensionless Coriolis frequency $f_c$ is defined as

(2.1)\begin{equation} \frac{Df_c}{\mathcal{U}_h}=2\varOmega\sin\phi, \end{equation}

where $\varOmega = 7.2921\times 10^{-5}$ rad s$^{-1}$ is the rotation rate of the Earth and $\phi$ is the latitude. Note that the dimensionless Coriolis frequency $f_c$ defined in (2.1) represents the ratio of Coriolis force to inertial force. This dimensionless parameter is in fact the inverse of the Rossby number $Ro$ that is commonly used in geophysical studies concerning flows in oceans and atmosphere (e.g. van der Laan et al. 2020).

The dimensionless form of the Reynolds-averaged Navier–Stokes equation reads as (Stull 2009)

(2.2)\begin{equation} \bar{u}_j\frac{\partial \bar{u}_i}{\partial x_j} =\varepsilon_{ij3}f_c\bar{u}_j-\frac{\partial \bar{p}}{\partial x_i}- \frac{\partial \overline{u_i'u_j'}}{\partial x_j}-\bar{f}_i, \end{equation}

where $u_i$ is the velocity component in the $x_i$ direction with $i = 1, 2, 3$ corresponding to the streamwise $x$, spanwise $y$ and vertical $z$ directions, respectively. Overbar denotes time averaging, and the turbulent fluctuating velocity is $u_i'=u_i-\bar {u}_i$. The permutation symbol is denoted by $\varepsilon _{ijk}$. Moreover, $\bar {f}_i$ is the time-averaged component of the turbine forces per unit volume, non-dimensionalised by $\rho \mathcal {U}_h^2/D$. Note that gravitational forces are neglected in (2.2) as only the streamwise and spanwise momentum directions are of interest in this study.

We separate the static pressure $p$ in (2.2) into that due to the background driving pressure gradient $p_g$ and that due to the presence of turbines $p_t$. The former is dictated by the force balance in the geostrophic layer on top of the ABL, while the latter is due to the pressure drop across the rotor disc and its recovery to the free-stream pressure downstream. The background driving pressure gradient $p_g$ can be written in terms of the geostrophic wind speed $\bar {u}_g$ (Stull 2009). This simplifies (2.2) to

(2.3)\begin{equation} \bar{u}_j\frac{\partial \bar{u}_i}{\partial x_j} ={-}\varepsilon_{ij3}f_c(\bar{u}_{gj}-\bar{u}_j)-\frac{\partial \bar{p}_t}{\partial x_i}- \frac{\partial \overline{u_i'u_j'}}{\partial x_j}-\bar{f}_i. \end{equation}

Different forms of spatial averaging such as surface or volumetric averaging have been frequently performed in prior studies on vegetation canopies or infinite wind farms (e.g. Calaf, Meneveau & Meyers 2010; Moltchanov, Bohbot-Raviv & Shavit 2011; Bai, Katz & Meneveau 2015; Goit & Meyers 2015). In this study, however, we perform lateral averaging (averaging only along the $y$ direction), because our aim is to determine how flow quantities at the hub-height level evolve along the streamwise direction $x$. Mathematically, this may be defined for an arbitrary variable $\psi$ as follows:

(2.4)\begin{equation} \langle \bar{\psi} \rangle (x,z)= \lim_{L\to\infty}\frac{1}{2L}\int_{{-}L}^L\bar{\psi}(x,y,z)\, {{\rm d}\kern0.05em y}, \end{equation}

where $\langle \rangle$ indicates lateral averaging and $[-L,L]$ is the lateral range over which the averaging is performed. The lateral fluctuation is defined as $\psi ''=\psi -\langle \psi \rangle$ and by definition $\langle \psi ''\rangle =0$. By performing the lateral average on (2.3), we obtain

(2.5)\begin{equation} \underbrace{\langle \overline{u_{j}}\rangle \frac{\partial \langle \overline{u_{i}}\rangle}{\partial x_{j}}}_{{A}} ={-}\underbrace{\varepsilon_{\mathit{ij3}}f_c (\langle \overline{u_{\mathit{gj}}}\rangle -\langle \overline{u_{j}}\rangle)\vphantom{\frac{\rangle\partial}{\partial}}}_{{C}} - \underbrace{{\frac{\partial \langle\overline{p_t}\rangle}{\partial x_i}}}_{{P}} - \underbrace{\frac{\partial \langle \overline{u_{i}' u_{j}'}\rangle}{\partial x_{j}}}_{{R}} - \underbrace{\frac{\partial \langle \overline{u_{i}}'' \overline{u_{j}}'' \rangle}{\partial x_{j}}}_{{D}}-\underbrace{\langle\, \bar{f}_i \rangle}_{{T}}. \end{equation}

The terms outlined in (2.5) are

1. $A$ Advection of momentum by mean flow.

2. $C$ Coriolis term.

3. $P$ Pressure gradient due to the presence of wind turbines.

4. $R$ Reynolds stress gradients.

5. $D$ Dispersive stress gradients.

6. $T$ Turbine forcing.

These terms are discussed in more detail in § 4. The dispersive stress term $\langle \overline {u_{i}}'' \overline {u_{j}}'' \rangle$ that is the product of spatial fluctuations in the lateral direction arises in (2.5) as a result of lateral averaging. It is also noteworthy that because of lateral averaging any terms including $\partial \langle \rangle /\partial x_j$ in (2.5) must be zero if $j=2$ (i.e. $\partial /\partial y=0$).

3. Numerical set-up: LES

Large-eddy simulations were performed using the open source software OpenFOAM (version 2.3.1) in conjunction with the Simulator fOr Wind Farm Applications (SOWFA) project libraries (Churchfield et al. 2012a) developed by the US National Renewable Energy Laboratory (NREL). The atmospheric solver used in SOWFA is called ABLSolver, which is a transient solver for turbulent flows of incompressible fluids and considers the Boussinesq approximation for buoyancy effects (Churchfield et al. 2012a).

A precursor–successor approach has been utilised to develop a conventionally neutral ABL flow for the simulation. The Coriolis force is calculated for $\phi =55.52^{\circ }$, which is a representative value for a wind farm in the North Sea (Hansen et al. 2012). In SOWFA, a prescribed streamwise velocity ($\mathcal {U}_h=8$ m s$^{-1}$) and wind direction ($\varphi =270^{\circ }$) at the turbine hub-height level can be achieved by adjusting the magnitude and direction of the driving pressure gradient. A capping inversion with a lapse rate of 0.05 K m$^{-1}$ is imposed at the top of the boundary layer covering the heights from $700$ to $800$ m. The height at the bottom of the capping inversion, denoted by $H$, is defined as the thickness of the ABL. The geostrophic layer above the capping inversion has a lapse rate of 0.003 K m$^{-1}$. The inclusion of the capping inversion helps to slow the vertical growth of the boundary layer with time in neutral conditions (Churchfield et al. 2012b). Due to the assumption of the fixed height of the capping inversion, however, the vertical displacement of the flow above the farm does not generate gravity waves in the capping inversion. This simplification may lead to errors in cases where gravity waves induce non-negligible pressure gradients at the turbine hub-height level (Allaerts & Meyers 2017, 2019; Lanzilao & Meyers 2023; Stipa et al. 2023). The precursor simulations are run without the turbines for a period of 10 h (36 000 s) to obtain a quasi-steady state. Next, the inlet conditions are recorded for a period of 9000 s to be fed into the successor simulation with turbines. The farm flow statistics are calculated for the last hour of the simulations. The convective terms are discretised using a second-order central difference scheme for the precursor and a local blend between linear (second-order) and upwind (first-order) schemes for the successor simulation depending on the cell size. This scheme uses 80 % linear and 20 % upwind in proximity of the turbines and 100 % linear in the rest of the domain. For temporal discretisation, an implicit second-order backward scheme is used. For the diffusion term, a Gauss linear second-order scheme is implemented using a non-orthogonality correction for surface-normal gradients. Subgrid-scale stresses are modelled using a one-equation turbulent-viscosity model (Yoshizawa 1986).

The turbine used in this study is NREL-5MW with a hub height of 90 m and a rotor diameter $D$ of 126 m (Jonkman et al. 2009). The turbines are modelled as an actuator disc with no rotation and a constant thrust coefficient of 0.776. This value of $C_T$ was found based on blade element momentum simulations of a turbine rotor with $\mathcal {U}_h$ of 8 m s$^{-1}$ and tip-speed ratio of 7.55 (Navarro Diaz et al. 2022). The body forces are spread across the rotor plane uniformly as axial forces. The equivalent inflow velocity is unknown for turbines that are subject to the wakes of upstream turbines, so a calibration table is used to relate the average velocity on the disc with the unperturbed inflow velocity (van der Laan et al. 2015).

The same domain size is used for both precursor and successor simulations. It extends 1000 m in the vertical direction. The first row of turbines is placed 15 rotor diameters (1890 m) downstream of the inlet, and the domain extends for 357 turbine diameters (approximately 45 km) after the last turbine row. A schematic of the computational domain is presented in figure 2. It is worth noting that the distance between the inlet and the first row of turbines is relatively short in these simulations. Our aim is to maximise the use of our computational resources to capture a very large extent of the farm wake. However, this may lead to an underestimation of the velocity slowdown in the upwind region caused by farm-scale blockage effects, as discussed in the recent parametric study by Lanzilao & Meyers (2023).

Figure 2. Schematic of the computational domain for the A0 case. Turbines are shown by black circles and the rotor diameter is denoted by $D$.

In the precursor simulations, grid cells are 21 m (i.e. $D/6$) long in the streamwise and lateral directions, and their height in the vertical direction grows with distance from the ground (from 2.5 to 60 m at the top of the domain). In the successor simulations including wind turbines, the mesh is refined in two steps. Each refinement halves the cell size. First, in a zone containing the wind farm and its downstream region, the mesh size is reduced to 10.5 m (i.e. $D/12$) in the streamwise and lateral directions. This refined region starts 1260 m (i.e. $10D$) upstream of the first turbine row to ensure that eddy structures are fully developed in the new refined mesh before reaching the wind turbines. In close proximity to the wind turbines, the mesh is further refined by a factor of two (i.e. cell size of 5.25 m or $D/24$ in the streamwise and lateral directions) to capture strong velocity gradients in this region.

Precursor simulations use cyclic boundary conditions at the inlet, outlet and sides. The nearest turbines to the sides are placed such that it resembles the infinite extent of the wind farm in the lateral direction. For instance, with a lateral spacing of $4D$, there is a $2D$ distance from each side as shown in figure 2. At the ground, a wall boundary condition with a prescribed roughness length based on the Schumann–Grotzbach formulation is implemented (Schumann 1975). At the domain top, a slip boundary condition is imposed for velocity and a fixed gradient for temperature. For successor simulations including wind turbines, the inlet uses the data from the precursor while a zero-gradient condition is applied at the outlet. The domain's sides, lower and upper parts have cyclic, wall and slip boundary conditions, respectively. In total, five simulations were performed to study the effect of wind farm layout, wind farm length, inter-turbine spacing and the incoming turbulence level on farm wake flows. The details of these simulations are summarised in table 1. In this table, $s_x$ and $s_y$ are, respectively, streamwise and spanwise inter-turbine spacing, normalised by the rotor diameter $D$. The surface roughness normalised by $D$ is shown by $z_{0,0}$. The number of turbine rows is denoted by $N$, and $u_*$ is the friction velocity normalised by $\mathcal {U}_h$. The incoming turbulence intensity at the hub height is shown by $I_0$, and $\Delta \varphi$ is the change in the incoming wind direction across the turbine rotor (i.e. from the bottom-tip height to the top-tip height). As seen in table 1, the first four cases (A0, S0, AS, AD) are subject to a smooth boundary layer with low surface roughness, whereas the incoming boundary layer in the AR case has a higher surface roughness. The two different inflow boundary-layer profiles are shown in figure 3. The instantaneous streamwise velocity field $u$ for a portion of the Aligned Baseline (A0) case is also shown in figure 4, where the highly turbulent nature of the atmospheric flow and low-speed wakes are clearly visible.

Table 1. Summary of LES for different semi-infinite wind farms. Here $s_x$ and $s_y$ are, respectively, streamwise and spanwise inter-turbine spacing, normalised by the rotor diameter $D$. The surface roughness normalised by $D$ is shown by $z_{0,0}$. The number of turbine rows is denoted by $N$ and $u_*$ is the friction velocity normalised by $\mathcal {U}_h$. The incoming turbulence intensity at the hub height is shown by $I_0$ and $\Delta \varphi$ is the change in the incoming wind direction across the turbine rotor (i.e. from the bottom-tip height to the top-tip height).

Figure 3. Spanwise-averaged vertical profiles of inflow characteristics obtained from precursor simulations. (a) The normalised streamwise velocity $\langle \bar {u}\rangle$, (b) the wind direction $\langle \varphi \rangle$, (c) the incoming turbulence intensity $I = \sigma _u/\mathcal {U}_h$, where $\sigma _u$ is the standard deviation of streamwise turbulent fluctuations, and (d) the horizontal turbulent shear stress defined as $\langle \overline {u_iu_j}\rangle _h=\sqrt {\langle \overline {u'w'}\rangle ^2+\langle \overline {u'v'}\rangle ^2}$. Horizontal dashed and dotted lines respectively indicate the turbine hub height and vertical positions of top/bottom blade tips.

Figure 4. Contours of instantaneous normalised streamwise velocity $u$ for the Aligned Baseline (A0) case. Turbines are shown by black circles.

4. Momentum budget analysis

In this section, the LES data for the Aligned Baseline (A0) case are employed to perform budget analysis on the laterally averaged momentum equations (2.5). This analysis determines dominant terms in the momentum equations both within and downwind of the wind farm. This serves as a basis for the development of the physics-based model later in § 5. Note that viscous terms in the momentum equations were found to be multiple orders of magnitude smaller than any other term, so they are neglected in this analysis. Moreover, regions immediately downstream and upstream of the turbine rows are removed due to steep flow gradients in these regions. This analysis investigates all terms in (2.5), except for the turbine forcing which is only relevant at the rotor disc.

4.1. Streamwise momentum equation

Writing (2.5) in the streamwise direction (i.e. $i=1$) and neglecting turbine forcing gives

(4.1)\begin{align} \underbrace{\langle\bar{u}\rangle\frac{\partial\langle\bar{u}\rangle}{\partial x}}_{{A1x}} + \underbrace{\langle\bar{w}\rangle \frac{\partial\langle\bar{u}\rangle}{\partial z}}_{{A2x}} &={-}\underbrace{f_c(\langle\overline{v_g}\rangle - \langle\bar{v}\rangle)\vphantom{\frac{\partial}{\partial}}}_{{Cx}} - \underbrace{\frac{\partial\langle\bar{u}''\bar{u}''\rangle}{\partial x}}_{{D1x}}\nonumber\\ &\quad - \underbrace{\frac{\partial\langle\bar{u}''\bar{w}''\rangle}{\partial z}}_{{D2x}} - \underbrace{\frac{\partial\langle\overline{u'u'}\rangle}{\partial x}}_{{R1x}} - \underbrace{\frac{\partial\langle\overline{u'w'}\rangle}{\partial z}}_{{R2x}} - \underbrace{\frac{\partial \langle \bar{p}_t\rangle}{\partial x}}_{{Px}}, \end{align}

where $\{u,v,w\}$ are respectively velocities in the streamwise, spanwise and vertical directions. The variation of all terms in (4.1) with respect to $x$ is illustrated in figure 5 until 150 rotor diameters downstream of the wind farm, beyond which limited change is observed in the flow quantities. As shown in figure 5 and as expected, the residual term is mostly negligible throughout the entire domain.

Figure 5. Non-dimensionalised streamwise momentum (4.1) budget at hub height for the Aligned Baseline (A0) case. Locations of turbine rows are denoted by vertical dotted lines. All variables are non-dimensionalised using a selection of $\mathcal {U}_h$ and $D$.

First, we start with the dominant streamwise advection term $A1x=\langle \bar {u}\rangle {\partial \langle \bar {u}\rangle }/{\partial x}$ representing the advection of streamwise momentum by the streamwise velocity. A positive value for $A1x$ means wake recovery/flow acceleration, and vice versa. Approaching the farm, $A1x$ becomes negative approximately $8D$ upstream of the farm. This is explained by the presence of an induction region preceding the farm that is caused by farm-scale blockage effects (Bleeg et al. 2018). Behind the first row of turbines, we observe that the flow acceleration is still suppressed by farm-scale blockage effects, and $A1x$ continues to be negative until about $3D$, where the maximum velocity deficit occurs (i.e. $A1x=0$). The induction entrance region of the wind farm can be also illustrated by positive vertical velocity $\langle \bar {w}\rangle$ shown in figure 6(b).

Figure 6. Variation of laterally averaged (a) streamwise $\langle \bar {u}\rangle$ and (b) spanwise $\langle \bar {v}\rangle$ and vertical $\langle \bar {w}\rangle$ velocities with $x$. All variables are non-dimensionalised using a selection of $\mathcal {U}_h$ and $D$. The locations of turbine rows are denoted by vertical dotted lines.

After the second turbine row, $A1x$ becomes immediately positive indicating that the maximum velocity deficit occurs much closer to the turbine, which is followed by flow acceleration (i.e. wake recovery). By inspection, the profile of vertical Reynolds stress gradient $R2x={\partial \langle \overline {u'w'}\rangle }/{\partial z}$ follows the profile of $A1x$, confirming that $R2x$ acts to replenish wake momentum. In other words, peak flow acceleration (i.e. maximum $A1x$) is observed, approximately where vertical momentum transport due to turbulence is also maximum. Note that terms on the right-hand side of (4.1) that are negative in figure 6 promote wake recovery and vice versa. The greater proportions of turbulent vertical momentum transport in later rows is also evident in figure 5, occurring due to increased flow shear from greater velocity deficits, as observed in figure 6(a). It is worth reminding ourselves that according to (4.1), the gradient of the Reynolds stress is responsible for wake recovery, as opposed to the common assumption that the mere presence of Reynolds stress promotes wake recovery (van der Laan, Baungaard & Kelly 2022).

Figure 5 also shows that the other advection term $A2x=\langle \bar {w}\rangle {\partial \langle \bar {u}\rangle }/{\partial z}$ has minimal impact on momentum transport within the domain. Moreover, although not discernible in figure 5, the Coriolis term $Cx=f_c(\langle \overline {v_g}\rangle - \langle \bar {v}\rangle )$ is negative across the domain as expected in the northern hemisphere, but like $A2x$, its value is negligible compared with the dominant terms.

The normal Reynolds stress gradient $R1x={\partial \langle \overline {u'u'}\rangle }/\partial x$ is illustrative of the rate of change of turbulence level (intensity) with $x$. Term $R1x$ in figure 5 indicates the turbulence level increases behind each turbine row, peaking about $3D$–$5D$ downstream (i.e. where $R1x=0$), before decreasing on the approach to the next row. This is in agreement with prior studies observing peak turbulence intensity occurring a few rotor diameters downstream of an individual turbine (Wu & Porté-Agel 2012). In the wake of the farm, $R1x$ quickly approaches zero as the turbulence decays to its background level.

Figure 5 shows that the turbine pressure gradient $Px={\partial \langle \bar {p}_t\rangle }/{\partial x}$ is significant within the entire farm. From actuator disc theory (Manwell, McGowan & Rogers 2010), we know that after the pressure increase upwind of turbines, there is a sudden pressure drop as the turbine extracts energy from the flow (not shown in figure 5). This is followed by a pressure increase as wake recovery occurs. Figure 5 shows the fast recovery of pressure downwind of each turbine row indicated by positive $P1x$. The value of $P1x$ (i.e. rate of pressure increase) decays with $x$ until it increases again due to the induction region of subsequent rows. It is worth noting that the variation of pressure is often neglected in wake models, but this figure and other recent studies (e.g. Bastankhah et al. 2021) highlights the importance of this term. According to figure 5, within the wind farm, this term is even comparable to other dominant terms (e.g. $A1x$ and $R2x$) in the momentum equation. The term $Px$, however, decays quickly in the wake of the wind farm as the pressure approaches its free-stream value.

Some of the most significant terms in figure 5 are the dispersive stress terms which are the product of deviations from the lateral averaging, and described as the tortuous streamlines induced by flow obstacles (Moltchanov et al. 2011). Dispersive stresses are correlated to obstacle density. Sparsely populated obstacle fields display increased dispersive stresses, due to greater disparity between flow over the obstacles and the mean flow (Moltchanov et al. 2011). For instance, dispersive stresses are expected to be greater in aligned wind farms (shown in figure 5) than in staggered ones (not shown here) although $D1x$ may be still considerable within a staggered wind farm. The term $D1x$ represents the amount of streamwise momentum transport caused by flow inhomogeneity. More precisely, it represents the rate of change in inhomogeneity with respect to the laterally averaged flow. According to figure 5, $D1x$ is negative within the wind farm indicating the homogeneity of the flow is increasing. This occurs as the wake recovers due to vertical momentum transport into the wake, reducing the magnitude of the spatial fluctuations of streamwise velocity (i.e. $\bar {u}''$). Accordingly, the location of minimum $D1x$ (i.e. maximum rate of approaching homogeneity) is correlated with where maximum vertical momentum transport $R2x$ and maximum wake recovery rate $A1x$ occur. As mixing increases so does the homogeneity, causing the reduction of the magnitude of $D1x$ displayed in figure 5. Despite its importance within the farm, $D1x$ decays sharply in the wake of the wind farm, where individual turbine wakes merge and form a holistic farm wake. This is discussed in more detail in § 6. Finally, we investigate the variation of $D2x={\partial \langle \bar {u}''\bar {w}''\rangle }/{\partial z}$. This term essentially quantifies the vertical transfer of streamwise momentum caused directly by the non-uniformity of the time-averaged wind farm flow field. Figure 5 shows that this term is clearly smaller than the other dispersive term $D1x$. It is also interesting to note that this term is mainly positive within the wind farm. This indicates that $D2x$ acts against wake recovery, in contrast to its turbulent counterpart $R2x$.

4.2. Spanwise momentum equation

Even though the terms of the spanwise momentum equation are of smaller magnitude than their streamwise counterparts, they are examined here due to their importance to the wake's trajectory. Writing (2.5) in the spanwise direction (i.e. $i=2$) yields

(4.2)\begin{align} \underbrace{\langle\bar{u}\rangle\frac{\partial\langle\bar{v}\rangle}{\partial x}}_{{A1y}} + \underbrace{\langle\bar{w}\rangle\frac{\partial\langle\bar{v}\rangle}{\partial z}}_{{A2y}} &={-}\underbrace{f_c(\langle\bar{u}\rangle-\langle\overline{u_g}\rangle) \vphantom{\frac{\partial}{\partial}}}_{{Cy}} - \underbrace{\frac{\partial \langle\bar{v}''\bar{u}''\rangle}{\partial x}}_{{D1y}} - \underbrace{\frac{\partial\langle\bar{v}''\bar{w}''\rangle}{\partial z}}_{{D2y}}\nonumber\\ &\quad -\underbrace{\frac{\partial\langle\overline{v'u'}\rangle}{\partial x}}_{{R1y}} - \underbrace{\frac{\partial\langle\overline{v'w'}\rangle}{\partial z}}_{{R2y}}. \end{align}

Variations of terms in (4.2) with the streamwise distance $x$ are shown in figure 7. Due to their small values, the terms in (4.2) are more prone to numerical errors, which may explain why their variations are more oscillatory and less smooth compared with their streamwise counterparts.

Figure 7. Non-dimensionalised spanwise momentum (4.2) budget at hub height for the Aligned Baseline (A0) case. Locations of turbine rows are denoted by vertical dotted lines. All variables are non-dimensionalised using a selection of $\mathcal {U}_h$ and $D$.

First, we start with the Coriolis term $Cy=f_c(\langle \bar {u}\rangle -\langle \overline {u_g}\rangle )$. The dominance of this term varies across the domain. From figure 7, $Cy$ is dominantly negative within the farm due to the streamwise flow deceleration, which according to (4.2) causes positive $A1y=\langle \bar {u}\rangle {\partial \langle \bar {v}\rangle }/{\partial x}$ and thereby positive $\langle \bar {v} \rangle$, as observed in figure 6(b). This deflects the wake to the left, which can be described as an anticlockwise flow rotation viewed from the top. In the wind farm wake, figure 7, however, shows that with flow acceleration the strength of the Coriolis force decays. This is where $R2y={\partial \langle \overline {v'w'}\rangle }/{\partial z}$ that is the vertical turbulent entrainment of veered momentum from above becomes more dominant. Consequently, this changes the sign of the advection term $A1y$ to negative as shown in figure 7, and therefore the far wake starts deflecting to the right (i.e. clockwise wake rotation viewed from top). In other words, the term $A2y$ turns the wind at the hub height towards the wind direction at higher altitudes. This ultimately leads to a negative spanwise velocity as depicted in figure 6(b). This interesting phenomenon is discussed in more detail in § 6. It is also noteworthy that the magnitudes of the second advection term $A2y$, the dispersive stress terms $D1y$ and $D2y$ and also $R1y$ are small, especially in the farm wake and can be neglected.

4.3. Approximate form of momentum equations

From the analysis in § 4.1, amongst others the dispersive stress ($D1x$) and pressure ($Px$) terms are evidently not negligible in the streamwise momentum equation, at least within the farm region. However, despite the evident importance of these two terms, the summation of the four terms $D1x$, $D2x$, $R1x$ and $Px$ – which are challenging to model – is rather small. This is illustrated by the dashed light green colour in figure 5. The combined value of these terms is negligible in the wake of the farm. Within the farm, the combined value is not negligible but smaller than the individual dispersive $D1x$ and pressure $Px$ terms. The term $D1x$ is negative, while $Px$ is positive; therefore, to some extent, they cancel each other out. For simplicity, we thus omit these terms from our model developed in § 5. Moreover, as discussed in § 4.2, it is apparent that in the spanwise direction the dominant terms are (i) the spanwise momentum advection by the streamwise velocity $A1y$, (ii) the Coriolis term $Cy$ and (iii) the vertical turbulent transport of veering wind $R2y$. Therefore, the approximate forms of the momentum equations including turbine forcing can be written as

(4.3a)$$\begin{gather} \langle \bar{u}\rangle \frac{\partial \langle \bar{u}\rangle}{\partial x} \approx{-}f_c (\langle \overline{v_{\mathit{g}}}\rangle -\langle \bar{v}\rangle)\vphantom{\frac{\rangle\partial}{\partial}} - \frac{\partial \langle \overline{u'w'}\rangle}{\partial z} -\langle\, \bar{f}_{x} \rangle, \end{gather}$$

(4.3b)$$\begin{gather}\langle \bar{u}\rangle \frac{\partial \langle \bar{v}\rangle}{\partial x} \approx{+}f_c (\langle \overline{u_{\mathit{g}}}\rangle -\langle \bar{u}\rangle)\vphantom{\frac{\rangle\partial}{\partial}} - \frac{\partial \langle \overline{v'w'}\rangle}{\partial z} -\langle\, \bar{f}_{y} \rangle. \end{gather}$$

In the following section, we simplify (4.3) to develop a system of ordinary differential equations that can be solved mathematically for a semi-infinite wind farm.

5. Derivation of physics-based fast-running farm wake model

5.1. Definition of semi-infinite wind farm

We start by assuming a semi-infinite wind farm which is infinite in the lateral direction but has a finite length in the streamwise direction. A right-handed Cartesian coordinate system $\{x,y,z\}$ aligned with the incoming wind at the turbine hub height $z_h$ is adopted such that $x$ is in the direction of the incoming wind, $y$ represents the horizontal direction normal to $x$ and $z$ measures the height from the ground. The total number of wind turbine rows is denoted by $N$. Wind turbines in the $n$th row are abbreviated to WT$_n$s, where the subscript $n=\{1,2,\ldots,N\}$ shows the row number labelled based on the streamwise position (i.e. ranging from $n=1$ for the first row to $n=N$ for the last row). The WT$_n$s are assumed to have the same values of thrust coefficient $C_{T,n}$ and yaw angle $\gamma _n$, which may, however, be different from $C_{T,m}$ and $\gamma _m$ if $n\neq m$. In other words, turbines in different rows may have different operating conditions. While the lateral spacing $s_y$ between turbines is assumed to be the same for all rows, the streamwise spacing between consecutive rows may be variable. The arbitrary streamwise positions of turbine rows is quite advantageous, as for instance the developed model can be even applied to a cluster of wind farms at once (not done in this study). Furthermore, turbine rows may have lateral offset with respect to each other as shown in figure 8. The lateral position of WT$_n$s is $y_n+ks_y$, where $k=-\infty \dots \infty$.

Figure 8. Schematic of a semi-infinite wind farm. Turbines in row $n$ are denoted by WT$_n$. The lateral spacing $s_y$ between turbines is assumed to be constant for the whole wind farm. However, the streamwise spacing can be variable, and moreover turbine rows may be laterally shifted with respect to each other.

5.2. Turbine force

The normalised aerodynamic force $\bar {f}$ exerted by wind turbines is given by

(5.1a)$$\begin{gather} \bar{f}_x ={-}\sum_{n=1}^N\sum_{k={-}\infty}^\infty\frac{1}{2} C_{T,n}\bar{u}_{h,n}^2\cos(\gamma_n)\delta(x-x_n)H(0.5^2-[(y-y_n-ks_y)^2+(z-z_h)^2]), \end{gather}$$

(5.1b)$$\begin{gather}\bar{f}_y={-}\sum_{n=1}^N\sum_{k={-}\infty}^\infty\frac{1}{2} C_{T,n}\bar{u}_{h,n}^2\sin(\gamma_n)\delta(x-x_n)H(0.5^2-[(y-y_n-ks_y)^2+(z-z_h)^2]), \end{gather}$$

where $\bar {u}_{h,n}$ is the local incoming hub-height velocity for WT$_n$s. In other words, it is the time-averaged velocity at the location of WT$_n$s’ rotor centre (i.e. $(x,y,z)=(x_n,y_n,z_h)$) in the absence of WT$_n$s. Here $\delta (x)$ is the Dirac delta function and $H(x)$ is the Heaviside step function, which is defined as $H(x)=0$ for $x\leq 0$ and $H(x)=1$ for $x>0$. We then apply the lateral averaging discussed in § 2 to obtain laterally averaged turbine forces at the hub height $z=z_h$ as follows:

(5.2a)$$\begin{gather} \langle\,\bar{f}_x\rangle ={-}\frac{1}{2s_y}\sum_{n=1}^N C_{T,n}\bar{u}_{h,n}^2\cos(\gamma_n)\delta(x-x_n), \end{gather}$$

(5.2b)$$\begin{gather}\langle\,\bar{f}_y\rangle ={-}\frac{1}{2s_y}\sum_{n=1}^N C_{T,n}\bar{u}_{h,n}^2\sin(\gamma_n)\delta(x-x_n). \end{gather}$$

5.3. Simplifying Reynolds shear stress terms in (4.3)

The Boussinesq eddy-viscosity hypothesis has been used in previous studies (Belcher, Jerram & Hunt 2003; Luzzatto-Fegiz & Colm-cille 2018) to model spatially averaged Reynolds stresses based on spatially averaged velocity gradients:

(5.3)\begin{equation} \langle \overline{u_i'u_j'}\rangle ={-}2\nu_t \left[ \underbrace{\frac{1}{2}\left(\frac{\partial \langle \overline{u_i}\rangle}{\partial x_j} + \frac{\partial \langle \overline{u_j}\rangle}{\partial x_i}\right)}_{{S_{ij}}}\right]\quad (\textrm{for}\ i\neq j), \end{equation}

where $S_{ij}$ is the dimensionless rate of strain tensor and $\nu _t$ is the turbulent viscosity non-dimensionalised by $\mathcal {U}_hD$. For $i=1$ and $j\neq 1$, $\partial u_i/\partial x_j$ is an order of magnitude larger than $\partial u_j/\partial x_i$ (Tennekes & Lumley 1972), so one can write

(5.4a,b)\begin{equation} \frac{\partial \langle \overline{u'w'}\rangle}{ \partial z} \approx{-}\nu_t \frac{\partial^2 \langle \bar{u}\rangle}{\partial z^2},\quad \frac{\partial \langle \overline{v'w'}\rangle}{\partial z} \approx{-}\nu_t \frac{\partial^2 \langle \bar{v}\rangle}{\partial z^2}. \end{equation}

We first assess the validity of the turbulent-viscosity hypothesis using our LES data. To do so, at a given streamwise position, we examine whether $-\langle \overline {u'w'}\rangle$ is linearly proportional to $\partial \langle \bar {u}\rangle /\partial z$, according to (5.3). Figure 9 shows values of $-\langle \overline {u'w'}\rangle$ and $\partial \langle \bar {u}\rangle /\partial z$ at different heights and streamwise positions. Results are shown for the A0 case, but they look qualitatively similar in other cases (not shown here). Data are plotted for a range of $z=0.25$ (shown in light blue) to $z=4.57$ (shown in magenta). Results in figure 9 are shown for four streamwise locations, with the first two in the farm region and the other two in the wake region: between WT$_3$ and WT$_4$ ($x=2.5s_x$) (figure 9a), between WT$_7$ and WT$_8$ ($x=6.5s_x$) (figure 9b), half of the farm length downstream in the wake ($x=1.5x_N$), where the farm length is $x_N$ (figure 9c) and an entire farm length downstream in the wake ($x=2x_N$) (figure 9d). Lines are fitted to the data at heights above the hub height. Figure 9(a,b) shows that within the farm region, the eddy-viscosity assumption seems to be a valid approximation for the entire vertical domain plotted in the figure, except for the region close to the ground. In the farm wake, however, this assumption is only valid at upper heights as shown in figure 9(c,d). Given that the height at which the eddy-viscosity assumption appears reasonable in the farm wake seems to grow with downstream distance, one possible explanation for this could involve the development of the secondary internal boundary layer (IBL) downwind of the wind farm due to the transition from rough to smooth terrain (see § 5.5 for more discussion on IBLs). However, further research is required to fully elucidate this phenomenon. Nevertheless, even in the farm wake, the turbulent-viscosity hypothesis is still valid at upper heights, where the top-down transport of energy by Reynolds shear stresses occurs, so we use the turbulent-viscosity assumption in this work to simplify the laterally averaged momentum equations.

Figure 9. Distribution of $-\langle \overline {u'w'}\rangle$ against $\partial \langle \bar {u}\rangle /\partial z$ at different heights for four different streamwise positions in the A0 case, where (a,b) are in the farm region and (c,d) are in the wake region, and $x_N=7s_x$ is the farm length.

The turbulent viscosity $\nu _t$ is then decomposed into two parts: one that is due to the ambient atmospheric turbulence denoted by $\nu _{t,0}$ and one shown by $\nu _{t,f}$ corresponding to the turbulence added by the wind farm, and it varies by $x$. In other words,

(5.5)\begin{equation} \nu_t(x)=\nu_{t,0}+\nu_{t,f}(x). \end{equation}

Specifying values of the ambient turbulent viscosity $\nu _{t,0}$ and the farm turbulent viscosity $\nu _{t,f}$ is deferred to § 5.5.

5.4. Mathematical solution of velocity deficit equations

Now, we develop and solve equations for the variation of laterally averaged velocity deficit in both streamwise $U_d$ and spanwise $V_d$ directions, defined as

(5.6)$$\begin{gather} U_d(x)=U_0- U(x), \end{gather}$$

(5.7)$$\begin{gather}V_d(x)=V_0-V(x). \end{gather}$$

In the above equation and hereafter, $U(x)=\langle \bar {u}\rangle (x,z=z_h)$, $V(x)=\langle \bar {v}\rangle (x,z=z_h)$. The subscript $0$ denotes the flow in the absence of the whole wind farm. Let us recall that the coordinate system is defined based on the incoming wind direction at the hub height (see § 5.1) and all velocities in this work are non-dimensionalised by $\mathcal {U}_h$, so $U_0=1$ and $V_0=0$. The latter means that $V_d=-V$ in this work. For generality, however, we write equations for $V_d$ so the model can still be used for a different coordinate system where $V_0\neq 0$.

The magnitude of the local velocity deficit caused by each individual turbine can be significant especially in the near-wake region (e.g. Zhang, Markfort & Porté-Agel 2012; Bastankhah & Porté-Agel 2017). However, laterally averaged values of the velocity deficit are fairly small, as is later shown in § 6. Therefore, we linearise (4.3) by replacing $\langle \bar {u}\rangle$ in the advection term with the incoming velocity $U_0=1$. Moreover, using the turbulent-viscosity hypothesis discussed in § 5.3 to simplify (4.3), we obtain

(5.8a)$$\begin{gather} \frac{\partial U}{\partial x} \approx{-}f_c (V_{\mathit{g}} -V) +\nu_t\frac{\partial^2 U}{\partial z^2} -\langle\, \bar{f}_{x} \rangle, \end{gather}$$

(5.8b)$$\begin{gather}\frac{\partial V}{\partial x} \approx{+}f_c (U_{\mathit{g}} -U) +\nu_t \frac{\partial^2 V}{\partial z^2} -\langle\, \bar{f}_{y} \rangle. \end{gather}$$

In the absence of the wind farm, (5.8) is simplified to

(5.9a)$$\begin{gather} 0 \approx{-}f_c (V_{\mathit{g}} -V_0)\vphantom{\frac{\rangle\partial}{\partial}} +\nu_{t,0} \frac{\partial^2 U_0}{\partial z^2}, \end{gather}$$

(5.9b)$$\begin{gather}0 \approx{+}f_c (U_{\mathit{g}} -U_0) +\nu_{t,0} \frac{\partial^2 V_0}{\partial z^2}. \end{gather}$$

It is worth noting that (5.9) was solved in Ekman (1905) for different altitudes to describe the well-known Ekman spiral. If we subtract (5.8) from (5.9) and also use the dimensional analysis of $\partial ^2 U_d/\partial z^2\propto U_d/l^2$, where $l\approx 1$ (i.e. length scale comparable to rotor diameter), we obtain

(5.10a)$$\begin{gather} \frac{\textrm{d} U_d}{\textrm{d} x} \approx{+} \underbrace{f_cV_d\vphantom{\frac{a}{b}}}_{Coriolis} -\underbrace{c_1 \nu_t U_d \vphantom{\frac{a}{b}}}_{recovery} +\underbrace{\langle\, \bar{f}_{x}\rangle \vphantom{\frac{a}{b}}}_{turbine\ forcing}+ \underbrace{C_{x}\vphantom{\frac{a}{b}}}_{shear}, \end{gather}$$

(5.10b)$$\begin{gather}\frac{\textrm{d} V_d}{\textrm{d} x} \approx{-} \underbrace{f_c U_d\vphantom{\frac{a}{b}}}_{Coriolis} -\underbrace{c_1 \nu_t V_d \vphantom{\frac{V_d}{Ro}}}_{recovery} +\underbrace{\langle\, \bar{f}_{y}\rangle \vphantom{\frac{V_d}{Ro}}}_{turbine\ forcing}+ \underbrace{C_{y}\vphantom{\frac{V_d}{Ro}}}_{veer}, \end{gather}$$

where $c_1$ is a constant. The term $C_{x}=-\nu _{t,f} \partial ^2 U_0/\partial z^2$ is related to the rate of incoming shear, so it is called shear in (5.10). Likewise, $C_{y}=-\nu _{t,f} \partial ^2 V_0/\partial z^2$ is related to the rate of incoming veer, and it is called veer in (5.10). Both $C_x$ and $C_y$ also depend on the amount of turbulence generated by the wind farm through $\nu _{t,f}$. Values of $C_x$ and $C_y$ are determined as a function of atmospheric and farm conditions later in § 5.7, and for now they are assumed to be known.

If we insert (5.2) into (5.10), the mathematical solution presented below for this system of ordinary differential equations can be obtained. The solution written in (5.11) is exact for a constant $\nu _{t,f}$. For a well-defined function of $\nu _{t,f}(x)$ with an arbitrary distribution, this provides an approximate solution with insignificant error with respect to the exact numerical solution (not shown here) of the ordinary differential equation system:

(5.11a)$$\begin{gather} U_d(x)\approx \frac{C_x}{c_1\nu_t }(1-\textrm{e}^{{-}c_1\int_{0}^{x} \nu_t\,\textrm{d} x})H(x)+\sum_{n=1}^N U_{d,n}(x), \end{gather}$$

(5.11b)$$\begin{gather}V_d(x)\approx \frac{C_y}{c_1\nu_t }(1-\textrm{e}^{{-}c_1\int_{0}^{x} \nu_t\,\textrm{d} x})H(x)+\sum_{n=1}^N V_{d,n}(x), \end{gather}$$

where $U_{d,n}(x)$ and $V_{d,n}(x)$ are, respectively, values of the laterally averaged streamwise and spanwise velocity deficit caused by WT$_n$s at $x$, and they are given by

(5.12a)$$\begin{gather} U_{d,n}(x)=\frac{C_{T,n}}{2s_y}(1-\eta_n U_d(x_n))^2\,\textrm{e}^{{-}c_1\int_{x_n}^{x} \nu_t\,\textrm{d} x}\cos(\gamma_n-f_c(x-x_n))H(x-x_n), \end{gather}$$

(5.12b)$$\begin{gather}V_{d,n}(x)=\frac{C_{T,n}}{2s_y}(1-\eta_n U_d(x_n))^2\,\textrm{e}^{{-}c_1\int_{x_n}^{x} \nu_t\,\textrm{d} x}\sin(\gamma_n-f_c(x-x_n))H(x-x_n). \end{gather}$$

In (5.12), $\bar {u}_{h,n}$ is replaced with

(5.13)\begin{equation} \bar{u}_{h,n}=1-\eta_n U_d(x_n), \end{equation}

where $U_d(x_n)=U_d(x=x_n)$, and we call $\eta$ the local deficit coefficient, which is defined as the ratio of the local velocity deficit experienced by wind turbines to the laterally averaged velocity deficit. This coefficient depends on the farm layout configuration, and it is determined by (5.26) developed later in § 5.6. To compute the ambient turbulent viscosity $\nu _{t,0}$ and farm turbulent viscosity $\nu _{t,f}$, (5.14) and (5.19) developed in § 5.5 are respectively used. Finally, values of $C_x$ and $C_y$ are estimated based on (5.30) in § 5.7. Once values of $C_x$, $C_y$, $\nu _{t,0}$ and $\nu _{t,f}$ are all known, a forward marching scheme in the streamwise direction is implemented using (5.11) to find the evolution of $U_d$ and $V_d$ with $x$. The forward marching scheme is stable and not sensitive to the streamwise resolution, which was tested (not shown here) for a range of values from $0.01D$ to $1D$. The solution in (5.11) is versatile as it can also be used for different distributions of $\nu _{t}$, $C_x$ and $C_y$ that can potentially be developed in future studies.

5.5. Estimation of turbulent viscosity $\nu _t=\nu _{t,0}+\nu _{t,f}$

In general, turbulent viscosity can be written as a product of a turbulence velocity scale $\hat {u}$ and a turbulence length scale (i.e. mixing length) $\hat {l}$ (Tennekes & Lumley 1972). Therefore, to estimate the ambient turbulent viscosity $\nu _{t,0}$ and the farm turbulent viscosity $\nu _{t,f}$, one needs to specify suitable values of $\hat {u}$ and $\hat {l}$ for each term.

For the ambient turbulent viscosity $\nu _{t,0}$, according to Prandtl's mixing-length hypothesis (Tennekes & Lumley 1972), we assume $\hat {u}_0=u_*$ and $\hat {l}_0=\kappa z_h$, where $\kappa \approx 0.41$ is the von Kármań constant. Therefore, the ambient turbulent viscosity is given by

(5.14)\begin{equation} \nu_{t,0}=\kappa u_*z_h. \end{equation}

It is worth noting that the assumption of $\hat {l}_0=\kappa z_h$ is expected to be valid only in the log-law region of the ABL (Pope 2000). For our simulations the top of the rotor disc is at a location $z/H = 0.22$ which is seemingly right at the outer limit of this assumption. More sophisticated models for the mixing length in ABLs can be used instead (e.g. van der Laan et al. 2020) but for simplicity we retain the simple log-law relationship for $\hat {l}_0$.

For the farm turbulent viscosity $\nu _{t,f}$, we use the height of the IBL that grows above the wind farm as the turbulence length scale. Several studies have shown that the IBL grows above wind farms following the classical Elliott's $x^{0.8}$ power law (e.g. Allaerts & Meyers 2017; Wu & Porté-Agel 2017). According to Wood (1982), the thickness of the IBL, $\delta$, due to a smooth to rough transition is given by

(5.15)\begin{equation} \frac{\delta}{z_{0,f}}=0.28\left(\frac{x}{z_{0,f}}\right)^{0.8}, \end{equation}

where $z_{0,f}$ is the equivalent roughness length of the wind farm. It is worth noting that the IBL may become separated from the ground surface as a new IBL starts developing due to the rough to smooth transition downwind of the wind farm (Oke 1976). For simplicity, we use the maximum height of the first IBL as the turbulence length scale and do not account for the second IBL development. In order to use (5.15), one needs to estimate the value of $z_{0,f}$ for the wind farm. For an infinite wind farm, several models have been already proposed in the literature to estimate $z_{0,f}$ (e.g. Frandsen 1992; Calaf et al. 2010; Yang et al. 2012; Abkar & Porté-Agel 2013). The one suggested by Frandsen (1992) for an infinite wind farm with uniformly distributed wind turbines states

(5.16)\begin{equation} z_{0,f}=z_h\exp\left(-\dfrac{\kappa}{\sqrt{\dfrac{1}{2}c_{ft} +\left[\dfrac{\kappa}{\ln(z_h/z_{0,0})}\right]^2}}\right), \end{equation}

where $c_{ft}={\rm \pi} C_T/(4s_xs_y)$, $s_x$ is the normalised streamwise spacing between turbine rows and as a reminder $z_{0,0}$ is the normalised surface roughness length in the absence of the wind farm. To maintain simplicity, we use the same relationship to estimate $z_{0,f}$ for the semi-infinite wind farm. To compute $c_{ft}$, we use the average streamwise inter-turbine spacing for $s_x$ (i.e. $s_x=\sum _{n=1}^{N-1}(x_{n+1}-x_{n})/(N-1)$), and the average value of thrust coefficient for $C_T$ (i.e. $C_T=\sum _{n=1}^{N}C_{T,n}/N$).

One can use (5.15) in conjunction with (5.16) to estimate the thickness of the IBL over the wind farm. This relationship is, however, only valid as long as $\delta$ is smaller than the ABL thickness $H$ (Wood 1982). It is known that the IBL growth is capped by the thermal inversion at the top of the ABL (Oke 1976), especially if the inversion layer has strong free-atmosphere stratification. In these cases, the inversion layer may act as a ‘lid’ on the top of the ABL and hinders the growth of the IBL. We therefore artificially limit the growth of the turbulence length scale $\hat {l}_f$ using the following relationship:

(5.17)\begin{equation} \hat{l}_f=\frac{\delta}{1+\delta/H}. \end{equation}

For small values of $x$, $\hat {l}_f\approx \delta$, while $\hat {l}_f\to H$ as the value of $x\to \infty$. Variation of $\hat {l}_f$ with $x$ based on (5.17) is shown in figure 10. For simplicity, we here assume that the value of $H$ is constant. It is, however, important to note that depending on the size of the wind farm and the level of thermal stratification in the inversion layer, the IBL may lead to the growth of the entire ABL, so in reality $H$ might change with $x$ (Allaerts & Meyers 2017; Wu & Porté-Agel 2017).

Figure 10. Variation of velocity scale $\hat {u}_f$, length scale $\hat {l}_f$ and farm turbulent viscosity $\nu _{t,f}$ based on the modelling approached elaborated in § 5.5. Variation of $\nu _{t,f}$ for the LES data (case S0) is also shown. In the figure, $\hat {u}_f$ is normalised by $\sqrt {c_{ft}}$, $\nu _{t,f}$ is normalised by $\nu _{t,f}(x=x_1+L_f)$ and $\hat {l}_f$ is normalised by $H$. Vertical dotted lines show the location of turbine rows, and the vertical dashed line shows $x_1+L_f$ where the turbulent viscosity $\nu _{t,f}$ is maximum.

Next, we determine the turbulence velocity scale $\hat {u}_f$. Within the wind farm, turbulence is mainly generated due to the shear caused by turbine forcing. Therefore, inspired by (Calaf et al. 2010), we use $\sqrt {c_{ft}}U_0$, where $U_0=1$, to estimate the turbulence velocity scale $\hat {u}_f$ within the wind farm. The turbulence generation is expected to peak at some distances downstream of the wind farm, which is then followed by a decay in turbulence generation due to wake recovery and reduction of flow shear (Stieren & Stevens 2022). The generated turbulence in the turbine wake usually peaks at around $5$ rotor diameters downstream (Chamorro & Porté-Agel 2009). The constant turbulence velocity scale $\hat {u}_f$ caused by turbine forcing is therefore assumed to be extended from $x_1$ to $x_N+5$ as shown in figure 10. Further downstream, $\hat {u}_f$ starts to decay due to the wake recovery. The wake of a semi-infinite wind farm can be modelled as a two-dimensional wake of a canopy of roughness elements in a turbulent boundary layer. According to Belcher et al. (2003), the velocity scale defined based on the maximum velocity deficit in this type of wake flows decays with $x^{-1}$. So in summary, we model the turbulence velocity scale $\hat {u}_f$ as

(5.18)\begin{equation} \hat{u}_f =\begin{cases} 0 & \text{if } x< x_1,\\ \sqrt{c_{ft}} & \text{if } x_1< x< x_1+L_f,\\ \sqrt{c_{ft}}\left(\dfrac{x-x_1}{L_f}\right)^{{-}1} & \text{if } x>x_1+L_f, \end{cases} \end{equation}

where $L_f=(x_N-x_1)+5$. Variations of $\hat {u}_f$ and $\nu _{t,f}$ are shown in figure 10, where the farm turbulent viscosity $\nu _{t,f}$ is given by

(5.19)\begin{equation} \nu_{t,f}=c_2\hat{u}_f\hat{l}_f, \end{equation}

with $c_2$ a constant. As can be seen in figure 10, the farm turbulent viscosity increases within the farm until it reaches its maximum value five rotor diameters downstream of the wind farm. It then decreases further downstream until it eventually becomes zero very far downstream. In general, this is in fairly good agreement with our LES data. As an example, the variation of $\nu _{t,f}$ for the Staggered Baseline (S0) case is shown in figure 10. Fairly similar variations of turbulent viscosity have been also reported in the literature for the wake of a single wind turbine (Scott et al. 2023).

5.6. Estimation of local deficit coefficient $\eta _n$

As discussed earlier, $\bar {u}_{h,n}=1-\eta _nU_d(x=x_n)$, where $\eta _n$ is the local deficit coefficient for WT$_n$s. A consequence of linearising the momentum equation (5.8) is that wake effects are linearly superposed (Bastankhah et al. 2021). Therefore, one can write

(5.20)\begin{equation} \eta_{n} =\frac{\displaystyle\sum_{m=1}^{n-1} \bar{u}_{d,m}(x_n,y_n,z_h)}{\displaystyle\sum_{m=1}^{n-1} U_{d,m}(x_n)}, \end{equation}

where $\bar {u}_{d,m}(x,y,z)$ is the local velocity deficit at $(x,y,z)$ caused by WT$_m$s and $U_{d,m}(x)$ is already defined in (5.12). To find $\bar {u}_{d,mn}$, one can use a Gaussian distribution $C\exp (-r^2/2\sigma ^2)$ to express the velocity deficit caused by each turbine, where $C$ is the wake-centre velocity deficit, $r$ is the distance from the wake centre and $\sigma$ is the characteristic wake width. Therefore, for a semi-infinite wind farm, we obtain

(5.21)\begin{equation} \bar{u}_{d,m}(x_n)=\lim_{K\to\infty}C_{mn}\sum_{k={-}K}^{K}\exp \left(-\frac{(y_n-y_m-ks)^2}{2\sigma_{mn}^2}\right), \end{equation}

where $C_{mn}$ is the maximum velocity deficit caused by each WT$_m$ at $x=x_n$, $k$ denotes the column number which varies from $-\infty$ to $\infty$ and $\sigma _{mn}$ is the width of WT$_m$ wakes at $x=x_n$. Solving (5.21) gives

(5.22)\begin{equation} \bar{u}_{d,m}(x_n)=\sqrt{2{\rm \pi}}C_{mn}\frac{\sigma_{mn}}{s_y}\vartheta_{3}\left[{\frac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left({-}2{\rm \pi}^2\frac{\sigma_{mn}^2}{s_y^2}\right)}\right], \end{equation}

where $\vartheta _{3}[z,q]$ is the Jacobi theta function defined as $1+2\sum _{l=1}^{\infty }q^{l^2}\cos (2lz)$ (Whittaker & Watson 2020). High-level programming languages (e.g. Python) may have a built-in function to compute the Jacobi theta function (mpmath development team 2023). The series describing the Jacobi theta function converges rather quickly, so as an approximation, one can alternatively compute the summation of the first few terms in the series. It is also noteworthy that $\vartheta _{3}[z\pm {\rm \pi},q]=\vartheta _{3}[z,q]$, so $y_n$ and $y_m$ in (5.22) can be the spanwise location of any turbines in WT$_n$ and WT$_m$, respectively.

By definition $U_{d,m}(x)=\langle \bar {u}\rangle _{d,m}(x)$, so we can write

(5.23)\begin{align} U_{d,m}(x_n)=\lim_{K\to\infty}\frac{C_{mn}}{Ms_y}\sum_{k={-}K}^{K}\int_{-\infty}^{+\infty}\exp \left(-\frac{(y-y_m-ks_y)^2}{2\sigma^2_{mn}}\right)\textrm{d}\kern0.05em y=\sqrt{2{\rm \pi}}C_{mn}\frac{\sigma_{mn}}{s_y}. \end{align}

Therefore, from (5.22) and (5.23), we conclude that

(5.24)\begin{equation} \bar{u}_{d,m}(x_n)=U_{d,m}(x_n)\vartheta_{3}\left[{\frac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left({-}2{\rm \pi}^2\frac{\sigma_{mn}^2}{s_y^2}\right)}\right]. \end{equation}

The wake width $\sigma _{mn}$ can be simplified by $k^{*}(x_n-x_m)+\sigma _0$, where the wake expansion $k^{*}$ is typically around 0.02–0.04 for offshore conditions and $\sigma _0$ is around 0.2–0.3 (Bastankhah & Porté-Agel 2014). For simplicity, we assume $k^{*}=0.025$ and $\sigma _0=0.25$ toreduce (5.24) to

(5.25)\begin{equation} \bar{u}_{d,m}(x_n)\approx U_{d,m}(x_n)\vartheta_{3}\left[{\frac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left(-\frac{(x_n-x_m+10)^2}{80s_y^2}\right)}\right]. \end{equation}

Therefore from (5.20) and (5.25), we conclude

(5.26)\begin{equation} \eta_{n} =\frac{\displaystyle\sum_{m=1}^{n-1} U_{d,m}(x_n)\vartheta_{3}\left[{\dfrac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left(-\dfrac{(x_n-x_m+10)^2}{80s_y^2}\right)}\right]} {\displaystyle\sum_{m=1}^{n-1}U_{d,m}(x_n)}, \end{equation}

where $n>1$ (for $n=1$, $\eta$ is zero) and $U_{d,m}(x)$ is defined in (5.12). Figure 11 shows the variation of the local deficit coefficient $\eta _n$ with the row number $n$ for both the Aligned Baseline (A0) and the Staggered Baseline (B0) wind farms. It is interesting to note that the value of $\eta _n$ is always greater than one for an aligned wind farm. This is expected as in this case turbines operate in full-waked conditions. Therefore, the local velocity deficit experienced by turbines is expected to always be larger than the laterally averaged velocity deficit. The maximum value of $\eta$ occurs at the second row where there is the maximum level of heterogeneity in the farm. Further downstream, due to the wake expansion and flow mixing, the value of $\eta _n$ decays and approaches a constant value of around two for this particular wind farm layout. For the staggered wind farm, however, we observe a very different behaviour. For the second row, the value of $\eta _n$ is almost zero. This is due to the fact that WT$_2$s are not affected by the wake of WT$_1$s due to the staggered layout of the wind farm. Wake interactions only occur from the third row where there is a sudden increase in the value of $\eta _n$. The value of $\eta _n$ for the staggered wind farm approaches one for WT$_4$ and downwind rows. This suggests a fairly uniform distribution of streamwise velocity deep inside a staggered wind farm. In figure 11(b), we can see how different distributions of the local velocity coefficient $\eta _n$ lead to different distribution of local hub-height velocity $\bar {u}_{h,n}$. The local hub-height velocity is clearly higher for the staggered wind farm layout which leads to more power generation. The trend shown in figure 11(b) is similar to the data reported in other works (e.g. Chamorro, Arndt & Sotiropoulos 2011; Stieren & Stevens 2022). We have only discussed aligned and staggered layouts here, but an important feature of our developed model is that it is generalisable, through variable local deficit coefficient $\eta _n$ developed in (5.26), to any conceivable wind farm layout provided the layout fulfils the requirements defined in § 5.1.

Figure 11. Variation of (a) local deficit coefficient $\eta _n$ and (b) local hub-height velocity $\bar {u}_{h,n}$ with row number $n$ for the Staggered Baseline (S0) and the Aligned Baseline (A0) cases based on (5.26).

5.7. Estimation of shear term $C_x$ and veer term $C_y$

As discussed earlier in (5.10), $C_{x}=-\nu _{t,f} \partial ^2 U_0/\partial z^2$ and $C_{y}=-\nu _{t,f} \partial ^2 V_0/\partial z^2$. Based on (5.9), $C_{x}$ and $C_{y}$ can be written as

(5.27a)$$\begin{gather} C_{x}={-}\frac{\nu_{t,f}}{\nu_{t,0}}f_c(V_g-V_0), \end{gather}$$

(5.27b)$$\begin{gather}C_{y}=\frac{\nu_{t,f}}{\nu_{t,0}}f_c(U_g-U_0). \end{gather}$$

As a first-order approximation, the vertical changes in the streamwise and spanwise velocities across the ABL are proportional to $G\cos \theta _0$ and $G\sin \theta _0$, respectively. Here, $G=\sqrt {U_g^2+V_g^2}$ is the geostrophic wind speed, and the cross-isobar angle $\theta _0$ is the angle between the wind direction on the ground surface and the geostrophic wind direction. One may thus write

(5.28a)$$\begin{gather} U_g-U_0=c_3G\cos\theta_0, \end{gather}$$

(5.28b)$$\begin{gather}V_g-V_0=c_3 G\sin\theta_0, \end{gather}$$

where $c_3$ is a constant. Values of $G$ and $\theta _0$ can be obtained from the widely used geostrophic drag law, which relates surface properties (e.g. $z_0$ and $u_*$) to geostrophic wind speed on top of the ABL (e.g. Blackadar & Tennekes 1968; Hess & Garratt 2002; van der Laan et al. 2020). For a neutrally stratified ABL, the geostrophic drag law reads as (Liu, Gadde & Stevens 2021)

(5.29a)$$\begin{gather} A=\ln\left(\frac{u_*}{z_{0,0}|f_c|}\right)-\frac{\kappa G}{u_*}\cos\theta_0, \end{gather}$$

(5.29b)$$\begin{gather}B={\mp} \frac{\kappa G}{u_*}\sin\theta_0, \end{gather}$$

where $A$ and $B$ are universal empirical constants, and values of $A = 1.8$ and $B= 4.5$ used in the Wind Atlas Analysis and Application Program (Floors, Troen & Kelly 2018) are adopted here. The minus sign on the right-hand side of the second equality in (5.29) relates to the northern hemisphere where $f_c$ is positive, whereas the positive sign is for the southern hemisphere where $f_c$ is negative (Liu et al. 2021). From (5.27), (5.28) and (5.29), we can therefore approximate the values of $C_x$ and $C_y$ as follows:

(5.30a)$$\begin{gather} C_{x}=c_3|f_c|\frac{\nu_{t,f}}{\nu_{t,0}}\frac{u_*}{\kappa}B, \end{gather}$$

(5.30b)$$\begin{gather}C_{y}=c_3f_c\frac{\nu_{t,f}}{\nu_{t,0}}\frac{u_*}{\kappa}\left[\ln\left(\frac{u_*}{z_{0,0}|f_c|}\right)-A\right]. \end{gather}$$

While for simplicity the conventional geostrophic drag law relationship is used here to estimate $C_x$ and $C_y$, recent works (e.g. Liu & Stevens 2022; Narasimhan, Gayme & Meneveau 2023) that describe the structure of Ekman boundary layer flows can be implemented in future works.

6. Results and discussion

In this section, we compare predictions of the model developed in § 5 with the LES data elaborated in § 3. The values of model coefficients, namely $c_1$, $c_2$ and $c_3$, used in this study are listed in table 2. It is worth recalling that $c_1$ is the coefficient that is multiplied to $\nu _t=\nu _{t,0}+\nu _{t,f}$ in the final solution, while $c_2$ is the coefficient within $\nu _{t,f}$. An increase of either $c_1$ or $c_2$ increases the wake recovery rate, but $c_2$ is the one that quantifies the importance of $\nu _{t,f}$ over $\nu _{t,0}$. On the other hand, $c_3$ is the coefficient of $C_x$ (i.e. shear term in the streamwise momentum equation) and $C_y$ (i.e. veer term in the spanwise momentum equation). Given that $C_x$ is fairly small compared with other terms in the streamwise momentum equation, the main impact of $c_3$ is on predictions of the spanwise velocity deficit. The reported coefficients were found manually based on comparing model outputs with the LES data, but a more systematic approach based on an optimisation algorithm might provide more suitable coefficients. We must also note that although the coefficients suggested in table 2 provide satisfactory predictions for several cases studied here, future research is indeed required to examine whether these values are universal.

Table 2. Empirical model coefficients used in this study.

First, we discuss the laterally averaged velocity deficit for the Aligned Baseline (A0) case shown in figure 12. For the streamwise direction, the figure shows a sudden jump in velocity deficit at each row, which is due to the turbine thrust force (i.e. $\langle \bar {f}_x\rangle$ in (5.10)). The model underpredicts the velocity deficit increase for the first row. This is likely due to the lack of modelling farm-scale blockage effects. The budget analysis in § 4 showed a considerable flow deceleration in the vicinity of WT$_1$s due to farm-scale blockage effects. Both LES data and model predictions suggest that the velocity deficit jump due to turbine thrust forcing is significantly reduced after the first row. This is mainly due to the fact that as shown in (5.2), the thrust force is proportional to the square of the local hub-height velocity, which is clearly lower for subsequent turbine rows.

Figure 12. Variation of laterally averaged (a) streamwise velocity deficit $U_d$ and (b) spanwise velocity deficit $V_d$ at the turbine hub height for the Aligned Baseline (A0) and the Staggered Baseline (S0) cases. The dashed lines show the LES data and the solid lines show predictions of the developed model. Vertical dotted lines show the locations of wind turbine rows.

After the last row of turbines, the velocity deficit diminishes rather rapidly in the farm near-wake region (e.g. for $x=50\unicode{x2013}100$). The primary reason for this fast recovery in the farm near-wake region is the large value of Reynolds stress gradient ($\partial \langle \overline {u'w'} \rangle \partial z$) as shown in figure 5. The dispersive stress gradient (i.e. ${\partial \langle \bar {u}''\bar {u}''\rangle }/{\partial x}$) is large immediately behind the wind farm, but it decays rapidly. As discussed previously, the dispersive stress quantifies the level of inhomogeneity in the flow. Right behind the wind farm, turbine wakes are not completely merged yet. This creates large velocity gradients over small length scales which in turn lead to higher turbulence generation. Further downstream, the turbine wakes merge and form a single farm wake. This is where the dispersive stress becomes negligible, and the gradient of Reynolds shear stress becomes the sole mechanism for the wake recovery. Merging of individual turbine wakes to form a single farm wake is evident in figure 13 that shows contours of the time-averaged streamwise velocity $\bar {u}$ at the hub height. This is conceptually similar to the transition of a wake array to a single wake that occurs downwind of multi-rotor turbines discussed in Bastankhah & Abkar (2019).

Figure 13. Contours of time-averaged streamwise velocity $\bar {u}$ on a horizontal plane at the turbine hub height for the Aligned Baseline (A0) case. Vertical dotted lines denote the location of wind turbine rows.

The developed model can capture the fast wake recovery in the farm near wake. In the momentum equation (5.10), the streamwise velocity deficit reduction is mainly caused by the recovery term $c_1\nu _tU_d$. In the farm near wake, the streamwise velocity deficit $U_d$ is still fairly large. Moreover, the farm turbulent viscosity $\nu _{t,f}$ has its maximum value immediately after the wind farm as shown in figure 10. Both of these contribute to a fast wake recovery in the farm near-wake region. In the farm far-wake region (e.g. for $x>150$), the velocity deficit clearly decays at a slower rate. According to the budget analysis in § 4, in this region, the advection term is in balance with the Reynolds stress gradient, whose effect is modelled by the recovery term $c_1\nu _tU_d$ in (5.10). However, both $\nu _t$ and $U_d$ are much smaller in this region which leads to a slower wake recovery rate. This explains the persistence of the wake over a very long distance. Figure 12 shows that the velocity deficit is still not negligible even after 20 km downstream of the wind farm (i.e. $x\approx 208$). It is also worth noting that for the LES cases studied here, the Coriolis term $f_cV_d$ is an order of magnitude smaller than other terms in the streamwise momentum equation (5.10). Moreover, while the shear term $C_x$ is not negligible, it is much smaller than the recovery term $c_1\nu _tU_d$ for this configuration. However, this is not case for the spanwise momentum equation as elaborated in the following.

The spanwise turbine forcing term $\langle \bar {f}_y\rangle$ defined in (5.2) is zero for the LES data given that the yaw angle $\gamma$ is zero. The Coriolis term in the spanwise momentum equation (5.10), however, is important because it is proportional to $U_d$, which has a considerable value especially within the farm. Therefore, the Coriolis term increases in the farm with an increase of streamwise velocity deficit. According to (5.10), the Coriolis term with $U_d>0$ leads to a negative $V_d$ (i.e. positive $V$) in the northern hemisphere, where $f_c>0$. This is described in the literature (e.g. van der Laan & Sørensen 2017) as an anticlockwise deflection based on a view from above. The initial anticlockwise deflection of the wake is shown in figure 12. Apart from the streamwise wake deficit, the farm-induced turbulence also impacts the distribution of spanwise velocity deficit. In particular, the increase of farm turbulent viscosity $\nu _{t,f}$ has two important effects. It increases both the recovery term $c_1\nu _t V_d$ and the veer term $C_y$ in (5.10). Both of these effects reduce the initial anticlockwise deflection of the wake such that after some downwind distances, the effect of the veer term $C_y$ becomes dominant, and the direction of the wake deflection changes to clockwise (i.e. $V_d>0$). The change in the direction of wake deflection was also reported in Gadde & Stevens (2019). Ultimately, very far downstream (not shown here), $C_y$ goes to zero as $\nu _{t,f}$ goes to zero, and therefore the spanwise velocity deficit is completely diminished by the recovery term. These conflicting behaviours of Coriolis and veer on the spanwise velocity coexist as the latter is the consequence of the former. However, depending on atmospheric conditions and simulation settings, their relative magnitude with respect to each other could be different at different downstream positions. For instance, the succeeding clockwise deflection observed in figure 12 can completely mask the initial anticlockwise deflection. This may happen either in the case of a strong wind veer which typically occurs in thermally stable boundary layers, or if the wind farm generates a high amount of turbulence. On the other hand, an anticlockwise deflection is mainly observed if both incoming wind veer and farm-generated turbulence are relatively small, or if only the farm near wake is of interest. This may explain why some prior works observed a clockwise deflection (e.g. Abkar et al. 2016; van der Laan & Sørensen 2017; Eriksson et al. 2019; Nouri, Vasel-Be-Hagh & Archer 2020), whereas others reported an anticlockwise deflection (e.g. Dörenkämper et al. 2015; Allaerts & Meyers 2017; Frank et al. 2023).

Overall, the agreement of the model predictions with the LES data is satisfactory for the streamwise velocity deficit $U_d$. The developed model can also successfully capture the overall trend for the spanwise velocity deficit $V_d$. In agreement with the LES data, it predicts the initial anticlockwise deflection followed by a clockwise deflection. However, we should note that figure 12 shows some differences in the magnitude of $V_d$ especially in the far wake. The value of $V_d$ is one order of magnitude smaller than $U_d$ for these LES data. While predicting $V_d$ with such small values in these cases can be challenging, the model's prediction for $V_d$ compared with the LES data remains within $1\,\%$ of $\mathcal {U}_h$, and the difference in wind direction $\varphi$ prediction is also within $1^{\circ }$ (not shown here). In this work, we only examined neutral boundary layers where wind veer is not strong. For thermally stable boundary layers, the cross-isobar angle $\theta _0$ in (5.28) is expected to be considerably larger (Peña, Gryning & Floors 2014), which in turn increases the value of $C_y$. It is therefore of great interest to extend the model to thermally stratified boundary layers in future works, where the deflection of the wake due to the wind veer is expected to be significantly larger. The discrepancy in $V_d$ observed between the LES data and the model can be also partly explained by the linearisation of momentum equations. In order to mathematically solve the equations, we replaced $U$ with $U_0=1$ in the advection term in (4.3). In the regions where the difference between $U$ and $U_0$ is not negligible (i.e. within the wind farm), this assumption leads to an underestimation of the velocity deficit. The error introduced by this assumption is expected to be more evident in spanwise velocity predictions given their small values.

Next, we discuss the effect of wind farm layout on the farm flow distribution. The variation of the velocity deficit for the Staggered Baseline (S0) case is shown in figure 12. The main notable difference in the streamwise velocity deficit between the two wind farm layouts (aligned versus staggered) is the fact that, from the second row of turbines, the jump in $U_d$ due to turbine forcing is clearly larger for the staggered wind farm, in agreement with prior studies (e.g. Stevens et al. 2016; Stieren & Stevens 2022). This leads to the maximum velocity deficit $U_d$ of 0.32 for the staggered wind farm, which is 28 % larger than that for the aligned wind farm. As discussed in § 5.6, turbines within the staggered wind farm experience a larger local velocity which according to (5.2) leads to a larger value of turbine thrust force. This highlights the importance of the local deficit coefficient $\eta$ discussed in § 5.6 and implemented in the model (5.11). Without this coefficient, model predictions will be the same for both cases of $A0$ and $S0$, which is clearly unrealistic according to the LES data. The enhanced turbulence mixing caused by large values of velocity deficit in the staggered wind farm accelerates the wake recovery downstream. At about $x=200$, the streamwise velocity deficit becomes approximately equal in the wake of both wind farms. The faster recovery of the wake for the staggered wind farm is captured in the model given that the recovery term in (5.10) depends on the value of $U_d$.

The LES data also show that the spanwise velocity deficit is larger for the staggered wind farm in comparison with the aligned wind farm. As discussed earlier, the Coriolis term in the spanwise momentum equation (5.10) directly depends on $U_d$. Therefore, one expects to observe a larger value of $V_d$ for the staggered wind farm. Due to the initial large anticlockwise wake deflection, the transition to the clockwise deflection occurs at a later downstream position for the staggered wind farm. However, the wake deflection for both wind farms eventually approaches the same value in the very far-wake region. Consistent with the LES data, the model predicts a larger negative peak of the spanwise velocity deficit in this case. The agreement is, however, less satisfactory for the staggered layout. The larger streamwise velocity deficit for the S0 case may exacerbate the error introduced by the linearisation of the momentum equations discussed earlier. Nonetheless, the results are still within a difference of $1\,\%\mathcal {U}_h$.

Next, we study the effect of surface roughness on the evolution of wind farm wakes. We compare the results for the two cases of the Aligned Baseline (A0) with $z_0=2\times 10^{-4}\ \textrm {m}/D$ and the Aligned Rough (AR) with $z_0=2\times 10^{-2} \textrm {m}/D$. Figure 14 compares the LES data for these two cases. Figure 14(a) shows that, as expected, the incoming streamwise turbulence intensity $\langle I\rangle$ is clearly larger for the case with the higher roughness. However, the difference in the turbulence level between the two cases becomes less clear within the wind farm, especially towards the end of the wind farm as shown in figure 14(a–c). In other words, these data suggest that the turbulence added by the wind farm is negatively proportional to the ambient turbulence level. This is consistent with the empirical relation of Crespo et al. (1996) for the added wake turbulence, as highlighted in Zehtabiyan-Rezaie & Abkar (2023). In addition, it is important to note that the higher turbulence level in the AR case promotes the wake recovery after each turbine row which in turn increases the incoming wind speed for the next turbine. This, however, increases the velocity deficit jump that occurs at the next turbine row according to (5.2). This is why despite the difference in the incoming turbulence level, the two cases have fairly similar velocity distribution within the wind farm as shown in figure 14(a). This suggests that the wind farm region is highly dominated by the turbulence generated by the wind farm, and it seems to be less dependent on the incoming turbulence level. However, in the farm wake where the turbulence added by the farm gradually diminishes, the impact of the ambient turbulence becomes more important. Figure 14(a) shows that the wake of the AR case recovers faster than that of the A0 case. Model predictions in comparison with the LES data for these two cases are depicted in figure 15, which overall shows a good agreement for the streamwise velocity deficit. The model predicts the further reduction of the velocity deficit in the wake of the wind farm for the AR case. For the spanwise velocity deficit, results seem to be fairly similar for the two cases. The only notable difference can be observed in the far-wake region of the farm where the spanwise deficit for the AR case is less than that for the A0 case. Similar to the streamwise velocity deficit, the smaller spanwise deficit for the AR case is due to the larger value of ambient turbulence and thereby faster wake recovery.

Figure 14. (a) Variation of laterally averaged streamwise turbulence intensity $\langle I \rangle$ and streamwise velocity $U$ with downwind distance $x$ for both cases of Aligned Baseline (A0) and Aligned Rough (AR). Contours of laterally averaged streamwise turbulence intensity $\langle I \rangle$ on a horizontal plane at the hub height for (b) the AR case and (c) the A0 case.

Figure 15. Variation of laterally averaged (a) streamwise velocity deficit $U_d$ and (b) spanwise velocity deficit $V_d$ at the turbine hub height for the Aligned Rough (AR) case. For comparison, the data for the Aligned Baseline (A0) case are also shown. The dashed lines show the LES data and the solid lines show predictions of the developed model. Vertical dotted lines show the locations of wind turbine rows.

Figures 16 and 17 show, respectively, the effect of wind farm length and turbine spacing on the wake evolution for both the LES and the model. Figure 16 shows that the streamwise velocity deficit is smaller for the short wind farm. This can be simply explained by the fact that there are fewer turbine rows in this wind farm and thereby less thrust forcing acting on the incoming wind. In figure 17, however, the number of turbine rows are the same, and they differ in both streamwise $s_x$ and spanwise $s_y$ inter-turbine spacing as shown in table 1. It is interesting to note that reducing $s_x$ and $s_y$ has multifaceted effects on the laterally averaged streamwise velocity deficit. The reduction of $s_y$ directly increases the velocity deficit, because by reducing $s_y$, turbine wakes occupy a larger potion of the flow field. On the other hand, with a reduction of $s_x$, the local incoming velocity for downwind turbines decreases as shown in figure 18. This reduces the amount of turbine thrust force, which is expected to decrease the total velocity deficit generated by wind turbines. However, the turbine spacing also affects the rate of wake recovery. Decreasing both $s_x$ and $s_y$ increases the turbulent farm velocity scale $\hat {u}_f$, which leads to a larger $\nu _{t,f}$ and thus faster wake recovery. Therefore, knowledge of the relative importance of each of these factors is needed for each case in order to predict the overall impact of changing the turbine spacing on the wake velocity deficit. Figure 17 shows that for this case the impact of $s_y$ change on the velocity deficit is initially dominant as the velocity deficit within the AD farm is much larger than that of the A0 farm. Further downstream, however, wakes of both wind farms experience a fairly similar level of streamwise velocity deficit. The spanwise velocity deficit is also approximately similar in both cases.

Figure 16. Variation of laterally averaged (a) streamwise velocity deficit $U_d$ and (b) spanwise velocity deficit $V_d$ at the turbine hub height for the Aligned Short (AS) case. For comparison, the data for the Aligned Baseline (A0) case are also shown. The dashed lines show the LES data and the solid lines show predictions of the developed model. Vertical dotted lines show the locations of wind turbine rows for the AS case.

Figure 17. Variation of laterally averaged (a) streamwise velocity deficit $U_d$ and (b) spanwise velocity deficit $V_d$ at the turbine hub height for the Aligned Dense (AD) case. For comparison, the data for the Aligned Baseline (A0) case are also shown. The dashed lines show the LES data and the solid lines show predictions of the developed model. Vertical dotted lines show the locations of wind turbine rows for the AD case.

Figure 18. Variation of the local hub-height velocity $\bar {u}_{h,n}$ as a function of turbine row number $n$ for the Aligned Dense (AD) and the Aligned Baseline (A0) cases based on the proposed model.

7. Summary

The aim of this work is to develop a new physics-based one-dimensional model to predict the variation of laterally averaged streamwise and spanwise velocities in the wake of a wind farm at the turbine hub-height level. Through a budget analysis based on the LES data of semi-infinite wind farms, dominant terms in the momentum equations were identified. This led to an approximate form of the momentum equations where the sum of the Coriolis force, the divergence of the Reynolds stresses and the turbine thrust force are in balance with the change in momentum by advection.

The linearised versions of the approximate form of the momentum equations in both the streamwise and spanwise directions were then mathematically solved to obtain the proposed model stated in (5.11). To derive this solution, the turbulent viscosity hypothesis was used to model the Reynolds shear stresses. The turbulent viscosity $\nu _t$ was decomposed into the ambient turbulent viscosity $\nu _{t,0}$ and the farm turbulent viscosity $\nu _{t,f}(x)$, where the latter changes with $x$. The dependency of $\nu _{t,f}$ on $x$ was modelled using a velocity scale proportional to the turbine forcing per unit area, and a length scale proportional to the thickness of the IBL $\delta$. The proposed model importantly accounts for the fact that wind farms with different layouts may generate noticeably different wakes. This is mainly due to the fact that the local incoming velocity experienced by wind turbines and thus the thrust force that is exerted by the wind turbines on the airflow depends on the farm layout. A geometric parameter called the local deficit coefficient $\eta$ was introduced to relate the local velocity deficit at the rotor-centre of wind turbines to the laterally averaged velocity deficit at the same streamwise position. Moreover, the gradients of the incoming wind shear and wind veer appeared in terms $C_x$ and $C_y$ of the final solution (5.11) and were estimated using the geostrophic drag law.

The model predictions are compared with LES data for five different cases to capture the response of the model to various changes in farm operating conditions. The Coriolis and shear terms in the governing equation (5.10) for the streamwise direction are relatively small. Therefore, the change in the streamwise velocity is mainly determined by how different operating conditions influence the turbine forcing term and the wake recovery term in (5.10). In general, a higher local incoming velocity increases the turbine forcing term, which leads to a higher velocity deficit. The wake recovery rate on the other hand is increased with an increase of either the turbulent viscosity or the velocity deficit. Therefore, the overall impact of changing a parameter such as inter-turbine spacing or incoming turbulence is not trivial, and these changes may lead to counteracting effects on the streamwise velocity deficit. Overall, our results showed that the proposed model is able to acceptably predict the variation of streamwise velocity deficit for the different cases studied here. For the spanwise velocity deficit, turbine forces were not present in the LES data (i.e. $\bar {f}_y=0$ as $\gamma =0$ for all turbines). However, in addition to the recovery term, the two Coriolis and veer terms in (5.10) are of great importance. Our results showed that the Coriolis term initially introduces an anticlockwise deflection (i.e. deflection to the left) in the northern hemisphere. Further downstream, the veer term becomes dominant and introduces a clockwise deflection (i.e. deflection to the right) in the northern hemisphere. The former is the direct effect of the Coriolis force, whereas the latter is present due to (i) an indirect effect of the Coriolis force through wind veer and (ii) the farm-generated turbulence. The total value and direction of the wake deflection due to the Coriolis force therefore depends on the streamwise location and more importantly on the relative magnitude of these two counteracting terms with respect to each other.

This work serves as the first study to model the wake of a semi-infinite wind farm. While only aligned and staggered layouts were modelled by our LESs, the developed model is capable of modelling different layouts. Therefore, future research can implement the model to study a wider range of layout configurations. Another interesting area of future research is to study the effect of yaw offset on the farm wake deflection. While this is not studied here, the effect of yaw angle is already incorporated in the model and can be used in future works. It is also of especial interest to study how atmospheric thermal stratification may affect the turbulent viscosity and also our estimation for the shear and veer terms. Finally, it is important to remind ourselves that by definition the model assumes a wind farm that extends infinitely in the lateral direction. Therefore, model predictions only resemble the flow in the wake centre of a wide finite farm where side effects are deemed to be insignificant. More research is thus essential to quantify the impact of lateral flow entrainment and side effects. Moreover, the dimensional analysis used to simplify the vertical gradient of velocities in the recovery term of (5.10) should be scrutinised in more detail. This simplification may imply that the cross-stream length scale associated with the vertical flow shear remains comparable to the rotor diameter. This may not be necessarily true especially in the far wake of a wind farm.