# Performance and wake characteristics of tidal turbines in an infinitely
large array

Pablo Ouro1,2<EMAIL_ADDRESS>Takafumi Nishino3  1Department of
Mechanical, Aerospace and Civil Engineering, The University of Manchester,
Manchester M13 9PL, UK 2Hydro-environmental Research Centre, School of
Engineering, Cardiff University, The Parade, Cardiff CF24 3AA, UK 3Department
of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, UK

###### Abstract

The efficiency of tidal-stream turbines in a large array depends on the
balance between negative effects of turbine-wake interactions and positive
effects of bypass-flow acceleration due to local blockage, both of which are
functions of the layout of turbines. In this study we investigate the
hydrodynamics of turbines in an infinitely large array with aligned or
staggered layouts for a range of streamwise and lateral turbine spacing.
First, we present a theoretical analysis based on an extension of the Linear
Momentum Actuator Disc Theory (LMADT) for perfectly aligned and staggered
layouts, employing a hybrid inviscid-viscous approach to account for the local
blockage effect within each row of turbines and the viscous (turbulent) wake
mixing behind each row in a coupled manner. We then present Large-Eddy
Simulation (LES) of open-channel flow for 28 different layouts of tidal
turbines using an Actuator-Line Method (ALM) with doubly periodic boundary
conditions. Both theoretical and LES results show that the efficiency of
turbines (or the power of turbines for a given bulk velocity) in an aligned
array decreases as we reduce the streamwise turbine spacing, whereas that in a
staggered array remains high and may even increase due to the positive local
blockage effect (causing the local flow velocity upstream of each turbine to
exceed the bulk velocity) if the lateral turbine spacing is sufficiently
small. The LES results further reveal that the amplitude of wake meandering
tends to decrease as we reduce the lateral turbine spacing, which leads to a
lower wake recovery rate in the near-wake region. These results will help to
understand and improve the efficiency of tidal turbines in future large tidal
arrays, even though the performance of real tidal arrays may depend not only
on turbine-to-turbine interactions within the array but also on macro-scale
interactions between the array and natural tidal currents, the latter of which
are outside the scope of this study.

###### keywords:

##  1 Introduction

Tidal stream energy is being developed rapidly as a promising complement to
intermittent wind and solar energy, owing to the high predictability of tides
(Adcock et al., 2021). To ensure tidal energy is successfully embedded into
the future net-zero carbon energy mix, we need to understand the flow physics
driving the efficiency of tidal turbine arrays when deployed at a large scale.
Due to the confined flow environment in which these turbines operate, we
cannot simply infer tidal array hydrodynamics from existing knowledge of wind
farm aerodynamics (Porte-Agel et al., 2020; Nishino & Dunstan, 2020). For
instance, relatively shallow waters may restrict the ability of turbine wakes
to expand vertically, potentially affecting the recovery of wakes in a large
turbine array.

Future tidal arrays will comprise multiple rows of turbines, whose design
requires to consider an appropriate spacing between turbines as well as their
impact on the tidal channel flow dynamics. These micro- (turbine to turbine)
and macro-scale (array to tidal channel) interactions need to be considered
when designing tidal arrays in order to find the optimal trade-off between
energy extraction and minimal changes to the tidal flow (De Dominicis et al.,
2018). Large tidal arrays deployed at a tidal channel lead to an added
resistance to its flow dynamics, which, if too large, can considerably
obstruct and modify the overall flow. Therefore, to optimise the design of
large tidal arrays, it is required to tune the operating conditions of
individual turbines for a given tidal channel (e.g. straight or variable-
section), tidal forcing and turbine density (Vennell, 2010, 2011; Vennell et
al., 2015).

At the micro-scale, the power production capability of tidal turbines is
driven by turbulent wake mixing and acceleration of bypass flow in-between
turbines. In relatively closely packed arrays, limiting the negative wake-
turbine interactions is often key to minimising power losses (Stallard et al.,
2013). However, local blockage due to a small lateral spacing between devices
(as well as a shallow water depth confining the flow) may lead to local flow
acceleration that enhances individual turbine power, as observed in
experimental tests, e.g. Stallard et al. (2013); Noble et al. (2020). The
effect of local blockage has also been investigated using the Linear Momentum
Actuator Disc Theory (LMADT) for a single lateral row of turbines (Garrett &
Cummins, 2007; Nishino & Willden, 2012, 2013; Vogel et al., 2016; Creed et
al., 2017) and two rows of turbines (Draper & Nishino, 2014). These studies
suggest that two staggered rows of turbines tend to be more efficient than two
perfectly aligned (or centred) rows of turbines, but less efficient than a
single row of the same total number of turbines with the same array width.

Whilst the above findings are important for the performance of arrays with a
small number of rows, future tidal arrays will require turbines to be deployed
in several rows to generate a sufficiently large amount of energy, which may
cause macro-scale flow interactions between the array and the tidal channel.
Vennell (2010, 2011) combined the LMADT with a simple theoretical tidal
channel flow model to analyse how the resistance and lateral spacing of
turbines within each row should be tuned for a given number of rows deployed
across a given tidal channel, to maximise the total power generation. However,
existing theoretical tidal array models based on the LMADT, including the two-
row model of Draper & Nishino (2014), do not fully account for the complex
effect of turbulent wake mixing. The Vennell-type array models assume that the
streamwise spacing between rows is large enough for individual turbine wakes
to be fully mixed after each row, whereas the two-row model of Draper &
Nishino (2014) assumes that the two rows are close enough to each other for
the effect of wake mixing to be negligible.

In narrow tidal channels or straits, turbines in a multi-row array may operate
in a fully-waked scenario (similar to perfectly aligned layout) during half of
the tidal cycle (e.g. ebb tide), but they may operate in partly-waked
conditions (akin to staggered layouts) during the other half of the cycle
(e.g. flood tide) (Garcia Novo & Kyozuka, 2019). This asymmetry between ebb
and flood tides further complicates the design of optimal array
configurations, requiring to understand the performance and hydrodynamics of a
given array design for various incident flow characteristics. Many existing
studies looking into tidal array optimisation have adopted low-fidelity flow
models, such as two-dimensional shallow water models (Culley et al., 2018),
analytical wake models, e.g. Gaussian models (Stansby & Stallard, 2015) and
aforementioned theoretical models based on the LMADT (Nishino & Willden, 2013;
Draper & Nishino, 2014), while high-fidelity simulations have been restricted
to relatively small arrays with a limited number of configurations, due to
their large computational expense (Afgan et al., 2013; Chawdhary et al., 2017;
Ouro et al., 2019c). However, as the performance of tidal devices in arrays is
driven by wake-turbine interactions as well as bathymetry-induced turbulence
(Stallard et al., 2013; Ouro & Stoesser, 2019), turbulence-resolving
approaches such as Large-Eddy Simulation (LES) are valuable to yield reliable
hydrodynamics results as well as to build more accurate low-order models that
can improve array optimisation tools.

In this paper, we investigate the flow characteristics and efficiency of
infinitely large tidal arrays with perfectly aligned and staggered
configurations, combining predictions from two contrasting approaches:
actuator disc theory and LES. Infinitely large arrays represent an asymptotic
case in which the flow passing the turbine rows is fully developed (i.e. flow
statistics become identical for all rows). This is similar to large wind farms
in which such flow conditions may be attained approximately after 10 to 15
rows, depending on atmospheric stability conditions (Porte-Agel et al., 2020;
Sharma et al., 2018). We consider a wide range of streamwise and lateral
turbine spacing to understand how the array efficiency can be maximised (from
the micro-scale perspective) by balancing the negative impact of turbine wakes
impinging downstream turbines and the positive local blockage effects. The
paper is structured as follows: in §2 we introduce an extended theoretical
model developed for periodic turbine arrays with perfectly aligned and
staggered configurations. Details of 28 LES runs comprising aligned and
staggered layouts are presented in §3, with results of flow characteristics
and hydrodynamic coefficients in §4. In §5 we provide further discussion on
the comparison between the predictions obtained from the theoretical analysis
and LES, followed by main conclusions in §6.

##  2 Theoretical analysis

We start with a simple theoretical analysis on the efficiency of an infinitely
large array of ideal turbines in a steady, uniform and vertically confined
flow. The analysis is based on the work of Draper & Nishino (2014), who
extended the LMADT for laterally confined flows (Garrett & Cummins, 2007;
Houlsby et al., 2008) to predict an upper limit to the efficiency of two
aligned or staggered rows of turbines. This two-row analysis is further
extended in this study to investigate an infinite number of aligned or
staggered rows of tidal turbines, following the idea of hybrid inviscid-
viscous approach recently proposed by Nishino & Draper (2019).

Figure 1 illustrates an example of a periodic staggered array of tidal
turbines. The key assumption employed in the hybrid inviscid-viscous approach
is that the streamwise extent of the region in which the expansion of the flow
through each turbine takes place is much shorter than the distance between
each row of turbines. With this assumption, we hypothetically divide the flow
field into two types of zones; namely the inviscid flow zones, which are
analysed using the LMADT neglecting the effect of viscous (or turbulent)
mixing, and the viscous flow zones, which are modelled separately to account
for the effect of mixing. This approach is also in line with the results of a
recent LES study of flow past a periodic array of actuator discs (West & Lele,
2020) showing that the effects of inviscid and viscous (turbulent) flow
processes are dominant, respectively, in the vicinity of the turbines and in
the rest of the flow field.

Figure 1: Schematic of the flow past a periodic staggered array of tidal
turbines, divided into the hypothetical inviscid (inv.) and viscous (visc.)
flow zones. The rectangular region enclosed by the magenta dashed line
corresponds to that depicted in figure 2.

Figure 2: Schematic of the quasi-1D theoretical model for a periodic staggered
array of tidal turbines, with examples of how cross-sectional flow patterns
may appear in a corresponding 3D flow problem. The three vertical thin black
lines indicate the locations of stations 1, 4 and 5 as well as the (non-
dimensional) cross-sectional average velocity,
$\psi=u_{\mathrm{av}}/u_{\mathrm{ref}}$, for the superposed plot of
$u/u_{\mathrm{ref}}$ at each station, whereas the thick black lines show the
profiles of $u/u_{\mathrm{ref}}$ at the three stations (calculated for the
case with $B=0.2$, $K=3$ and $m=0.7$). Reproduced from Nishino & Draper (2019)
with minor modifications.

###  2.1 Staggered rows of actuator discs

A schematic of the theoretical model for the staggered case is shown in figure
2. Here we consider a straight local flow passage containing only one-half of
an actuator disc due to the periodic and symmetric nature of the array as
shown in figure 1. The cross-sectional area of the flow passage is $A/B$,
where $A$ is the half-disc area and $B$ is the area blockage ratio. Although
the figure is depicted in a two-dimensional manner, this is still a quasi-one-
dimensional flow model as we do not consider any variation of flow quantities
(velocity and pressure) over the cross-section of each streamtube in the
inviscid zone. Assuming that the flow is incompressible, the average velocity
over the cross-section of the entire flow passage, $u_{\mathrm{av}}$, does not
change in the streamwise direction.

