2.1 Channel-scale conservation equations

In the absence of turbines, the axial velocity profile is assumed to be uniform in the vertical and cross-span directions. The mass and momentum conservation equations across the small axial element of length $\text{d}x$ shown in figure 1, which must be satisfied under all conditions, are:

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle W_{c}h(x)U(x)=W_{c}h(x+\text{d}x)U(x+\text{d}x)\;\Longrightarrow \;h(x)\frac{\text{d}U}{\text{d}x}=-U(x)\frac{\text{d}h}{\text{d}x}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & & \displaystyle W_{c}\int _{0}^{h(x)}[p_{a}+\unicode[STIX]{x1D70C}g(h(x)-y)]\,\text{d}y-W_{c}\int _{0}^{h(x+\text{d}x)}[p_{a}+\unicode[STIX]{x1D70C}g(h(x+\text{d}x)-y)]\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle \qquad -\,W_{c}p_{a}(h(x)-h(x+\text{d}x))-W_{c}E_{\unicode[STIX]{x0394}h}\unicode[STIX]{x1D6FF}_{xe}\,\text{d}x-W_{c}\unicode[STIX]{x1D70F}_{w}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =W_{c}\unicode[STIX]{x1D70C}U(x+\text{d}x)^{2}h(x+\text{d}x)-W_{c}\unicode[STIX]{x1D70C}U(x)^{2}h(x).\end{eqnarray}$$

In the momentum equation, the first three terms are the pressure force terms, $E_{\unicode[STIX]{x0394}h}\unicode[STIX]{x1D6FF}_{xe}$ is the kinetic energy loss term, $\unicode[STIX]{x1D70F}_{w}$ is the seabed friction force per unit area and the two terms on the right-hand side are the inertial terms. The kinetic energy loss term ( $E_{\unicode[STIX]{x0394}h}\unicode[STIX]{x1D6FF}_{xe}$ ) and seabed friction ( $\unicode[STIX]{x1D70F}_{w}$ ) are also included in the models given by Garrett & Cummins (Reference Garrett and Cummins2005) and Bryden & Couch (Reference Bryden and Couch2007) and together they make up the inherent channel drag.

To understand the origin of the kinetic energy loss term, consider the fact that flow enters the channel from a very large basin. It, therefore, has a very small velocity before the inlet. The flow then gains kinetic energy, which is proportional to the square of the channel flow velocity at inlet ( $U_{1}$ ). This gain, which is given as $W_{c}\int _{x_{1}}^{x_{5}}E_{\unicode[STIX]{x0394}h}\unicode[STIX]{x1D6FF}_{xe}\,\text{d}x\approx W_{c}(h_{1}/2)\unicode[STIX]{x1D70C}U_{1}^{2}$ , must come at the cost of tidal head. Here, the delta function, $\unicode[STIX]{x1D6FF}_{xe}$ , signifies that we account for this term at the channel exit (location 5).

The seabed friction term is commonly taken to be proportional to the square of the flow velocity as $\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D70C}f_{1}U(x)^{2}$ , where $f_{1}$ is a non-dimensional friction coefficient which is assumed to be constant throughout the channel. The channel modelled here has a constant width and the tidal head drop ( $\unicode[STIX]{x0394}h$ ) is usually small. Consequently, $U(x)$ does not change much throughout the channel (see (2.3)). So we further simplify the seabed friction to $\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D70C}f_{1}U_{1}^{2}$ .

The models for the two terms of the inherent channel drag are inserted into (2.3) along with (2.2). The resulting equation is then integrated to give a relationship between the undisturbed channel flow velocity ( $U_{1}$ ), the tidal head term ( $\unicode[STIX]{x0394}h$ ) and the drag terms ( $f=f_{1}L_{c}+h_{1}/2$ ):

(2.4) $$\begin{eqnarray}U_{1}^{2}=\frac{g\unicode[STIX]{x0394}h{\displaystyle \frac{h_{1}+h_{5}}{2}}}{{\displaystyle \frac{h_{1}}{h_{5}}}\unicode[STIX]{x0394}h+f}\;\Longrightarrow \;Fr^{2}=\frac{U_{1}^{2}}{gh_{1}}=\frac{\unicode[STIX]{x0394}H(1+H_{5})}{C_{f}},\end{eqnarray}$$

where $C_{f}=2[\unicode[STIX]{x0394}h/h_{5}+f/h_{1}]$ is the inherent channel drag coefficient and $\unicode[STIX]{x0394}H=\unicode[STIX]{x0394}h/h_{1}$ and $H_{5}=h_{5}/h_{1}$ are the non-dimensional tidal head drop and channel outlet water column height, respectively.

In the presence of a turbine array at one location (say $x_{T}$ ) in the channel, the momentum equation (2.3) is modified such that the thrust imparted by turbines ( $T_{A}$ ) is added to the equation as shown by the bold term below:

(2.5) $$\begin{eqnarray}\displaystyle & & \displaystyle W_{c}\int _{0}^{h(x)}[p_{a}+\unicode[STIX]{x1D70C}g(h(x)-y)]\,\text{d}y-W_{c}\int _{0}^{h(x+\text{d}x)}[p_{a}+\unicode[STIX]{x1D70C}g(h(x+\text{d}x)-y)]\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle \qquad -\,W_{c}p_{a}(h(x)-h(x+\text{d}x))-W_{c}E_{\unicode[STIX]{x0394}h}\unicode[STIX]{x1D6FF}_{xe}\,\text{d}x-W_{c}\unicode[STIX]{x1D70F}_{w}\,\text{d}x-\boldsymbol{T}_{\boldsymbol{A}}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{\boldsymbol{T}})\,\mathbf{d}\boldsymbol{x}\nonumber\\ \displaystyle & & \displaystyle \quad =W_{c}\unicode[STIX]{x1D70C}U(x+\text{d}x)^{2}h(x+\text{d}x)-W_{c}\unicode[STIX]{x1D70C}U(x)^{2}h(x).\end{eqnarray}$$

The cross-stream variation in the axial flow velocity caused by the array will affect the friction term in the channel $W_{c}\unicode[STIX]{x1D70F}_{w}\,\text{d}x$ . This effect is ignored here for now, but it will be considered in § 2.3, after the array-scale dynamics is introduced. Integrating this modified momentum equation along the length of the channel gives:

(2.6) $$\begin{eqnarray}U_{1nd}^{2}=\frac{U_{1d}^{2}}{U_{1}^{2}}=1-\frac{T_{A}}{{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}U_{1}^{2}nA_{d}}\frac{nA_{d}}{W_{c}h_{1}}\frac{1}{C_{f}},\end{eqnarray}$$

where $U_{1d}$ is the reduced upstream channel flow velocity and $U_{1nd}$ is its non-dimensional form. The term $T_{A}/((\unicode[STIX]{x1D70C}U_{1}^{2}nA_{d})/2)$ is a non-dimensional form of the thrust imparted by the turbines and will be referred to as the global thrust coefficient ( $C_{TG}$ ). Equation (2.6) can then be rewritten as:

(2.7) $$\begin{eqnarray}U_{1nd}^{2}=1-C_{TG}B_{g}/C_{f}.\end{eqnarray}$$

This equation shows that the channel flow velocity at the inlet reduces as the thrust imparted by the turbines, $C_{TG}$ , increases. The channel-scale conservation equation (2.7) is our first model equation and contains two unknown variables: $U_{1nd}$ and  $C_{TG}$ .

2.1.1 Estimation of $C_{f}$

The inherent channel drag coefficient ( $C_{f}$ ) is a channel parameter which governs the effect of added drag from the turbines. In our model, we vary $C_{f}$ continuously, but here we give a brief discussion of its physical meaning and of its magnitude in typical marine flows. By assuming that the non-dimensional tidal head drop is small ( $\unicode[STIX]{x0394}H\ll 1$ – this will be true even in shallow channels), the inherent channel drag coefficient is approximated as:

(2.8) $$\begin{eqnarray}C_{f}=2\frac{\unicode[STIX]{x0394}H}{H_{5}}+2\frac{f_{1}L_{c}+h_{1}/2}{h_{1}}\approx 1+2\,\frac{f_{1}L_{c}}{h_{1}},\end{eqnarray}$$

where $f_{1}$ depends up on the roughness of the seabed. This estimation of $C_{f}$ is similar to that of Bryden & Couch (Reference Bryden and Couch2007), who use a different friction model. They give the coefficient as $C_{f}=1+2gn_{m}^{2}L_{c}/R_{d}^{(4/3)}$ , where $R_{d}$ is the hydraulic radius of the channel and $n_{m}$ is the Manning seabed friction coefficient. For the channel considered in Bryden & Couch (Reference Bryden and Couch2007) (where $L_{c}=4000~\text{m}$ , $W_{c}=500~\text{m}$ , $h_{1}=40~\text{m}$ and $n_{m}=0.025~\text{sm}^{-1/3}$ ), the value of $C_{f}$ is approximately 1.44 (or $1/C_{f}\approx 0.70$ ).

Equation (2.7) shows that small values of $C_{f}$ are desirable, because smaller $C_{f}$ means larger $U_{1}$ for a given $\unicode[STIX]{x0394}H$ and therefore greater available kinetic energy for a given tidal head. For channels with small $C_{f}$ , the drag added by the turbines is comparable to the inherent channel drag. This implies that the added drag reduces the channel flow velocity at inlet, i.e. $U_{1nd}<1$ . For those channels with large $C_{f}$ , the drag added by the turbines is negligible compared to the inherent channel drag. This implies that the added drag has no effect on the upstream channel flow velocity, i.e. $U_{1nd}\approx 1$ .

2.2 Array-scale conservation equations

Figure 2 shows the (a) front, (b) top and (c) side views of the turbine array. The front view shows that $n$ turbines are arranged symmetrically in the centre of the channel. They are regularly spaced with distance $L_{a}$ between each other, and each turbine is represented as an actuator disc of diameter $L_{d}$ (area $A_{d}=\unicode[STIX]{x03C0}L_{d}^{2}/4$ ). The height of the centre of turbines from the bottom of the channel is $h_{T}$ , but this term will not appear in the final equations because the upstream velocity profile is assumed to be uniform in the vertical direction. The array inlet and outlet are marked as $1A$ and $5A$ , respectively. We assume that the tidal array is close to the channel inlet such that axial locations 1 in figure 1 and  $1A$ in figure 2 coincide. We also assume that the channel is long enough that location  $5A$ is upstream of location 5, which means that all the wake mixing behind the array is completed before the channel outlet.

