2.1. Unconfined plate

The current model is based on the work of Taylor (Reference Taylor1944), in which the wake is assumed to be described by potential flow theory at the initial expansion region near the plate. The porous plate is then represented by a distribution of sources of strength $m$ and width $L_p$ (see figure 1). The sources induce the streamwise and crosswise velocities

(2.1)\begin{equation} \hat{u} = \frac{m}{2{\rm \pi}} \left( \arctan{\frac{y+L_p/2}{x}} - \arctan{\frac{y-L_p/2}{x}} \right) \end{equation}

and

(2.2)\begin{equation} \hat{v} = \frac{m}{4{\rm \pi}} \left( \ln{\frac{(y+L_p/2)^2 + x^2}{(y-L_p/2)^2 + x^2}} \right), \end{equation}

respectively. The flow velocity may then be computed by superposing the above velocities to the unperturbed flow $\boldsymbol {U} = (U_0,0 )$ which, for this flow case, is equal to the free-stream velocity $U_\infty$. Using this flow representation, Taylor (Reference Taylor1944) obtained the following prediction for the drag coefficient of an unconfined porous plate:

(2.3)\begin{equation} C_D = 4 u^* (1-u^*), \end{equation}

where $u^* = u/U_\infty$, with $u$ being the velocity through the porous plate. The above prediction is identical to the one obtained using the Rankine–Froude theory (Hansen Reference Hansen2015), which constitutes the basis of most blockage correction models for porous plates and wind turbines. Figure 1 shows the well-known failure of the Rankine–Froude–Taylor theory when porosity becomes small, which limits the applicability of blockage correction models to high-porosity regimes (Apelt & West Reference Apelt and West1975; Ross & Polagye Reference Ross and Polagye2020a; Miller et al. Reference Miller, Duvvuri and Hultmark2021). Appendix A presents a relation that connects $u^*$ with $\beta =A_p/A$ (where $A_p$ and $A$ are the porous and gross areas of the plate, respectively), proposed by Taylor & Davies (Reference Taylor and Davies1944).

Figure 1.

(a) Comparison of the predictions of Taylor (Reference Taylor1944) and Steiros & Hultmark (Reference Steiros and Hultmark2018) with drag measurements of plates of varying porosity $\beta$. In certain low-porosity cases, a splitter plate is used to suppress vortex shedding. Data are taken from Steiros & Hultmark (Reference Steiros and Hultmark2018). (b) Superposition of a distribution of sources of strength $m$ on a free stream, corresponding to the Taylor model.

The above limitation can be corrected by refining the Taylor model so that (i) continuity is imposed across the plate (Koo & James Reference Koo and James1973), and (ii) the wake pressure can assume lower than ambient values, which are calculated implicitly using momentum theory (Steiros & Hultmark Reference Steiros and Hultmark2018). This representation then yields the revised prediction

(2.4)\begin{equation} C_D = \frac{4}{3}\,\frac{(1-u^*)(2+u^*)}{2-u^*}, \end{equation}

which, as shown in figure 1, is in much better agreement with experimental data for small plate porosities, as long as vortex shedding is absent or suppressed via the use of a splitter plate. A blockage correction model based on the refined theory of Steiros & Hultmark (Reference Steiros and Hultmark2018) could thus offer increased accuracy in high-solidity regimes.

2.2. Confined plate kinematics

Following Koo & James (Reference Koo and James1973), the effect of channel walls in the above flow representation is introduced via the method of mirror images. An infinite number of plates are positioned side by side, with the distance between their centres being equal to the channel width $L_c$ (see figure 2). The induced velocities from each plate are then

(2.5)\begin{equation} \hat{u}_n = \frac{m}{2{\rm \pi}} \left( \arctan{\frac{y+L_p/2-n L_c}{x}} - \arctan{\frac{y-L_p/2-n L_c}{x}} \right) \end{equation}

and

(2.6)\begin{equation} \hat{v}_n = \frac{m}{4{\rm \pi}} \left( \ln{\frac{(y+L_p/2 -n L_c)^2 + x^2}{(y-L_p/2 -n L_c )^2 + x^2}} \right), \end{equation}

where $n \in \{\ldots, -2, -1, 0, 1, 2,\ldots \}$ and $n=0$ is the original plate. To calculate the velocity of the fluid, we need to sum the contribution of all plates, and on top of that add the unperturbed velocity of the stream if no plates were present, $(U_0,0)$.

Figure 2.

Schematic of the flow model. An infinite number of plates, positioned side by side at a distance $L_c$, represent the flow around a porous plate inside a channel.