Following the common notations used in the LMADT (Houlsby et al., 2008) we
define four stations within the inviscid zone: station 1 is at the inlet of
the inviscid zone where the core-flow streamtube (that encompasses the
actuator disc) starts to expand, stations 2 and 3 are immediately upstream and
downstream of the disc, respectively, and station 4 is at the outlet of the
inviscid zone where the pressure is equalised between the core- and bypass-
flow streamtubes. In addition, we describe the outlet of the viscous zone as
station 5 (which is station 1 for the next row of discs). The pressure at
stations 1 to 5 is denoted by $p_{1}$ to $p_{5}$, respectively, whereas the
velocity is described using the velocity coefficients $\alpha$ and $\beta$
with subscripts as in figure 2. Each velocity coefficient represents the ratio
of the velocity there to a reference velocity, $u_{\mathrm{ref}}$. In the
following we take the core-flow velocity at station 1 as $u_{\mathrm{ref}}$
(i.e. $\alpha_{1}=1$) for convenience.

We now consider the conservation of mass, momentum and energy in the inviscid
flow zone in a similar manner to the work of Draper & Nishino (2014). First,
since the Bernoulli equation must be satisfied to conserve energy in each of
the three bypass-flow streamtubes between stations 1 and 4, we obtain

$\displaystyle p_{1}-p_{4}=\frac{1}{2}\rho
u_{\mathrm{ref}}^{2}\left(\beta_{4a}^{2}-1\right),$ (1) $\displaystyle
p_{1}-p_{4}=\frac{1}{2}\rho
u_{\mathrm{ref}}^{2}\left(\beta_{4b}^{2}-\beta_{5}^{2}\right),$ (2)
$\displaystyle p_{1}-p_{4}=\frac{1}{2}\rho
u_{\mathrm{ref}}^{2}\left(\alpha_{8}^{2}-\alpha_{5}^{2}\right),$ (3)

where $\rho$ is the fluid density. Substituting (2) and (3), respectively,
into (1) gives

$\displaystyle\beta_{4b}=\left(\beta_{4a}^{2}+\beta_{5}^{2}-1\right)^{1/2},$
(4)
$\displaystyle\alpha_{8}=\left(\beta_{4a}^{2}+\alpha_{5}^{2}-1\right)^{1/2}.$
(5)

As the Bernoulli equation must also be satisfied in the upstream part (between
stations 1 and 2) and downstream part (between stations 3 and 4) of the core-
flow streamtube, we also obtain

$\displaystyle p_{2}-p_{3}=p_{1}-p_{4}+\frac{1}{2}\rho
u_{\mathrm{ref}}^{2}\left(1-\alpha_{4}^{2}\right),$ (6)

which, together with (1), leads to an expression for the half-disc thrust,
$T$, as

$\displaystyle T=(p_{2}-p_{3})A=\frac{1}{2}\rho
u_{\mathrm{ref}}^{2}A\left(\beta_{4a}^{2}-\alpha_{4}^{2}\right),$ (7)

and thus an expression for the half-disc power, $P$, as

$\displaystyle P=T\alpha_{2}u_{\mathrm{ref}}=\frac{1}{2}\rho
u_{\mathrm{ref}}^{3}A\alpha_{2}\left(\beta_{4a}^{2}-\alpha_{4}^{2}\right).$
(8)

Hence, to obtain $P$ for a given $\alpha_{2}$, for example, we need to know
how $\alpha_{4}$ and $\beta_{4a}$ depend on $\alpha_{2}$. By considering the
conservation of mass in the ‘main’ bypass-flow streamtube between station 1
(where the velocity is $\beta_{5}u_{\mathrm{ref}}$) and station 4 (where the
velocity is $\beta_{4b}u_{\mathrm{ref}}$) we obtain

$\displaystyle\beta_{5}u_{\mathrm{ref}}A\left(\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}-\gamma\right)=\beta_{4b}u_{\mathrm{ref}}A\left(\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}-\frac{\gamma-\alpha_{2}}{\beta_{4a}}-\gamma\right),$
(9)

which leads to

$\displaystyle\alpha_{4}=\frac{\alpha_{2}}{\dfrac{1}{B}-\dfrac{\beta_{4b}}{\beta_{4a}}\left(\dfrac{\gamma-\alpha_{2}}{\beta_{4b}-\beta_{5}}\right)-\gamma},$
(10)

where

$\displaystyle\gamma=\frac{\alpha_{2}\alpha_{5}}{\alpha_{4}\alpha_{8}}$ (11)

is the area ratio for the wake flow upstream of the disc as depicted in figure
2. Note that this wake originates from the disc located two rows upstream (and
its velocity increases from $\alpha_{5}u_{\mathrm{ref}}$ to
$\alpha_{8}u_{\mathrm{ref}}$ when it passes through the row immediately
upstream) due to the staggered configuration. As for $\beta_{4a}$ we consider
the conservation of momentum for the entire flow passage to obtain

$\begin{split}&(p_{1}-p_{4})\frac{A}{B}-T\\\ &=\rho
u_{\mathrm{ref}}^{2}A\left[\alpha_{2}(\alpha_{4}-1)+(\gamma-\alpha_{2})(\beta_{4a}-1)+\gamma\alpha_{8}(\alpha_{8}-\alpha_{5})+\beta_{5}\left(\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}-\gamma\right)\left(\beta_{4b}-\beta_{5}\right)\right],\end{split}$
(12)

which, together with (1) and (7), leads to

$\begin{split}\left(1-B\right)\beta_{4a}^{2}-2B\left(\gamma-\alpha_{2}\right)\beta_{4a}&-2B\left[\alpha_{2}\alpha_{4}-\frac{1}{2}\alpha_{4}^{2}+\gamma\left(\alpha_{8}^{2}-\alpha_{5}\alpha_{8}-1\right)\right]\\\
&-2\beta_{5}\left(\beta_{4b}-\beta_{5}\right)\left[1-B\left(\frac{\alpha_{2}}{\alpha_{4}}+\gamma\right)\right]-1=0.\end{split}$
(13)

The above set of equations for the inviscid zone is not closed as it includes
$\alpha_{5}$ and $\beta_{5}$, which need to be modelled considering the effect
of mixing in the viscous zone. There are several possible ways to model the
effect of mixing but here we employ a very simple approach using a single non-
dimensional input parameter called the mixing factor, $m$, which represents
the completeness of mixing in the viscous zone. Specifically, the value of $m$
(between 0 and 1) describes how much the velocity of flow at each cross-
sectional position returns to the cross-sectional average velocity of the
entire flow passage, $u_{\mathrm{av}}$, as the flow passes through the viscous
zone. By applying this to the wake flow of the disc we obtain

$\displaystyle\alpha_{5}=m\psi+\left(1-m\right)\alpha_{4},$ (14)

where $\psi=u_{\mathrm{av}}/u_{\mathrm{ref}}$ is the non-dimensional cross-
sectional average velocity and this can be calculated from the velocity
profile at station 1, for example, as

$\displaystyle\psi=B\left[\gamma\left(\alpha_{8}+1\right)+\beta_{5}\left(\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}-\gamma\right)\right].$
(15)

For the bypass flow, however, the difficulty is that the number of flow
passages at station 4 does not agree with that at station 1, since the
actuator disc creates the additional narrow bypass streamtube that has
$\beta_{4a}$ at station 4. To obtain a closed set of equations for this
periodic flow problem within the framework of quasi-one-dimensional modelling,
here we assume that the narrow bypass flow immediately outside of the wake is
‘merged’ or fully mixed with the main bypass flow regardless of the value of
$m$. This assumption seems reasonable since the difference between
$\beta_{4a}$ and $\beta_{4b}$ tends to be small as depicted in figure 2, and
eventually leads to

$\displaystyle\beta_{5}=m\psi+\left(1-m\right)\beta_{4m},$ (16)

where $\beta_{4m}$ is the ‘area-weighted’ average of $\beta_{4a}$ and
$\beta_{4b}$, i.e.,

$\displaystyle\beta_{4m}=\frac{\left(\gamma-\alpha_{2}\right)+\beta_{5}\left(\dfrac{1}{B}-\dfrac{\alpha_{2}}{\alpha_{4}}-\gamma\right)}{\dfrac{\left(\gamma-\alpha_{2}\right)}{\beta_{4a}}+\dfrac{\beta_{5}}{\beta_{4b}}\left(\dfrac{1}{B}-\dfrac{\alpha_{2}}{\alpha_{4}}-\gamma\right)}.$
(17)

It should be noted that $\alpha_{9}=m\psi+\left(1-m\right)\alpha_{8}$ is
required together with (14) to (17) to conserve the total mass in the viscous
zone, but this is automatically satisfied as we enforce
$\alpha_{9}=\alpha_{1}=1$ in this analysis due to the periodicity of the flow.

Finally, the thrust and power coefficients of the disc are expressed (for a
given cross-sectional average velocity of the entire flow,
$u_{\mathrm{av}}=\psi u_{\mathrm{ref}}$) as

$\displaystyle C_{T}=\frac{T}{\frac{1}{2}\rho\left(\psi
u_{\mathrm{ref}}\right)^{2}A}=\frac{\beta_{4a}^{2}-\alpha_{4}^{2}}{\psi^{2}},$
(18) $\displaystyle C_{P}=\frac{P}{\frac{1}{2}\rho\left(\psi
u_{\mathrm{ref}}\right)^{3}A}=\frac{\alpha_{2}\left(\beta_{4a}^{2}-\alpha_{4}^{2}\right)}{\psi^{3}}.$
(19)

For convenience, we also define the resistance coefficient (or local thrust
coefficient) of the disc as

$\displaystyle
K=\frac{T}{\frac{1}{2}\rho\left(\alpha_{2}u_{\mathrm{ref}}\right)^{2}A},$ (20)

from which and (7), we obtain

$\displaystyle\alpha_{2}=\left(\frac{\beta_{4a}^{2}-\alpha_{4}^{2}}{K}\right)^{1/2}.$
(21)

In summary, the above theoretical model for an infinite number of staggered
rows of identical ideal turbines consists of three non-dimensional input
parameters ($B$, $K$, $m$), ten non-dimensional unknowns to be determined
($\alpha_{2}$, $\alpha_{4}$, $\alpha_{5}$, $\alpha_{8}$, $\beta_{4a}$,
$\beta_{4b}$, $\beta_{4m}$, $\beta_{5}$, $\gamma$, $\psi$) and a set of ten
equations to be solved numerically: (4), (5), (10), (11), (13), (14), (15),
(16), (17) and (21).

