Effect of blockage on wind turbine power and wake development

Recent work by Lanzilao and Meyers (2024) has shown that wind-farm blockage introduces an unfavourable pressure gradient in front of the farm and a favourable pressure gradient in the farm, which are strongly correlated with the nonlocal efficiency and wake efficiency, respectively. In particular, the favourable pressure gradient in the farm increases the farm wake efficiency, defined as the average farm power normalized by the average front-row power. Here, we investigate the impact of blockage on wake development and the power of wind turbines using an idealized large-eddy simulation setup in which blockage conditions are artificially introduced using a rigid lid, in addition to using neutral stratification and no wind veer. We simulate both infinite and finite single turbine rows, as well as a setup with two staggered rows. Blockage strength is adjusted by varying the boundary layer height (H) and turbine spacing (S). We find that blockage strongly affects near-wake behaviour, altering Froude momentum theory, by introducing a favourable pressure difference (Δp_NW) across the turbine row. The same setup also leads to an unfavourable pressure difference (Δp_FW) in the far wake, which simply follows from the rigid-lid conditions and the change in momentum flux due to wake recovery. A strong positive correlation of −Δp_NW with both the power coefficient (C_P) and thrust coefficient (C_T) is observed. Specifically, as S and H decrease, −Δp_NW, C_P, and C_T increase. At the same time, a lower induction is observed at the rotor disc, and a lower wake deficit, in the near wake. The reduction of near-wake velocity deficit as a result of blockage also translates into lower deficits and wake widths in the far wake. When scaling the far-wake development with the initial far-wake deficit and width, we do not see a direct effect of the adverse pressure gradient on the wake recovery. However, we do see a profound effect of H on the wake spreading, with higher boundary layers leading to faster spreading. This relates to the fact that the wake can more freely expand vertically in high-boundary layer cases into a larger region of high-speed flow than for shallow boundary layers. Finally, we introduce a simplified Froude momentum balance to parameterize the relation between blockage, pressure drop, and near-wake properties and compare it to the large-eddy simulation results.

1 Introduction

In large wind farms, significant flow slowdown occurs upstream of the farm, called blockage. Measurements conducted by Bleeg et al. (2018) before and after wind-farm commissioning revealed significant reductions of wind speed upstream of the farm (up to 4 % for 10 rotor diameters). Around the same time, Allaerts and Meyers (2017, 2018) described the excitation of gravity waves in their wind-farm large-eddy simulations (LESs), also observing significant slow down of wind speed in front of the farm. Since then, various other studies have also discussed the excitation of gravity waves by wind farms and related blockage effects (Wu and Porté-Agel, 2017; Maas, 2023; Stipa et al., 2024). Recently, Lanzilao and Meyers (2024) performed LESs of a fixed 1.6 GW wind farm in different atmospheric stratified conditions to investigate effects on wind-farm efficiency, blockage, and related gravity-wave excitation. They found that the upstream slowdown associated with blockage mainly originates from the vertical displacement of the capping inversion at the top the boundary layer and the associated hydrostatic pressure induced by the increased column of cold (higher density) air below the inversion. For the atmosphere conditions considered in their simulations, purely hydrodynamic blockage effects (i.e. associated with Bernoulli's law and induction by the turbines; see also Segalini and Dahlberg, 2020) were at least an order of magnitude smaller (Lanzilao and Meyers, 2022, 2024).

Similar to Allaerts and Meyers (2017, 2018), Lanzilao and Meyers (2024) identified a strong unfavourable pressure gradient upstream of the farm, but they also identified a favourable pressure gradient within the wind farm. They defined the nonlocal efficiency ${\mathit{\eta }}_{\mathrm{nl}}=P_{\mathrm{1}}/P_{\mathrm{\infty }}$ as the ratio between the power of a freestanding turbine P_∞ and the average turbine power of the first row of a wind farm P₁, and they defined the wake efficiency ${\mathit{\eta }}_{\mathrm{w}}=P_{\mathrm{tot}}/\left(N_{\mathrm{t}}{P}_{\mathrm{1}}\right)$ as the ratio between the average turbine power in the wind farm ($P_{\mathrm{tot}}/N_{\mathrm{t}}$) and the average turbine power of the first row. A strong negative correlation was found between the unfavourable upstream pressure gradient and η_nl and a strong positive correlation between the favourable pressure gradient within the farm and η_w. Additionally, for wind farms in the conventionally neutral boundary layer, depending on the atmospheric conditions, the beneficial effects of the favourable pressure gradient can offset the negative effects of the upstream unfavourable pressure gradient, sometimes leading to a larger farm efficiency than in a similar fully neutral case without free-atmosphere stratification.

Figure 1 Sketch (side view) of a single turbine row (infinitely wide) inside an idealized atmospheric boundary layer with the rigid-lid top condition.

In the current paper, our aim is to better understand the relationship between a favourable pressure gradient and improved efficiency using a new set of large-eddy simulations. To this end, we strongly simplify the simulation setup by replacing blockage induced by free-atmosphere stratification with blockage induced by a rigid lid at the top of the boundary layer. Such a rigid-lid condition could be the equivalent of an infinitely strong capping inversion. Although such a condition cannot exist, because the strongest density jump conceivable over the capping inversion should be well below the difference between density at 1 atm and vacuum, rigid-lid blockage nonetheless shows some similarities with blockage induced by free-atmosphere stratification (Smith, 2024). The advantage of using a rigid-lid condition is that the relationship between blockage, pressure gradient, and efficiency is easier to quantify. This is illustrated in Fig. 1 for a single, infinitely wide row of turbines. Looking at a control volume around the near wake and induction region of a turbine, it is directly clear that a favourable pressure difference should exist (Δp_NW<0) – presuming negligible friction at the ground, this is a direct result of Newton's second law, the presence of the turbine thrust force, and continuity (such that outflowing momentum is larger than inflowing momentum). Similarly, an unfavourable pressure difference (Δp_FW>0) should apply in the far wake. Moreover, applying the conservation of momentum in the streamwise direction, considering the entire domain shown in Fig. 1 for an infinite row, that is, from the beginning of the induction zone to the end of the far wake, as wide as the turbine spacing S and as high as the boundary layer height H, we find:

$$\begin{array}{}\text{(1)}& HS(\mathrm{\Delta }{p}_{\text{NW}}+\mathrm{\Delta }{p}_{\text{FW}}+\mathrm{\Delta }{p}_{\mathrm{bg}})+F_{\mathrm{T}}+F_{\text{fric}}=\mathrm{0},\end{array}$$

where F_T (>0) is the turbine thrust, F_fric (>0) is the friction from the ground, and Δp_bg (<0) is the background pressure difference in the absence of the turbines. We remark that the velocity profile at the end of the far wake has recovered back to the velocity profile at the start of the induction zone; therefore, no change in momentum is present in Eq. (1). We assume that F_fric≈Δp_bg; thus, Δp_NW and Δp_FW are dynamic pressure perturbations, superimposed on the background pressure. This simplifies Eq. (1) to $-\mathrm{\Delta }{p}_{\text{NW}}\approx \mathrm{\Delta }{p}_{\text{FW}}+F_{\mathrm{T}}/HS$, and thus $-\mathrm{\Delta }{p}_{\text{NW}}>\mathrm{\Delta }{p}_{\text{FW}}$, from which we hypothesize that the effect of blockage on efficiency is mostly related to changes in the near-waking behaviour.

The effects of domain blockage by rigid boundaries on turbine performance have been studied before, mostly in the context of blockage corrections in wind tunnel experiments of single model wind turbines (Mikkelsen and Sørensen, 2002; Werle, 2010; Segalini and Inghels, 2014). In these studies, the focus was on estimating free-flow conditions, thus eliminating the effects of wind tunnel walls. These corrections were built by extending Froude momentum theory (Froude, 1889) to account for the pressure gradients included by any (small) domain blockage present in the wind tunnel. Here, we will use similar ideas to develop a simple model that correlates favourable pressure gradients induced by blockage with near-wake induction and turbine power extraction, but now to quantify the effect of blockage rather than to exclude it. A similar approach was already used by Garrett and Cummins (2007) for tidal turbines in constant water level channels, and some elements were also used by Nishino and Willden (2012, 2013) for modelling the blockage effect of a finite array of tidal turbines in a wide channel cross-section.

The setup of a row of closely spaced turbines that we are studying in the current work has also been explored in the past by McTavish et al. (2013, 2015) and Strickland and Stevens (2020, 2022). They reported so called “in-field” blockage, leading to increased power production of the turbines, which they attributed to mutual favourable interactions between turbines. However, by performing a careful domain sensitivity study, Bleeg and Montavon (2022) later showed that these beneficial effects disappear when the domain size is sufficiently large. Thus, in these earlier reports, power increase due to “in-field” blockage was effectively a result of domain blockage and consequently may be related to the increased wake efficiency observed by Lanzilao and Meyers (2024) under favourable pressure gradient conditions.

Looking at the sketch in Fig. 1 and referring to the discussion above, an unfavourable pressure gradient is expected in the far wake, albeit smaller than the favourable pressure gradient in the near wake. The effects of pressure gradients on wake development have been studied by Liu et al. (2002) and Shamsoddin and Porté-Agel (2017) for plane wakes and by Shamsoddin and Porté-Agel (2018) for axisymmetric wakes, and later more general conditions were considered by Dar and Porté-Agel (2022). All of these studies conclude that unfavourable pressure gradients lead to slower wake recovery and favourable pressure gradients lead to faster wake recovery. However, we note already that for all simulations considered in the current paper, we found the unfavourable pressure gradients in the far wake to be much smaller than the values appearing in the above studies such that wake recovery is unaffected by them.

In the current work, we set up a range of large-eddy simulations of a neutral rigid-lid pressure-driven boundary layer, in which we represent both infinite and finite single turbine rows, as well as a setup with staggered rows. Blockage strength is adjusted by varying the boundary layer height (H) and turbine spacing (S). Simulation results are further compared to a simple model that extends classical Froude momentum theory, parameterizing the effects of pressure gradients on wind turbine power, thrust, and induction. Such a model can straightforwardly be incorporated into engineering blockage and wake models such as WAYVE (Allaerts and Meyers, 2019; Devesse et al., 2022; Stipa et al., 2024; Devesse et al., 2024a, b), in which large-scale pressure gradients coming from free-atmosphere stratification and gravity wave feedback are available. The latter is, however, a topic of future research and not in the scope of the current paper.

The article is structured as follows. The setup of the LESs is elaborated in Sect. 2. The extension of Froude momentum theory is discussed in Sect. 3. Next, Sect. 4 presents the LES results for the far- and near-wake analysis, including model validation. Lastly, conclusions are provided in Sect. 5.

2 Methodology

2.1 Governing equations

We consider the filtered Navier–Stokes equations for a neutral pressure-driven boundary layer, given by

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_i}}& {\displaystyle }=\mathrm{0},\text{(3)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial t}}+{\stackrel{\mathrm{̃}}{u}}_j{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{p}}^{*}}{\partial x_i}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}{p}_{\mathrm{\infty }}}{\mathrm{d}{x}_{\mathrm{1}}}}{\mathit{\delta }}_{\mathrm{1}i}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}^{\mathrm{sgs}}}{\partial x_j}}+f_i,\end{array}$$