The flow that passes through each turbine is called the device-scale core flow and the flow that passes between the turbines is called the device-scale bypass flow. Similarly, the flow that passes through the array is called the array-scale core flow (this flow includes the device-scale core and bypass flows) and the flow that passes between the array and the channel boundaries is called the array-scale bypass flow. The wake mixing regions between $4L$ and  $5L$ (device-scale far-wake mixing) and between $4A$ and  $5A$ (array-scale far-wake mixing) are shaded in the figure. It is assumed in the model that wake mixing is limited to the shaded regions and that the device-scale and array-scale wake mixing processes happen in separate regions. In reality, however, much of the array-scale wake mixing in the vertical direction is likely to happen between $4L$ and  $5L$ (array-scale near-wake region) as opposed to between $4A$ and  $5A$ . These assumptions, of separate regions of wake mixing for the device and array scales and of the absence of wake mixing in the near-wake region, are limitations of the present model and their implications are discussed in § 4.

At $1A$ , both the velocity ( $U_{1nd}$ ) and channel height ( $h_{1}$ ) are uniform across the whole channel width. The locations immediately upstream and downstream of the turbines are marked as  $U$ and  $D$ , respectively. Neither the velocity nor the channel height are uniform at these locations. The average values of the array-scale dynamical variables between $U$ and  $D$ are denoted by the subscript  $2A$ . The location downstream of turbines where the channel height ( $h_{4}$ ) is uniform across the channel width but the flow velocity is still varying is marked as  $4A$ . The location far downstream where the channel height ( $h_{5}$ ) and flow velocity are both uniform is marked as  $5A$ .

At first, we neglect seabed friction at the array scale and derive the conservation equations. Mass conservation between $1A$ and  $4A$ gives:

(2.9) $$\begin{eqnarray}\displaystyle U_{1d}h_{1}W_{c} & = & \displaystyle U_{1d}\unicode[STIX]{x1D6FD}_{4A}(h_{4}W_{c}-A_{4Ac})+U_{1d}\unicode[STIX]{x1D6FC}_{4A}A_{4Ac}\nonumber\\ \displaystyle & & \displaystyle \;\Longrightarrow \;1=H_{4}\unicode[STIX]{x1D6FD}_{4A}+\unicode[STIX]{x1D6FC}_{2A}B_{a}\left(1-\frac{\unicode[STIX]{x1D6FD}_{4A}}{\unicode[STIX]{x1D6FC}_{4A}}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{2A}$ is the ratio of the core flow velocity at $2A$ to the upstream flow velocity at $1A$ , $\unicode[STIX]{x1D6FC}_{4A}$ and $\unicode[STIX]{x1D6FD}_{4A}$ are the ratios of the core and bypass flow velocities, respectively, at $4A$ to the upstream flow velocity at $1A$ . Finally, $A_{4Ac}=nL_{a}L_{d}(\unicode[STIX]{x1D6FC}_{2A}/\unicode[STIX]{x1D6FC}_{4A})$ is the array core flow area at  $4A$ , and $H_{4}=h_{4}/h_{1}$ .

Considering energy conservation in the array core flow region between $1A$ and  $4A$ gives:

(2.10) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2}\unicode[STIX]{x1D70C}U_{1d}^{3}\unicode[STIX]{x1D6FC}_{4A}^{3}A_{4Ac}+U_{1d}\unicode[STIX]{x1D6FC}_{4A}\frac{A_{4Ac}}{L_{d4}}\int _{h_{T}-L_{d4}/2}^{h_{T}+L_{d4}/2}[p_{a}+\unicode[STIX]{x1D70C}g(h_{4}-y)]\,\text{d}y+P_{A}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{2}\unicode[STIX]{x1D70C}U_{1d}^{3}A_{1Ac}+U_{1d}\frac{A_{1Ac}}{L_{d1}}\int _{h_{T}-L_{d1}/2}^{h_{T}+L_{d1}/2}[p_{a}+\unicode[STIX]{x1D70C}g(h_{1}-y)]\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle \qquad \;\Longrightarrow \;U_{1nd}^{2}(1-\unicode[STIX]{x1D6FC}_{4A}^{2})+\frac{2}{Fr^{2}}(1-H_{4})=B_{l}C_{TG},\end{eqnarray}$$

where $A_{1Ac}=nL_{a}L_{d}\unicode[STIX]{x1D6FC}_{2A}$ is the array core flow area at $1A$ , $L_{d1}$ and $L_{d4}$ are the vertical extents of the array core flow regions at $1A$ and  $4A$ , respectively, and $P_{A}$ is the power extracted from the flow by the array of turbines (given by $T_{A}U_{1d}\unicode[STIX]{x1D6FC}_{2A}$ ). Note that $P_{A}$ is the sum of the useful power extracted by the turbines ( $nP_{D}=T_{A}U_{1d}\unicode[STIX]{x1D6FC}_{2A}\unicode[STIX]{x1D6FC}_{2L}$ , see § 2.4) and the power dissipated in all the wake mixing regions between $4L$ and  $5L$ .

Energy conservation for the bypass flow can be satisfied by using Bernoulli’s equation, which gives:

(2.11) $$\begin{eqnarray}1-H_{4}=\frac{Fr^{2}}{2}U_{1nd}^{2}(\unicode[STIX]{x1D6FD}_{4A}^{2}-1).\end{eqnarray}$$

The momentum conservation equation for the array-scale flow between $1A$ and $4A$ is given by the relation:

(2.12) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}U_{1d}^{2}h_{1}W_{c}-\unicode[STIX]{x1D70C}U_{1d}^{2}\unicode[STIX]{x1D6FD}_{4A}^{2}(h_{4}W_{c}-A_{4Ac})-\unicode[STIX]{x1D70C}U_{1d}^{2}\unicode[STIX]{x1D6FC}_{4A}^{2}A_{4Ac}=p_{a}(h_{1}-h_{4})W_{c}\nonumber\\ \displaystyle & & \displaystyle \quad +\,T_{A}+W_{c}\int _{0}^{h_{4}}[p_{a}+\unicode[STIX]{x1D70C}gy]\,\text{d}y-W_{c}\int _{0}^{h_{1}}[p_{a}+\unicode[STIX]{x1D70C}gy]\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle \quad \;\Longrightarrow \;2U_{1nd}^{2}[1-\unicode[STIX]{x1D6FD}_{4A}+\unicode[STIX]{x1D6FC}_{2A}B_{a}(\unicode[STIX]{x1D6FD}_{4A}-\unicode[STIX]{x1D6FC}_{4A})]=C_{TG}B_{g}-\frac{1}{Fr^{2}}(1-H_{4}^{2}).\end{eqnarray}$$

The array-scale conservation equations, therefore, give us four additional equations, ((2.9)–(2.12)), and four additional unknown variables, which are $H_{4}$ , $\unicode[STIX]{x1D6FC}_{2A}$ , $\unicode[STIX]{x1D6FC}_{4A}$ and $\unicode[STIX]{x1D6FD}_{4A}$ . Together with the equation from the channel scale (2.7), there are now five equations and six unknown variables.

2.3 Channel and array-scale friction terms

The installation of a tidal array in a channel causes an increase in the total channel drag, which is referred to as added drag in this paper. This added drag comes from two sources. The first is directly, from the drag imparted by the turbines on the flow, its effect appears as $C_{TG}$ in (2.7). The second is indirectly, from the cross-stream axial flow variations caused by the tidal array. The friction force in a channel depends on the flow velocity and here it is modelled to be proportional to the square of the velocity as explained in § 2.1. Consequently, the cross-stream axial flow variations, if significant, cause an increase in the friction force in the channel. This second source of added drag is non-negligible only when the array length is comparable to the channel length because only then the cross-stream axial flow variations persist for a significant proportion of the channel length to affect the total friction force. When the array length is long, such that it is comparable to the channel length, then another effect that can come into play is the friction at the array scale. In this subsection, we will include the effects of (i) the tidal array-induced axial flow variations on the friction force in the channel and (ii) the friction force on the array-scale dynamics.

Figure 3 shows a top view schematic of the channel-scale flow with the tidal array represented as a box with blockage area ratio $B_{a}$ and length $L_{15A}$ . As mentioned before, the change in the water column height along the axial direction is usually small, so we assume that the average velocity remains constant at every axial location and take it to be $U_{1d}$ everywhere for the purpose of estimating the friction force. We make a further approximation by using average array bypass and core flow velocities throughout the array length ( $L_{15A}$ ). The maximum array bypass flow velocity factor is $\unicode[STIX]{x1D6FD}_{4A}$ , which is achieved at location  $4A$ , and the minimum is 1, which is achieved at locations  $1A$ and  $5A$ . The average array bypass flow from  $1A$ to  $5A$ ( $\unicode[STIX]{x1D6FD}_{b}$ ), therefore, is modelled approximately as: $\unicode[STIX]{x1D6FD}_{b}=0.5(1+\unicode[STIX]{x1D6FD}_{4A})$ . In order to satisfy the conservation of mass, the average core flow from $1A$ to $5A$ ( $\unicode[STIX]{x1D6FC}_{c}$ ) is given by the relation below.

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{c}B_{a}+0.5(\unicode[STIX]{x1D6FD}_{4A}+1)(1-B_{a})=1.\end{eqnarray}$$

Figure 3.

Top view schematic of the channel-scale flow with the tidal array represented as a box with blockage area ratio  $B_{a}$ .

The array length ( $L_{15A}$ ) varies with the local turbine arrangement ( $B_{l}$ ). On the one hand, for $B_{l}\approx 0$ , i.e. when devices are very far apart from each other, there is very little interaction between the devices and the array-scale wake length will be close to the wake length for a single device. On the other hand, for $B_{l}\approx 1$ , i.e. when devices are almost touching, the flow will see the devices as one large device and the array-scale wake length will be scaled accordingly. We assume a simple linear relation for $L_{15A}$ at intermediate values of $B_{l}$ as:

(2.14) $$\begin{eqnarray}L_{15A}=L_{Array}\left(\frac{(n-1)B_{l}+1}{n}\right),\end{eqnarray}$$

where $L_{Array}$ is the maximum length of the array that is achieved when $B_{l}=1$ and $L_{Array}$ itself depends on the global blockage. It will increase monotonically with $B_{g}$ for a given channel, and we assume this increase will also be linear. In this paper, the channel length is always given in terms of the ratio $L_{A/C}=L_{Array}/L_{c}$ at a given $B_{g}$ , and this ratio ( $L_{A/C}$ ) changes linearly with $B_{g}$ . For example, if the channel length is given as $L_{A/C}=0.5$ at $B_{g}=0.2$ , then for this channel $L_{A/C}=0.25$ at $B_{g}=0.1$ .

Based on the above assumptions, the friction force across the channel is given as:

(2.15) $$\begin{eqnarray}\displaystyle W_{c}\int _{0}^{L_{c}}\unicode[STIX]{x1D70F}_{w}\,\text{d}x & = & \displaystyle W_{c}\unicode[STIX]{x1D70C}f_{1}U_{1d}^{2}(L_{c}-L_{15A})+W_{c}B_{a}\unicode[STIX]{x1D70C}f_{1}U_{1d}^{2}\unicode[STIX]{x1D6FC}_{c}^{2}L_{15A}\nonumber\\ \displaystyle & & \displaystyle +\,W_{c}(1-B_{a})\unicode[STIX]{x1D70C}f_{1}U_{1d}^{2}\left(\frac{1+\unicode[STIX]{x1D6FD}_{4A}}{2}\right)^{2}L_{15A}\nonumber\\ \displaystyle & = & \displaystyle W_{c}\unicode[STIX]{x1D70C}f_{1}U_{1d}^{2}\left[L_{c}+L_{15A}\left(B_{a}\unicode[STIX]{x1D6FC}_{c}^{2}+(1-B_{a})\left(\frac{1+\unicode[STIX]{x1D6FD}_{4A}}{2}\right)^{2}-1\right)\right].\qquad\end{eqnarray}$$

This friction term is then inserted in (2.5), and the resulting equation after integration and non-dimensionalisation is given as:

(2.16) $$\begin{eqnarray}U_{1nd}^{2}=\frac{C_{f}}{C_{f}+C_{fm}}-\frac{C_{TG}B_{g}}{C_{f}+C_{fm}},\end{eqnarray}$$

where $C_{fm}=(2f_{1}L_{15A}/h_{1})(B_{a}\unicode[STIX]{x1D6FC}_{c}^{2}+(1-B_{a})((1+\unicode[STIX]{x1D6FD}_{4A})/2)^{2}-1)$ is the non-dimensional tidal array-induced increase in the channel friction force. A higher value of $C_{fm}$ means a higher reduction in the incoming flow velocity $U_{1nd}$ . It can be seen that $C_{fm}$ is higher when the bypass flow coefficient $\unicode[STIX]{x1D6FD}_{4A}$ and the array length $L_{15A}$ are higher; both of these terms increase with $B_{l}$ and  $B_{g}$ .

In order to include the friction terms in the array-scale conservation equations, we also need an estimation of the length from $1A$ to $4A$ ( $L_{14A}$ ) as a fraction of $L_{15A}$ . The length from $2A$ to $4A$ , which is the near-wake length, is known to be 2–4 times the turbine diameter (Mehta et al. Reference Mehta, van Zuijlen, Koren, Holierhoek and Bijl2014), and the length from $2A$ to  $5A$ , which includes the near as well as the far-wake lengths, is known to be 6–10 times the turbine diameter (Mehta et al. Reference Mehta, van Zuijlen, Koren, Holierhoek and Bijl2014). The length from $1A$ to $2A$ is usually taken to be twice the turbine diameter (Porté-Agel et al. Reference Porté-Agel, Wu, Lu and Conzemius2011). Based on these numbers, $L_{14A}$ is 0.5–0.75 times $L_{15A}$ . In this paper, we assume $L_{14A}=0.7L_{15A}$ . The friction force from $1A$ to $4A$ is then given as:

(2.17) $$\begin{eqnarray}W_{c}\int _{0}^{L_{14A}}\unicode[STIX]{x1D70F}_{w}\,\text{d}x=W_{c}\unicode[STIX]{x1D70C}f_{1}U_{1d}^{2}L_{14A}\left(B_{a}\unicode[STIX]{x1D6FC}_{c}^{2}+(1-B_{a})\left(\frac{1+\unicode[STIX]{x1D6FD}_{4A}}{2}\right)^{2}\right).\end{eqnarray}$$

This term is substituted into (2.12) (momentum conservation between $1A$ and  $4A$ ), which becomes:

(2.18) $$\begin{eqnarray}2U_{1nd}^{2}[1-\unicode[STIX]{x1D6FD}_{4A}+\unicode[STIX]{x1D6FC}_{2A}B_{a}(\unicode[STIX]{x1D6FD}_{4A}-\unicode[STIX]{x1D6FC}_{4A})]=U_{1nd}^{2}C_{fm2}+C_{TG}B_{g}-\frac{1}{Fr^{2}}(1-H_{4}^{2}),\end{eqnarray}$$

where $C_{fm2}=(2f_{1}L_{14A}/h_{1})(B_{a}\unicode[STIX]{x1D6FC}_{c}^{2}+(1-B_{a})((1+\unicode[STIX]{x1D6FD}_{4A})/2)^{2})$ is the non-dimensional friction force at the array scale. This term increases with $\unicode[STIX]{x1D6FD}_{4A}$ , $L_{14A}$ and  $f_{1}$ .