###  2.2 Aligned rows of actuator discs

Compared to the staggered case, the theoretical model for the aligned case
becomes simpler as only two bypass-flow streamtubes (instead of three) need to
be considered in the inviscid zone. This is because the wake encounters the
disc immediately downstream (instead of two rows downstream). In the
following, we again consider the conservation of mass, momentum and energy in
the inviscid zone, and the simplified mixing process in the viscous zone, to
derive a similar but smaller set of equations for the aligned case.

By applying the Bernoulli equation to the two bypass-flow streamtubes and the
core-flow streamtube in the same manner as in the staggered case, we obtain
(1), (2) and (6), which lead to (4), (7) and (8) for $\beta_{4b}$, $T$ and
$P$, respectively. Note that these equations are identical for both aligned
and staggered cases, whereas (3) and (5) are only for the staggered case
(since the third bypass-flow streamtube does not exist in the aligned case).
Next, from the conservation of mass in the ‘main’ bypass-flow streamtube, we
obtain

$\displaystyle\beta_{5}u_{\mathrm{ref}}A\left(\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}\right)=\beta_{4b}u_{\mathrm{ref}}A\left[\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}-\frac{1}{\beta_{4a}}\left(\frac{\alpha_{2}}{\alpha_{4}}-\alpha_{2}\right)\right],$
(22)

which leads to an expression for $\alpha_{4}$ as

$\displaystyle\alpha_{4}=\frac{\alpha_{2}\left(1+\dfrac{1}{\beta_{4a}}-\dfrac{\beta_{5}}{\beta_{4b}}\right)}{\dfrac{1}{B}\left(1-\dfrac{\beta_{5}}{\beta_{4b}}\right)+\dfrac{\alpha_{2}}{\beta_{4a}}}.$
(23)

For $\beta_{4a}$ we consider the conservation of momentum for the entire flow
to obtain

$\begin{split}&(p_{1}-p_{4})\frac{A}{B}-T\\\ &=\rho
u_{\mathrm{ref}}^{2}A\left[\alpha_{2}(\alpha_{4}-1)+\left(\frac{\alpha_{2}}{\alpha_{4}}-\alpha_{2}\right)(\beta_{4a}-1)+\beta_{5}\left(\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}\right)\left(\beta_{4b}-\beta_{5}\right)\right],\end{split}$
(24)

which, together with (1) and (7), leads to

$\begin{split}\left(1-B\right)\beta_{4a}^{2}-2B\left(\frac{\alpha_{2}}{\alpha_{4}}-\alpha_{2}\right)\beta_{4a}&-2B\left(\alpha_{2}\alpha_{4}-\frac{1}{2}\alpha_{4}^{2}-\frac{\alpha_{2}}{\alpha_{4}}\right)\\\
&-2\beta_{5}\left(\beta_{4b}-\beta_{5}\right)\left(1-B\frac{\alpha_{2}}{\alpha_{4}}\right)-1=0.\end{split}$
(25)

For $\beta_{5}$ we need to consider the effect of mixing in the viscous zone.
By employing the same approach as in the staggered case, we obtain (16) for
the aligned case as well. Note, however, that $\psi$ and $\beta_{4m}$ are
different for the aligned case, which are

$\displaystyle\psi=B\left[\frac{\alpha_{2}}{\alpha_{4}}+\beta_{5}\left(\frac{1}{B}-\frac{\alpha_{2}}{\alpha_{4}}\right)\right],$
(26)

and

$\displaystyle\beta_{4m}=\frac{\left(\dfrac{\alpha_{2}}{\alpha_{4}}-\alpha_{2}\right)+\beta_{5}\left(\dfrac{1}{B}-\dfrac{\alpha_{2}}{\alpha_{4}}\right)}{\dfrac{1}{\beta_{4a}}\left(\dfrac{\alpha_{2}}{\alpha_{4}}-\alpha_{2}\right)+\dfrac{\beta_{5}}{\beta_{4b}}\left(\dfrac{1}{B}-\dfrac{\alpha_{2}}{\alpha_{4}}\right)}.$
(27)

We also need (14) together with the above equations to conserve the total mass
in the viscous zone, but this is automatically satisfied since here we enforce
$\alpha_{5}=\alpha_{1}=1$ due to the periodicity of the flow. Finally,
$C_{T}$, $C_{P}$ and $K$ are all defined as in the staggered case, yielding
the same equations (18) to (20) and thus (21) for $\alpha_{2}$.

In summary, the theoretical model for the aligned case consists of three non-
dimensional input parameters ($B$, $K$, $m$), seven non-dimensional unknowns
($\alpha_{2}$, $\alpha_{4}$, $\beta_{4a}$, $\beta_{4b}$, $\beta_{4m}$,
$\beta_{5}$, $\psi$) and a set of seven equations to be solved numerically:
(4), (16), (21), (23), (25), (26) and (27).

###  2.3 Example solutions

Here we present some example solutions of the above theoretical model, first
for the aligned case and then for the staggered case. As this is a three
parameter ($B$-$K$-$m$) problem, we start with fixing the disc resistance
coefficient, $K$, at its optimal value for the complete mixing case, i.e.,
$m=1$. When $m=1$, our problem (for both aligned and staggered cases) reduces
to the two-parameter ($B$-$K$) problem of Garrett & Cummins (2007) (hereafter
referred to as GC07); hence, this optimal value for $m=1$ is known to be
$K=2(1+B)^{3}/(1-B)^{2}$, which we refer to as $K_{\mathrm{GC07}}$.

Figure 3 shows how the performance of aligned rows of such ‘suboptimal’
actuator discs (with $K=K_{\mathrm{GC07}}$) changes with $m$, at two different
blockage ratios, $B=0.05$ and $0.2$. As can be expected intuitively, the power
coefficient $C_{P}$ of aligned discs decreases monotonically as the mixing
factor $m$ decreases (regardless of the blockage ratio). This agrees with the
common observation that the power of aligned rows of turbines tends to
decrease as we reduce the streamwise spacing between the rows, noting that the
mixing tends to be less complete (i.e., $m$ tends to be small) when the
spacing is small.

We also present in figure 3 the values of
$\alpha_{2}(\beta_{4a}^{2}-\alpha_{4}^{2})$ and $1/\psi^{3}$, to explain why
$C_{P}$ decreases as $m$ decreases. As can be seen from (19), $C_{P}$ is
equivalent to the product of these two values; the former is a different power
coefficient defined using the velocity upstream of the disc
($u_{\mathrm{ref}}$) instead of the cross-sectional average velocity
($u_{\mathrm{av}}$), and the latter represents the effect of the difference
between $u_{\mathrm{ref}}$ and $u_{\mathrm{av}}$. Here we can see that the
value of $\alpha_{2}(\beta_{4a}^{2}-\alpha_{4}^{2})$ actually increases as $m$
decreases. This is essentially because the local blockage effect is enhanced
(or the ‘effective’ blockage ratio increases) when the upstream core flow is
slower than the upstream bypass flow (Draper et al., 2016). This enhancement
of the local blockage effect caused by the incomplete wake mixing of the
upstream disc, however, is not strong enough to compensate for the ‘loss’ of
power possessed by the upstream core flow compared to the ‘original’ power
possessed by the cross-sectionally averaged flow (i.e., the increase rate of
$\alpha_{2}(\beta_{4a}^{2}-\alpha_{4}^{2})$ is smaller than the decrease rate
of $1/\psi^{3}$); therefore, $C_{P}$ eventually decreases as $m$ decreases.

Figure 3: Theoretical prediction of the effect of mixing on the performance of
aligned rows of actuator discs at two different blockage ratios: ($a$)
$B=0.05$ and ($b$) $B=0.2$. The disc resistance coefficient is fixed at the
optimal value for the case with complete mixing ($m=1$), i.e.,
$K=2(1+B)^{3}/(1-B)^{2}$.

Figure 4: Theoretical $C_{P}$ versus $K$ for the aligned case at ($a$)
$B=0.05$ and ($b$) $B=0.2$; and ($c$) the maximum $C_{P}$ obtained by varying
$K$ for a given set of $B$ and $m$.

Although the above analysis was for ‘suboptimal’ discs with
$K=K_{\mathrm{GC07}}$, the same relationship between $C_{P}$ and $m$ is
generally found for aligned discs with any given $K$ values. Figures 4($a$)
and ($b$) show $C_{P}$ versus $K$ curves for five different $m$ values, again
at $B=0.05$ and $0.2$, respectively, demonstrating the trend. It can also be
seen that the optimal $K$ value (to maximise $C_{P}$) tends to decrease as $m$
decreases, although the $C_{P}$ obtained at $K=K_{\mathrm{GC07}}$ is still
close to the maximum $C_{P}$ for a given $m$ (especially when the blockage
ratio is high). Figure 4($c$) shows how the maximum $C_{P}$, or
$C_{P\mathrm{max}}$, changes with $m$. As can be expected, for aligned rows of
actuator discs, $C_{P\mathrm{max}}$ also decreases monotonically as $m$
decreases.

Next, we focus on staggered rows of actuator discs. Similarly to the aligned
case, we start with the performance of ‘suboptimal’ discs with
$K=K_{\mathrm{GC07}}$, which is shown in figures 5($a$) and ($b$) for $B=0.05$
and $B=0.2$, respectively. At $B=0.05$, the value of $1/\psi^{3}$ slightly
increases (meaning that the upstream core flow velocity becomes slightly
higher than the cross-sectional average velocity) as $m$ decreases from 1 to
about 0.95, and then decreases as $m$ further decreases. In contrast, the
value of $\alpha_{2}(\beta_{4a}^{2}-\alpha_{4}^{2})$ first decreases slightly
and then increases as $m$ decreases from 1; this is again due to changes in
the ‘effective’ blockage ratio, i.e., the local blockage effect is enhanced
(or diminished) when the upstream core flow is surrounded by a faster (or
slower) bypass flow. However, the increase (or decrease) rate of
$\alpha_{2}(\beta_{4a}^{2}-\alpha_{4}^{2})$ tends to be smaller than the
decrease (or increase) rate of $1/\psi^{3}$, and hence $C_{P}$ follows the
trend of $1/\psi^{3}$, i.e., $C_{P}$ slightly increases first and then
decreases as $m$ decreases from 1. This initial increase in $1/\psi^{3}$ (and
thus in $C_{P}$) is more clearly seen at $B=0.2$, demonstrating that the
beneficial effect of the staggered layout (allowing the velocity upstream of
the disc to become higher than the cross-sectional average velocity) is
enhanced at a higher blockage ratio.