where the horizontal and vertical directions are represented by indices $i=\mathrm{1},\mathrm{2},$ and 3, with $(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})=(x,y,z)$. The filtered velocity components are denoted by ${\stackrel{\mathrm{̃}}{u}}_i$ for the three-dimensional flow field, with $({\stackrel{\mathrm{̃}}{u}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{u}}_{\mathrm{2}},{\stackrel{\mathrm{̃}}{u}}_{\mathrm{3}})=(\stackrel{\mathrm{̃}}{u},\stackrel{\mathrm{̃}}{v},\stackrel{\mathrm{̃}}{w})$. The filtered modified pressure is defined as ${\stackrel{\mathrm{̃}}{p}}^{*}=\stackrel{\mathrm{̃}}{p}-p_{\mathrm{\infty }}+{\mathit{\rho }}_{\mathrm{0}}{\mathit{\tau }}_{kk}/\mathrm{3}$, where p_∞ represents the mean background pressure and ${\mathit{\tau }}_{kk}/\mathrm{3}$ denotes the trace of the subgrid-scale stress tensor ${\mathit{\tau }}_{ij}=\stackrel{\mathrm{̃}}{u_i{u}_j}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j$, and where a subgrid-scale model is used to model the anisotropic component of the residual-stress tensor ${\mathit{\tau }}_{ij}^{\mathrm{sgs}}={\mathit{\tau }}_{ij}-{\mathit{\delta }}_{ij}{\mathit{\tau }}_{kk}/\mathrm{3}$. Note that the direct effect of viscosity on resolved scales in the LES is negligible; all the dissipation is handled by the subgrid-scale stresses. The forces (per unit of density) f_i exerted by the wind turbines on the flow are modelled using a non-rotating actuator disc model (Goit and Meyers, 2015) with a Shapiro correction factor (Shapiro et al., 2019) to avoid over prediction of turbine power on typical LES grid resolutions.

The governing Eqs. (2)–(3) are solved using the SP-Wind solver, an in-house software developed at KU Leuven (Calaf et al., 2010; Goit and Meyers, 2015; Munters et al., 2016b; Allaerts and Meyers, 2017; Lanzilao and Meyers, 2024). The equations are integrated over time using a standard fourth-order Runge–Kutta scheme, with the time step determined by a Courant–Friedrichs–Lewy (CFL) number of 0.4. Discretization in the horizontal direction is performed using a Fourier pseudo-spectral method, employing the $\mathrm{3}/\mathrm{2}$ dealiasing rule. In the vertical direction, an energy-preserving fourth-order finite difference scheme is applied (Verstappen and Veldman, 2003). Continuity is enforced by direct solving the Poisson equation at each stage of the Runge–Kutta method. The influence of subgrid-scale motions on the resolved flow is captured by the Smagorinsky model in combination with Mason and Thomson (1992) wall damping, using a Smagorinsky length ${\mathrm{\ell }}_{\mathrm{s}}^{-n}=(C_{\mathrm{s}}\mathrm{\Delta }{)}^{-n}+[\mathit{\kappa }(z+z_{\mathrm{0}}){]}^{-n}$, with C_s=0.14, n=1, $\mathrm{\Delta }=(\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z{)}^{\mathrm{1}/\mathrm{3}}$ the grid spacing, and κ=0.41 the Von Kármán constant. This is consistent with prior studies using SP-Wind (Meyers, 2011; Allaerts and Meyers, 2017; Lanzilao and Meyers, 2024). We also refer to Calaf et al. (2010), Lignarolo et al. (2016), Martínez-Tossas et al. (2018), and Sood et al. (2022) for code benchmarking and validation.

Figure 2 The plan view of the (a) infinite row, (b) finite row, and (c) infinite staggered rows main computational domains.

Table 1 Overview of the wind-farm simulation cases and their computational domain dimensions $L_x\times L_y\times H$, the turbine spacing $S/D$, the row averaged power coefficient 〈C_P〉_row, and the thrust coefficient 〈C_T〉_row (values of both rows are shown for the infinite staggered row simulations).

2.2 Simulation setup

Figure 2 provides an overview of the different simulations setups considered in the current study, and a list of all cases is provided in Table 1. All turbines correspond to the International Energy Agency (IEA) 15 MW offshore reference wind turbine (Gaertner et al., 2020), with hub height z_h=150 m and diameter D=240 m. Conditions are such that the turbines operate below their rated power. Consequently, a constant disc-based thrust coefficient, $C_{\mathrm{t}}^{\prime }=\mathrm{1.44}$ (Calaf et al., 2010), is used as an input parameter for the simulations.

First of all, a set of single infinite row simulations is considered (Fig. 2a). We select a domain width L_y=9.6 km with spanwise periodic boundary conditions, and by varying the number of turbines in the domain, we can change the turbine spacing S. Secondly, we also consider a finite row setup consisting of seven turbines (Fig. 2b). Here, we take L_y=50 km, so that $L_y/L_y^{\mathrm{f}}\approx \mathrm{7}$ for the case with widest spacing (S=5), following Lanzilao and Meyers (2024), who recommend $L_y/L_y^{\mathrm{f}}\ge \mathrm{6}$ to avoid artificial effects from spanwise boundaries. Finally (Fig. 2c), two infinitely wide staggered rows of turbines are considered, with $S_x/D=\mathrm{5}$, also using a domain width of L_y=9.6 km.