It can be checked readily that the addition of the image plates does not alter the axial velocity through the plate $u$, which retains its unconfined value (see Steiros & Hultmark Reference Steiros and Hultmark2018):

(2.7)\begin{equation} u = U_0 - \frac{m}{2}. \end{equation}

Far upstream or downstream from the plate, however, the induced velocities of the images add up to a non-negligible axial velocity contribution. In particular, the velocity at the centreline $y=0$ becomes

(2.8)\begin{equation} u_{x \rightarrow \pm\infty} = U_0 \pm \lim_{x\rightarrow \pm\infty}{ \sum_{n={-}\infty}^{\infty} \hat{u}_n (x,0 )}, \end{equation}

where the minus corresponds to upstream, and plus to downstream. For large distances from the plate $|x|$, we have $|\hat {u}_{n+1} (x,0 ) - \hat {u}_{n} (x,0 ) | \ll |\hat {u}_{n} (x,0 )|$, and thus $\sum _{n=-\infty }^{\infty } \hat {u}_n (x,0 )$ can to a good approximation be considered a continuous function of the continuous variable $\tilde {n}$. The series can therefore be replaced by an integral, i.e. (note that (2.7) is also used)

(2.9)\begin{equation} u_{x \rightarrow \pm \infty} = U_0 \pm \lim_{x\rightarrow \pm \infty}{ \int_{-\infty}^{\infty} \hat{u}_{\tilde{n}} (x,0 ) \,{\rm d} \tilde{n} } = U_0 \pm (U_0-u) b, \end{equation}

where $b = {L_p}/{L_c}$. This yields an expression for the free-stream velocity in the channel $U_\infty = u_{x \rightarrow - \infty }$, i.e.

(2.10)\begin{equation} \frac{U_\infty}{U_0} = 1 - \left(1-\frac{u}{U_0} \right)b = \frac{1-b}{1-u^*b}, \end{equation}

where $u^* = u/U_\infty$. Expression (2.10) shows that the free-stream velocity depends on the level of blockage. For instance, when a plate is fully blocking the channel ($b=1$), the free-stream velocity will be equal to the velocity through the plate, i.e. $U_\infty = u$, even if the ‘nominal’ channel velocity $U_0$ is larger than that.

The wake velocity can be obtained from $u_{x \rightarrow \infty }$ after applying a correction to it. In particular, as Koo & James (Reference Koo and James1973) first noted, the source representation produces a discontinuity at the plate location $(x,y)=(0,0)$, as the velocity immediately upstream of the plate $u^{-}=U_0 - m/2$ is different to the velocity immediately downstream of the plate $u^{+}=U_0 + m/2$, violating the conservation of mass and momentum. To impose continuity across the plate, Koo & James (Reference Koo and James1973) proposed a rescaling of the velocities of the whole wake by a factor $E = u^{-}/u^{+}$. Using (2.7) and (2.10), $E$ becomes in our case

(2.11)\begin{equation} E = \frac{U_0 - m/2}{U_0 + m/2} = \frac{u^*}{2U_0/U_\infty-u^*} = \frac{u^*(1-b)}{2 - u^* - u^*b} , \end{equation}

and thus the far wake velocity after the rescaling becomes

(2.12)\begin{equation} \frac{U_w}{U_\infty} = \frac{Eu_{x \rightarrow \infty} }{U_\infty} = E \left ( \frac{U_0}{U_\infty} + \left(\frac{U_0}{U_\infty} - u^* \right) b \right) = \frac{u^*(1-2 u^*b +b)}{2 - u^* - u^*b} . \end{equation}

This equation shows that by increasing the tunnel blockage $b=L_p/L_c$ (i.e. by narrowing the channel), the wake velocities are also increased. Then conservation of mass dictates that the wake would be squeezed, as expected.

Finally, we consider the square crosswise velocity $v_1^2$ at the neighbourhood of the plate, i.e. at $x\to 0$. This varies along $y$, but for simplicity we consider it constant, and equal to its average value $\overline {v_1^2}$. Appendix B shows that the latter can be approximated by

(2.13)\begin{equation} \overline{v_1^2} \approx \frac{(U_\infty-u)^2}{3}. \end{equation}

2.3. Conservation laws

Given the above description of the flow kinematics, we may apply conservation of mass, momentum and energy in the control volume shown in figure 2 to calculate the drag and wake pressure coefficients of the confined plate, given the normalised velocity through the plate $u^*=u/U_\infty$ and the blockage ratio $b=L_p/L_c$. Conservation of mass within the wake and channel regions yields