Figure 5($c$) compares $C_{P}$ versus $B$ curves for the complete mixing case
($m=1$) and for the ‘optimal mixing’ case ($m=m_{C_{P}}$, where $m_{C_{P}}$ is
the value of $m$ that maximises $C_{P}$), again for staggered rows of
‘suboptimal’ discs ($K=K_{\mathrm{GC07}}$). Also plotted together are
$m_{C_{P}}$ and $m_{\psi}$, the latter of which represents the value of $m$
that minimises $\psi$ (and thus maximises $1/\psi^{3}$). It can be seen how
the increase in $C_{P}$ due to the beneficial effect of the staggered layout
is enhanced, and how the optimal level of mixing decreases, as we increase the
blockage ratio. It is also worth noting that the difference between
$m_{C_{P}}$ and $m_{\psi}$ is small, especially at low blockage ratios. This
reflects the dominant influence of $1/\psi^{3}$ on $C_{P}$ as described
earlier in figures 5($a$) and ($b$).

The $C_{P}$ values presented above are for ‘suboptimal’ discs (with
$K=K_{\mathrm{GC07}}$) and can therefore be increased further by adjusting $K$
for a given $m$ (or for a given streamwise distance between the rows).
However, this additional increase in $C_{P}$ achieved by varying $K$ is rather
small for the staggered case. Figures 6($a$) and ($b$) show $C_{P}$ versus $K$
curves for five different $m$ values, again at $B=0.05$ and $0.2$,
respectively. When $B$ is small, $K_{\mathrm{GC07}}$ tends to be already close
to the optimal $K$ value for a given $m$; hence, the benefit of further
adjusting $K$ is small. When $B$ is large, $K_{\mathrm{GC07}}$ is not so close
to the optimal $K$ for a given $m$, but the $C_{P}$ versus $K$ curve tends to
have a flatter peak; hence again, the additional gain in $C_{P}$ by adjusting
$K$ is small. Figure 6($c$) shows the maximum $C_{P}$ (achieved by adjusting
$K$) versus $m$ curves for five different blockage ratios. The curves for
$B=0.05$ and $0.2$ are in fact very similar to the $C_{P}$ versus $m$ curves
for $K=K_{\mathrm{GC07}}$ plotted earlier in figures 5($a$) and ($b$). At
$B=0.2$, for example, the highest $C_{P}$ value achieved with
$K=K_{\mathrm{GC07}}$ is only 1.1$\%$ lower than that achieved by optimising
$K$.

Figure 5: Theoretical prediction of the performance of staggered rows of
actuator discs: the effect of mixing at ($a$) $B=0.05$ and ($b$) $B=0.2$; and
($c$) comparison of $C_{P}$ versus $B$ between the complete mixing case
($m=1$) and optimal mixing case ($m=m_{C_{P}}$). Plotted together are
$m_{C_{P}}$ and $m_{\psi}$, which are the values of $m$ to maximise $C_{P}$
and to minimise $\psi$, respectively. The disc resistance coefficient is fixed
at the optimal value for $m=1$, i.e., $K=2(1+B)^{3}/(1-B)^{2}$.

Figure 6: Theoretical $C_{P}$ versus $K$ for the staggered case at ($a$)
$B=0.05$ and ($b$) $B=0.2$; and ($c$) the maximum $C_{P}$ obtained by varying
$K$ for a given set of $B$ and $m$.

In summary, these example solutions of the simple three-parameter
($B$-$K$-$m$) theoretical model suggest the following: (i) For an infinitely
large staggered rows of tidal turbines, there is an optimal streamwise turbine
spacing (for a given cross-sectional blockage ratio or a given lateral turbine
spacing) to maximise the power of each turbine relative to the power of cross-
sectional average flow. This is an outcome of the combined effects of local
blockage and wake mixing, i.e., this optimum exists because the local flow
velocity upstream of each disc can become higher than the cross-sectional
average velocity (and this does help to increase the relative power despite
reducing the ‘effective’ blockage ratio) when the streamwise turbine spacing
is reasonably (not excessively) large for the wake mixing behind each disc to
be largely (but not entirely) completed. (ii) The optimal turbine resistance
to maximise this relative power also depends on both lateral and streamwise
turbine spacing, but the optimal value for a single row, namely
$K=2(1+B)^{3}/(1-B)^{2}$, is expected to yield close to the maximum relative
power. (iii) For an infinitely large aligned rows of turbines, however, this
relative power is maximised simply when the streamwise turbine spacing is
large enough for the wake mixing to be entirely completed; hence, the optimal
turbine resistance becomes identical to that for a single row.

Finally, it should be remembered that, in the real world, the array of
turbines cannot be infinitely large and that the cross-sectional average
velocity (which was considered as a fixed parameter in the above analysis)
will depend on how the flow resistance caused by the entire array alters the
natural tidal-channel-scale momentum balance (see, e.g., Vennell (2012);
Vennell et al. (2015); Gupta & Young (2017)). Therefore, the true optimal
turbine resistance (to maximise the power of turbines in a given tidal
channel) will differ from that discussed in the above analysis.

##  3 Large-eddy simulation set-up

The theoretical model presented above has been developed on the assumption
that the flow around each turbine is inviscid and steady. However, real
turbine wakes are unsteady with coherent flow structures developed at a rotor
level, e.g. tip-vortices, and further downstream, e.g. wake meandering. Thus,
we now investigate how the theoretical predictions of the performance of
infinitely large tidal arrays compare to high-fidelity numerical predictions
using Large-Eddy Simulation (LES). We adopt the well-validated in-house LES
code Digital Offshore Farms Simulator (DOFAS) in which turbine blades are
represented using an Actuator Line Method (ALM) (Ouro et al., 2019c) and the
flow solver is fully parallelised using the Message Passing Interface (MPI),
providing a great computational scalability and performance (Ouro et al.,
2019a). Details of the flow solver are provided in the Appendix.

We intentionally set the flow conditions to be as similar as possible to the
conditions considered in the theoretical analysis, whilst ensuring that the
physical dimensions of the turbines are close to those found in real tidal
stream turbines. The 1MW DEEPGen IV tidal stream turbine design from the
ReDAPT project is adopted, and the details of the blade hydrodynamic data used
for the ALM are available in Scarlett et al. (2019). The turbines have a
diameter ($D$) of 12 m, rotating at a constant speed that corresponds to a
tip-speed ratio of 4.0 (which is known to be optimal for the case of single-
turbine operation), and include 10 m long ($0.8D$) nacelles. For convenience,
we set the cross-sectional average velocity $U_{0}$ to 2.0 m s-1 and consider
this as our reference velocity, which yields a rotational speed of
$\Omega=1.35$ rad s-1 and a full-revolution period of $T$ = 4.724 s.

The computational domain is 432 m long ($L$), 144 m wide ($W$) and 24 m high
($H$), equivalent to $36D\times 12D\times 2D$, which is close to $6\pi H\times
2\pi H\times H$ commonly used in turbulent channel-flow simulations. Turbines
are vertically centred at mid-water depth, i.e. $z=H/2$, to reduce vertical
asymmetry effects that may complicate the comparison between LMADT and LES. A
fixed time step $\Delta t$ is set to 0.045 s together with a uniform spatial
resolution $\Delta x$ equal to 0.25 m, yielding a total of 48 mesh elements
across the rotor diameter which is similar to the resolution adopted in other
LES-ALM studies, e.g. Churchfield et al. (2013); Yang & Sotiropoulos (2019);
Foti et al. (2019). Hence, the domain is divided into 1692 $\times$ 576
$\times$ 96 grid cells over the three spatial directions with a total of about
93.5 $\times$ 106 elements. Simulations are run using 864 processors on three
High-Performance Computing facilities, namely Supercomputing Wales, GW4
Isambard, and ARCHER.

In the present LES, we represent infinite turbine arrays by imposing periodic
boundary conditions in the streamwise and transverse directions. The flow is
driven by the pressure gradient term $\Pi_{i}$ which enforces the mass flux
across the entire domain to be constant, providing the cross-sectional
averaged velocity is equal to $U_{0}$. Fixing the bulk velocity enables a
direct comparison between the LES and the theoretical analysis in §2, unlike
fixing $\Pi_{i}$ as in other infinitely-large wind farm simulations (Calaf et
al., 2010; Sharma et al., 2018). Note that, to simulate a real large tidal
array, both of the bulk velocity and $\Pi_{i}$ would need to vary depending on
the macro-scale momentum balance over the entire tidal channel, which is
outside the scope of the present study. Shear-free conditions are adopted at
the top boundary to represent a free surface, whilst the bottom shear stress
is calculated using wall functions for a hydrodynamically smooth wall,
similarly to Kang et al. (2014). A representation of the computational domain
is presented in figure 7 together with instantaneous flow structures
visualised using Q-criterion iso-surfaces for one of the configurations
studied.

Figure 7: Dimensions of the computational domain and instantaneous flow
structures generated in the ST-9x6 case, visualised using iso-surfaces of the
Q-criterion.

An initial precursor simulation was performed for over 60 eddy turnover time
units ($t_{e}$), corresponding to nearly 27 flow-through ($FT=L/U_{0}$), in
order to generate fully-developed turbulent flow conditions to be used as
initial flow field for each array simulation. We perform 28 infinitely-large
array simulations, whose details are presented in table 1 providing values of
normalised streamwise and transverse separation between turbines ($S_{x}/D$,
$S_{y}/D$), local blockage ratio $B$ (relating the turbine’s projected area to
the open-channel cross-section and number of turbines per row, i.e.
$B=n_{y}\pi D^{2}/4HW$), physical simulated time in terms of flow-through
time, and the number of computed revolutions. Also presented in this table are
the key results of each simulation, namely the thrust coefficient ($C_{T}$),
power coefficient ($C_{P}$), their relative difference ($\Delta C_{T}$,
$\Delta C_{P}$) from the AL-36x12 configuration as a reference case, and the
resistance coefficient ($K$). The arrays adopted in the LES are perfectly
aligned and staggered layouts labelled as AL and ST followed by the device
spacing $S_{x}$x$S_{y}$, e.g. case ST-9x6 corresponds to a staggered array
with $S_{x}/D$ = 9 and $S_{y}/D$ = 6 as in figure 7.