The simulation domain $L_x\times L_y\times H$ is further selected as follows. All simulations use L_x=20 km, with an upstream region L_ind=4 km in front of the turbines and a fringe region length L_fr=2.6 km. Different boundary layer heights are considered, i.e. H=350, 500, and 700 m. All simulations use wall stress boundary conditions at the bottom (with a surface roughness $z_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ m, which is a typical offshore value (Taylor and Yelland, 2001)) and symmetry conditions at the top. Boundary conditions in the horizontal direction are periodic, as a result of our pseudo-spectral discretization method. To break the periodicity in the streamwise direction and prescribe an inflow condition, we use a fringe-region technique to drive the main domain by turbulent, fully developed, statistically steady flow fields obtained from a concurrent precursor simulation (Stevens et al., 2014). Using a precursor simulation ensures that turbulent inflow conditions are generated directly by the Navier–Stokes equations, thereby enhancing the accuracy and realism of the inflow turbulence (Munters et al., 2016a).

To prevent persistent spanwise locking of large-scale streamwise turbulent structures, a shifted periodic boundary condition is applied within the fringe region, bypassing the need for excessive streamwise domain lengths (Munters et al., 2016b).

The precursor simulation domain $L_x^p\times L_y^p\times H$ is defined as follows. Three different precursor simulations are performed for the infinite row cases, i.e. $L_x^p=\mathrm{10}$ km and $L_y^p=\mathrm{9.6}$ with a height of H=350, 500, and 700 m. One additional precursor simulation is performed for the finite row cases, with $L_x^p=L_y^p=\mathrm{10}$ km and H=500 m. However, because the fringe region spans the full width and height of the main domain, SP-Wind requires matching heights and widths between the precursor and main domain when they run concurrently. To achieve this, the tiling technique of Sanchez Gomez et al. (2023) is used to extend the precursor flow fields in the y direction from 10 to 50 km.

All precursor simulations are driven by a constant pressure gradient $\mathrm{d}{p}_{\mathrm{\infty }}/\mathrm{d}{x}_{\mathrm{1}}=-u_{\mathit{\tau }}^{\mathrm{2}}/H$, with a friction velocity u_τ≈0.275 m s⁻¹, typical for offshore conditions (Lanzilao and Meyers, 2024). Fixing u_τ yields approximately the same hub-height wind speed for a given z₀ and varying H simulations. Finally, all simulations and domains use the same grid resolution, i.e. Δx=40 m, Δy=24 m, and Δz=7.93 m in the streamwise, spanwise, and vertical directions, respectively. This is consistent with resolutions used earlier in, e.g., Allaerts and Meyers (2017, 2018); Lanzilao and Meyers (2024).

Precursor simulations are initialized using a log profile $u=(u_{\mathit{\tau }}/\mathit{\kappa })\ln (z/z_{\mathrm{0}})$ in combination with random noise and simulated for 20 000 s, such that a fully developed, statistically stationary, pressure-driven boundary layer can develop. Subsequently, the precursor domain and main domain are run concurrently for another 3000 s for the main domain to fully develop given the precursor inflow. Finally, the precursor and main domain are concurrently progressed for an average time T_av of 10 000 s.

Figure 3 Top view of the general setup of the control volume around a turbine to define blockage effects.

3 Froude momentum theory with blockage

We extend the Froude momentum theory to parameterize the relation between blockage, pressure drop, and near-wake properties for a single turbine present within a row or farm. The model is an extension of classic axial momentum theory developed by Rankine (1865) and Froude (1889). A similar approach was, e.g., used by Werle (2010) and Segalini and Inghels (2014) to derive blockage corrections for wind tunnels or by Garrett and Cummins (2007) in the context of tidal turbines, although we start from a slightly more general formulation that does not define inlet and outlet areas of the control volume in advance (see Fig. 3).

Consider a control volume around a single turbine within a row or farm (Fig. 3) that is based on a streamtube, extending from the start of the induction zone until the end of the near wake. We choose the area A₁=SH, which corresponds to the “available” inflow for the turbine. Note that the outlet area A₂ is a priori unknown but expected to be larger than A₁. In the case of an infinite row of turbines, A₂=A₁ is obtained. Moreover, similarly to Garrett and Cummins (2007) and Werle (2010), A₁ and A₂ do not need to be cylindrical, as the momentum theory is a one-dimensional analysis that does not intrinsically define the cross-sectional shapes. We refer to Fig. A1 in Appendix A for a visualization of the control volume around turbines in an infinite and a finite row, based on the streamlines calculated from the LES. We further define the disc area $A_{\mathrm{d}}=\mathit{\pi }{D}^{\mathrm{2}}/\mathrm{4}$ and the wake area A_w, which is not known a priori.

We presume uniform inflow U_in into the control volume and uniform pressures p_in and p_in+Δp_NW over the inlet and outlet, respectively. Furthermore, the disc velocity $U_{\mathrm{d}}=U_{\text{in}}(\mathrm{1}-a)$, with a being the axial induction factor. This leads to a thrust force $F_{\mathrm{T}}=\mathrm{1}/\mathrm{2}\mathit{\rho }{C}_{\mathrm{T}}^{\prime }{U}_{\text{in}}^{\mathrm{2}}(\mathrm{1}-a{)}^{\mathrm{2}}{A}_{\mathrm{d}}$ and, in the absence of drag or mechanical losses, also to $P=\mathrm{1}/\mathrm{2}\mathit{\rho }{C}_{\mathrm{T}}^{\prime }{U}_{\text{in}}^{\mathrm{3}}(\mathrm{1}-a{)}^{\mathrm{3}}{A}_{\mathrm{d}}$.

First, applying the principle of conservation of mass on the streamtube that passes through the rotor and the entire streamtube, we obtain

$$\begin{array}{}\text{(4)}& {\displaystyle }{A}_{\mathrm{d}}{U}_{\mathrm{in}}(\mathrm{1}-a)& {\displaystyle }=A_{\mathrm{w}}{U}_{\mathrm{w}},\text{(5)}& {\displaystyle }{A}_{\mathrm{w}}{U}_{\mathrm{w}}+(A_{\mathrm{2}}-A_{\mathrm{w}})U_{\mathrm{s}}& {\displaystyle }=A_{\mathrm{1}}{U}_{\mathrm{in}}.\end{array}$$

Secondly, we apply the principle of conservation of momentum in the streamwise direction over the entire streamtube. Similar to Eq. (1), we presume that the friction at the ground is balanced by the background pressure difference present in the absence of a turbine (F_fric≈Δp_bg); hence, to the first order, Δp_NW contains the pressure perturbations due to the presence of the wind turbines only. A similar strategy was used by, e.g., Kirby et al. (2023). As a result, the momentum balance equation is

$$\begin{array}{}\text{(6)}& \begin{array}{rl}-& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{T}}^{\prime }{U}_{\text{in}}^{\mathrm{2}}(\mathrm{1}-a{)}^{\mathrm{2}}{A}_{\mathrm{d}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathit{\alpha }\mathrm{\Delta }{p}_{\mathrm{NW}}(A_{\mathrm{1}}+A_{\mathrm{2}})\\ & =(A_{\mathrm{2}}-A_{\mathrm{w}})U_{\mathrm{s}}^{\mathrm{2}}+A_{\mathrm{w}}{U}_{\mathrm{w}}^{\mathrm{2}}-A_{\mathrm{1}}{U}_{\text{in}}^{\mathrm{2}},\end{array}\end{array}$$

where the density ρ is a known constant value. Here, we have further presumed that the average pressure on the mantle of the streamtube corresponds approximately to $p_{\mathrm{side}}\approx p_{\text{in}}+\mathit{\alpha }\mathrm{\Delta }{p}_{\mathrm{NW}}$. Thirdly, using $F_{\mathrm{T}}=(p^{+}-p^{-})A_{\mathrm{d}}$, and eliminating p⁺ and p⁻ using Bernoulli for streamlines in front of and behind the disc, respectively (similar to what is done in derivations of the classical Betz theory), we arrive at

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{T}}^{\prime }{U}_{\text{in}}^{\mathrm{2}}(\mathrm{1}-a{)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_{\mathrm{in}}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_{\mathrm{w}}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\Delta }{p}_{\mathrm{NW}}.\end{array}$$

