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

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 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 perform large-eddy simulation (LES) of open-channel flow for 28 layouts of tidal turbines using an actuator line method 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 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.

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 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 perform large-eddy simulation (LES) of open-channel flow for 28 layouts of tidal turbines using an actuator line method 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 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.

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 (Nishino & Dunstan 2020;Porté-Agel, Bastankhah & Shamsoddin 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 array. Effects of vertical mixing between turbine wakes and the atmospheric flow aloft, known to play a key role in large wind farms (Cal et al. 2010;Camp & Cal 2016;Bossuyt, Meneveau & Meyers 2017), are also absent in tidal arrays.
Future tidal arrays will comprise multiple rows of turbines, whose design requires us 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, Wolf & O'Hara 2018). Large tidal arrays deployed at a tidal channel would 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(Vennell , 2011Vennell 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 . 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) and 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 2012Vogel, Houlsby & Willden 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 (2010Vennell ( , 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.

A30-2
Tidal turbines in an infinitely large array 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 a 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 tide directions further complicates the design of optimal array configurations, requiring us 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, Ramírez & Harrold 2019c). However, as the performance of tidal devices in arrays is driven by wake-turbine interactions as well as bathymetry-induced turbulence , 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. An infinitely large array represents 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 the atmospheric stability conditions (Bossuyt et al. 2017;Sharma et al. 2018;Porté-Agel et al. 2020). 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 the 28 LES runs comprising aligned and staggered layouts are presented in § 3, with results of flow characteristics, hydrodynamic coefficients and turbine loading unsteadiness 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.

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, Draper & Oldfield 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 Inv.
Inv. Inv. Inv. Visc. Visc. Visc. rows of tidal turbines, following the idea of the 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.
Although we refer to this periodic flow problem as an 'infinitely large' array problem, the analysis presented below is still relevant to finite-size tidal arrays to be built in the future. The key conditions that would need to be satisfied (for the analysis to be directly relevant) are: (i) the number of turbine rows is large enough for the flow through the array to reach (approximately) a 'fully developed' state; and (ii) the array occupies the entire width of a given tidal channel so that the incoming flow will not expand laterally to bypass the whole array. For the first condition, the number of rows required could be relatively small for tidal arrays (compared with that for wind farms) as the flow is vertically confined, whereas the second condition means that the number of turbines in each row could also be relatively small if the channel is narrow. Hence, this 'infinitely large' array study could be relevant to even medium-sized tidal arrays.
In this study we consider that the mass flow rate through the array is constant and not affected by the resistance caused by the array. In reality, the mass flow rate would be affected by the array if the number of rows is large enough for the array-induced resistance to become non-negligible compared with the natural resistance of the channel. Such a scenario could be studied if the analysis presented below is combined with, for example, a channel-scale momentum balance model (Vennell 2010). The analysis may also be combined with a shallow channel flow model (Creed et al. 2017) to study the array performance in a laterally unconfined (or less confined) situation, where the effect of bed friction becomes important as it affects how the incoming flow expands laterally to bypass the whole array. However, these channel-scale effects are outside the scope of Station 1 Station 4 Station 5 the present study, as our main focus here is on the effect of turbine layout on 'local' power characteristics of turbines within a fully developed region of an array (i.e. power of turbines for a given mass flow rate).
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 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 925 A30-5 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 α and β with subscripts as in figure 2. Each velocity coefficient represents the ratio of the velocity there to a reference velocity, u ref . In the following we take the core-flow velocity at station 1 as u ref (i.e. α 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 where ρ is the fluid density. Substituting (2.2) and (2.3), respectively, into (2.1) gives 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 which, together with (2.1), leads to an expression for the half-disc thrust, T, as and, thus, an expression for the half-disc power, P, as 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 α 5 u ref to α 8 u ref when it passes through the row immediately upstream) due to the staggered configuration. As for β 4a , we consider the conservation of momentum for the entire flow passage to obtain which, together with (2.1) and (2.7), leads to (2.13) The above set of equations for the inviscid zone is not closed as it includes α 5 and β 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 av , as the flow passes through the viscous zone. By applying this to the wake flow of the disc we obtain where ψ = u av /u ref is the non-dimensional cross-sectional average velocity and this can be calculated from the velocity profile at station 1, for example, as (2.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 β 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 β 4a and β 4b tends to be small as depicted in figure 2, and eventually leads to where β 4m is the 'area-weighted' average of β 4a and β 4b , i.e. (2.17) It should be noted that α 9 = mψ + (1 − m)α 8 is required together with (2.14) to (2.17) to conserve the total mass in the viscous zone, but this is automatically satisfied as we enforce α 9 = α 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 av = ψu ref ) as For convenience, we also define the resistance coefficient (or local thrust coefficient) of the disc as from which and (2.7), we obtain 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 (α 2 , α 4 , α 5 , α 8 , β 4a , β 4b , β 4m , β 5 , γ , ψ) and a set of ten equations to be solved numerically: (2.4), (2.5), (2.10), (2.11), (2.13), (2.14), (2.15), (2.16), (2.17) and (2.21).