(2.14)\begin{gather} L_p u = L_w U_w , \end{gather}

(2.15)\begin{gather}L_c U_\infty = L_w U_w + U_e (L_c-L_w ) . \end{gather}

Conservation of momentum in the channel yields

(2.16)\begin{equation} D - p_\infty L_c + p_w L_w + p_e (L_c - L_w ) = \rho L_c U_\infty^2 - \rho L_w U_w^2 - \rho (L_c - L_w ) U_e^2 , \end{equation}

where $D= (p_1-p_2 ) L_p$ is the drag per unit length. We note that after the end of the control volume, turbulent stresses are assumed to emerge, balancing the vertical pressure gradient (see Tennekes & Lumley Reference Tennekes and Lumley1972). In the inviscid control volume, on the other hand, turbulent stresses are neglected, even though in reality they will exist in small magnitudes. We may thus expect a slight imbalance of the pressure gradient in the crosswise momentum equation.

Applying Bernoulli's equation upstream, downstream and around the plate, we obtain

(2.17)\begin{gather} p_\infty + 0.5\rho U_\infty^2 = p_1 + 0.5\rho (u^2+v_1^2 ) , \end{gather}

(2.18)\begin{gather}p_2 + 0.5\rho (u^2 + v_2^2 ) = p_w + 0.5\rho U_w^2 , \end{gather}

(2.19)\begin{gather}p_\infty + 0.5\rho U_\infty^2 = p_e + 0.5\rho U_e^2 , \end{gather}

where $v_1^2$ is given by (2.13), $v_2=E v_1$, and $U_w$ comes from (2.12). By solving the system (2.14)–(2.19), the following analytical expressions of the flow quantities of interest can be derived, given as functions of $u^*$ and $b$:

(2.20)\begin{gather} C_D = \frac{4 (u^*b-1)(1-u^*)}{(1-b)(2-u^*-u^*b)} \left( \frac{1-u^*}{3} - \frac{1-2u^*b+b}{1-b}\right), \end{gather}

(2.21)\begin{gather}u_e^* = \frac{U_e}{U_\infty} = \frac{ 1-2u^*b + b }{1 - b}, \end{gather}

(2.22)\begin{gather}\frac{L_w}{L_c} = \frac{\left(2-u^* - u^*b \right)b}{1-2u^*b + b}, \end{gather}

(2.23)\begin{gather}Cp_w ={-} C_D + \frac{4 (1-u^*)(u^*b-1)}{(2-u^*-u^*b)^2} \left( \frac{(1-u^*)^2}{3} - (1-u^{*2}b)\right), \end{gather}

which reduce to the unconfined plate expressions of Steiros & Hultmark (Reference Steiros and Hultmark2018) for $b=0$.

2.4. Blockage correction method

The above theoretical framework can be used to convert the flow quantities from one level of blockage (i.e. experimental conditions) to another (typically unconfined conditions). Similar to most existing blockage correction models (see e.g. Barnsley & Wellicome Reference Barnsley and Wellicome1990; Mikkelsen & Sørensen Reference Mikkelsen and Sørensen2002; Houlsby et al. Reference Houlsby, Draper and Oldfield2008; Werle Reference Werle2010), we achieve that following the rationale of Glauert (Reference Glauert1935), i.e. we seek the ‘corrected’ free-stream velocity $U_\infty ^c$ that would produce the same values of thrust $T$, through-velocity $u$, and angular velocity $\omega$ (if turbines are considered) of the blocked model, but in unconfined conditions, i.e.

(2.24a–c)\begin{equation} T^c = T^m,\quad u^c = u^m,\quad \omega^c = \omega^m, \end{equation}

where the superscripts $c$ and $m$ denote corrected/unconfined data and measured/confined data, respectively. Using the actuator disk relation for power $P= T u$, the above also yield $P^c=P^m$. Expressed in non-dimensional form, the correction model consists of a set of expressions for the ‘corrected’ flow quantities as functions of the confined flow data and the free-stream velocity ratio $R_V={U_\infty ^m}/{U_\infty ^c}$:

(2.25)\begin{gather} \lambda^{c} = \lambda^{m} R_V , \end{gather}

(2.26)\begin{gather}C_T^c = C_T^m R_V^2 , \end{gather}

(2.27)\begin{gather}C_P^c = C_P^m R_V^3 , \end{gather}