Finally, using Bernoulli for streamlines that do not pass through the rotor area, we also find

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_{\text{in}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\Delta }{p}_{\mathrm{NW}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_{\mathrm{s}}^{\mathrm{2}}.\end{array}$$

The above leads to a set of five model equations (Eqs. 4–8) with, in principle, four known input variables $A_{\mathrm{1}},A_{\mathrm{d}},C_{\mathrm{t}}^{\prime }$, and U_in, which leaves six unknown variables: $a,U_{\mathrm{w}},U_{\mathrm{s}},A_{\mathrm{w}},A_{\mathrm{2}}$, and Δp_NW. Thus, we lack one equation to arrive at a closed system.

There are two situations that lead to a simple closure of the system. First of all, imposing Δp_NW=0 leads to the classical Betz–Joukowsky theory for a single isolated turbine. The other case corresponds to an infinite row of turbines, in which case $A_{\mathrm{2}}=A_{\mathrm{1}}=SH$ (see Appendix A, Fig. A1a). This leads to the wind tunnel blockage corrections discussed earlier (e.g. Werle, 2010; Segalini and Inghels, 2014; note that a converging or diverging wind tunnel would require a known A₂≠A₁).

In the more general setting of a finite row of turbines, neither A₂ nor Δp_NW is known a priori. At the sides of a finite row, turbines clearly have more space to expand sideways, and possibly a lower Δp_NW applies than in the centre of the row, where the expansion of the inflow area A₁ is hindered by the surrounding turbines (see Appendix A, Fig. A1a). Thus, in such a system, a coupling through a larger pressure system can be expected, with a pressure gradient not only in the streamwise direction but also in the spanwise direction. An additional relation for the pressure system may be, e.g., obtained from an atmospheric perturbation model (Allaerts and Meyers, 2019; Devesse et al., 2022; Stipa et al., 2024; Devesse et al., 2024a, b; see also the open-source model WAYVE), which incidentally also parameterizes the more complex relation between capping inversion displacement and gravity wave excitation in wind farms. It is, however, not in the scope of the current work to develop such a coupling. Instead, we will evaluate the pressure system arising in our large-eddy simulations in more detail and use Δp_NW measured from the simulations as an additional input to evaluate the system above, comparing its output in terms of power and thrust with those obtained directly from LES.

4 Results and discussion

4.1 Velocity deficit and wake recovery

4.1.1 Induction zone and near-wake analysis

Here, we look in more detail into the wake development for the different simulation cases. In Fig. 4, we present the streamwise centreline velocity (U_c), scaled by the free-stream velocity at hub height (U_∞) and averaged over the turbines in each row, for all cases, providing a general overview of the differences and possible similarities in wake development. In Fig. 4a, two regions of interest are depicted, i.e. the induction zone and near wake on the one hand and the far wake on the other hand. We use the minimum of the velocity to mark the end of the near-wake region, which can depend significantly on the boundary layer height H as, e.g., seen in Fig. 4a. The start of the far-wake region is selected as $x/D=\mathrm{6}$ for all cases. As we will show below (see Fig. 6a and b), the far-wake velocity profile can be fitted very well with a Gaussian function, which we can use to quantify the wake recovery behaviour in detail. Upstream of our selected starting point for the far wake, such a fit does not work very well, as the wake profile is transitioning from a top-hat to a Gaussian shape, potentially requiring more advanced fitting shapes.

Figure 4 Streamwise centreline velocity U_c, scaled with the free-stream velocity U_∞, averaged over all turbines within the row for (a) single infinite row cases with H=500 m, (b) single infinite row cases with S=5D, (c) finite row cases, and (d) staggered infinite row cases.

First, looking at the near wake in Fig. 4, we observe that the maximum velocity deficit in the near wake decreases with S and H, i.e. with increasing blockage, for all cases. Consequently, we also expect the axial induction factor (a) to depend significantly on blockage. In Fig. 5, we show the axial induction factor for all cases, as obtained from the LES, as a function of the inverse geometrical blockage ratio $(A/A_{\mathrm{d}}=\mathrm{4}SH/\mathit{\pi }{D}^{\mathrm{2}})$, i.e. cases with strong blockage have a low inverse geometrical blockage ratio. We note that an alternative geometric blockage parameter based on the cross-section of the domain, i.e. $A_{\mathrm{domain}}/A_{\mathrm{farm}}=\mathrm{4}{L}_{y}H/N\mathit{\pi }{D}^{\mathrm{2}}$, with L_y being the domain width, as, e.g., used for analysis of the tidal channels (see, e.g., Nishino and Willden, 2012), is not meaningful here, as, for our finite wind-farm simulations, we are considering L_y→∞. We evaluate $a=\mathrm{1}-U_{\mathrm{d}}/U_{\text{in}}$, with U_d being the disc averaged (and time averaged) velocity and U_in the turbine inflow velocity at hub height. For all infinite row cases, we also average over the different turbines in the row (as they have, in principle, the same induction factor); for the finite row cases, we show the induction of the individual turbines. For all cases, except the staggered row cases, U_in is simply the far upstream inflow velocity at hub height U_∞, obtained from the concurrent precursor simulation. It is calculated as the streamwise velocity averaged along the streamwise line passing through the turbine hub location. For the staggered cases, particularly the second row, we observe, however, a significant acceleration of the flow upstream of the second-row turbines, which passes in between the turbines of the first row (see Fig. 4d). Therefore, for the second row, U_in is defined as the maximum velocity upstream of the turbine, which is located at the end of the near wake of the turbines in the first row (see Fig. 4d). Looking in particular at the comparison between the axial induction factor in the first and second row of the staggered row case and the single infinite row case, for the same boundary layer height (and thus the same inverse geometrical blockage ratio $A/A_{\mathrm{d}}=\mathrm{4}SH/\mathit{\pi }{D}^{\mathrm{2}}$), we observe that they are nearly equal and within averaging uncertainty, indicating that our choice of U_in works well. Later (see Sect. 4), we will show that this definition also leads to good performance of the simple model (Eqs. 4–8) when compared to LES data.

The error bars in Fig. 5 are constructed using moving block bootstrapping (MBB). Time averaging was performed over a time interval T_av=10 000 s, sampling every 2 s. The MBB method splits the original time series with n data samples into $N_{\mathrm{b}}=n-L+\mathrm{1}$ overlapping blocks containing L samples. From this pool of N_b blocks, a new time series is assembled by randomly choosing $K=n/L$ blocks with replacement, and then the mean of the new time series is calculated. This process is repeated B times, resulting in a distribution of means. Finally, the 2.5th and 97.5th percentiles of the distribution of means mark the 95 % confidence interval. The MBB approach is defined by selecting the number of bootstrapping runs B and the size of the blocks L. The preliminary sensitivity study showed convergence, for L=20, to a robust value that does not significantly change for longer block lengths. A similar sensitivity study showed that B=1500 iterations is sufficient for this purpose.

Looking further at Fig. 5a, we observe that the axial induction factor converges to a≈0.264 for a high inverse geometrical blockage ratio ($A/A_{\mathrm{d}}$ high). This is inline with the expected Betz–Joukowsky value of $a=C_{\mathrm{T}}^{\prime }/(\mathrm{4}+C_{\mathrm{T}}^{\prime })=\mathrm{0.265}$. The plot also shows that the finite row cases have an overall higher induction than their infinite row counterparts for the same inverse geometrical blockage ratio, implying that less blockage is present in the finite row cases. This will be explained further in Sect. 4.2.2. Overall, at low inverse geometrical blockage ratios ($A/A_{\mathrm{d}}$ low), we see a significant deviation of the induction factor towards values lower than the expected Betz–Joukowsky value. This difference further translates towards the far wake, as seen before in Fig. 4. It is further interesting to look at the induced near-wake velocity. In Fig. 5b, we show $(\mathrm{1}-U_{\mathrm{w}}/U_{\text{in}})/\left(\mathrm{2}a\right)$, where U_w is evaluated as the minimum wake velocity (see Fig. 4) at hub height. Note that the expected Betz–Joukowsky value of this ratio corresponds to 1. As can be seen in Fig. 5b, the induced near-wake velocity approaches the Betz–Joukowsky value, for high inverse geometrical blockage, but some differences remain. In particular, the onset of wake recovery can reduce the maximum wake deficit expected from the pure inviscid solution predicted by Betz–Joukowsky. Next to that, subtle effects related to the presence of shear may also play a role. When looking at decreasing inverse geometric blockage ratios, we observe that not only a lower induction occurs at the turbine disc but also an even lower wake deficit in the near wake (i.e. $(\mathrm{1}-U_{\mathrm{w}}/U_{\mathrm{in}})<\mathrm{2}a$).

Figure 5 LES (a) axial induction factor a and (b) induced near-wake velocity scaled with the axial induction obtained, as a function of the inverse geometrical blockage ratio $A/A_{\mathrm{d}}=\mathrm{4}SH/\left(\mathit{\pi }{D}^{\mathrm{2}}\right)$. The row averages are shown for the infinite row cases, while individual turbine values are shown for the finite row cases. Error bars represent the 95 % confidence intervals, obtained using moving block bootstrapping.

4.1.2 Far-wake analysis

We first discuss the evolution of the far wake as seen in a horizontal plane at hub height. Here, we observe that a classical Gaussian shape function provides good fits along the downstream direction. This is in line with Liu et al. (2002) and Shamsoddin and Porté-Agel (2017, 2018), who observed that the turbulent wake profiles retained a Gaussian shape under non-zero pressure gradient conditions. Thus, we use a Gaussian profile (Pope, 2000):

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{U_{\mathrm{s}}\left(x\right)-u(x,y)}{U_{\mathrm{s}}\left(x\right)}}=C\left(x\right)e^{-y^{\mathrm{2}}/\left(\mathrm{2}{\mathit{\delta }}_y\right(x{)}^{\mathrm{2}})},\end{array}$$