Case | $S_{x}/D$ | $S_{y}/D$ | $B$ | Time ($FT$) | Revs | $C_{T}$ | $C_{P}$ | $\Delta C_{T}$ | $\Delta C_{P}$ | $K$
---|---|---|---|---|---|---|---|---|---|---
AL-36x12 | 36 | 12 | 0.033 | 27.1 | 1257 | 0.767 | 0.287 | 1.000 | 1.000 | 2.55
AL-18x12 | 18 | 12 | 0.033 | 27.1 | 1257 | 0.766 | 0.277 | 0.999 | 0.965 | 2.74
AL-9x12 | 9 | 12 | 0.033 | 27.1 | 1257 | 0.751 | 0.264 | 0.979 | 0.922 | 2.80
AL-6x12 | 6 | 12 | 0.033 | 20.6 | 957 | 0.737 | 0.252 | 0.961 | 0.879 | 2.82
ST-18x12 | 18 | 12 | 0.033 | 25.0 | 1160 | 0.777 | 0.292 | 1.014 | 1.021 | 2.29
ST-9x12 | 9 | 12 | 0.033 | 25.0 | 1160 | 0.774 | 0.291 | 1.009 | 1.016 | 2.18
ST-6x12 | 6 | 12 | 0.033 | 25.0 | 1160 | 0.774 | 0.292 | 1.009 | 1.019 | 2.07
AL-36x6 | 36 | 6 | 0.065 | 24.6 | 1141 | 0.786 | 0.301 | 1.025 | 1.050 | 2.20
AL-18x6 | 18 | 6 | 0.065 | 26.6 | 1233 | 0.785 | 0.302 | 1.024 | 1.054 | 2.11
AL-9x6 | 9 | 6 | 0.065 | 20.9 | 972 | 0.772 | 0.293 | 1.006 | 1.021 | 2.12
AL-6x6 | 6 | 6 | 0.065 | 20.8 | 967 | 0.764 | 0.285 | 0.996 | 0.996 | 2.02
ST-18x6 | 18 | 6 | 0.065 | 24.8 | 1151 | 0.784 | 0.300 | 1.022 | 1.045 | 2.28
ST-9x6 | 9 | 6 | 0.065 | 23.4 | 1088 | 0.783 | 0.299 | 1.020 | 1.045 | 2.25
ST-6x6 | 6 | 6 | 0.065 | 25.0 | 1160 | 0.777 | 0.293 | 1.013 | 1.024 | 2.27
AL-36x4 | 36 | 4 | 0.098 | 24.7 | 1146 | 0.787 | 0.307 | 1.026 | 1.073 | 1.92
AL-18x4 | 18 | 4 | 0.098 | 25.0 | 1160 | 0.787 | 0.310 | 1.026 | 1.082 | 1.89
AL-9x4 | 9 | 4 | 0.098 | 23.3 | 1083 | 0.777 | 0.300 | 1.013 | 1.046 | 1.88
AL-6x4 | 6 | 4 | 0.098 | 23.4 | 1088 | 0.765 | 0.285 | 0.997 | 0.995 | 2.08
ST-18x4 | 18 | 4 | 0.098 | 23.1 | 1074 | 0.787 | 0.310 | 1.027 | 1.082 | 1.93
ST-9x4 | 9 | 4 | 0.098 | 23.4 | 1088 | 0.781 | 0.314 | 1.018 | 1.094 | 1.82
ST-6x4 | 6 | 4 | 0.098 | 19.9 | 924 | 0.780 | 0.312 | 1.017 | 1.089 | 1.78
AL-36x3 | 36 | 3 | 0.131 | 27.1 | 1257 | 0.796 | 0.316 | 1.038 | 1.102 | 2.09
AL-18x3 | 18 | 3 | 0.131 | 22.7 | 1054 | 0.796 | 0.316 | 1.038 | 1.104 | 1.90
AL-9x3 | 9 | 3 | 0.131 | 23.8 | 1102 | 0.780 | 0.300 | 1.017 | 1.047 | 1.99
AL-6x3 | 6 | 3 | 0.131 | 22.2 | 1030 | 0.767 | 0.285 | 1.000 | 0.996 | 2.19
ST-18x3 | 18 | 3 | 0.131 | 26.9 | 1246 | 0.785 | 0.319 | 1.024 | 1.112 | 1.58
ST-9x3 | 9 | 3 | 0.131 | 22.0 | 1020 | 0.787 | 0.320 | 1.027 | 1.116 | 1.64
ST-6x3 | 6 | 3 | 0.131 | 23.3 | 1083 | 0.788 | 0.320 | 1.028 | 1.116 | 1.64
Table 1: Details of the LES cases with streamwise and spanwise spacing
($S_{x}/D$, $S_{y}/D$), local blockage ratio ($B$), simulated physical time in
terms of flow-through times, number of revolutions, hydrodynamic coefficients
($C_{T}$ and $C_{P}$), their variation compared to the reference case
AL-36x12, and resistance coefficient ($K$).

##  4 Large-eddy simulation results

First, we focus on the array layouts with $S_{x}/D$ = 9 in order to elucidate
how the flow field varies with $S_{y}/D$. The time-averaged and instantaneous
flow fields are presented in §4.1, followed by some turbulence statistics in
§4.2. We then present the wake centre-line velocities for all 28 cases in §4.3
to analyse the effects of both $S_{x}/D$ and $S_{y}/D$ on the wake recovery,
and finally in §4.4 the variations of the power and thrust coefficients are
presented and compared with the theoretical predictions. Hereafter, the time-
averaged value of any variable is denoted as $\langle\cdot\rangle$ and the
instantaneous fluctuation value is represented as $(\cdot)^{\prime}$ obtained
from the Reynolds decomposition.

###  4.1 Streamwise velocity field

We present in figure 8 contours of the normalised time-averaged streamwise
velocity ($\langle u\rangle/U_{0}$) comparing aligned and staggered layouts
with the same streamwise spacing of $S_{x}/D=9$ and lateral spacing of
$S_{y}/D$ = 6, 4 and 3. For aligned arrays, reducing the lateral separation
between devices in the same row leads to an increased local blockage that
induces larger flow acceleration in the bypass flow between them, which is
well observed for AL-9x3 with the bypass-flow velocity being approximately 20%
higher than the bulk velocity. In perfectly staggered arrays, the wake
generated behind each turbine is mostly recovered when reaching the following
row at a downstream distance of $S_{x}$. Then, due to the lateral blockage
caused by the turbines located on both sides, the mostly-recovered wake
accelerates further and then impinges the turbine located in the next row,
i.e. at 2$S_{x}$ downstream of the turbine that originally generated the wake.
For layouts with low lateral blockage, i.e. $S_{y}/D$ = 6 and 4, there is a
wider lateral spread of the wakes, which is well-observed in figure 8 for the
aligned cases. However, in staggered configurations this lateral wake
expansion appears limited compared to the corresponding aligned cases.

Figure 8: Contours of time-averaged streamwise velocity $\langle
u\rangle/U_{0}$ at hub height for aligned and staggered cases with $S_{x}/D$ =
9 and lateral spacing $S_{y}/D$ = 6 (top), 4 (mid) and 3 (bottom).

The different spreading rates of the time-averaged turbine wakes are partly
explained by their instantaneous behaviour. In figure 9 we present contours of
instantaneous streamwise velocities ($u/U_{0}$) for the previous cases shown
in figure 8 with $S_{x}/D$ = 9 and changing $S_{y}$, which reveals the
pronounced meandering nature of the wakes with large oscillation amplitudes
(Foti et al., 2019). In aligned configurations, the meandering of the wake is
affected by the lateral distance between turbines in the same row, i.e.,
increasing $S_{y}$ leads to a larger meandering amplitude. This is presumably
because a higher lateral blockage can constrain the spanwise wake motion. It
can also be seen that there is a larger wake meandering amplitude in aligned
arrays compared to their staggered counterparts. The lateral wake motion due
to wake meandering seems almost negligible for ST-9x3, explaining the narrow
time-averaged wakes shown earlier in figure 8.

Figure 9: Contours of instantaneous streamwise velocity $u/U_{0}$ at hub
height for aligned and staggered cases with $S_{x}/D$ = 9 and lateral spacing
$S_{y}/D$ = 6 (top), 4 (mid) and 3 (bottom). Same legend as 8.

###  4.2 Turbulence statistics

We further analyse the flow field in figures 10 and 11 with contours of
turbulence intensity in the streamwise ($\sigma_{u}/U_{0}$, with
$\sigma_{u}=\langle u^{\prime}u^{\prime}\rangle^{0.5}$) and transverse
directions ($\sigma_{v}/U_{0}$, with $\sigma_{v}=\langle
v^{\prime}v^{\prime}\rangle^{0.5}$), again for the cases with $S_{x}/D$ = 9.
These second-order statistics indicate that having turbines in an aligned
layout leads to a notably stronger flow unsteadiness both inside and outside
the wake of each turbine compared to the staggered cases. This agrees with the
earlier observation that the wake meandering is stronger in the aligned cases.
In the near-wake region, in which tip vortices are still coherent and maintain
a shear layer between the core flow and the bypass flow (Ouro et al., 2019c),
both streamwise and spanwise turbulence intensities are substantially higher
in the aligned cases than in the staggered cases. Low turbulence intensity
regions in the staggered cases are evident in the flow bypassing each turbine,
coinciding with the regions in which the wake of an upstream turbine is almost
fully recovered as shown earlier in figures 8 and 9. The results also show
that increasing the blockage ratio $B$ (or reducing the distance $S_{y}$
between turbines in the same row) reduces the turbulence intensity
irrespective of the overall turbine arrangement, as a consequence of
constraining the meandering motion.

Figure 10: Contours of streamwise turbulence intensity ($\sigma_{u}/U_{0}$) at
hub height for aligned and staggered cases with $S_{x}/D$ = 9 and lateral
spacing $S_{y}/D$ = 6 (top), 4 (mid) and 3 (bottom).

Figure 11: Contours of spanwise turbulence intensity ($\sigma_{v}/U_{0}$) at
hub height for aligned and staggered cases with $S_{x}/D$ = 9 and lateral
spacing $S_{y}/D$ = 6 (top), 4 (mid) and 3 (bottom).

As can be expected from the above results, the wake meandering also notably
affects the distribution of turbulent momentum exchange, which is described in
figure 12 with contours of the Reynolds shear stress $-\langle
u^{\prime}v^{\prime}\rangle/U_{0}^{2}$ at hub height. For each array layout,
narrowing the lateral spacing reduces the level of Reynolds shear stress,
indicating that there is a lower level of momentum exchange between the wakes
and the surrounding bypass flows. Staggered layouts consistently attain a
lower shear stress level compared to the aligned layouts irrespective of the
number of turbines deployed.

Figure 12: Contours of Reynolds shear stress ($-\langle
u^{\prime}v^{\prime}\rangle/U_{0}^{2}$) obtained at hub height for aligned and
staggered cases with $S_{x}/D$ = 9 and lateral spacing $S_{y}/D$ = 6 (top), 4
(mid) and 3 (bottom).

###  4.3 Centre-line velocities