Aligned rows of actuator discs
Compared with 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 925 A30-8 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 (2.1), (2.2) and (2.6), which lead to (2.4), (2.7) and (2.8) for β 4b , T and P, respectively. Note that these equations are identical for both aligned and staggered cases, whereas (2.3) and (2.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 which leads to an expression for α 4 as (2.23) For β 4a , we consider the conservation of momentum for the entire flow to obtain (2.24) which, together with (2.1) and (2.7), leads to For β 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 (2.16) for the aligned case as well. Note, however, that ψ and β 4m are different for the aligned case, which are We also need (2.14) together with the above equations to conserve the total mass in the viscous zone, but this is automatically satisfied since here we enforce α 5 = α 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 (2.18) to (2.20) and, thus, (2.21) for α 2 .

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 GC07 . Figure 3 shows how the performance of aligned rows of such 'suboptimal' actuator discs (with K = K 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 α 2 (β 2 4a − α 2 4 ) and 1/ψ 3 , to explain why C P decreases as m decreases. As can be seen from (2.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 ref ) instead of the cross-sectional average velocity (u av ), and the latter represents the effect of the difference between u ref and u av . Here we can see that the value of α 2 (β 2 4a − α 2 4 ) 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 with the 'original' power possessed by the 925 A30-10 cross-sectionally averaged flow (i.e. the increase rate of α 2 (β 2 4a − α 2 4 ) is smaller than the decrease rate of 1/ψ 3 ); therefore, C P eventually decreases as m decreases.
Although the above analysis was for 'suboptimal' discs with K = K GC07 , the same relationship between C P and m is generally found for aligned discs with any given K values. Figures 4(a) and 4(b) show C P vs 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 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 Pmax , changes with m. As can be expected, for aligned rows of actuator discs, C Pmax also decreases monotonically as m decreases. Figure 5 shows C P vs C T curves for the same five different m values as in figures 4(a) and 4(b), again at B = 0.05 and 0.2. The same trend can be confirmed here as well.
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 GC07 , which is shown in figures 6(a) and 6(b) for B = 0.05 and 0.2, respectively. At B = 0.05, the value of 1/ψ 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 α 2 (β 2 4a − α 2 4 ) first decreases slightly and then increases as m decreases from 1; this is again due to changes in the  Figure 6. 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 vs B between the complete mixing case (m = 1) and optimal mixing case (m = m C P ). Plotted together are m C P and m ψ , which are the values of m to maximise C P and to minimise ψ, respectively. The disc resistance coefficient is fixed at the optimal value for m = 1, i.e.
'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 α 2 (β 2 4a − α 2 4 ) tends to be smaller than the decrease (or increase) rate of 1/ψ 3 , and, hence, C P follows the trend of 1/ψ 3 , i.e. C P slightly increases first and then decreases as m decreases from 1. This initial increase in 1/ψ 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 6(c) compares C P vs 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 GC07 ). Also plotted together are m C P and m ψ , the latter of which represents the value of m that minimises ψ (and, thus, maximises 1/ψ 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 ψ is small, especially at low blockage ratios. This reflects the dominant influence of 1/ψ 3 on C P as described earlier in figures 6(a) and 6(b).
The C P values presented above are for 'suboptimal' discs (with K = K 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 7(a) and 7(b) show C P vs K curves for five different m values, again at B = 0.05 and 0.2, respectively. When B is small, K 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 GC07 is not so close to the optimal K for a given m, but the C P vs K curve tends to have a flatter peak; hence again, the additional gain in C P by adjusting K is small. Figure 7(c) shows the maximum C P (achieved by adjusting K) vs 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 vs m curves for K = K GC07 plotted earlier in figures 6(a) and 6(b). At B = 0.2, for example, the highest C P value achieved with K = K GC07 is only 1.1 % lower than that achieved by optimising K. The impact of m on C P can also be confirmed from the C P vs C T curves in figure 8, plotted for the same sets of m and B as in figures 7(a) and 7(b). In summary, these example solutions of the simple three-parameter (B-K-m) theoretical model suggest the following. (i) For an infinitely large staggered array 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 array of tidal 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 a real tidal array to be built in the future, the cross-sectional average velocity (which was considered as a fixed parameter in the above analysis) would 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). As the additional resistance caused by the array would reduce the mass flow rate through the array, the turbine resistance required to maximise the power for a given 'natural' mass flow rate (observed in the 'natural' channel without the array) would be lower than the optimal resistance discussed above. It should also be noted that, in the real world, tidal currents are oscillatory and, thus, it would be beneficial to vary the turbine resistance during the tidal cycle (Vennell & Adcock 2014).

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 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 Appendix A.
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 1 MW DeepGen IV tidal stream turbine design from the ReDAPT project is adopted, and the details of the blade hydrodynamic data used in 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 Ω = 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 × 12D × 2D, which is close to 6πH × 2πH × 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. The corresponding Froude number (Fr = U 0 / √ gH, with g denoting the gravitational acceleration) is equal to 0.13, which, together with the given rotor size relative to the water depth, is deemed small enough to ignore any water surface deformation in this study (Houlsby & Vogel 2017;Yan et al. 2017). A fixed time step Δt is set to 0.045 s together with a uniform spatial resolution Δ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, Li & Moriarty (2013), Foti et al. (2019) and Yang & Sotiropoulos (2019b). Hence, the domain is divided into 1692 × 576 × 96 grid cells over the three spatial directions with a total of about 93.5 × 10 6 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 streamwise pressure gradient term Π 1 that 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 Π 1 that would vary U 0 as in previous infinitely large wind farm simulations (Calaf, Meneveau & Meyers 2010;Yang, Kang & Sotiropoulos 2012;Stevens, Gayme & Meneveau 2014). Note that, to simulate a real large tidal array, both of the bulk velocity and Π 1 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 rigid-lid 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, Yang & Sotiropoulos (2014). A representation of the computational domain is presented in figure 9 together with instantaneous flow structures visualised using Q-criterion isosurfaces for one of the configurations studied. 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 an 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 (n y ), i.e. B = n y πD 2 /4HW), physical simulated time in terms of flow-through time and the number of computed turbine revolutions. Also presented in this table are the results of the thrust coefficient (C T ) and power coefficient (C P ) obtained from averaging the values for all the turbines simulated in each case, their relative difference (ΔC T , ΔC P ) with the AL-36 × 12 configuration as a reference case (a single turbine is considered), and the resistance coefficient (K) calculated using the velocity at the rotor (U D ). As we adopt the ALM, the value of U D is estimated as the rotor-averaged velocity integrated two grid cells upstream of the turbine rotors and excluding the area occupied by the nacelles. The arrays adopted in the LES are perfectly aligned and staggered layouts, labelled as AL and ST followed by the device spacing, e.g. the case ST-9 × 6 corresponds to a staggered array with S x /D = 9 and S y /D = 6 as in figure 9.
In our simulations the streamwise forces from the turbines (thrust) and the bed surface friction contribute to the overall resistance to the periodic mean flow, which is quantified from the time-averaged streamwise pressure gradient in terms of the friction velocity, i.e. u * = √ Π 1 H, (Yang et al. 2012). In table 1 we report the resulting 'overall' friction coefficient (C f = 2u 2 * /U 2 0 ) as this corresponds to the one often used in shallow water models for representing the flow resistance due to tidal turbines and bottom  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 (FT) times, number of revolutions, hydrodynamic coefficients (C T and C P ) with their variation compared with the reference case AL-36 × 12, resistance coefficient (K), velocity at the rotor (U D /U 0 ), and 'overall' friction coefficient (C f ).
surface stress. Our results indicate that the overall resistance increases approximately linearly with the number of turbines deployed, regardless of whether the turbines are aligned or staggered. In the precursor simulation the friction velocity only accounts for the temporal and spatially averaged bed friction, yielding a value of C f = 0.0009.

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 velocity fields are presented in § 4.1, followed by some turbulence statistics in § 4.2. We then present the wake-centreline 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 · and the instantaneous fluctuation value is represented as (·) obtained following the Reynolds decomposition.

Streamwise velocity field
We present in figure 10 contours of the normalised time-averaged streamwise velocity ( u /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-9 × 3 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 2S 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 10 for the aligned cases. However, in staggered configurations this lateral wake expansion appears limited compared with their aligned counterparts.
The different spreading rates of the time-averaged turbine wakes are partly explained by their instantaneous behaviour. In figure 11 we present contours of instantaneous streamwise velocities (u/U 0 ) for the previous cases shown in figure 10 with S x /D = 9 and changing S y , which reveals a pronounced meandering nature of the wakes with large oscillation amplitudes (Foti et al. 2019). In both aligned and staggered configurations, the wake meandering is influenced by the lateral distance between turbines in each row, i.e. reducing S y leads to a smaller meandering amplitude. This trend seems clearer in the staggered cases, e.g. the wake meandering almost disappears in the staggered cases with S y /D = 4 and S y /D = 3, explaining the narrow time-averaged wakes shown earlier in figure 10. This is presumably because, in the aligned cases, the flow going through the lateral gaps between turbines creates high-speed streaks, resulting in a larger velocity difference between the wake (core) and surrounding (bypass) flows (Stevens et al. 2014). Consequently, there is a stronger instability in the shear layer of turbine wakes that enhances the wake meandering (Yang & Sotiropoulos 2019a) more in aligned arrays compared with their staggered counterparts. It is worth noting that the vertical confinement and shorter lateral spacing (commonly expected for future tidal arrays) both lead to a larger velocity difference in aligned arrays. Hence, tidal turbine wakes in a large array could be more sensitive to their arrangement than wind turbine wakes in a large wind farm.

Turbulence statistics
We further analyse the flow field in figures 12 and 13 with contours of turbulence intensity in the streamwise (σ u /U 0 , with σ u = u u 0.5 ) and transverse directions (σ v /U 0 , with σ v = v v 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 with 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 10 and 11. 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 formation of high-momentum streaks and, thus, meandering motion. As can be expected from the above results, the resulting unsteady meandering wake also notably affects the distribution of turbulent momentum exchange, which is described in figure 14 with contours of the Reynolds shear stress (− u v /U 2 0 ) at hub height. For each array layout, narrowing the lateral spacing reduces the magnitude 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 lower magnitudes of shear stresses compared with the aligned layouts irrespective of the number of turbines deployed.

Centreline velocities
Next, we quantitatively compare the mean streamwise velocity ( u /U 0 ) distribution along the wake centreline in figure 15 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 as a shaded area the location of the intermediate row of turbines. When the streamwise spacing is relatively small (S x /D = 9 and 6), the values of u /U 0 behind the turbine nacelles increase faster for the lower lateral blockage cases, i.e. S y /D = 6 and 12.  Figure 13. Contours of spanwise turbulence intensity (σ v /U 0 ) at hub height for aligned and staggered cases with S x /D = 9 and lateral spacing S y /D = 6 (a,b), 4 (c,d) and 3 (e, f ).
This quicker wake recovery in the lower blockage cases results from a larger entrainment of ambient flow into the wake, which can be enhanced by a stronger meandering behaviour as seen in the distribution of Reynolds shear stresses in figure 14.
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 u /U 0 are slightly over unity shortly upstream of the turbines for configurations with S y /D ≤ 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 u /U 0 of about 0.9 to 0.95 at one rotor-diameter upstream of the turbines (i.e. at x/D = 8 and 17). It is clearly seen that the wakes in the cases with a higher blockage (AL-9 × 4 and AL-9 × 3) 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 cases with S x /D = 6 in which 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.
In staggered arrays with S x /D = 18, the centreline mean velocity notably increases at around x/D ≈ 18 due to the local blockage caused by the turbines in the intermediate second row. Such flow acceleration is more noticeable when the lateral spacing S y /D is small, e.g. for ST-18 × 3 and ST-18 × 4, the maximum velocity is about 1.1U 0 whereas, for ST-18 × 6 or ST-18 × 12, it is approx. 1.02U 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 larger 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 ≥ S x ) the centreline 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 configuration, due to stronger turbulent mixing enhanced by the wake meandering as shown earlier.

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 16 with error bars indicating the standard deviation of array-averaged values. For aligned layouts with a given S y /D, both thrust and power   Figure 15. Distribution of time-averaged streamwise velocities ( u /U 0 ) along the wake centreline at the hub height for aligned (a,b,d, f ) and staggered (c,e,g) configurations with S x /D = 36 (a), 18 (b,c), 9 (d,e) and 6 ( f,g), and S y /D = 12 (black), 6 (blue), 4 (green) and 3 (orange). (e) Figure 16. Results of C T (a-d) and C P (e-h) obtained from the LES for the aligned ( ) and staggered ( ) layouts.
coefficients tend to have their maximum values at the largest streamwise spacing of 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 (Ali et al. 2018), as shown earlier in figure 15. In contrast, the performance of turbines in staggered configurations appears less sensitive to S x /D. Comparing arrays with a 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 attained at 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 from earlier wind farm studies, staggered configurations are in general more efficient than the aligned counterparts (Tian, Ozbay & Hu 2018).

Turbine loading unsteadiness
Next, we analyse the temporal fluctuations of thrust and power from their root-mean-square (r.m.s.) values. The results averaged over all turbines comprising the array are presented in figure 17, which shows that decreasing S x /D leads to an increase in the fluctuations of both thrust and power, as expected from an 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 11. 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.
To further link the fluctuation of turbine forces with the wake meandering, we present in figure 18 the power spectral density of the power signal from a single turbine comparing the same sets of staggered and aligned layouts discussed earlier in figures 10 to 14. The spectra show distinct peaks at the blade passing frequency f p = 3f 0 , where f 0 = Ω/(2π) is the rotational frequency, and its harmonics (6f 0 and 9f 0 ) which indicate turbine-induced high-frequency changes in the power generation (Ouro et al. 2019c). Wake meandering is, however, a low-frequency phenomenon that is identified by the peaks in these spectra at f wm = 0.08f 0 , i.e. its undulating motion has a period of about 12 times that of the turbines' rotation. These peaks are well observed in the aligned cases irrespective of S y , whilst in the staggered cases the signal at those low frequencies reduces with lower S y and, for S y /D = 4 and 3, the peak due to the wake meandering disappears. This agrees well with our earlier observation of the disappearance of wake meandering from the instantaneous velocity fields presented in figure 11.
The power spectra from the aligned cases appear to follow a −7/3 decay after the peak at f wm irrespective of S y , and this is also observed for the staggered array with S y /D = 6. Narrowing the lateral spacing for the ST-9 × 4 and ST-9 × 3 cases leads to three major changes: a progressive transition in the decay slope from −7/3 to −5/3, a decrease in the energy associated to the low-frequency (inertial) range, and the absence of a distinct peak at f wm . Earlier studies with tidal turbines reported a slope of −5/3 from the spectrum of the power output signal (e.g. Deskos et al. 2020) resulting from the ambient turbulence governing the power fluctuation. However, our power spectra in figure 18 suggest that low-frequency changes in the power output of turbines in a large tidal array could be strongly affected by the wake meandering, depending on the turbine spacing and array layout. The Strouhal number (St = f wm D/U 0 ) associated with the wake meandering frequency is found to be approximately 0.11, which is similar to the values reported for wind turbines, e.g. 0.15 obtained by Foti et al. (2016).
The observed variability in the decay slope for different array layouts is presumably because the dynamics of flow through a large tidal array is mostly driven by the turbine-induced changes (unlike wind farms, where low-frequency dynamics are also impacted by large-scale flow structures in the atmospheric boundary layer; see, e.g. Stevens, Gayme & Meneveau 2016;Bossuyt et al. 2017). Further investigations into these spectra for cases with irregular turbine layouts and lower vertical confinements (larger water depths) would be helpful in future studies.

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 ∝ 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 ∝ S x ). It should be borne in mind that the theoretical model allows us to predict overall trends and upper bound estimates of the array efficiency with almost negligible computational cost, whereas the LES-ALM allows us to resolve the details of a 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 19. 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-36 × 12), 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 19 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 ≤ 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, e.g. using LES-computed wake-centreline velocities for a wider range of array configurations. 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 overpredicted by the theoretical model compared with 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).

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 LMADT, and LES with an ALM (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 ≤ S x /D ≤ 36) and lateral (3 ≤ S y /D ≤ 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, whilst its frequency is found to be independent of the spacing adopted. The main consequence of this change in wake dynamics is that a lower amplitude of wake meandering means there is 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 centreline 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 us 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. The power spectral density distribution of the turbine power signal is found to follow a −7/3 slope for the aligned arrays with a peak at the wake meandering frequency, whilst this progressively changes to a −5/3 slope in the staggered arrays as the lateral spacing is reduced and the wake meandering effects disappear.
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 925 A30-27 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 because the local flow velocity upstream of each turbine can exceed the bulk velocity only when the streamwise spacing is reasonably (not excessively) large and, thus, the mixing is largely (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.
where u i = (u, v, w) T is the vector of spatially filtered velocities, the coordinates vector is x i = (x, y, z) T , ρ denotes the fluid density, p is the relative pressure, ν is the kinematic viscosity of the fluid, and f t is a source term resulting from the ALM and immersed boundary forcing used for the representation of turbine rotors and nacelles, respectively. Here Π 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 ν t is calculated using the wall-adapting local eddy-viscosity 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, McSherry & Stoesser 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;, geophysical flows (Ouro et al. 2017b), rough open-channel flows (Stoesser 2010;Bomminayuni & Stoesser 2011;Nikora et al. 2019) and fluid-structure interaction (Kara, Stoesser & McSherry 2015;Ouro, Muhawenimana & Wilson 2019b). In the present work the IBM is adopted to represent the turbine nacelles, using the φ 4 delta function for the interpolation procedures. Turbine rotors are represented using an actuator line model 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, Ω, 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 torque Q, and, thus, determine the generated power P = QΩ. Eventually, the coefficients of thrust (C T ) and power (C P ) are computed as After the computation of the 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, Sørensen & Mikkelsen 2005), where r L denotes the radial distance between the marker L and the considered cell face i, and ε is the interpolation stencil set to 3.0Δx i . A Prandtl-type tip-loss correction is adopted