where U_s(x) is the streamwise velocity far from the wake centre, u(x,y) is the wake velocity, C(x) is the normalized wake velocity deficit, and δ_y(x) is the horizontal wake width. We fit this function to the LES data, thus extracting U_s(x), C(x), and δ_y(x) as a result at each downstream position. Figure 6a and b show indeed that the horizontal wake profile, averaged across turbines in the row, has a Gaussian shape in the specified downstream range. Similar profiles are also found for all the other simulation cases.

Figure 6 Horizontal wake velocity profiles, averaged over the turbines within the row, at different positions downstream in the far wake. The circles show LES results for the cases (a) Inf-H500-S2.5 and (b) Inf-H500-S5, and the orange lines represent the classical Gaussian shape function from Eq. (9). The R² value denotes the coefficient of determination.

Figure 7 Evolution of the far-wake properties for all simulations, averaged over the turbines within the row. (a–c) Normalized wake velocity deficit C(x) and (d–f) wake width δ(x), scaled with their respective value at the start of the far wake C₀=C(x₀) and ${\mathit{\delta }}_{y,\mathrm{0}}={\mathit{\delta }}_{y,\mathrm{0}}\left(x_{\mathrm{0}}\right)$, for various H, plotted against the downstream distance $(x-x_{\mathrm{0}})/{\mathit{\delta }}_{y,\mathrm{0}}$, with x₀=6D being the starting location of the far wake.

Figure 7a–c show the evolution of C(x), and Fig. 7d–f show the evolution of δ_y(x), averaged over the different turbines in each row, scaled with their values at the start of the far wake for all simulations. The downstream distance is reformulated as $(x-x_{\mathrm{0}})/{\mathit{\delta }}_{y,\mathrm{0}}$, with x₀ being the starting location of the far wake. Figure 7b and e show that the far-wake development of the two infinite staggered rows, as well as the single finite row, follows the same trend as that of the single infinite row. As H increases, C(x) shows a slightly steeper decline, while δ_y(x) spreads faster. In particular, the difference between H=350 m and the other two boundary layer heights is clearly visible. This can be understood by looking at the evolution of the vertical velocity profiles (in the wake centre) shown in Fig. 8 for infinite row cases with S=5D. In particular for cases where H=350 m, it is observed that vertical spreading of the wake is limited by the presence of the rigid lid, while the maximum wake deficit seems minimally impacted. Note that the turbine tip height (270 m) is not much lower than the boundary layer height in this case. Thus, for this case, the turbine far wake behaves less as an axisymmetric and more as a planar wake, which is known to have a lower spreading rate (Pope, 2000). We remark that the rigid lid reduces the overall wake spreading only at low H, while the wake recovery seems only slightly affected. We suggest that this imbalance is due to the unusual shape of the wake in the vertical direction for H=350 m.

Returning to Fig. 7, we further observe that for constant H, the C(x) and δ_y(x) curves align closely across varying S values and thus varying adverse pressure strengths. This contrasts with the findings of Liu et al. (2002) and Shamsoddin and Porté-Agel (2018), which suggest that a stronger adverse pressure gradient will slow down the wake deficit recovery and enhance wake spreading. However, the adverse pressure gradient in the far wake of our simulations (normalized by $\mathit{\rho }{U}_{\mathrm{in}}^{\mathrm{2}}$), turns out to be at least an order of magnitude smaller than the one in Liu et al. (2002) and Shamsoddin and Porté-Agel (2018), which they obtained as a result of diverging domain boundaries. Finally, looking back at Fig. 6a, which shows the turbine wake for the infinite row case Inf-H500-S2.5, which is periodic in the y direction, we see no uniform flow at the sides of the wake, starting from x=15D. This implies that for cases with S=2.5D (and also S=5D for the two staggered rows), neighbouring wakes start touching around this downstream location. This also drastically changes the wake spreading rate, as wakes essentially start to merge. We thus conclude that blockage has negligible direct effect on the far-wake recovery. Nevertheless, Fig. 4 shows a positive correlation between blockage and the near-wake deficit. Therefore, we suggest that blockage enhances wake efficiency, with the impact occurring primarily in the near wake.

Figure 8 Vertical wake velocity profiles at different downstream positions, averaged over the turbines within the row for the infinite row turbine cases with S=5D and H=350, 500, and 700 m.

4.2 Evaluation of the near-wake blockage model

We now focus on the effect of blockage on the power and thrust and evaluate the simple near-wake blockage model developed based on momentum theory (Eqs. 4–8) against the LES data. To this end, we evaluate from the large-eddy simulations:

$$\begin{array}{}\text{(10)}& {\displaystyle }{C}_{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{F{U}_{\mathrm{d}}}{\frac{\mathrm{1}}{\mathrm{2}}{A}_{\mathrm{d}}{U}_{\mathrm{in}}^{\mathrm{3}}}}=C_{\mathrm{T}}^{\prime }{\left({\displaystyle \frac{U_{\mathrm{d}}}{U_{\mathrm{in}}}}\right)}^{\mathrm{3}}=C_{\mathrm{T}}^{\prime }(\mathrm{1}-a{)}^{\mathrm{3}},\text{(11)}& {\displaystyle }{C}_{\mathrm{T}}& {\displaystyle }={\displaystyle \frac{F}{\frac{\mathrm{1}}{\mathrm{2}}{A}_{\mathrm{d}}{U}_{\mathrm{in}}^{\mathrm{2}}}}=C_{\mathrm{t}}^{\prime }{\left({\displaystyle \frac{U_{\mathrm{d}}}{U_{\mathrm{in}}}}\right)}^{\mathrm{2}}=C_{\mathrm{t}}^{\prime }(\mathrm{1}-a{)}^{\mathrm{2}},\end{array}$$

where a is evaluated as discussed above (see Sect. 4.1.1 and Fig. 5).

4.2.1 Infinite row

For the infinite row cases, all input variables ($A_{\mathrm{1}}=A_{\mathrm{2}}=SH,A_{\mathrm{d}},C_{\mathrm{t}}^{\prime }$, and U_in) required for the blockage model are known. U_in is obtained from the concurrent precursor as the streamwise velocity averaged over the streamwise line through the turbine hub location, except in the case of the second staggered row, where U_in is derived from the wind-farm simulation as the maximum velocity upstream of the turbine. We also briefly compare with calculating C_P and C_T in the LES using a rotor-disc-averaged velocity at rotor height in the precursor domain. In the results (not shown here), this produced similar C_P and C_T, with less than a 1 % difference for the highest blockage case. With $C_{\mathrm{T}}^{\prime }$ fixed, C_P and C_T depend on the ratio $U_{\mathrm{d}}/U_{\text{in}}$ (i.e. the induction factor), which depends solely on the inverse geometrical blockage ratio. Thus, for varying U_in but constant SH, the model returns a single C_P and C_T value.

Figure 9a and c show 〈C_P〉_row and 〈C_T〉_row as a function of $S/D$, scaled with C_P,0 and C_T,0, which are the C_P and C_T values of a single freestanding turbine without blockage. To this end, we used the turbine in case Inf-H700-S40, as the geometric blockage ratio $A_{\mathrm{d}}/A=\mathrm{0.0067}$, which can be regarded as negligible (Segalini and Inghels, 2014). As can be observed in the figure, the model agrees well with the LES data, although it tends to overpredict both thrust and power, particularly at high blockage ratios. We further see that at a turbine spacing of $S/D=\mathrm{5}$, there is a notable increase in 〈C_P〉_row of 11 %, 8 %, and 5 % for H=350, 500, and 700 m, respectively, compared to the single freestanding turbine.