This second additional friction term ( $C_{fm2}$ ) also affects the array-scale core and bypass flow energy equations. The effect of friction is to create shear in the entire flow such that energy gets dissipated everywhere. This dissipation is given as the friction force multiplied by the velocity flux. This means that the dissipation in the array-scale core and bypass flow will be given by multiplying the friction force, $C_{fm2}$ , by $U_{1d}\unicode[STIX]{x1D6FC}_{2A}B_{a}W_{c}h_{1}$ and $U_{1d}W_{c}h_{1}(1-\unicode[STIX]{x1D6FC}_{2A}B_{a})$ , respectively. The energy equations in the array-scale core and bypass flows between $1A$ and  $4A$ , which were given by (2.10) and (2.11), respectively, are now given as:

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle U_{1nd}^{2}(1-\unicode[STIX]{x1D6FC}_{4A}^{2})+\frac{2}{Fr^{2}}(1-H_{4})=B_{l}C_{TG}+U_{1nd}^{2}C_{fm2}, & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle 1-H_{4}=\frac{Fr^{2}}{2}U_{1nd}^{2}[(\unicode[STIX]{x1D6FD}_{4A}^{2}-1)+C_{fm2}], & \displaystyle\end{eqnarray}$$

where (2.20) is equivalent to the Darcy–Weisbach equation, which is a phenomenological equation used for estimating head loss in a pipe or a channel with wall friction. This equation also shows that $C_{fm2}$ causes a reduction in $H_{4}$ , and therefore makes the flow through the turbines faster and thus increases the power generation. It seems counter-intuitive that friction is causing an increase in power generation. However, it must be noted that for a given tidal head, the increase in channel friction will also cause a reduction in the incoming flow velocity (2.16) and so the overall effect will be a reduction in power generation.

Lastly, we modify the definition of the global power coefficient ( $C_{PG}$ ), given in (2.1), to:

(2.21) $$\begin{eqnarray}C_{PG}=\frac{P_{D}}{{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}U_{5A}^{3}A_{d}}=\frac{nP_{D}}{{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}U_{5A}^{3}nA_{d}},\end{eqnarray}$$

where $U_{5A}$ is the flow velocity through location $5A$ when turbines are not present in the channel. This change is to account for the increase in the flow through the turbines because of their placement in the channel.

The introduction of these two effects – the tidal array-induced increase in the channel friction and the friction at the array scale – therefore gives us an additional variable $\unicode[STIX]{x1D6FC}_{c}$ and an additional equation (2.13). Further to this, the channel and array-scale momentum equations, and the array-scale core and bypass flow energy equations are replaced by (2.16), (2.18), (2.19) and (2.20), respectively.

2.4 Device-scale conservation equations

The device scale is between locations $1L$ and $5L$ as marked in figure 2. We assume that the length of the device-scale region is small enough for the effects of friction to be ignored. The flow through the individual devices is separated by the dashed lines and comprises of the device-scale core and bypass flows (shown in figure 2 b). The dimensions of the device-scale flow width and vertical height at locations $1L$ and $4L$ are given in terms of the array-scale area factors ( $m_{1},k_{1},m_{4}$ and $k_{4}$ ) (shown in figure 2 b,c). These factors are modelled empirically in § 2.5. At the upstream location, $1L$ , the velocity ( $U_{1L}=m_{1}k_{1}U_{1d}\unicode[STIX]{x1D6FC}_{2A}$ ) and channel height ( $h_{1L}$ ) are both uniform across the local passage. Locations  $U$ and  $D$ are the positions immediately upstream and downstream of the turbines, respectively. The average values of the device-scale variables at the turbine location have subscript  $2L$ . The location downstream of the turbines where the channel height ( $h_{4L}$ ) is uniform across the local passage width ( $L_{d}/k_{4}$ ) but the flow velocity is still varying is marked  $4L$ . Location  $5L$ is further downstream in the channel after the device-scale wake mixing is completed. The channel height and flow velocity are both uniform across the local passage width at this location.

Mass conservation between $1L$ and $4L$ gives:

(2.22) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{4L}=\frac{1-\unicode[STIX]{x1D6FC}_{2L}B_{l}}{1-{\displaystyle \frac{\unicode[STIX]{x1D6FC}_{2L}}{\unicode[STIX]{x1D6FC}_{4L}}}B_{l}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{2L}$ is the local core flow velocity coefficient at the array location and $\unicode[STIX]{x1D6FC}_{4L}$ and $\unicode[STIX]{x1D6FD}_{4L}$ are the local core and bypass flow velocity coefficients at  $4L$ , respectively.

Energy conservation in the local core flow between $1L$ and $4L$ gives:

(2.23) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2}\unicode[STIX]{x1D70C}(m_{4}k_{4}U_{1d}\unicode[STIX]{x1D6FC}_{2A}\unicode[STIX]{x1D6FC}_{4L})^{3}A_{4Lc}-\frac{1}{2}\unicode[STIX]{x1D70C}(m_{1}k_{1}U_{1d}\unicode[STIX]{x1D6FC}_{2A})^{3}A_{1Lc}+P_{D}\nonumber\\ \displaystyle & & \displaystyle \quad =m_{1}k_{1}U_{1d}\unicode[STIX]{x1D6FC}_{2A}A_{1Lc}p_{avg1L}-m_{4}k_{4}U_{1d}\unicode[STIX]{x1D6FC}_{2A}\unicode[STIX]{x1D6FC}_{4L}A_{4Lc}p_{avg4L}\nonumber\\ \displaystyle & & \displaystyle \qquad \;\Longrightarrow \;U_{1nd}^{2}\unicode[STIX]{x1D6FC}_{2A}^{2}[(m_{1}k_{1})^{2}-(m_{4}k_{4})^{2}\unicode[STIX]{x1D6FC}_{4L}^{2}]+\frac{2}{Fr^{2}}(H_{1L}-H_{4L})=C_{TG},\end{eqnarray}$$

