## Abstract

At tidal energy sites, large arrays of hundreds of turbines will be required to generate economically significant amounts of energy. Owing to wake effects within the array, the placement of turbines within will be vital to capturing the maximum energy from the resource. This study presents preliminary results using Gerris, an adaptive mesh flow solver, to investigate the flow through four different arrays of 15 turbines each. The goal is to optimize the position of turbines within an array in an idealized channel. The turbines are represented as areas of increased bottom friction in an adaptive mesh model so that the flow and power capture in tidally reversing flow through large arrays can be studied. The effect of oscillating tides is studied, with interesting dynamics generated as the tidal current reverses direction, forcing turbulent flow through the array. The energy removed from the flow by each of the four arrays is compared over a tidal cycle. A staggered array is found to extract 54 per cent more energy than a non-staggered array. Furthermore, an array positioned to one side of the channel is found to remove a similar amount of energy compared with an array in the centre of the channel.

## 1. Introduction

To capture significant energy, tidal turbines need to be built in large arrays within regions such as tidal channels and straits where the local bathymetry focuses the flow. The layout of these arrays of turbines can potentially significantly change the amount of energy captured from the flow. This study presents preliminary results of work using numerical models of different array layouts to predict the flow field and the energy capture by each layout.

Previously, Harrison *et al*. [1], Lee *et al*. [2], Batten & Bahaj [3] and MacLeod *et al*. [4] have modelled flow around arrays of tidal turbines using some form of the Reynolds-averaged Navier–Stokes equations with the turbine represented as some form of actuator disc or frozen rotor. As they used the Reynolds-averaged approach, they focused on average quantities rather than the time-dependent dynamical structure of the vortex wake which this study focuses on. In this study, the time-dependent structure is modelled using a large eddy simulation method, similar to that used by McCombes *et al*. [5] who modelled the flow past a single tidal turbine using a vorticity–velocity formulation of the Navier–Stokes equations with no sub-grid-scale turbulence model. In McCombes *et al*.’s [5] model, the artificial diffusive mixing acts as viscosity.

Harrison *et al*.’s [1,6] models of tidal turbines in unsteady flow show that ambient turbulence levels in the free stream greatly affect the wake structure behind a single turbine and in multiple turbine arrays.

An alternative approach is actuator disc theory, a simple analytical model of a turbine where the turbine is represented by an infinitely thin actuator disc which removes along stream velocity from the flow [7]. The expansion of a stream tube of fluid as it passes through the disc is modelled. This model predicts a step drop in the pressure across the disc and a gradual decrease in the velocity upstream and downstream of the disc with half the along stream velocity loss occurring upstream of the disc and half downstream. Actuator disc theory also leads to the Lanchester–Betz limit for the optimal thrust coefficient for a turbine of *C*_{T}=0.89 and optimal power coefficient of *C*_{P}=0.59 [8,9].

Vennell’s [10] one-dimensional analytical model for an array of tidal turbines in a channel shows that to get the most energy from a channel, each row must block as much of the channel as possible. The model is extended in the study of Vennell [11] to show that each turbine in an array must be tuned in concert to have the optimum drag coefficient within the scope of the whole array.

Draper *et al*. [12] represented an array of turbines in a channel as regions of increased, depth-averaged bottom friction in a shallow water equation solver. They found good agreement between their numerical model and their one-dimensional analytical model. Similarly, Blunden & Bahaj [13] present an analytical model where turbines are represented by distributed roughness elements in a boundary layer comparable to flow through vegetation in a water column of finite depth. This study uses a similar method of representing the turbines applied to different array layouts.

Actuator disc and distributed roughness approaches have been used to numerically model array effects in wind farms [14,15] but there are some significant differences in the fluid dynamics between tidal energy and wind energy. In a tidal channel, the free surface varies which is the driving force for the current and hence energy extraction. Also, tidal turbines are likely to be roughly 20 m in diameter [16] and positioned in channels up to 50 m deep. This means that the free water surface will be at most 1.5 diameters above the top of the turbine, allowing very little undisturbed flow above the turbine to mix the turbulent wake into the free stream flow compared with wind turbine arrays where there are kilometres of atmosphere above the turbines.

## 2. Method: model set-up

Here, we used Gerris to solve the linear shallow water equations around four different possible layouts for turbine arrays. The total power captured by each array was compared in order to find the array that captured the most energy over a tidal cycle.