Finally, Fig. 9b and d show $\langle C_{\mathrm{P}}{\rangle }_{\mathrm{row}}/C_{\mathrm{P},\mathrm{0}}$ and $\langle C_{\mathrm{T}}{\rangle }_{\mathrm{row}}/C_{\mathrm{T},\mathrm{0}}$ results from Fig. 9a and c as a function of the inverse geometrical blockage ratio $(A/A_{\mathrm{d}})$. As mentioned above, the model output in this case depends only on the blockage ratio. We observe that the different LES results also collapse well onto a single curve when this scaling is used.

Figure 9 (a) Row averaged scaled power coefficient $\langle C_{\mathrm{P}}{\rangle }_{\mathrm{row}}/C_{\mathrm{P},\mathrm{0}}$ as a function of spacing $S/D$ and (b) as a function of inverse geometric blockage ratio $A/A_{\mathrm{d}}$ for varying boundary layer heights H. (c) Row averaged scaled thrust coefficient $\langle C_{\mathrm{T}}{\rangle }_{\mathrm{row}}/C_{\mathrm{T},\mathrm{0}}$ as a function of $S/D$ and (d) as a function of $A/A_{\mathrm{d}}$ for varying H. Data shown for all infinite row cases (single and staggered). C_P,0 and C_T,0 are the C_P and C_T of a single freestanding turbine with no blockage, i.e. case Inf-H700-S40. The error bars are bootstrap 95 % confidence intervals.

4.2.2 Finite row

For the finite row cases, the known input variables are A₁=SH, A_d, $C_{\mathrm{t}}^{\prime }$, and U_in, which is obtained from the concurrent precursor as the streamwise velocity averaged over the streamwise line through the turbine hub location. For this case, the area A₂ is not known a priori, nor is the favourable pressure gradient known. This would require an additional closure relation that related the wind-farm thrust to the larger pressure system around the farm.

Figure 10 (a) Pressure drop at the end of the near wake (−Δp_nw(y)) and (b) pressure drop plus the lateral Reynolds stress ($\stackrel{\mathrm{‾}}{v^{\prime }{v}^{\prime }}\left(y\right))$, averaged over H, as a function of the spanwise direction, where y=0 km is the location of the centre turbine in the row.

To appraise this large-scale pressure system, we evaluate for the finite row cases in Fig. 10a the pressure drop −Δp_nw(y), containing only the pressure perturbations superimposed on the background pressure gradient (see Eq. 3). −Δp_nw(y) is calculated as the difference between the domain inlet and the end of the near wake (defined as the location of maximum wake deficit), averaged over H, and shown as a function of the spanwise direction y, where y=0 is located at the row centre. It is observed that the pressure drop reaches a maximum at the row centre, decreasing for larger absolute values of y until a (near-zero) pressure difference is reached far from the turbine row. As expected, the wakes of turbines at the edge of the row have more space for lateral expansion and thus experience a lower $\left|\mathrm{\Delta }{p}_{\mathrm{NW}}\right|$ than those at the centre, where expansion is constrained by neighbouring turbines and the boundary layer height. This explains the overall lower axial induction found in the finite row cases compared to the equivalent infinite row cases, as seen in Fig. 5. All turbines in a finite row have fewer neighbouring turbines, resulting in a less constrained flow around the turbines. This can be considered as the collective blockage effect. Turbines within an infinite row are thus an extreme example of the centre turbines in a finite row. Additional simulations are needed to quantify the transition of a finite row into an infinite row, but this is beyond the scope of the current study. Looking in more detail at the pressure drop around the turbines, we see small local variations. These are partially related to the normal Reynolds stresses. For a conventional axisymmetric wake, it is well understood that $(\mathrm{1}/\mathit{\rho })\partial p/\partial r+\partial \stackrel{\mathrm{‾}}{u_r^{\prime }{u}_r^{\prime }}/\partial r=\mathrm{0}$, with $\stackrel{\mathrm{‾}}{u_r^{\prime }{u}_r^{\prime }}$ being the radial normal stress (Pope, 2000). Thus, in Fig. 10b, we have plotted $\mathrm{\Delta }[p/\mathit{\rho }+\stackrel{\mathrm{‾}}{v^{\prime }{v}^{\prime }}]$ as a function of y (and averaged over the boundary layer height). Note that the spanwise fluctuation v^′ is, strictly speaking, only in the radial direction in a horizontal plane, but as appreciated in the figure, the local variations in the pressure drop (in particular in the wake centres) disappear.

To close the system of Eqs. (4)–(8), we will use Δp_NW as observed in the LES and compare the thrust and power output of the near-wake blockage model with that of the LES. To obtain Δp_NW, we simply average Δp_nw(y) over the interval $[y_h-S/\mathrm{2},y_h+S/\mathrm{2}]$ for each turbine in the row. We note that a pressure closure at the farm level could, e.g., be provided by typical gravity wave models such as, e.g., WAYVE (Allaerts and Meyers, 2019; Devesse et al., 2022; Stipa et al., 2024; Devesse et al., 2024a, b), in which a coupling between the larger atmospheric flow around the farm and local wake models is foreseen. However, here, we aim specifically at assessing the relevance of a simple near-wake model that accounts for a favourable pressure drop. A full interaction with the effects of gravity-wave feedback is a topic for future research.

In addition to closing the system with the pressure gradient from LES, we also evaluate the two-scale model proposed by Nishino and Willden (2012, 2013) for half-open channel flows (formulated in the context of tidal turbines). At the turbine scale (i.e. the near-wake model), they use relations similar to Eqs. (4)–(8), while at the farm scale, they use classical Froude momentum theory without a pressure gradient to obtain a relation for the overall farm wake expansion. By further assuming that all turbines in a row produce the same power, they arrive at a closed system of equations. Although, the drawback of the approach is that there is no spanwise variation in turbine power output, it does allow for an overall assessment of farm power. We refer to Nishino and Willden (2012, 2013) for more details about this model and its implementation.

Figure 11a and b present the C_P and C_T for each turbine of cases Fin-H500-S2.5 and Fin-H500-S5, together with the predicted values from the extended Froude momentum theory model and the model developed by Nishino and Willden (2012). We notice that 〈C_P〉_row is higher for the case with smaller S, i.e. stronger blockage, which was also observed for the infinite row cases. Additionally, the C_P value increases towards the row centre, with $C_{\mathrm{P},\mathrm{centre}}>\langle C_{\mathrm{P}}{\rangle }_{\mathrm{row}}>C_{\mathrm{P},\mathrm{side}}$, which was also found in the Reynolds-averaged Navier–Stokes results of Nishino and Draper (2015). The same findings are observed for the C_T data in Fig. 11b. We appreciate from the figures that, for the high blockage case (S2.5), both models overpredict the power. For the lower blockage case (S5), Nishino and Willden (2012, 2013)'s model is reasonably close to the row average predicted by the LES, whereas our model manages to predict also the shape of the power and thrust along the row rather well.

Figure 11 (a) C_P and (b) C_T for each turbine in the finite row cases, as found by the large-eddy simulations, the extended Froude momentum theory model, and the model from Nishino and Willden (2012). The error bars are bootstrap 95 % confidence intervals. Horizontal dashed lines mark the row averaged values of the corresponding LES cases.

Finally, in Fig. 12, the power and thrust coefficients C_P and C_T are provided for all simulations, assembling the results form Figs. 9 and 11 into one figure. Overall, the model agrees well with the LES data, although it tends to overpredict both thrust and power, in particular at high blockage ratios, and for both the infinite and finite row cases, the error increases. We observe smaller increases in C_P and C_T for the finite row cases, compared to their infinite row counterparts, as less surrounding turbines are present to constrain the inflow from expanding, thus resulting in weaker pressure drops and blockage.

Figure 12 (a) The scaled power coefficient $C_{\mathrm{P},\mathrm{LES}}/C_{\mathrm{P},\mathrm{0}}$ and (b) the scaled thrust coefficient $C_{\mathrm{T},\mathrm{LES}}/C_{\mathrm{T},\mathrm{0}}$, as found by the extended Froude momentum theory model (x axis) and the large-eddy simulations (y axis). The squares, circles, and triangles represent the row averaged values of the single infinite row cases with respect to H. The stars represent the values of the individual turbines for the finite row cases. The diamonds represent the row averaged values of the infinite staggered row cases. The R² value denotes the coefficient of determination.