where $P_{D}$ is the useful power extracted from the flow by a single turbine and $H_{1L}=h_{1L}/h_{1}$ and $H_{4L}=h_{4L}/h_{1}$ are the non-dimensional heights across the device scale at $1L$ and $4L$ , respectively. $A_{1Lc}$ and $A_{4Lc}$ are local core flow areas at $1L$ and  $4L$ , respectively, and are given as:

(2.24a,b ) $$\begin{eqnarray}A_{1Lc}=\frac{A_{d}}{m_{1}k_{1}}\unicode[STIX]{x1D6FC}_{2L},\quad A_{4Lc}=\frac{A_{d}}{m_{4}k_{4}}\frac{\unicode[STIX]{x1D6FC}_{2L}}{\unicode[STIX]{x1D6FC}_{4L}}.\end{eqnarray}$$

The terms $p_{avg1L}$ and $p_{avg4L}$ are the average pressures in the local flow passage at $1L$ and $4L$ , respectively, and are approximated by the hydrostatic pressure at the midpoint of the core flow at $1L$ and $4L$ as:

(2.25a,b ) $$\begin{eqnarray}p_{avg1L}=[p_{a}+\unicode[STIX]{x1D70C}g(h_{1L}-h_{T})],\quad p_{avg4L}=[p_{a}+\unicode[STIX]{x1D70C}g(h_{4L}-h_{T})].\end{eqnarray}$$

Energy conservation in the local bypass flow between $1L$ and $4L$ gives:

(2.26) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2}\unicode[STIX]{x1D70C}(m_{4}k_{4}U_{1d}\unicode[STIX]{x1D6FC}_{2A}\unicode[STIX]{x1D6FD}_{4L})^{3}\left(\frac{A_{l}}{m_{4}k_{4}}-A_{4Lc}\right)-\frac{1}{2}\unicode[STIX]{x1D70C}(m_{1}k_{1}U_{1d}\unicode[STIX]{x1D6FC}_{2A})^{3}\left(\frac{A_{l}}{m_{1}k_{1}}-A_{1Lc}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =U_{1d}\unicode[STIX]{x1D6FC}_{2A}\left[m_{1}k_{1}\left(\frac{A_{l}}{m_{1}k_{1}}-A_{1Lc}\right)p_{avg1L}-m_{4}k_{4}\unicode[STIX]{x1D6FC}_{4L}\left(\frac{A_{l}}{m_{4}k_{4}}-A_{4Lc}\right)p_{avg4L}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad \;\Longrightarrow \;\frac{2}{Fr^{2}}(H_{1L}-H_{4L})=U_{1nd}^{2}\unicode[STIX]{x1D6FC}_{2A}^{2}[(m_{4}k_{4}\unicode[STIX]{x1D6FD}_{4L})^{2}-(m_{1}k_{1})^{2}].\end{eqnarray}$$

The term $H_{1L}-H_{4L}$ is eliminated by combining the device-scale core flow and bypass flow energy equations ((2.23) and (2.26)). The resulting equation is:

(2.27) $$\begin{eqnarray}U_{1nd}^{2}\unicode[STIX]{x1D6FC}_{2A}^{2}(m_{4}k_{4})^{2}[\unicode[STIX]{x1D6FD}_{4L}^{2}-\unicode[STIX]{x1D6FC}_{4L}^{2}]=C_{TG}.\end{eqnarray}$$

Momentum conservation for the local passage between $1L$ and $4L$ gives:

(2.28) $$\begin{eqnarray}\displaystyle T_{D} & = & \displaystyle p_{avg1L}\frac{A_{l}}{m_{1}k_{1}}-p_{avg4L}\frac{A_{l}}{m_{4}k_{4}}+p^{\ast }A_{14L}^{\ast }+\unicode[STIX]{x1D70C}(m_{1}k_{1}U_{1d}\unicode[STIX]{x1D6FC}_{2A})^{2}\frac{A_{l}}{m_{1}k_{1}}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D70C}(m_{4}k_{4}U_{1d}\unicode[STIX]{x1D6FC}_{2A}\unicode[STIX]{x1D6FC}_{4L})^{2}A_{4Lc}-\unicode[STIX]{x1D70C}(m_{4}k_{4}U_{1d}\unicode[STIX]{x1D6FC}_{2A}\unicode[STIX]{x1D6FD}_{4L})^{2}\left(\frac{A_{l}}{m_{4}k_{4}}-A_{4Lc}\right)\nonumber\\ \displaystyle & & \displaystyle \;\Longrightarrow \;\frac{(m_{4}k_{4})^{2}}{2}(\unicode[STIX]{x1D6FD}_{4L}^{2}-\unicode[STIX]{x1D6FC}_{4L}^{2})=\frac{m_{1}k_{1}}{Bl}-m_{4}k_{4}\unicode[STIX]{x1D6FC}_{4L}\unicode[STIX]{x1D6FC}_{2L}-\frac{m_{4}k_{4}\unicode[STIX]{x1D6FD}_{4L}^{2}}{Bl}\nonumber\\ \displaystyle & & \displaystyle +\,m_{4}k_{4}\unicode[STIX]{x1D6FD}_{4L}^{2}\frac{\unicode[STIX]{x1D6FC}_{2L}}{\unicode[STIX]{x1D6FC}_{4L}}+\frac{1}{2Blm_{4}k_{4}}[(m_{4}k_{4}\unicode[STIX]{x1D6FD}_{4L})^{2}-(m_{1}k_{1})^{2}].\end{eqnarray}$$

The term $p^{\ast }A_{14L}^{\ast }$ is the pressure force term arising from the non-zero axial projection of the lateral surface area between $1L$ and  $4L$ . This term is modelled in the same way as in Nishino & Willden (Reference Nishino and Willden2013):

(2.29) $$\begin{eqnarray}p_{avg1L}\frac{A_{l}}{m_{1}k_{1}}+p^{\ast }A_{14L}^{\ast }=p_{avg1L}\frac{A_{l}}{m_{4}k_{4}}.\end{eqnarray}$$

The resulting equation then is written as:

(2.30) $$\begin{eqnarray}\displaystyle \frac{(m_{4}k_{4})^{2}}{2}(\unicode[STIX]{x1D6FD}_{4L}^{2}-\unicode[STIX]{x1D6FC}_{4L}^{2}) & = & \displaystyle \frac{m_{1}k_{1}}{B_{l}}-m_{4}k_{4}\unicode[STIX]{x1D6FC}_{4L}\unicode[STIX]{x1D6FC}_{2L}-\frac{m_{4}k_{4}\unicode[STIX]{x1D6FD}_{4L}^{2}}{B_{l}}+m_{4}k_{4}\unicode[STIX]{x1D6FD}_{4L}^{2}\frac{\unicode[STIX]{x1D6FC}_{2L}}{\unicode[STIX]{x1D6FC}_{4L}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{2B_{l}m_{4}k_{4}}[(m_{4}k_{4}\unicode[STIX]{x1D6FD}_{4L})^{2}-(m_{1}k_{1})^{2}].\end{eqnarray}$$

The device-scale conservation equations therefore give us 3 additional model equations, these are (2.22), (2.27) and (2.30). The new variables introduced are $\unicode[STIX]{x1D6FC}_{2L}$ , $\unicode[STIX]{x1D6FC}_{4L}$ and $\unicode[STIX]{x1D6FD}_{4L}$ , as well as the products $m_{1}k_{1}$ and $m_{4}k_{4}$ which are the array-scale area factors.

Having considered the conservation equations over the channel, array and device scales, we have nine equations and twelve unknown variables. In order to close the equations, we now derive three additional equations based on empirical models. These equations govern the drag of the turbines and the array-scale area factors.

2.5 Empirically modelled equations necessary for closure

The turbines are modelled empirically in terms of a local induction factor ( $a$ ), which is defined as the ratio of the deceleration in the local core flow velocity between $1L$ and $2L$ to the local core flow velocity at  $1L$ :

(2.31) $$\begin{eqnarray}a=\frac{m_{1}k_{1}\unicode[STIX]{x1D6FC}_{2A}U_{1d}-\unicode[STIX]{x1D6FC}_{2A}\unicode[STIX]{x1D6FC}_{2L}U_{1d}}{m_{1}k_{1}\unicode[STIX]{x1D6FC}_{2A}U_{1d}}=1-\frac{\unicode[STIX]{x1D6FC}_{2L}}{m_{1}k_{1}}.\end{eqnarray}$$

A higher value of induction factor implies that more drag is imparted on the flow by the turbines. This parameter will be used as the turbine tuning parameter in this paper.

The other empirically modelled quantities are the array-scale area factors $m_{1}$ , $k_{1}$ , $m_{4}$ and $k_{4}$ . These factors quantify the expansion in the array core flow between $1L$ and $4L$ . Because these factors only appear as products, $m_{1}k_{1}$ and $m_{4}k_{4}$ , in the equations, they are modelled together as:

(2.32) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle m_{1}k_{1}=\left[1+\left(\frac{1}{n}\right)(\unicode[STIX]{x1D6FC}_{2A}-1)\right]^{-1},\\ \displaystyle m_{4}k_{4}=\left[1+\left(\frac{1}{n}\right)\left(\frac{\unicode[STIX]{x1D6FC}_{2A}}{\unicode[STIX]{x1D6FC}_{4A}}-1\right)\right]^{-1}.\end{array}\right\}\end{eqnarray}$$

This model for the array-scale expansion factors is given by Nishino & Willden (Reference Nishino and Willden2013). For $n\rightarrow \infty$ , the array scale is much bigger than the device scale, and therefore, $m_{1}k_{1}$ and $m_{4}k_{4}$ both approach 1. For $n=1$ , the array scale and device scale are the same, and therefore, $m_{1}k_{1}=1/\unicode[STIX]{x1D6FC}_{2A}$ and $m_{4}k_{4}=\unicode[STIX]{x1D6FC}_{4A}/\unicode[STIX]{x1D6FC}_{2A}$ . For intermediate values of $n$ , equation (2.32) gives empirical power laws for these area factors.

The final list of 12 equations in terms of 12 unknown variables is given in appendix A. The channel parameters – the inherent drag coefficient ( $C_{f}$ ), the channel length ( $L_{c}$ ) and the Froude number ( $Fr$ ) – are constant for a given tidal channel. The design parameters – the global blockage ( $B_{g}$ ), the local blockage ( $B_{l}$ ) and the induction factor ( $a$ ) – can be varied in order to maximise the power extracted per turbine ( $C_{PG}$ ). This final set of equations is solved using a Newton–Raphson method. Python is used for the numerical calculations.