Gerris is an adaptive mesh, finite-volume flow solver [17] which has been used to solve the shallow water equations with a linear free surface approximation. By adapting the mesh in space and time, the solver is able to represent the wide range of scales present in the model, from the kilometre-wide channel to metre-wide vortices in the wake structure and adapt to the vorticity as the wake travels downstream. This adaptivity is a key feature in modelling the individual wakes within an array of turbines, as it reduces the computation time to roughly 10 per cent of that of a simulation with a fixed mesh.

### (a) Adaptive grid

An initial grid was specified with enough resolution near the centre to resolve the turbines and a large enough cell size near the boundaries to produce smooth input flow and satisfy the CFL criteria for a maximum time step of 5 s, d*L*=d*t*(*gH*)^{0.5}=111 m as shown in figure 1. The grid adapts based on vorticity level at each later time step as shown in figure 1 so that regions with higher vorticity have a smaller cell size down to a minimum specified size. When vorticity reduces at a later time, the mesh adapts to a larger cell size. Cells near the boundaries are kept large to reduce the impact of high vorticity on the boundary forcing.

Simulations were run over a range of minimum cell size from 4 to 0.5 m to test for convergence of the power extracted from the flow by one turbine. The energy starts to converge as the minimum cell size reaches 0.5 m as shown in figure 2. This convergence is evidence that the simulation is consistent. However, the simulations with coarser mesh show much lower vorticity in the wake owing to increased numerical diffusion. A minimum size of 1 m was chosen for the following simulations as a good compromise of simulation length and complexity.

The time step also adapts based on vorticity so that the time step reduces when there is higher total vorticity in the simulation.

### (b) Tidal channel is 3 km long by 1 km wide by 50 m deep

These simulations are of an array of turbines in a tidal channel. The channel is 3 km long, 1 km wide and 50 m deep. A uniform bottom friction of 2.5×10^{−3} is applied across the domain. A sinusoindally varying flow of 2 m s^{−1} amplitude is driven by applying a cosine function for the free surface elevation that is 180^{°} out of phase at each end using a Flather boundary condition [18].

### (c) Boundary conditions establish a pressure gradient

Flather boundary conditions are a form of radiative boundary condition where an expected free surface elevation and velocity are prescribed and any difference between the expected and calculated free surface is allowed to propagate out of the domain at the long wave phase speed. Flather boundary conditions are less prone to reflections into the domain than clamped boundary conditions [18].

The fluid is initially at rest with a high free surface on the left-hand side of the domain and a low free surface on the right-hand side so that the fluid slowly starts moving left to right then reverses to flow right to left. This results in a smooth sinusoidally varying flow with minimal reflections from the upstream and downstream boundaries that dissipate in the first hundred seconds of the simulation. The boundary conditions on the side walls are zero normal derivatives for pressure and zero across stream flow representing a vertical-sided channel with no friction on the impermeable side walls.

### (d) No ambient turbulence or sub-grid-scale turbulence

As a preliminary step, this study does not include any representation of ambient turbulence in the channel. However, a turbulent wake does flow back through the array in the second half of the tidal period, owing to the periodic nature of the oscillating flow. Thus, the first half of the simulation represents flow in a uniform, non-turbulent flow, whereas the second half demonstrates the effect of a turbulent wake flowing back through an array, still with no ambient turbulence.

There is no sub-grid turbulence model; instead Gerris relies on numerical truncation error terms in the advection scheme to dampen the unresolved turbulence at scales below the minimum mesh size. Turbulence is also minimally damped by a background bottom friction coefficient of 0.0025, applied uniformly in the domain. Hence, while turbulent structures are dissipated numerically and by the bottom friction, the structures potentially survive longer than they would in a scheme with a sub-grid-scale turbulence closure scheme.

### (e) Reduced tidal period

The period of tidal oscillation in these simulations is 1.24 h. This is 10 per cent of the M2 tidal period and was chosen to give easy scalability to an M2 tide while keeping the total simulation time for a full tidal cycle reasonable and still covering the full range of tidal velocities expected. This reduced period produces a shorter tidal excursion than a full tidal period would, so the vortex structures associated with peak flow travel through the farm earlier in the phase and with less dissipation than in a full tidal cycle.

### (f) Turbines are rectangles of increased bottom friction