4.3 Discussion

To locate the origin of the model error, we perform a streamwise momentum budget analysis. The budget equation is derived by time-averaging the streamwise momentum equation and then integrating it over the control volume Ω, which is the same as the model's control volume. In the streamwise direction, the control volume extends from the computational domain entrance to the end of the near wake (defined as the location with maximum wake deficit), x₁=0 to x₂=x_NW. The vertical dimension of the control volume coincides with the vertical boundaries of the computational domain, that is, from z₁=0 to z₂=L_z. Therefore, there is no momentum flux across those faces. The lateral extension of the control volume follows the streamlines originating from the domain entrance at locations $y=y_{\mathrm{hub}}\pm S/\mathrm{2}$. Therefore, $(\overline{u}{n}_x+\overline{v}{n}_y)=\mathrm{0}$ and thus no advection of momentum across those faces (see Appendix A, Fig. A1, for a visualization of the streamlines). We denote $x_{\mathrm{1}},x_{\mathrm{2}},s_{\mathrm{1}},s_{\mathrm{2}},z_{\mathrm{1}},$ and z₂ as the boundaries of the control volume, where s represents the streamline coordinate. The y-z and x-y boundary faces are denoted by A_x and Ψ_z, respectively, while the s-z streamline boundary faces are denoted by Γ_s. By applying the divergence theorem to the streamwise momentum equation, the divergence operator is removed, allowing the conversion of the volume integral into a surface integral. The resulting equation is

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Delta }{\mathcal{M}}_{A_x}}{\underbrace{-{\left[\underset{A_x}{\int }\overline{u}\overline{u}\mathrm{d}{A}_x\right]}_{x_{\mathrm{1}}}^{x_{\mathrm{2}}}}}\underset{\mathrm{\Delta }{\mathcal{F}}_{A_x}}{\underbrace{-{\left[\underset{A_x}{\int }\left({\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}}^r+{\mathit{\tau }}_{xx}^{\text{sgs}}\right)\mathrm{d}{A}_x\right]}_{x_{\mathrm{1}}}^{x_{\mathrm{2}}}}}\\ {\displaystyle }& {\displaystyle }+\underset{\mathrm{\Delta }{\mathcal{F}}_{{\mathrm{\Gamma }}_s}}{\underbrace{\sum _{j=\mathrm{1}}^{\mathrm{2}}{\left[\underset{{\mathrm{\Gamma }}_s}{\int }\left(\left({\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}}^r+{\mathit{\tau }}_{xx}^{\text{sgs}}\right)n_x+\left({\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}}^r+{\mathit{\tau }}_{xy}^{\text{sgs}}\right)n_y\right)\mathrm{d}{\mathrm{\Gamma }}_s\right]}_{/}^{s_j}}}\\ {\displaystyle }& {\displaystyle }+\underset{\mathrm{\Delta }{\mathcal{F}}_{\mathrm{fric}}^{*}}{\underbrace{{\left[\underset{{\mathrm{\Psi }}_z}{\int }{\mathit{\tau }}_{xz}{n}_z\mathrm{d}{\mathrm{\Psi }}_z\right]}_{/}^{z_{\mathrm{1}}}-\underset{\mathrm{\Omega }}{\int }{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_o}}{\displaystyle \frac{\partial {\overline{p}}_{\mathrm{\infty }}}{\partial x}}\mathrm{d}\mathrm{\Omega }}}\\ {\displaystyle }& {\displaystyle }\underset{{\mathcal{P}}_x^{*}}{\underbrace{-\underset{\mathrm{\Omega }}{\int }{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_o}}{\displaystyle \frac{\partial {\overline{p}}^{*}}{\partial x}}\mathrm{d}\mathrm{\Omega }}}+\underset{{\mathcal{F}}_{\mathrm{t},x}}{\underbrace{\underset{\mathrm{\Omega }}{\int }{\overline{f}}_x\mathrm{d}\mathrm{\Omega }}}=\mathrm{0}\end{array},\end{array}$$

where $u^{\prime }=u-\overline{u}$ and

$$\begin{array}{}\text{(13)}& {\displaystyle }\begin{array}{rl}{\displaystyle }\underset{{\mathcal{P}}_x^{*}}{\underbrace{\underset{\mathrm{\Omega }}{\int }{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_o}}{\displaystyle \frac{\partial {\overline{p}}^{*}}{\partial x}}\mathrm{d}\mathrm{\Omega }}}& {\displaystyle }=\underset{{\mathcal{P}}_{A_x}^{*}}{\underbrace{{\left[\underset{A_x}{\int }{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_o}}{\overline{p}}^{*}\mathrm{d}{A}_x\right]}_{x_{\mathrm{1}}}^{x_{\mathrm{2}}}}}\\ {\displaystyle }& {\displaystyle }+\underset{{\mathcal{P}}_{{\mathrm{\Gamma }}_s}^{*}}{\underbrace{\sum _{j=\mathrm{1}}^{\mathrm{2}}{\left[\underset{{\mathrm{\Gamma }}_s}{\int }{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_o}}{\overline{p}}^{*}{n}_x\mathrm{d}{\mathrm{\Gamma }}_s\right]}^{s_j}}}.\end{array}\end{array}$$

The term $\mathrm{\Delta }{\mathcal{M}}_{A_x}$ describes the advection of streamwise momentum by the mean streamwise velocity, while $\mathrm{\Delta }{F}_{A_x}$ accounts for the total (both resolved and modelled) turbulent transport of streamwise momentum along the streamwise direction. Further, $\mathrm{\Delta }{\mathcal{F}}_{\mathrm{fric}}^{*}$ captures the additional ground friction, not balanced by the background pressure gradient, due to the side flow acceleration induced by the turbines. Lastly, ${\mathcal{P}}_x^{*}$ accounts for pressure gradients caused by wind-farm effects, and ℱ_t,x stands for the turbine thrust force. By the sign convention used, the first six terms in Eq. (12) are positive when the flow or turbulence brings more momentum into the control volume than out and negative when the opposite occurs. For instance, $\mathrm{\Delta }{\mathcal{M}}_{A_x}$ is positive if the mean streamwise momentum integrated over $A_{x_{\mathrm{1}}}$ is greater than that over $A_{x_{\mathrm{2}}}$. An increase in momentum along x would make this term negative.

The terms in the above equation are calculated and compared in Fig. 13 for the LESs and the model. We remark that, in the model, the terms $\mathrm{\Delta }{\mathcal{F}}_{A_x},\mathrm{\Delta }{\mathcal{F}}_{{\mathrm{\Gamma }}_s}$, and $\mathrm{\Delta }{\mathcal{F}}_{\mathrm{fric}}^{*}$ are not represented and therefore equal to 0. Figure 13a shows the comparison of the terms between the model and the LES for cases Inf-H350-S2.5 and Inf-H500-S2.5, scaled with the momentum flux going into the domain ${\mathcal{M}}_{A_{\mathrm{1}}}$. It is clear that the model is capable of representing the dominant terms ($\mathrm{\Delta }{\mathcal{M}}_{A_x},{\mathcal{P}}_x^{*}$, and ℱ_t,x) in the momentum equation. More specifically, as also shown in Fig. 13b, the model can accurately predict the thrust yet overpredicts the magnitude of $\mathrm{\Delta }{\mathcal{M}}_{A_x}$ and ${\mathcal{P}}_x^{*}$. This error increases for higher blockage cases. We identify three main causes for the overprediction of the model. First, the overprediction may partly be due to the model not accounting for the smaller remaining terms $\mathrm{\Delta }{\mathcal{F}}_{A_x},\mathrm{\Delta }{\mathcal{F}}_{{\mathrm{\Gamma }}_s}$, and $\mathrm{\Delta }{\mathcal{F}}_{\mathrm{fric}}^{*}$. For the infinite row cases, $\mathrm{\Delta }{\mathcal{F}}_{{\mathrm{\Gamma }}_s}$ is expected to be 0 because the control volume is not expanding (see Fig. A1a in Appendix A); therefore, the in- and outflowing turbulent momentum should be equal at the side boundaries. The term $\mathrm{\Delta }{\mathcal{F}}_{A_x}$ could be added to the pressure term ${\mathcal{P}}_x^{*}$ if $\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}$ has no direct effect on the thrust. Moreover, $\mathrm{\Delta }{\mathcal{F}}_{\mathrm{fric}}^{*}$ is stated to be 0 in the model, as we assumed that the ground friction is fully balanced by the background pressure gradient. The LES shows that the ground friction increases with the blockage strength and therefore $\mathrm{\Delta }{\mathcal{F}}_{\mathrm{fric}}^{*}$ is not 0, although the term is small. Secondly, the overprediction may be because of the poor representation of the near-wake velocity profile. The model assumes a top-hat function, whereas the LES already exhibits a smoothed top-hat profile. Therefore, the model will overpredict the streamwise momentum flux ${\mathcal{M}}_{A_{\mathrm{2}}}$ at the end of the near wake. Doubling the grid resolution of the LES in all three directions, as done for case Inf-H350-S2.5 in Fig. 13a and b, produces similar results. This suggests that within the current range of grid resolutions, the filter width of the actuator disc model has an insignificant effect on the near-wake shape. Thus, we conclude that the smoothed near-wake velocity profile is predominantly a result of turbulent mixing in the near wake. In future work, the model could be improved by defining the near-wake velocity profile with a shape function and parameterizing the near-wake energy as a function of the turbine thrust and ambient turbulence. Lastly, the overprediction could be due to the Shapiro correction. The correction was developed for a single turbine with no blockage. LES without the Shapiro correction shows that the correction introduces additional uncertainty for the highest blockage case (an ≈8 %, 6 %, and 4 % increase for $\mathrm{\Delta }{\mathcal{M}}_{A_x},\mathrm{\Delta }{\mathcal{P}}_x^{*}$, and Δℱ_t,x compared to the case without the Shapiro correction). In a future study, the turbines could be simulated using the actuator line model instead of the actuator disc model; this would circumvent the need for the correction.