Next, we quantitatively compare the mean streamwise velocity ($\langle
u\rangle/U_{0}$) distribution along the wake centre-line in figure 13 for
every configuration simulated (see table 1). Note that for $S_{x}/D$ = 18, 9
and 6 we plot the velocity distribution over two rows of turbines, including a
shaded area representing the location at which turbines from the second row
are located. When the streamwise spacing is relatively small ($S_{x}/D$ = 9
and 6), the values of $\langle u\rangle/U_{0}$ behind the turbine nacelles
increase faster for the lower lateral blockage cases, i.e. $S_{y}/D$ = 6 and
12. This faster wake recovery in lower blockage cases results from a larger
entrainment of ambient flow into the wake enhanced by the stronger meandering
behaviour, as seen earlier in the distribution of Reynolds shear stresses in
figure 12.

Although lower blockage ratios (or larger $S_{y}/D$ values) tend to result in
higher wake recovery rates in the near-wake region, this trend does not hold
in the far-wake region. In the aligned cases with $S_{x}/D$ = 18, the values
of $\langle u\rangle/U_{0}$ are slightly over unity shortly upstream of the
turbines for configurations with $S_{y}/D\leq$ 6, whereas for $S_{y}/D$ = 12
the wake velocity is not fully recovered to the bulk velocity value despite
the high wake recovery rate in the near-wake region. Further decreasing
$S_{x}/D$ to 9 in aligned cases reduces the amount of wake recovery for every
$S_{y}/D$, resulting in $\langle u\rangle/U_{0}$ of about 0.9 to 0.95 at one
rotor-diameter upstream of the turbines ($x/D$ = 8 and 17). It is clearly seen
that the wakes in higher blockage cases (AL-9x4 and AL-9x3) have a lower
recovery rate in the near-wake region but a higher recovery rate in the far-
wake region. The same trend can be observed in the smallest streamwise spacing
cases ($S_{x}/D$ = 6), where the streamwise velocities upstream of the
turbines are even lower due to the lack of space for the wake to recover.
Overall, for aligned arrays we may conclude that $S_{x}/D$ is the main design
parameter that determines how much the wake recovers before approaching the
turbine downstream, with $S_{y}/D$ contributing to a lower extent by changing
the turbulent wake characteristics.

Figure 13: Distribution of time-averaged streamwise velocities ($\langle
u\rangle/U_{0}$) along the wake centre-line at the hub height for aligned
(left) and staggered (right) configurations with $S_{x}/D$ = 36 (top), 18
(mid-top), 9 (mid-bottom) and 6 (bottom), and $S_{y}/D$ = 12 (black), 6
(blue), 4 (green) and 3 (orange).

In staggered arrays with $S_{x}/D$ = 18, the centre-line mean velocity notably
increases at around $x/D$ = 18, which coincides with the location of the
second row shifted laterally with respect to the first row of turbines. This
flow acceleration is due to the local blockage caused by the turbines in the
second row, and thus is more noticeable when the lateral spacing $S_{y}/D$ is
small. For example, for ST-18x3 and ST-18x4 the maximum velocity is about
1.1$U_{0}$, whereas for ST-18x6 or ST-18x12 it is about 1.02$U_{0}$. This
agrees qualitatively with the theoretical analysis presented earlier in §2.3,
i.e., the local flow velocity upstream of each turbine may exceed the cross-
sectional average velocity in staggered arrays. Reducing the streamwise
spacing $S_{x}/D$ makes this additional local flow acceleration less
noticeable, but the rate at which the wake velocity recovers still varies with
$S_{y}/D$, being higher for the higher blockage cases with $S_{y}/D$ = 4 and
3. Another distinct feature of time-averaged wakes in staggered arrays is that
after passing through the streamwise location of the second row (i.e., $x\geq
S_{x}$) the centre-line mean velocity remains fairly constant until about two
rotor-diameters upstream of the next turbine.

Finally, comparing aligned and staggered cases for a given streamwise spacing,
we can confirm that the velocity recovery rate in the near-wake region is
higher in the former configurations than in the latter. This is due to
stronger wake mixing enhanced by the wake meandering as shown earlier.

###  4.4 Array efficiency

As the impact of array layout and turbine spacing on the flow field has been
confirmed, we now focus on the thrust and power coefficients of turbines to
analyse their efficiency in the simulated infinitely large tidal arrays. For
each configuration, time-averaged hydrodynamic coefficients are further
averaged for all turbines comprising the array, and the results are presented
in figure 14 with error bars indicating the standard deviation of array-
averaged values. For aligned layouts with a given $S_{y}/D$, both thrust and
power coefficients tend to have their maximum values at $S_{x}/D$ = 36, but
the results at $S_{x}/D$ = 18 are similar as this streamwise spacing is still
large enough for the wakes to almost fully recover. However, when $S_{x}/D$ is
further reduced, both $C_{T}$ and $C_{P}$ drop considerably in line with the
decrease in the mean streamwise velocity upstream of each turbine, as shown
earlier in figure 13. In contrast, the performance of turbines in staggered
configurations appears less sensitive to $S_{x}/D$.

Figure 14: Results of $C_{T}$ (top) and $C_{P}$ (bottom) obtained from the LES
for the aligned ($\square$) and staggered ($\bigcirc$) layouts.

Comparing arrays with large lateral spacing ($S_{y}/D$ = 12), staggered
configurations consistently provide higher $C_{T}$ and $C_{P}$ values than the
aligned counterparts, with larger differences for shorter streamwise spacing.
For cases with $S_{y}/D$ = 6 and 4, staggered configurations still tend to
provide higher $C_{T}$ and $C_{P}$ than the aligned counterparts, although the
differences are negligibly small at $S_{x}/D$ = 18. Further reducing the
lateral spacing to $S_{y}/D$ = 3 leads to higher $C_{T}$ in the aligned case
than in the staggered case at $S_{x}/D$ = 18, but again $C_{P}$ values are
consistently higher in the staggered cases. These results suggest that, as
expected, staggered configurations are in general more efficient than the
aligned counterparts.

Next, we analyse the temporal fluctuations of thrust and power from their
root-mean-square (rms) values. The results averaged over all turbines
comprising the array are presented in figure 15, which shows that decreasing
$S_{x}/D$ leads to an increase in the fluctuations of both thrust and power,
as expected from enhanced unsteadiness of the approaching flow. Adopting a
large lateral spacing also increases the load fluctuations in comparison to
cases with the same streamwise spacing and a smaller lateral spacing. This is
in line with the strength of wake oscillations increasing with $S_{y}/D$, as
shown earlier in figure 9. It is also worth noting that turbines in aligned
configurations tend to have larger load variations than the staggered
counterparts, especially when decreasing $S_{x}/D$. These results suggest that
staggered configurations not only enhance the array’s efficiency (i.e., mean
power coefficient) but can also reduce the load fluctuations and thus the
long-term fatigue damage of turbine rotors.

Figure 15: Results of root-mean-squared temporal fluctuations of $C_{T}$
(left) and $C_{P}$ (right) obtained from the LES of the considered array
layouts.

##  5 Discussion

The layout of turbines in a large tidal array has important implications to
its power generation capability. This is because the flow through a large
array is characterised by complex wake-turbine interactions involving the
effects of local blockage and wake mixing, both of which are functions of the
layout of turbines. In this study we aimed to understand the impact of turbine
layout for infinitely large tidal arrays, using a new theoretical model based
on the LMADT and high-fidelity LES-ALM, focusing on how the array efficiency
changes depending on the turbine resistance, local blockage ratio within each
row of turbines ($B$ in the theoretical analysis, which is a function of
$S_{y}$ in the LES, i.e. $B\propto S_{y}$) and the completeness of wake mixing
between each row ($m$ in the theoretical analysis, which is a function of
$S_{x}$ in the LES, i.e. $m\propto S_{x}$). It should be borne in mind that
the theoretical model allows to predict overall trends and upper bound
estimates of the array efficiency with almost negligible computational cost,
whereas the LES-ALM allows to resolve the details of complex turbulent flow
field within the array but at a high computational expense. Thus, these are
two highly contrasting approaches to this fluid flow problem.

Before comparing and further discussing the theoretical and LES-ALM results,
we emphasise two key differences between these two approaches: (i) In the LES-
ALM the operating point of turbines (i.e., rotational speed) is kept constant
for all cases, resulting in slightly different $K$ values as presented in
table 1. No further tuning of the rotational speed to obtain a constant $K$
(to match the theoretical analysis) is performed due to the high computational
cost that would be required for it. Alternatively, a fairer comparison with
the theoretical analysis could be made using LES with an Actuator Disc Model
(ADM), but we adopted an ALM in this study as it allows a more realistic
representation of the turbulent flow field within the array. Thus,
consideration needs to be taken when comparing the theoretical results for a
constant $K$ value with scattered LES results with slightly different $K$
values. (ii) The current theoretical model does not provide an explicit
relationship between the mixing factor $m$ and the array configuration,
whereas the LES-ALM automatically predicts the wake recovery rate for a given
configuration by resolving the turbulent flow field. To make a direct
comparison, here the mixing factor in the theoretical analysis is assumed to
be $m=1-(S_{x}/D)^{-1}$, which is arguably the simplest model to relate $m$ to
$S_{x}/D$ without knowing any detailed characteristics of turbine wake mixing
for a given array configuration a priori.

We now compare array efficiency predictions between the theoretical analysis
(assuming $K$ = 2) and LES-ALM in figure 16. Since the theoretical $C_{P}$
values are for ideal turbines (or actuator discs) and they are not directly
comparable to $C_{P}$ values for real rotors, here we normalise $C_{P}$ with a
reference power coefficient $C_{P_{ref}}$, which is defined for the
theoretical analysis and LES-ALM separately. For the LES-ALM, $C_{P_{ref}}$ is
the power coefficient obtained for the sparsest array (AL-36x12), in which
turbine-to-turbine interactions are deemed negligible, whereas for the
theoretical analysis, this is the power coefficient for $B=0.033$, $K=2$ and
$m=1$. Hence, $C_{P}/C_{P_{ref}}$ plotted in figure 16 represents the change
rate of $C_{P}$ due to the effect of different turbine layouts. Overall, there
is a qualitative agreement in $C_{P}/C_{P_{ref}}$ between the two approaches.
In aligned arrays, the efficiency monotonically decreases with $S_{x}/D$
(except when $1-(S_{x}/D)^{-1}$ is close to unity and $S_{y}/D\leq$ 6) as
turbines increasingly operate in the wake of upstream turbines. The rate of
decrease in $C_{P}/C_{P_{ref}}$ is different between the theoretical and LES
results, largely due to the simple relationship between $m$ and $S_{x}/D$
assumed to make this comparison. It is therefore expected that the agreement
would improve if the relationship between $m$ and $S_{x}/D$ is modelled
appropriately in future work. For staggered arrays, again the theoretical
predictions agree qualitatively with the LES, showing that $C_{P}/C_{P_{ref}}$
is insensitive to the streamwise spacing at least within the range of
conditions considered here. However, the effect of lateral spacing $S_{y}/D$
is slightly over-predicted by the theoretical model compared to the LES.