The turbines are represented by 20 m across stream×5 m along stream rectangular boxes of increased depth-averaged friction, *C*_{d}. This drag coefficient represents the total drag owing to the turbine and the structure on the flow.

To find the optimum drag coefficient, simulations were run on an individual turbine for a quarter of a tidal cycle, capturing the full range of velocity magnitude, over a range of drag coefficients from 0.05 to 50. The maximum energy extracted from the flow over a quarter tidal cycle was for a turbine with *C*_{d}=12 so this drag coefficient was used for all subsequent simulations.

This peak is not easy to compare with the Lanchester–Betz peak power coefficient of 0.59 as it differs conceptually in two ways. The drag coefficient used in this study is a depth-averaged value, averaged across the water column above, through and below the turbine so it represents an average of regions of no drag and the region of the turbine. Also, the force applied by this friction term at a certain cell is calculated using the velocity in that cell as opposed to the Lanchester–Betz drag force calculation which uses a free stream velocity upstream of the turbine.

The depth-averaged drag coefficient per unit width that is used in this study (*C*_{d}=12) is divided by the depth (50 m) in the model to get a depth-averaged drag per unit width. If the turbine is assumed to occupy 20 m of the 50 m water column, then the drag per unit area of an equivalent three-dimensional turbine is (12/50)×(50/20)=0.6.

The upstream peak velocity that would be used for a Lanchester–Betz type calculation is 2 m s^{−1}. The average velocity across cells inside the turbine at the same peak flow is 1.1 m s^{−1}. So the power coefficient at peak flow in the Lanchester–Betz model is lower than the power coefficient to extract the same power in this simulation by a factor of (1.1/2)^{3}=0.17. However, the drag coefficient used in this study is not the maximum for peak flow. The drag coefficient in this study maximizes the power output over a quarter of a tidal cycle.

### (g) Single turbine wake shown for comparison

Initially, simulations were run using a single turbine in the centre of the channel. The pressure and velocity fields through the centreline of the turbine are investigated and compared with actuator disc theory. Velocity and vorticity fields around the turbine are plotted for comparison with other techniques.

### (h) Four different array layouts compared

Subsequently, different configurations of 15-turbine arrays were simulated to compare the power from four different array configurations: centred, offset, larger spacing and staggered. Each array consists of 15 turbines configured in five rows of turbines, each row containing three equally spaced turbines. Rows are 10 diameters apart and turbines in each row are 7.5 diameters away from each other.

The centred array is centred in the middle of the channel and each row of turbines is directly downstream of the upstream row as shown in figure 4*b*. The offset array is the same arrangement but positioned to the side of the channel to represent an arrangement that allows shipping traffic to pass on one side of the turbine array as shown if figure 4*c*. The larger array is the same as the centred array but rows are 20 diameters apart, twice the spacing of the centred array as shown in figure 6*a*. The staggered array is also centred in the channel with the same row spacing as the centred array but the second and fourth rows are staggered so that the turbines are halfway between the turbines in the adjacent rows and the bottom turbine is below the bottom turbine in the adjacent rows as shown in figure 6*b*.

The simulations are run over one tidal cycle and the total instantaneous power removed from the flow by each array is calculated using equation (2.1). This power is summed over a tidal cycle to give the total energy captured by each array. This energy and power are normalized by the power and energy in the centred array as this study only focuses on the relative energy captured by different arrays. 2.1

## 3. Results

The preliminary results in figures 3–6 depict the flow moving left to right, in the first half of the tidal cycle where there is no turbulence in the upstream flow. The structure of the wakes is more easily identifiable in the non-turbulent flow than in figure 7 where the flow has reversed and flows back through the array. However, the second half of the tidal cycle is closer to realistic flow in a turbulent channel, with wake turbulence from the first half of the tidal period moving back through the array. There is still no ambient background turbulence, only the structure due to the wakes. The extra mixing due to the wake turbulence flowing back through the array shortens the shear region of each turbine’s wake as shown in figure 7, and increases the total power captured by the array relative to the first phase of the tidal period as shown in figure 8.

### (a) Single turbine wake

The flow around a single turbine is presented in figures 3, 4*a* and 5 for comparison with other predictions of wake effects at the array scale rather than as a prediction of the wake expected from a single turbine in isolation. This wake structure does not capture the full dynamics seen in single turbine wake models such as Draper *et al*. [12] or in Myers & Bahaj’s [19] tank studies but does show similar structure at the array scale.