Finally, we look at the finite row case Fin-H500-S2.5 in Fig. 13c and d. The conclusions drawn for the infinite row cases also hold for the finite row cases. We remark that now the model error on ${\mathcal{P}}_x^{*}$ is drastically smaller because Δp_NW obtained from the LES is directly used as an input to the model. However, the error on this term is not zero because the model assumption of the side pressure $p_{\mathrm{side}}=p_{\mathrm{in}}+\mathit{\alpha }\mathrm{\Delta }{p}_{\mathrm{NW}}$, with α=0.5, is not entirely correct. The LES shows that p_side shifts closer to p_in, i.e. α<0.5; however, α depends on the blockage strength. Nevertheless, the contribution of p_side in the streamwise momentum balance is significantly smaller than the other terms due to the small streamwise projection of the side surfaces.

Figure 13 Comparison of the momentum sources and sinks as calculated by Eq. (12) between the developed model and the LES cases. (a) Row averaged results of the Inf-H350-S2.5 standard case, both with and without the Shapiro correction and with a finer grid, and of the Inf-H500-S2.5 standard case, both with and without the Shapiro correction. (c) Side and centre turbines of Fin-H500-S2.5. (b) and (d) show the corresponding model errors.

5 Conclusions

This study set out to analyse the effect of blockage on the wake development behind turbines and the turbine power. Recent research by Lanzilao and Meyers (2024) discovered a strong positive correlation between the favourable pressure gradient in the farm and the wake efficiency η_w. This implies that the favourable pressure gradient, induced by blockage, enhances the wake recovery mechanism. We performed 17 LESs consisting of infinite and finite single turbine rows, as well as two staggered turbine rows, with constant $C_{\mathrm{T}}^{\prime }$, in an idealized atmospheric boundary layer setting. Blockage conditions were artificially introduced using a rigid lid, inducing a favourable pressure difference (Δp_NW<0) over the turbine row and an adverse pressure difference (Δp_FW>0) in the far wake. The blockage strength was adjusted by varying the turbine spacing ($S/D=\mathrm{2.5},\mathrm{5},\mathrm{10},\mathrm{20},\mathrm{40}$) and boundary layer height ($H=\mathrm{350},\mathrm{500},\mathrm{700}$ m).

A strong positive correlation was identified between Δp_NW and both the power coefficient (C_P) and thrust coefficient (C_T). Specifically, as S and H decrease, −Δp_NW, C_P, and C_T increase. Simultaneously, the rotor disc experiences lower induction, and the near wake shows a reduced wake deficit. Specifically, the infinite row cases, at realistic spacings $S/D=\mathrm{5}$, already show a significant increase in C_P of 11 %, 8 %, and 5 % for H=350, 500, and 700 m, respectively, compared to a single freestanding turbine without blockage. Smaller increases are found for the finite row cases, as less surrounding turbines are present to constrain the inflow from expanding, thus resulting in weaker pressure drops and blockage. In this case, power and thrust are distributed and are maximum at the centre of the row, and −Δp_NW increases towards the centre of the row.

The reduction in near-wake velocity deficit due to blockage also leads to smaller velocity deficits and narrower wake widths in the far wake. However, blockage has a negligible direct effect on far-wake development when scaling the far-wake deficit and width with their initial far-wake values. We note that the adverse pressure gradient in the far wake that we observed in our study is at least an order of magnitude smaller than those found by Liu et al. (2002) and Shamsoddin and Porté-Agel (2018) in diverging channels (for which case adverse effects on wake recovery were noticed). We thus conclude that the blockage enhances the wake efficiency, an effect that originates primarily from the lower near-wake deficit. However, we do see a profound effect of H on the wake spreading, with higher boundary layers leading to faster spreading. This relates to the fact that the wake can more freely expand vertically in high boundary layer cases into a larger region of high-speed flow than for shallow boundary layers.

We further developed an analytical model consisting of five equations to predict the blockage effect on near-wake properties (similar to earlier models to correct for wind tunnel blockage in experiments – see Mikkelsen and Sørensen, 2002; Werle, 2010; Segalini and Inghels, 2014). Based on theoretical Froude momentum theory applied to a confined domain, the model uses five known input variables, i.e. rotor diameter (D), turbine spacing (S), boundary layer height (H), disc-based thrust coefficient ($C_{\mathrm{T}}^{\prime }$), and inflow velocity (U_in), but requires one additional closure relation. In the infinite row case, this is given by A₁=A₂ (no expansion possible at the farm level), whereas in the finite row case, the near-wake pressure drop needs to be provided. To this end, we provided the pressure drop as measured in the large-eddy simulations. Overall, we found a good agreement between LES power and thrust predictions and the simple model based on momentum theory, indicating that the latter is a suitable candidate to improve turbine power prediction under blocking conditions. However, the model does underperform for the very high blockage cases. A momentum budget analysis identified the sources of the model overprediction, which can be improved primarily by parameterizing the near-wake velocity profile using a shape function, on the one hand, and possibly by parameterizing the ground friction, on the other hand. Furthermore, more detailed validation using LES in combination with an actuator line model, so that uncertainties related to the Shapiro correction can be excluded, is also of interest. In a practical implementation, such a model needs to be coupled to a farm-scale model that provides input on the large pressure distribution in the farm. Atmospheric perturbation models (Allaerts and Meyers, 2019; Devesse et al., 2022; Stipa et al., 2024; Devesse et al., 2024b; see also the open-source model WAYVE) that explicitly model the pressure feedback coming from gravity waves are an interesting application for this. This is an ongoing topic of further research.

Appendix A: Control volume streamlines

This appendix aims to visualize the control volume around the turbines, as used in both the developed model and the momentum budget analysis. The control volume is bounded by streamlines derived from the time-averaged velocity fields of the LES, extending from the domain inlet to the end of the near-wake region. Figure A1a and b illustrate these streamlines at hub height for the infinite and finite row cases. For clarity, only four turbines are shown in each case. In the finite row case, these include the side turbine up to the row-centred turbine. We remark that the assumption A₁=A₂ holds for the infinite row case, while A₂>A₁ is valid for the finite row case.

Figure A1 Plan view at hub height of four turbines in (a) case Inf-H500-S2.5 and (b) case Fin-H500-S2.5. The solid blue lines are the streamlines calculated from the LES. The blue dashed lines mark the beginning and end of the control volumes. The vertical black lines represent the turbines. The turbine at the centre of the finite row has coordinates $(x/D,y/D)=(\mathrm{0},\mathrm{0})$, and at the edge, $(x/D,y/D)=(\mathrm{0},-\mathrm{7.5})$.