Figure 16: Comparison of normalised power coefficient ($C_{P}$/$C_{P_{ref}}$)
predicted with the theory (assuming $K=2$ and $m=1-(S_{x}/D)^{-1}$) and
computed from the LES.

The above comparison of $C_{P}/C_{P_{ref}}$ suggests that the new theoretical
tidal array model is promising as it seems to capture the combined effects of
local blockage and wake mixing qualitatively correctly, despite not accounting
for some key transient flow phenomena, such as blade tip vortices and wake
meandering, which are well-captured in the LES-ALM. As the efficiency of
turbines in a large array is determined mainly by the time-averaged flow field
within the array, further improvements of the theoretical model could be made
in future studies using LES. In particular, further results of LES-ALM for a
wider range of parameters would allow us to empirically model the mixing
factor $m$ as a function of both $S_{x}/D$ and $S_{y}/D$, where it would be
important to account for the dependency of wake meandering and other transient
flow phenomena (that collectively determine the wake recovery rate) on the
layout of turbines. It should be remembered, however, that the efficiency of
real tidal arrays would depend not only on micro-scale flow interactions
within the array, namely the turbine-to-turbine interactions studied here, but
also on macro-scale flow interactions outside the array (Vennell, 2012;
Vennell et al., 2015; Gupta & Young, 2017).

##  6 Conclusions

This paper has investigated the performance of tidal stream turbines in an
infinitely large array, using two different approaches: a quasi-one-
dimensional theoretical model based on the Linear Momentum Actuator Disc
Theory (LMADT), and Large-Eddy Simulation with an Actuator Line Method (LES-
ALM). Two different types of turbine layouts were considered in this study,
i.e., perfectly aligned and staggered layouts. For the LES-ALM, 28 different
arrays with various streamwise ($6\leq S_{x}/D\leq 36$) and lateral ($3\leq
S_{y}/D\leq 12$) turbine spacing were considered. For the theoretical
analysis, a hybrid inviscid-viscous approach was employed to model an
infinitely large array using only three input parameters, namely the local
blockage ratio $B$, disc resistance coefficient $K$ and wake mixing factor
$m$.

Our LES results have shown that the lateral spacing has a pronounced effect on
the characteristics of wake meandering. In particular, the amplitude of wake
meandering is found to decrease as $S_{y}/D$ is reduced. The main consequence
of this change in wake dynamics is that a lower amplitude of wake meandering
leads to less entrainment of momentum from the surrounding bypass flow into
the wake, and thus a lower wake recovery rate. However, this negative effect
of small lateral spacing on the wake recovery rate is observed mainly in the
near-wake region only, and the completeness of wake recovery (i.e., how much
the wake velocity is recovered before the wake approaches the next turbine)
tends to depend more on the streamwise spacing, especially for aligned arrays.
We have also confirmed from our LES results that, in staggered arrays, the
wake experiences an additional acceleration when it passes through the
laterally-shifted row of turbines immediately downstream. This additional
acceleration is due to the effect of local blockage, which is enhanced when
the lateral spacing is small. When $S_{y}/D$ is sufficiently small, the
centre-line wake velocity is found to even exceed the bulk velocity, resulting
in a high power of turbines for a fixed bulk velocity. Resolving the turbulent
flow field with LES has also allowed to study the temporal fluctuations of
turbine loads in a large tidal array, showing that these are approximately
twice larger in aligned arrays than in staggered arrays for a given turbine
spacing.

Whilst the LES results have revealed the complexity of turbulent flow
phenomena that collectively determine the performance of turbines in a large
tidal array, the simple theoretical model has captured the basic trend of
turbine performance in both aligned and staggered arrays qualitatively
correctly. In particular, the theoretical model suggests that there is an
optimal streamwise spacing to maximise the performance of turbines (or the
power of turbines for a fixed bulk velocity) in a large staggered array. This
optimum exists as the local flow velocity upstream of each turbine can exceed
the bulk velocity only when the streamwise spacing is reasonably (not
excessively) large for the mixing between locally faster and slower flows to
be largely (but not entirely) completed within that streamwise distance.
However, both theoretical and LES results show that, at least within the range
of conditions tested, the effect of streamwise spacing is less than that of
lateral spacing, i.e., the performance of turbines in staggered arrays depends
more on $S_{y}/D$ than on $S_{x}/D$. We have also observed some quantitative
differences between the theoretical predictions and the LES. One of the main
causes of differences is that, to make a comparison between the two
approaches, we have assumed a simple relationship between the mixing factor
$m$ (representing the completeness of mixing after each row of turbines) and
the streamwise spacing between rows. Further LES results for a wider range of
parameters would be helpful in future studies to develop an empirical model of
$m$ as a function of both streamwise and lateral spacing.

The results obtained in this study will help understand and improve the
performance of tidal turbines in future large tidal arrays. It should be borne
in mind, however, that the performance of real tidal arrays may depend not
only on micro-scale (turbine-to-turbine) flow interactions within the array
but also on macro-scale flow interactions between the array and tidal-channel
flow, the latter of which was outside the scope of this study.

## Declaration of interests

The authors report no conflict of interest

## Acknowledgements

This research has been partially funded by the UK’s Engineering and Physical
Sciences Research Council (EPSRC) (grant number EP/R51150X/1). The first
author would like to acknowledge the support of the Supercomputing Wales
project, which is partially funded by the European Regional Development Fund
(ERDF) via the Welsh Government, the Isambard project funded by the EPSRC
(EP/P020224/1), the GW4 alliance, the Met Office, Cray and Arm. This work also
used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk).
The second author would like to thank Dr Scott Draper for useful discussions
on LMADT.

## Appendix A: Large-eddy simulation code DOFAS

We used the Digital Offshore Farms Simulator (DOFAS) (Ouro et al., 2019c), an
in-house LES code fully parallelised with Message Passing Interface (MPI) that
also features a hybrid MPI/OpenMP scheme to maximise its computational
performance (Ouro et al., 2019a). In DOFAS, the spatial domain is divided into
rectangular sub-domains and discretised using Cartesian grids with staggered
storage of velocities, i.e. velocity components are computed at the cell faces
whilst pressure and scalar values are calculated at the cell centres. This
scheme allows an even subdivision of the computational region into sub-domains
to effectively perform the simulations. The governing equations resolved in
DOFAS are the incompressible spatially-filtered Navier-Stokes equations:

$\displaystyle\frac{\partial u_{i}}{\partial x_{i}}=0,$ (28)
$\displaystyle\frac{\partial u_{i}}{\partial t}+u_{j}\frac{\partial
u_{i}}{\partial x_{j}}=-\frac{1}{\rho}\frac{\partial p}{\partial
x_{i}}+(\nu+\nu_{t})\frac{\partial^{2}u_{i}}{\partial
x_{j}^{2}}+f_{t}+\Pi_{i},$ (29)

where $u_{i}$ = $(u,v,w)^{T}$ is the vector of spatially-filtered velocities,
the coordinates vector is $x_{i}$ = $(x,y,z)^{T}$, $\rho$ denotes the fluid
density, $p$ is the relative pressure, $\nu$ is the kinematic viscosity of the
fluid, and $f_{t}$ is a source term resulting from the actuator line method
and immersed boundary forcing used for the representation of turbine rotors
and nacelles, respectively. $\Pi_{i}$ is a source term representing the
driving pressure gradient responsible for keeping a constant flow rate when
periodic boundary conditions are used in the streamwise direction. The eddy-
viscosity $\nu_{t}$ is calculated using the Wall-Adapting Local Eddy-viscosity
(WALE) sub-grid scale model from Nicoud & Ducros (1999).

The velocity field is spatially discretised using a fourth-order central
differences scheme. Simulations are advanced in time using a fractional step
method, with a three-step low-storage Runge-Kutta scheme to obtain the non-
solenoidal velocity field by explicitly computing the convection and diffusion
terms, which is then corrected after the Poisson pressure equation is solved
using an efficient multigrid solver (Cevheri et al., 2016).

For the representation of solid bodies, DOFAS adopts a discrete direct-forcing
Immersed Boundary Method (IBM) with pointwise interpolating delta functions,
which has been validated in studies including tidal turbines (Ouro et al.,
2017a; Ouro & Stoesser, 2019), geophysical flows (Ouro et al., 2017b), rough
open-channel flows (Stoesser, 2010; Bomminayuni & Stoesser, 2011; Nikora et
al., 2019), and fluid-structure interaction (Kara et al., 2015; Ouro et al.,
2019b). In the present work, the IBM is adopted to represent the turbine
nacelles, using the $\phi_{4}$ delta function for the interpolation
procedures. Turbine rotors are represented using an Actuator Line Model (ALM)
validated in Ouro et al. (2019c) which discretises the blades into a set of
$N_{L}$ points, evenly spaced as a function of the mesh resolution. The ALM
has been proven to provide an adequate description of the wake dynamics for
wind and tidal turbines (Breton et al., 2017).

In this study we set turbine rotors to have a constant rotational speed,
$\Omega$, and use prescribed lift and drag coefficients of hydrofoils
tabulated for a range of angles-of-attack to obtain the lift and drag forces
at every point comprising the turbine blades. From the drag and lift forces we
calculate the thrust force $T$ and tangential force $Q$, and thus determine
the generated power $P=Q\Omega$. Eventually, the coefficients of thrust
($C_{T}$) and power ($C_{P}$) are computed as

$\displaystyle C_{T}=\frac{T}{\frac{1}{2}\rho\pi(D/2)^{2}U_{0}^{2}}$ (30)
$\displaystyle C_{P}=\frac{P}{\frac{1}{2}\rho\pi(D/2)^{2}U_{0}^{3}}$ (31)

After the computation of hydrodynamic force from ALM at each time step, the
force exerted from every Lagrangian point comprising the turbine rotors (as
well as the force computed from IBM for nacelles) is transferred back to the
fluid grid to correct the Eulerian velocity field. This interpolation
procedure is performed using an isotropic Gaussian projection (Shen et al.,
2005),

$f_{L_{ALM}}(x_{i})=\frac{1}{\varepsilon^{3}\pi^{3/2}}\text{exp}\left(-\frac{r_{L}^{2}}{\varepsilon^{2}}\right)$
(32)

where $r_{L}$ denotes the radial distance between the marker $L$ and the
considered cell face $i$, and $\varepsilon$ is the interpolation stencil set
to 3.0$\Delta x_{i}$. A Prandtl-type tip-loss correction is adopted to correct
the ALM forcing near the blade tip as a function of the number of blades and
the tip-speed ratio (Shen et al., 2005).