where $\lambda = \omega R/U_\infty$ is the turbine tip speed ratio, $R$ is the turbine radius, and $C_T = 2T/(\rho U_\infty ^2 {\rm \pi}R^2)$, $C_P = 2P/(\rho U_\infty ^3 {\rm \pi}R^2)$ are the turbine thrust and power coefficients, respectively. The blockage correction thus rests on the determination of the free-stream velocity ratio $R_V$. The latter is calculated using potential flow theory, i.e. the Rankine–Froude theory in existing models, and the framework of §§ 2.2 and 2.3 in the current work. The necessary inputs are the blockage ratio and one confined flow property, typically the thrust/drag coefficient, even though other quantities can also be used (e.g. flow-through velocity, wake width). Summarising, the blockage correction process rests on the calculation of $R_V = {U_\infty ^m}/{U_\infty ^c} = u^{*,c}/u^{*,m}$ and consists of the following steps.

1. 1. Given the measured thrust coefficient $C_T^m$ and blockage ratio $b_m$, solve (2.20) to obtain $u^{*,m}$.

2. 2. Solve (2.26) for $u^{*,c}$, where $C_T^c$ is given by (2.20) and $b=b_c$ is the target blockage ratio (typically zero).

3. 3. Compute the free-stream velocity ratio $R_V = u^{*,c}/u^{*,m}$.

4. 4. Compute the ‘blockage-corrected’ rotor/plate performance using (2.25)–(2.27).

The above described blockage correction method is expected to work for two-dimensional plates (strips) of arbitrary porosity, as long as vortex shedding is suppressed or absent. When considering turbines, the additional complexity of three-dimensionality appears, i.e. horizontal axis turbines are circular, while vertical axis turbines have a finite aspect ratio. Therefore, the current blockage model can be thought strictly appropriate for correcting rows of closely spaced turbines (Draper et al. Reference Draper, Nishino, Adcock and Taylor2016). However, in § 4.3 we show that the current model produces accurate corrections for the thrust and power of single horizontal axis wind turbines of any porosity, while it yields results almost identical to those from conventional correction models at high effective porosities, despite the conceptual differences of two- and three-dimensional flows.

3.1. Experimental set-up

To validate the predictions of (2.20)–(2.22), experiments were carried out in the aerodynamics flume of the Aeronautics Department at Imperial College London. The facility has a test section that is $L_c=0.6$ m wide with variable water height, which can achieve flow speed up to $0.8$ m s$^{-1}$, with free-stream turbulence intensities below 1.5 % (Cicolin et al. Reference Cicolin, Buxton, Assi and Bearman2021). The experiments were performed at the speed range from 0.2 to 0.5 m s$^{-1}$, leading to a variation in the Reynolds number based on the plate width, $Re = U_\infty L_p / \nu$, between $9 \times 10^3$ and $5.5 \times 10^4$; all porosity/blockage combinations (see below) were tested at two different Reynolds numbers (different by a factor of $\simeq 1.5\unicode{x2013}1.7$). The variation in the measured quantities was found to be negligible.

Force measurements were made using two force transducers (see figure 3), one attached to the bottom of the plates (i.e. submerged in the water), and the other attached to the top of the plates (outside the water). Each transducer (model ATI Mini 40) measures forces up to 80 N with resolution 0.005 N. The sum of the values from the two transducers yielded the total drag $D$ acting on the plate. The drag coefficient was then calculated as $C_D = 2D/(\rho U_\infty ^2 A)$, with $A$ the gross wetted area of the plate (i.e. without subtracting the perforations). Particle image velocimetry (PIV) measurements were acquired at two different positions. The illumination source was a high-speed Litron LDY 304 Nd:YLF laser, which operated at constant frequency 1 kHz. An optical system split the laser beam into two components. One illuminated a small field of view upstream, intended to measure the incoming flow speed, while the second laser sheet illuminated the wake of the plate, as illustrated in figure 3. The images were captured and recorded by two cameras, model Phantom v641, with a maximum resolution $2460\times 1600$ px and internal storage memory for 5700 images at the highest resolution. The laser pulse and each camera operated at different acquisition rates, and their synchronisation was performed through an external pulse generator (Stanford DG645).

Figure 3.

Schematic of the experimental apparatus and measurement techniques.