The wake behind a single turbine in this representation forms three distinct regions as seen in figures 3 and 4*a*. First, there is a shear region immediately downstream of the turbine where the slower flow directly behind the turbine is separate from the faster, free stream flow to the sides of the wake. The two flows are separated by bands of high vorticity corresponding to the two sides of the turbine. Downstream of the shear region is a separated region of high localized vorticity where the clear lines of shear mix into vortical structures of a Von Karman vortex street. Further downstream, the vorticity is more mixed into the free stream flow where strong vorticity structure persists but vortices are more diffuse and the wake is wider. Some vortex pairs are present in this mixed region, meandering up to a few turbine diameters away from the centreline of the wake. In figure 3*a*, vorticity is eventually dissipated numerically by a region of larger cell size before hitting the downstream boundary.

Along stream velocity and free surface elevation along the centreline of the turbine at maximum velocity, *t*=*T*/4, are plotted in figure 5. The small region shown is all within the shear region shown in figure 4*a*. Further downstream, the flow is more disturbed, with sharp changes in pressure and velocity matching the regions of high vorticity in figure 4*a*. The along stream velocity slows gradually across a region upstream and downstream of the turbine, with more slowing occurring across the region closest to the turbine which is consistent with actuator disc theory. There is a small region of deviation from this gross pattern within the 5 m of the turbine itself where flow is slowed more than in the infinitely thin plane of an actuator disc. The free surface elevation is also consistent with actuator disc theory, showing a steady increase in pressure upstream of the turbine, a sharp drop across the region of the turbine and a steady increase downstream of the turbine, recovering to a level a little lower than upstream.

### (b) Centred array

Figure 4*b* shows the array layout (black boxes) and vorticity field around the centred array at maximum tidal flow, *t*=*T*/4. The wakes in the first row of turbines look similar to the single turbine case until the second row. The second row of turbines is in the separated region of the first row’s wake so the power available to the second row is much reduced, relative to the first row, and the turbines will be subject to high shear loading on the blades which could damage the turbines. The wake from the second and subsequent rows is quite different from the first row as it lacks the early shear region. The flow that hits the downstream turbines is sufficiently disturbed that the wakes form straight into separated regions with strong vortexes immediately behind the turbines. The last three rows are in an even more mixed region of the wake where the region of vorticity structure is wider and more diffuse but they are still subject to strong gradients and reduced average flow which results in lower power available to the turbine and strong shear forces. The total energy extracted from the flow over the tidal cycle is 65 per cent of the energy removed by the staggered array as seen in figure 8 where this centred array has the lowest power of all four configurations, so this is clearly not an optimal layout for an array of tidal turbines but illustrates the effect of placing turbines in lines.

### (c) Offset array

Figure 4*c* shows the array layout (black boxes) and vorticity field around the offset array at maximum tidal flow of 2 m s^{−1}. The wake structure is very similar to that of the centred array, but is constrained a little by being near the wall and spreads a little further laterally into the undisturbed flow than the wake for the centred array. The total energy removed from the flow over a tidal cycle is 4 per cent more than the energy removed by the centred array. Figure 8 shows that the extra power is extracted from the flow on the second phase of the tide, relative to the centred array. This difference in energy is small and may not be significant; however, it does suggest that positioning an array near the side of a channel to allow shipping traffic to travel on one side does not change the power removed from the flow by the array when the flow is uniform across the channel.

### (d) Larger spacing

Figure 6*a* shows the array layout (black boxes) and vorticity field around the larger-spaced array at maximum tidal flow of 2 m s^{−1}. The second row of turbines is in the separated region of the first row’s wake so there is more power available to the second and subsequent rows than in the centred array. Hence the total energy removed from the flow by this array is 31 per cent more than that removed by the centred array. However, the downstream turbines are still subject to large shear forces in the turbulent wake of upstream turbines and the power available to downstream turbines is not as high as that for the first row.

### (e) Staggered array

Figure 6*b* and electronic supplementary material, movie S1, show the array layout (black boxes) and vorticity field at maximum tidal flow of 2 m s^{−1} for the staggered array where the second and third rows are staggered so that the turbines lie between the turbines in the other rows. This array removes the most energy of all the arrays tested, 54 per cent more energy than the centred array.