## References

  * Adcock et al. (2021) Adcock, T.A.A., Draper, S., Willden, R.H.J. & Vogel, C.R. 2021 The fluid mechanics of tidal stream energy conversion. Annual Review of Fluid Mechanics 53, 287–310.
  * Afgan et al. (2013) Afgan, I., McNaughton, J., Rolfo, S., Apsley, D.D., Stallard, T. & Stansby, P. 2013 Turbulent flow and loading on a tidal stream turbine by LES and RANS. International Journal of Heat and Fluid Flow 43, 96–108.
  * Bomminayuni & Stoesser (2011) Bomminayuni, S. & Stoesser, T. 2011 Turbulence statistics in an open-channel flow over a rough bed. Journal of Hydraulic Engineering 137 (11), 1347–1358.
  * Breton et al. (2017) Breton, S.P., Sumner, J., Sorensen, J.N., Hansen, K.S., Sarmast, S. & Ivanell, S. 2017 A survey of modelling methods for high-fidelity wind farm simulations using large eddy simulation. Phil. Trans. R. Soc. A 375, 20160097.
  * Calaf et al. (2010) Calaf, M., Meneveau, C. & Meyers, J. 2010 Large eddy simulation study of fully developed wind-turbine array boundary layers. Physics of Fluids 22, 015100.
  * Cevheri et al. (2016) Cevheri, M., McSherry, R. & Stoesser, T. 2016 A local mesh refinement approach for large-eddy simulations of turbulent flows. International Journal for Numerical Methods in Fluids 82, 261–285.
  * Chawdhary et al. (2017) Chawdhary, S., Hill, C., Yang, X., Guala, M., Corren, D., Colby, J. & Sotiropoulos, F. 2017 Wake characteristics of a TriFrame of axial-flow hydrokinetic turbines. Renewable Energy 109, 332–345.
  * Churchfield et al. (2013) Churchfield, M.J., Li, Y. & Moriarty, P.J. 2013 A large-eddy simulation study of wake propagation and power production in an array of tidal-current turbines. Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 371, 20120421\.
  * Creed et al. (2017) Creed, M.J., Draper, S., Nishino, T. & Borthwick, A.G.L. 2017 Flow through a very porous obstacle in a shallow channel. Proceedings of the Royal Society A 473, 20160672\.
  * Culley et al. (2018) Culley, D.M., Funke, S.W., Kramer, S.C. & Piggott, M.D. 2018 Integration of cost modelling within the micro-siting design optimisation of tidal turbine arrays. Renewable Energy 85, 215–227.
  * De Dominicis et al. (2018) De Dominicis, M., Wolf, J. & M.R., O’Hara 2018 Comparative effects of climate change and tidal stream energy extraction in a shelf sea. Journal of Geophysical Research: Oceans 123, 5041–5067.
  * Draper & Nishino (2014) Draper, S. & Nishino, T. 2014 Centred and staggered arrangements of tidal turbines. Journal of Fluid Mechanics 739, 72–93.
  * Draper et al. (2016) Draper, S., Nishino, T., Adcock, T.A.A. & Taylor, P.H. 2016 Performance of an ideal turbine in an inviscid shear flow. Journal of Fluid Mechanics 796, 86–112.
  * Foti et al. (2019) Foti, D., Yang, X., Shen, L. & Sotiropoulos, F. 2019 Effect of wind turbine nacelle on turbine wake dynamics in large wind farms. Journal of Fluid Mechanics 869, 1–26.
  * Garcia Novo & Kyozuka (2019) Garcia Novo, P. & Kyozuka, Y. 2019 Analysis of turbulence and extreme current velocity values in a tidal channel. Journal of Marine Science and Technology 24, 659–672.
  * Garrett & Cummins (2007) Garrett, C. & Cummins, P. 2007 The efficiency of a turbine in a tidal channel. Journal of Fluid Mechanics 588, 243–251.
  * Gupta & Young (2017) Gupta, V. & Young, A.M. 2017 A one-dimensional model for tidal array design based on three-scale dynamics. Journal of Fluid Mechanics 825, 651–676.
  * Houlsby et al. (2008) Houlsby, G.T., Draper, S. & Oldfield, M.L.G. 2008 Application of linear momentum actuator disc theory to open channel flow. Tech. Rep. OUEL 2296/08. Dept. Engineering Science, University of Oxford, UK.
  * Kang et al. (2014) Kang, S., Yang, X. & Sotiropoulos, F. 2014 On the onset of wake meandering for an axial flow turbine in a turbulent open channel flow. Journal of Fluid Mechanics 744, 376–403.
  * Kara et al. (2015) Kara, M.C., Stoesser, T. & McSherry, R. 2015 Calculation of fluid–structure interaction: methods, refinements, applications. Proceedings of the ICE - Engineering and Computational Mechanics 168, 59–78.
  * Nicoud & Ducros (1999) Nicoud, F & Ducros, F 1999 Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion 62 (3), 183–200.
  * Nikora et al. (2019) Nikora, V.I., Stoesser, T., Cameron, S. M., Stewart, M, Papadopoulos, K, Ouro, P., McSherry, R., Zampiron, A., Marusic, I. & Falconer, R.A. 2019 Friction factor decomposition for rough-wall flows: theoretical background and application to open-channel flows. Journal of Fluid Mechanics 872, 626–664.
  * Nishino & Draper (2019) Nishino, T. & Draper, S. 2019 Theoretical prediction of the efficiency of very large turbine arrays: combined effects of local blockage and wake mixing. In 7th Oxford Tidal Energy Workshop, pp. 31–32.
  * Nishino & Dunstan (2020) Nishino, T. & Dunstan, T.D. 2020 Two-scale momentum theory for time-dependent modelling of large wind farms. Journal of Fluid Mechanics 894, A2.
  * Nishino & Willden (2012) Nishino, T. & Willden, R.H.J. 2012 The efficiency of an array of tidal turbines partially blocking a wide channel. Journal of Fluid Mechanics 708, 596–606.
  * Nishino & Willden (2013) Nishino, T. & Willden, R.H.J. 2013 Two-scale dynamics of flow past a partial cross-stream array of tidal turbines. Journal of Fluid Mechanics 730, 220–244.
  * Noble et al. (2020) Noble, D.R., Draycott, S., Nambiar, A., Sellar, B.G., Steynor, J. & Kiprakis, A. 2020 Experimental assessment of flow, performance, and loads for tidal turbines in a closely-spaced array. Energies 13, 1977.
  * Ouro et al. (2019a) Ouro, P., Fraga, B., Lopez-Novoa, U. & Stoesser, T. 2019a Scalability of an Eulerian-Lagrangian large-eddy simulation solver with hybrid MPI/OpenMP parallelisation. Computers & Fluids 179, 123–136.
  * Ouro et al. (2017a) Ouro, P., Harrold, M., Stoesser, T. & Bromley, P. 2017a Hydrodynamic loadings on a horizontal axis tidal turbine prototype. Journal of Fluids and Structures 71, 78–95.
  * Ouro et al. (2019b) Ouro, P., Muhawenimana, V. & Wilson, C.A.M.E. 2019b Asymmetric wake of a horizontal cylinder in close proximity to a solid boundary for Reynolds numbers in the subcritical turbulence regime. Physical Review Fluids 4, 104604\.
  * Ouro et al. (2019c) Ouro, P., Ramírez, L. & Harrold, M. 2019c Analysis of array spacing on tidal stream turbine farm performance using Large-Eddy Simulation. Journal of Fluids and Structures 91, 102732.
  * Ouro & Stoesser (2019) Ouro, P. & Stoesser, T. 2019 Impact of Environmental Turbulence on the Performance and Loadings of a Tidal Stream Turbine. Flow, Turbulence and Combustion 102 (3), 613–639.
  * Ouro et al. (2017b) Ouro, P., Wilson, C.A.M.E., Evans, P. & Angeloudis, A. 2017b Large-eddy simulation of shallow turbulent wakes behind a conical island. Physics of Fluids 29 (12), 126601.
  * Porte-Agel et al. (2020) Porte-Agel, F., Bastankhah, M. & Shamsoddin, S. 2020 Wind-turbine and wind-farm flows: a review. Boundary-Layer Meteorology 174, 1–59.
  * Scarlett et al. (2019) Scarlett, G.T., Sellar, B., van den Bremer, T. & Viola, I.M. 2019 Unsteady hydrodynamics of a full-scale tidal turbine operating in large wave conditions. Renewable Energy 143, 199–213.
  * Sharma et al. (2018) Sharma, V., Cortina, G., Margairaz, F., Parlange, M. B. & Calaf, M. 2018 Evolution of flow characteristics through finite-sized wind farms and influence of turbine arrangement. Renewable Energy 115, 1196–1208.
  * Shen et al. (2005) Shen, W.Z., Sørensen, J.N. & Mikkelsen, R. 2005 Tip loss correction for actuator/Navier–Stokes computations. Journal of Solar Energy Engineering 127, 209–213.
  * Stallard et al. (2013) Stallard, T., Collings, R, Feng, T & Whelan, J 2013 Interactions between tidal turbine wakes: experimental study of a group of three-bladed rotors. Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 371, 20120159\.
  * Stansby & Stallard (2015) Stansby, P. & Stallard, T. 2015 Fast optimisation of tidal stream turbine positions for power generation in small arrays with low blockage based on superposition of self-similar far-wake velocity deficit profiles. Renewable Energy 92, 366–375.
  * Stoesser (2010) Stoesser, T. 2010 Physically realistic roughness closure scheme to simulate turbulent channel flow over rough beds within the framework of LES. Journal of Hydraulic Engineering 136, 812–819.
  * Vennell (2010) Vennell, R. 2010 Tuning turbines in a tidal channel. Journal of Fluid Mechanics 663, 253–267.
  * Vennell (2011) Vennell, R. 2011 Tuning tidal turbines in-concert to maximise farm efficiency. Journal of Fluid Mechanics 671, 587–604.
  * Vennell (2012) Vennell, R. 2012 The energetics of large tidal turbine arrays. Renewable Energy 48, 210–219.
  * Vennell et al. (2015) Vennell, R., Funke, S.W., Draper, S., Stevens, C. & Divett, T. 2015 The energetics of large tidal turbine arrays. Renewable and Sustainable Energy Reviews 41, 454–472.
  * Vogel et al. (2016) Vogel, C.R., Houlsby, G.T. & Willden, R.H.J. 2016 Effect of free surface deformation on the extractable power of a finite width turbine array. Renewable Energy 88, 317–324.
  * West & Lele (2020) West, J.R. & Lele, S.K. 2020 Wind turbine performance in very large wind farms: Betz analysis revisited. Energies 13, 1078.
  * Yang & Sotiropoulos (2019) Yang, X. & Sotiropoulos, F. 2019 Wake characteristics of a utility-scale wind turbine under coherent inflow structures and different operating conditions. Physical Review Fluids 4, 024604.