In total, 16 plates of 5 mm thickness aluminum were tested, where circular holes were water jetted leading to porosity ratios $\beta$ of 0 %, 15 %, 30 % and $50\,\%$. The plates spanned the whole height of the flume (with the exception of a small area underneath the plates where the force sensor was positioned). Each porosity ratio was applied to four different plate widths $L_p$, leading to four blockage ratios for each porosity case (see table 1). For the $\beta$ cases of $0\,\%$ and $15\,\%$, a 450 mm long splitter plate was positioned in the wake of the plates, in order to suppress vortex shedding, which otherwise appears in low porosities (Steiros, Bempedelis & Ding Reference Steiros, Bempedelis and Ding2021). Acquisition started after the flow was stabilised, typically after two minutes from the moment the flume speed was set up. Force and PIV measurements were acquired simultaneously. The forces were acquired at constant frequency 100 Hz, whereas the PIV acquisition varied from 4 to 20 Hz, depending on the temporal scale of the model. The acquisition period was chosen to be sufficient to obtain at least 100 cycles of vortex shedding, considering the reference case of the solid plate. The field of view was constant in all experiments, covering an area $340\times 210$ mm$^2$, as shown in figure 4. The wake width $L_w$ was calculated as the maximum distance from the centreline reached by the streamline of the mean flow starting from the tip of the plate. The velocity between the channel wall and the wake border $U_e$ was calculated as the asymptotic value of the streamwise velocity away from the plate.

Table 1.

Details of experimental configurations. Channel width is $L_c=600$ mm, $L_p$ is the plate width, and $b=L_p/L_c$ is the blockage ratio.

Figure 4.

Mean velocity fields for two plates of $\beta =50\,\%$ at different blockage ratios obtained using PIV: (a) $b=11.67\,\%$, (b) $b=25\,\%$.

3.2. Numerical simulations

Large-eddy simulations (LES) of the flow around model wind turbines were performed using Xcompact3d, a finite-difference framework that uses sixth-order compact schemes to solve the incompressible Navier–Stokes equations (Laizet & Lamballais Reference Laizet and Lamballais2009; Laizet & Li Reference Laizet and Li2011; Bartholomew et al. Reference Bartholomew, Deskos, Frantz, Schuch, Lamballais and Laizet2020; Deskos, Laizet & Palacios Reference Deskos, Laizet and Palacios2020). In the present work, the wind turbines were modelled by the actuator line method, the standard Smagorinsky model was used to model the effects of the unresolved fluid motions, and a third-order Adams–Bashforth method was used for time integration.

The considered configurations were based loosely on the experiments of the ‘Blind Test 1’ (BT1) workshop (Krogstad & Eriksen Reference Krogstad and Eriksen2013). Xcompact3d has been used previously to study this particular problem, showing very good agreement with the experimental results (Deskos, Laizet & Piggott Reference Deskos, Laizet and Piggott2019; Bempedelis & Steiros Reference Bempedelis and Steiros2022). In this work, we considered the high-solidity six-bladed variant of the BT1 turbine. Two square channels, a small one of two turbine diameters width, and a large one of five diameters width, allowed the simulation of blockage ratios of $\simeq 19.6\,\%$ and $\simeq 3.1\,\%$, respectively. The turbines were positioned at the centre of the square channel, creating close to isotropic blockage. The channel length in both cases was 12 turbine diameters. Slip conditions were used at the lateral domain boundaries, and laminar inflow ($U_\infty =10$ m s$^{-1}$, yielding $Re=U_\infty D / \nu = 5.95 \times 10^5$) was considered at the inlet. (Very low turbulence levels, $I \simeq 0.3\,\%$, were present in the experiments; Krogstad & Eriksen Reference Krogstad and Eriksen2013.) A uniform grid consisting of 32 fluid mesh points per turbine radius was used, resulting in meshes of size $\simeq 13$ and $\simeq 79 \times 10^6$ points, respectively; such a resolution has been shown to be sufficient for the purposes of this study through comparisons with experimental data (Deskos et al. Reference Deskos, Laizet and Piggott2019; Bempedelis & Steiros Reference Bempedelis and Steiros2022). Each blade/actuator line was discretised with 52 elements, and the tower and nacelle were omitted. After an initialisation period, time series of thrust $T$ and power $P$ were extracted over $\simeq 64\unicode{x2013}128$ rotations, depending on the angular velocity of the turbine $\omega$. The thrust and power coefficients were subsequently calculated as $C_T = 2T/(\rho U_\infty ^2 {\rm \pi}R^2)$ and $C_P = 2P/(\rho U_\infty ^3 {\rm \pi}R^2)$, respectively, with $R$ being the turbine radius. An example view of the mean streamwise velocity of the same wind turbine in the two different channels is shown in figure 5.

Figure 5.

Mean streamwise velocity contours at a horizontal cut through the turbine centreline, $\lambda =6$, for two channels of different sizes: (a) $b=3.1\,\%$, (b) $b=19.6\,\%$.