The second row of turbines here sits in a region of flow that has accelerated around the first row of turbines and actually removes more energy from the flow than the first row of turbines do as shown in table 1. The third row is spaced the same distance downstream from the first row as each row is spaced in the array with double spacing between rows. However, in the staggered array, the third row possibly also benefits from the same squeezed flow effect as the second row experiences from the first row. This is not easily discernible in table 1 because the power from rows three, four and five is reduced significantly by being in the wake of the first half of the array. By the fourth and fifth rows, the wake from the upstream turbines in the offset row is mixing with the wake from two rows upstream so the vorticity field is more mixed and less distinctly in three lines. This is clearly the best arrangement of the four arrays tested by a significant fraction.

### (f) Power comparison

Figure 8 shows the instantaneous power at each point in the tidal cycle for each of the arrays, with a 200 s moving window filter applied. The filter is applied to highlight the difference between average power for different array layouts while still showing some of the variation in time. Interestingly, more power is removed from the flow on the second phase of the tide than during the first phase.

This difference in power in each phase is due to the difference in incident flow in each phase and to two effects that result from the difference in flow. In the first half of the tidal cycle, the flow upstream of the array is linear and smooth compared with the incident flow in the second half of the period when the turbulent wake moves back through the farm on the reversing tide as shown in figure 7. The two effects that result from the difference in upstream turbulence are: firstly, the shorter wake recovery in the higher turbulence increases power available to downstream turbines in the second phase; secondly, the cubic relationship between velocity and power means that the turbulent flow field that moves back through the arrays in the second half of the tide contains more energy than the smooth upstream flow that hits the array in the first phase.

Figure 7 shows the vorticity field a little after the tide has reversed and shows the regions of high vorticity moving back through the array. In practice, the second half of the tidal cycle is more representative of real flow where the disturbed wake will slosh back and forward through the array on each changing tide. Figures 7 and 8 also show a limitation of the reduced tidal period which produces a shorter tidal excursion than a full tidal period would. Consequently, the vortex structures associated with peak flow are less dissipated than they would be in reality when they flow back through the farm.

## 4. Conclusions

Four different array arrangements for tidal turbines in a tidal channel have been modelled. The best arrangement of these is a staggered array where alternating rows of turbines are positioned directly between the upstream turbines. This arrangement captures 54 per cent more energy than a grid arrangement of turbines. Also, very little difference has been found between an array in the middle of a channel compared with an array positioned to one side of the channel. Doubling the spacing between rows of turbines increases the energy removed from the flow by 31 per cent but not as much as staggering alternate rows.

While representing a sophisticated three-dimensional turbine as a simple rectangle of increased drag friction will not capture the rotating helical wake structure of tip vortices or the full dynamics of a turbine, it does capture the wake behaviour predicted by actuator disc theory at the array scale. Further validation of the single turbine wake against experimental results is planned.

This representation of turbines as rectangular friction elements in the shallow water equations is limited and does not show the same flow structure as a more accurate representation of a single turbine, such as [12] or [19]. However, the goal of this paper is to represent the turbines and wakes in the simplest realistic way at an array scale so that the flow and power capture in tidally reversing flow through large arrays can be studied. As such, these initial results show that a simple representation of a turbine array in an adaptive mesh model shows promise. With some refinement, this technique should provide a useful method to further investigate flow fields through turbine arrays.

Further studies should look at more than one tidal cycle as the first half of simulation shows a different power capture from the second half of the tidal cycle.

The effect of background turbulence from complex bathymetry has been ignored here. The difference between the power in the first half of the tidal cycle where the flow upstream of the arrays is laminar and smooth and the power in the second half of the tidal cycle, where the flow upstream of the turbine is highly turbulent, demonstrates that potentially more power will be available in a turbulent flow than a smooth flow. However, this will come at the expense of higher shear stresses on the turbines. Furthermore, tidal streams are non-uniform turbulent, so the power and flow behaviour in the first phase of this study’s tide is only useful for comparison with other results in non-turbulent, uniform flow. The second phase of the tide is more realistic as the flow passing through the turbines is turbulent. Further work should attempt to include ambient turbulence to make the model more representative of real channel flows.

This study has only investigated the effect of a tidal current consisting of a simple sinusoidal tide. Complex tidal flow composed of a mixture of higher harmonics and other tidal constituents may change this result.

## Footnotes

One contribution of 14 to a Theme Issue ‘New research in tidal current energy’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.