## Abstract

This paper presents our initial work in performing large-eddy simulations of tidal turbine array flows. First, a horizontally periodic precursor simulation is performed to create turbulent flow data. Then those data are used as inflow into a tidal turbine array two rows deep and infinitely wide. The turbines are modelled using rotating actuator lines, and the finite-volume method is used to solve the governing equations. In studying the wakes created by the turbines, we observed that the vertical shear of the inflow combined with wake rotation causes lateral wake asymmetry. Also, various turbine configurations are simulated, and the total power production relative to isolated turbines is examined. We found that staggering consecutive rows of turbines in the simulated configurations allows the greatest efficiency using the least downstream row spacing. Counter-rotating consecutive downstream turbines in a non-staggered array shows a small benefit. This work has identified areas for improvement. For example, using a larger precursor domain would better capture elongated turbulent structures, and including salinity and temperature equations would account for density stratification and its effect on turbulence. Additionally, the wall shear stress modelling could be improved, and more array configurations could be examined.

## 1. Introduction

Interest is growing in the large-scale deployment of marine hydrokinetic (MHK) energy conversion devices. One such device is the horizontal-axis MHK turbine that harnesses energy from tidal currents, ocean currents or in-stream river currents. In this paper, when turbine is mentioned, we refer to the horizontal-axis MHK turbine, unless it is specified otherwise. At present, companies such as Verdant Power and Marine Current Turbines (MCT) are well on their way to developing commercially viable turbines. For example, Verdant has deployed six turbines, with rotors that are 5 m in diameter (*D*), in an array as a demonstration project in the East River in New York City (RITE project) [1], and MCT has deployed a 1.2 MW dual-rotor system in Strangford Lough in Northern Ireland (SeaGen project) [2].

To make cost-efficient use of tidal regions, turbines will be arranged in arrays of multiple turbines [3]. These devices are much like horizontal-axis wind turbines, creating wakes as a by-product of energy extraction. If a downstream turbine is within the wake of an upstream turbine, the energy production of that turbine will be affected, and it will probably encounter increased turbulence intensity. Researchers and engineers can use analytical methods to obtain computationally inexpensive yet reasonable estimates of energy production [4,5] and of turbine performance characteristics, including the effects of wake interaction [6,7]. However, high-fidelity methods to analyse ambient turbulence and wake effects are just starting to be used. Hence, there is a limited understanding of how these wakes affect the performance and mechanical loading of downstream turbines. It also is unclear how wakes interact with the turbulence in the ambient tidal flow. Among other things, turbulence in the tidal flow is affected by the roughness of the lower surface of the tidal channel, and by density gradients in the flow caused by temperature and salinity stratification, which change throughout a tidal cycle. An array of turbines can be designed to take advantage of the bi-directionality of the ebb and flow of tides. There is a predominant streamwise flow direction, so the turbines can be closely spaced in the cross-stream direction; the low-energy wakes are unlikely to propagate in that direction.

A necessary step to better understanding the flow physics is to take wake and inflow measurements in an array of turbines; however, this is a costly endeavour. Computational fluid dynamics (CFD) tools complement such measurements. CFD can give insight into the physics of the flow and highlight phenomena that should be verified by experimental measurement. CFD also gives a complete picture of the flow field. Using a combination of CFD and experimental measurements will lead to an increased understanding of tidal turbine array and tidal channel physics. This knowledge is necessary to create lower-order engineering tools that developers can use to plan efficient tidal turbine arrays.

The focus of this paper is to present our initial work, using CFD, to simulate flow through turbine arrays. Our current simulations are not meant to capture all of the physics of tidal flows—this is a work in progress. The effects of density stratification, tidal channel topography, tidal-cycle flow variation or free-surface effects are not included. We have simplified the tidal channel flow simulation by making it periodic in the cross-stream direction during the turbine array simulations and periodic in both horizontal directions during the precursor turbulence spin-up simulation. Furthermore, the flow is approximated to have constant density, and is driven by a body force to a time-invariant, turbine hub-height average speed. Rotating actuator lines impose forces on the flow to mimic the effect of the turbines. The type of CFD used is large-eddy simulation (LES), in which the larger turbulent scales are resolved rather than modelled. This work is meant to develop a framework for generating inflow tidal channel turbulence data, feeding those into a tidal turbine array simulation, and modelling the effects of the turbines on the flow and on each other. Further detail will be incrementally added in the future. In setting up this framework, we studied the effect of array layout on overall power production and the characteristics of the turbine wakes.

To our knowledge, there are no other published works in which LES is used to simulate tidal turbines. However, a number of researchers have used CFD, based on the Reynolds-averaged Navier–Stokes (RANS) equations, for this purpose. Others have used LES to simulate the tidal flow without turbines. We also are not aware of any tidal turbine work in which an actuator line model for the turbines has been used, although that model is used for wind turbine work.

MacLeod *et al.* [8] use RANS-based CFD, in which the turbine is modelled as a disc through which there is a pressure jump. That group focused on the wakes generated by one turbine or tandem turbines. They explored the effect of ambient turbulence on wake recovery and found that, as expected, the wake recovers more rapidly with downstream distance as the ambient turbulence level is increased.

Gant & Stallard [9] used steady and unsteady RANS-based CFD to study a single turbine in a tidal channel flow. Like Macleod *et al.* [8], Gant and Stallard model the turbine as a disc through which pressure decreases. This study attempted to provide unsteady turbulent inflow generated by imposing fluctuations onto a power-law mean velocity profile. The fluctuations were derived from either a von Kármán spectral approach or a ‘synthetic eddy method.’ Gant and Stallard [9] found that the unsteady flow imposed at the upstream boundary decays with downstream distance, as is to be expected with RANS. Nonetheless, they found that wake recovery with downstream distance was enhanced by using the unsteady inflow. This result highlights the point that proper inflow modelling and characterization are important in accurately simulating tidal turbine arrays.

Sun [10] used RANS-Based CFD to study near-wake fluctuations and compared them with experimental measurements. Additionally, Sun *et al.* [11] used this CFD method to examine the flow past a small turbine in a small channel. They modelled the turbine as a porous disc and the upper surface of the channel as a free surface using the volume-of-fluid method. They performed both two- and three-dimensional calculations to show that the turbine has a considerable effect on the free-surface elevation, which, in turn, causes a local acceleration of fluid above the turbine. This free-surface disturbance, they argue, has a significant effect on the way in which the turbine wake develops with downstream distance.

Bai *et al.* [12] used RANS with the Fluent CFD solver to investigate the effects of different array layouts. Harrison *et al.* [13] also studied array effects by using an actuator disc method to study the far wake velocities and turbulence of the array, and those researchers showed agreement with an actuator disc experimental test.

With the exception of Gant & Stallard [9], these previous studies focus more on the effects of placing turbines in the flow than on the turbulent inflow. Li *et al.* [14] examined the inflow by performing LES to simulate the tidal-cycle variation of tidal flow, without the influence of turbines. Their work includes the effects of salinity and temperature stratification and a cyclically varying force that drives the flow. It shows how the mean velocity, salinity and the nature of the turbulence vary during a tidal cycle. These effects are important in understanding the performance of tidal turbine arrays and how turbine wakes develop, and they will be included in our future simulations. Lower-order tools used to plan tidal turbine arrays should be redesigned or modified to incorporate such effects, also.

The effects of tidal channel topography could be important in simulating tidal turbine arrays, but we are unaware of research that explores this in the context of array simulations. Researchers are, however, using various types of CFD to study how channel topography, in the absence of turbines, affects the flow. Examples include the work of Zhao [15], in which LES was used to compute the flow through a tidal channel on the central California coast, and that of Booij [16], in which RANS and LES are compared when computing flow in curved tidal channels and rivers.

## 2. Large-eddy simulation methodology

LES is the CFD method used in this study. The larger turbulent scales are directly resolved by solving the spatially filtered Navier–Stokes equations, whereas the effects of the remaining, more isotropic, smaller scales are modelled with a subfilter-scale (SFS) turbulence model. LES provides more flow physics detail, and it places less reliance on turbulence modelling, than the more computationally efficient RANS type of CFD, in which all turbulent scales are modelled. Unlike direct numerical simulation (DNS), in which all turbulent scales are resolved, LES remains computationally tractable at the high Reynolds number of this flow. The larger scales are not resolved near the lower tidal surface; instead, we use a wall model.

The turbine geometry is not resolved (i.e. a grid is not fitted around the blades and the support structure). Rather, we use actuator lines in which the forces created by the turbine are modelled and applied to the flow field as a body force.

The simulations are performed in two stages. In the first stage, the precursor simulation, the turbulence is allowed to develop without the influence of the turbines. The tidal flow is approximated as infinite in both of the horizontal directions through the use of periodic boundary conditions. As the simulation progresses in time, the turbulence reaches a quasi-steady state. Then, data are taken from a streamwise boundary plane and saved at each time step over some time interval. The data are used as the inflow boundary condition for the second stage of the simulations, in which a turbine array is introduced. Streamwise periodicity is no longer assumed. Turbulent inflow data generated during the precursor stage enter the domain, are perturbed by the turbines, and then the flow containing the turbine wakes exits the domain without cycling back through.

### (a) Resolved-scale governing equations

With LES, the incompressible Navier–Stokes equations are spatially filtered to arrive at the resolved-scale (large-eddy scale) equations of motion. The filtered continuity equation is
2.1where the overbar denotes filtering, and is the resolved-scale velocity vector, which is the instantaneous velocity vector, *u*_{j}, minus the SFS velocity vector, *u*′_{j}. (Please note that a finite-volume formulation is used and a filter is not explicitly applied, as is often done with pseudo-spectral calculations. Rather, the finite-volume formulation itself acts as a box filter.) The filtered momentum equation is
2.2where is a modified pressure consisting of the resolved pressure, normalized by the constant density, *ρ*, plus a contribution from the normal SFS stresses, *R*_{kk}. The tensor is the deviatoric part of the full density-normalized SFS stress tensor, *R*_{ij}, which arises from spatially filtering the nonlinear convective term; and and are the density-normalized forces that drive the tidal flow and are imposed by the actuator line representation of the turbine, respectively. In reality, is a cyclically varying force that causes tides to ebb and flow. In our simulations, adjusts with time to drive the flow to a time-invariant mean value of 1.9 m s^{−1} at the turbine hub height. Because this is a high-Reynolds-number flow, the viscous stresses only become significant near the wall. However, we use a wall model that includes the total viscous and SFS stresses in that region. Therefore, equation (2.2) does not explicitly include a viscous term. It is important to note that we neglected the effects of buoyancy, which can suppress or enhance turbulence production, and we also neglected the effects of Coriolis forces caused by the Earth's rotation. Buoyancy effects arise from density variations caused by both temperature and salinity stratification.

### (b) Subfilter-scale model

The deviatoric part of the SFS stress, the quantity upon which the divergence operator is acting in the second term of the right-hand side of equation (2.2), is defined as
2.3where *ν*^{SFS} is the SFS viscosity, and
2.4is the resolved strain-rate tensor. The SFS viscosity in equation (2.3) is found using the Smagorinsky model [17]
2.5where *C*_{s} is the model constant set to 0.125, and Δ is the filter length scale, which is defined as Δ=(*ΔxΔyΔz*)^{1/3}, where *Δx*,*Δy* and *Δz* are the grid spacings in the streamwise, cross-stream and vertical directions, respectively.

### (c) Boundary/initial conditions and computational domain

In this study, we approximate a tidal flow through a wide channel. For example, channels within Washington's Puget Sound have been identified as potentially excellent sites for deployment of MHK energy devices [18]. Some of those channels are a few kilometres wide and tens of metres deep.

Because we are simulating a large channel, in the precursor stage of the simulation when inflow turbulence is developed without the influence of the turbines, the flow is assumed to be horizontally infinite and of uniform depth through the use of periodic boundary conditions on velocity, pressure and SFS stresses in the streamwise and cross-stream directions. The simulation is initialized with divergence-free organized perturbations near the lower surface superimposed upon a logarithmic mean velocity profile.

In the second stage of the simulation, we simulate the tidal turbine array. Periodic boundary conditions on velocity, pressure and SFS stresses are applied in the cross-stream direction. Space- and time-varying velocity data saved from the precursor simulation are used as an upstream boundary condition. On the downstream boundary, the normal gradient of velocity is zero, and the resulting velocity flux through that boundary is adjusted to maintain global continuity. SFS stresses on the upstream and downstream boundaries are based on velocity gradients computed using one-sided differences. The gradient of pressure normal to the upstream and downstream boundaries is zero. The array simulations are initialized with the precursor flow field that was saved when we began to collect precursor boundary data.

The upper boundary remains the same during the precursor and array simulations. It is approximated as an impenetrable, no-stress lid, instead of simulating a free surface. The velocity normal to the boundary, the gradient of velocity parallel to the boundary, all components of the SFS stress tensor and the gradient of pressure normal to the boundary are all zero.

The lower boundary remains the same during the precursor and array simulations, as well. The normal gradient of pressure is zero. Because we are modelling the lower boundary as a rough surface, the velocity normal to the surface is zero, and there is no need to specify the horizontal velocity. The no-slip condition is not appropriate on a rough surface in which the roughness elements are not resolved because the mean velocity goes to zero at the aerodynamic roughness height, *z*_{0}. Instead, the total of the viscous and SFS stresses on the lower boundary is modelled and applied following the work of Moeng [19]. On this boundary, the total stress tensor only has two finite independent components, and they are modelled as
2.6where the subscript *i* is 1 or 2 referring to the *x*- or *y*-directions, respectively; the subscript ‘1/2’ refers to values at the centres of the cells adjacent to the surface (values at one-half the grid cell height); angle brackets denote horizontal averaging; and is the magnitude of the resolved horizontal velocity. The symbol *u*_{*} is the friction velocity, which is related to the mean surface stress as
2.7To solve equation (2.6), *u*_{*} is needed. Assuming that the first grid level cell centres lie within the surface layer of the tidal boundary layer, friction velocity is estimated using the log law for rough walls,
2.8where *z*_{1/2} is one-half the height of the grid cells adjacent to the lower surface, *z*_{0} is the aerodynamic roughness height and *κ* is the von Kármán constant. In this study, we have taken *z*_{0} as 4 cm and *κ* as 0.41.

The domain measures 240×80×70 m in the streamwise, cross-stream and vertical directions, respectively. There are 410×180×160 grid cells in each of these directions, yielding a grid resolution of roughly 0.5 m in all directions. The domain is longer in the streamwise direction than the cross-stream in an attempt to capture the elongated streamwise turbulent structures in the boundary layer, and to allow for a small array of turbines that would be less closely spaced in the streamwise direction than in the cross-stream.

### (d) Turbine model

The geometry of the turbine is not resolved because that would require such small grid spacing near the surface of the structure as to render the computation infeasible. The wall model that we use to avoid that issue for the lower surface of the channel is not designed to model the turbulence near the turbine structure surfaces. Instead, the turbine is represented using Sørensen & Shen's [20] actuator line model that imposes body forces on the flow equal and opposite to the lift and drag forces experienced by the turbine.

Each turbine blade is represented as an actuator line. Each actuator line is divided into a number of equally spaced segments. In this study, there are 40 segments per line. The blade aerofoil type, chord, twist and local flow velocity are known at the centre point of each actuator segment. With that information and the rotor speed, the velocity magnitude, *V* _{mag}, and local flow angle, *α*_{local}, can be computed. The local flow angle is different from the angle of attack, *α*. The angle of attack is the angle between the aerofoil chord line and the free stream, whereas the local flow angle is the angle between the aerofoil chord line and the local flow. When an aerofoil creates lift, there is an upwash in front of the aerofoil caused by the bound circulation, which causes the flow to curve upwards. This means that *α*<*α*_{local}. We currently do not include a correction for this, and we assume *α*=*α*_{local}. In the future, we will incorporate a correction like that discussed in Shen *et al.* [21]. Employing lift and drag look-up tables for the rotor aerofoils, the magnitude of lift, *L*, and drag, *D*, normalized by density at each actuator segment can be computed using
2.9and
2.10where *C*_{l} and *C*_{d} are the lift and drag coefficients, respectively; *c* is the chord length; and *w* is the actuator segment width. Knowing that lift and drag forces act perpendicular and parallel to the oncoming wind vector, respectively, one uses those forces to form the total force vector experienced by the turbine, .

Each actuator segment has its own value of , which is a point force that cannot be directly applied to the flow field; rather it must be projected onto the flow field volume. The volume integral of the volume-projected force must be the same as the original actuator segment force. Following Sørensen & Shen [20], the projected force at a grid cell located at (*x*,*y*,*z*) is related to the actuator segment forces using a Gaussian projection as follows:
2.11The summation is calculated over all *N* actuator segments of the turbine. The location of the *j*th segment is (*x*_{j},*y*_{j},*z*_{j}), where is the magnitude of the vector between (*x*,*y*,*z*) and (*x*_{j},*y*_{j},*z*_{j}), *ε* controls the width of the Gaussian, and the negative sign accounts for the fact that the force that the turbine exerts on the flow field is equal and opposite to the force it experiences from the flow. The value of *ε* is equal to twice the cube root of the volume of the cells local to the turbine. (The cells local to the turbine all have the same volume so *ε* is a constant.) This value of *ε* is roughly the minimum value at which the force is smoothed enough to avoid spurious oscillations in the resulting velocity field.

The turbine being modelled is a hypothetical 550 kW two-bladed machine meant to operate in flow speeds ranging from 0.5 to 3 m s^{−1}. It has a 20 m rotor diameter with NACA 63 series sections. This design has been extensively analysed by researchers at the National Renewable Energy Laboratory (NREL), and is outlined in Lawson *et al.* [22].

### (e) Numerical method

Our solver was created using the Open Field Operations and Manipulations (OpenFOAM) CFD toolbox [23]. The OpenFOAM software is a set of C++ libraries meant for solving partial differential equations. The governing equations are solved using the finite-volume method. Although the solver is unstructured, our grid is composed of hexahedral elements and could be described in a structured way. All variables, except SFS viscosity and stress, are cell-centred and collocated on the grid. SFS viscosity and stress are located on cell faces. To avoid the pressure–velocity decoupling that occurs with collocated, incompressible solvers, the velocity fluxes at the finite-volume faces are constructed using an interpolation similar to that in Rhie & Chow [24]. All other interpolations from cell centres to faces are linear, which is similar to second-order central finite differencing. Time advancement uses the pressure-implicit splitting operation (PISO) algorithm [25], which is a semi-implicit predictor–corrector scheme. We say it is ‘semi-implicit’ because we do not treat all the terms implicitly. The SFS stresses, turbine forces and velocity flux within the convective term of the momentum equation are treated explicitly to avoid iterating during each time step. The implicit terms are integrated in time using Crank–Nicholson discretization. One predictor followed by three correctors are used. The momentum equation (equation (2.2)) is solved directly. However, to enforce the continuity equation (equation (2.1)), the divergence of the discrete momentum equation is used, which results in an elliptic equation for pressure. The pressure equation is solved using an iterative diagonal incomplete Cholesky preconditioned conjugate gradient linear system solver. The momentum equation is solved using an iterative diagonal incomplete LU preconditioned biconjugate gradient linear system solver. The code is parallelized using the Message-Passing Interface.

## 3. Cases simulated

### (a) Precursor simulation

The precursor simulation is run for 6000 s, which corresponds to seven large-eddy turnover times, where a large-eddy turnover time is defined as *τ*_{*}=*d*/*u*_{*} and *d* is the channel depth of 70 m. Although Moeng & Sullivan [26] state that six large-eddy turnover times are sufficient to reach a quasi-steady state in planetary boundary layer flows, we need to perform a more definitive study to determine the time at which the quasi-steady state occurs for this type of flow. Beginning at 6000 s, planes of velocity boundary data in the streamwise direction are collected for an additional 500 s to be used as inflow for the turbine array simulations. The volume flow field at 6000 s is used as the initial flow field for the array simulations.

### (b) Tidal turbine array simulations

Ten different turbine array cases are computed using the same inflow and initialization data from the precursor simulation in all cases. Figure 1 shows the plan view of the turbine array, in a non-staggered or staggered configuration. The downstream spacing, *S*_{1}, and the cross-stream spacing between rotor discs, *S*_{2}, are variable. The turbines are labelled T0, T1, T2, T3, and T4. Outlined in table 1, each case has different streamwise and cross-stream turbine spacings and a different turbine rotor rotation sense (cw, clockwise; ccw, counter-clockwise); consecutive rows have stagger or no stagger.

## 4. Results

### (a) Precursor simulation

During the precursor simulation, statistics are gathered between 5500 and 6500 s. Figure 2*a* shows the mean horizontal velocity profile of the precursor tidal channel flow. The turbine rotor disc lies between the dotted lines. It shows that the mean velocity at hub height (*z*/*d*=0.5) is 1.9 m s^{−1}.

Resolved velocity variance profiles of the precursor flow are shown in figure 2*b*. They follow the expected progression in which the streamwise velocity variance is greatest in magnitude, followed by the cross-stream, and finally the vertical. The streamwise velocity variance at the point nearest the wall should be greater than the values farther from the wall. There should not be a peak in streamwise velocity variance, as observed here, because wall-modelled LES does not resolve down to where the actual physically occurring peak occurs. This problem is linked to the log-layer mismatch problem, and remedies are discussed in Sullivan *et al.* [27] and Wei & Brasseur [28].

Velocity spectra are also useful in assessing how well the precursor flow has been simulated. One-dimensional spectra of horizontal and vertical velocity taken in the streamwise direction at various heights are shown in figure 3. The horizontal velocity spectra exhibits a cascade region of −5/3 slope (grey line) followed by the LES filter cut-off at a wavenumber of roughly 2, which is as expected. The vertical velocity contains a smaller −5/3 slope region. None of the spectra contain a peak, indicating that the domain size we chose is too small. The largest resolvable scales on this domain are the most energetic, but the domain should be able to resolve scales that are even larger, but less energetic.

The resolved scales can be visualized by examining contours of velocity on horizontal planes through the flow. Figure 4 shows contours of instantaneous streamwise and vertical resolved velocity fluctuations on horizontal planes at two different heights. Both the streamwise and vertical velocity fluctuations are larger nearer the lower surface of the tidal channel, which agrees with the velocity variances being largest near the surface in figure 2*b*. Elongated regions of lower- and higher-speed streamwise velocity are seen near the surface, and the domain is not long enough to accommodate an entire structure. At hub height, these elongated structures are not present, but large pockets of lower- and higher-speed regions exist. Here, the grid is not sufficiently wide to accommodate these pockets. The scale of the vertical velocity structures near the wall is much smaller than that of the streamwise velocity, and the domain seems large enough to accommodate them. This condition is reflected in the vertical velocity spectrum at the lowest level in figure 3*b*. It shows that a peak would probably occur if the domain were enlarged slightly. At hub height, the vertical velocity structures are larger, but not as large as those of the streamwise velocity. Again, it appears that the domain should be larger to properly capture these large structures.

### (b) Turbine array simulations

An overview of a typical turbine array simulation flow field is shown in figure 5, which depicts an instance of the flow field from case J. The red areas are positive isosurfaces of the second invariant of the velocity gradient tensor, sometimes called the ‘Q criterion.’ This quantity is positive in regions of rotating flow, such as vortices. The tip and root vortices are seen spiralling behind the first two turbines. The vortex systems appear to break down about 1*D* downstream. The turbulence created in the wake of the first row of turbines aids in breaking down the vortices of the second row of turbines, visible 4*D* downstream of the first row. Note that the flow is periodic in the cross-stream direction, so only half of turbines T2 and T4 are seen in the second row. Streamwise velocity is volume rendered in figure 5, meaning that both colour and opacity values are given to each ‘voxel’ of the flow field. The opacity is set such that only the lower-speed fluid, in the turbine wakes and near the tidal channel's lower surface, is visible in blue to light aqua. Note that the flow speeds in the wakes are significantly lower than the 1.9 m s^{−1} inflow hub-height average speed.

Figure 6*a* shows contours of resolved streamwise velocity in a horizontal plane at hub height from case B. The wake behind the first row of turbines is clearly visible with flow speeds reduced to roughly 1.2 m s^{−1}. Note that because we are not modelling the nacelle, there is a region of higher-speed fluid through the centre of the rotor. To maintain continuity, the flow speed increases between the rotors. This suggests that a staggered configuration may be advantageous. The low-speed wakes of the first row of turbines impinge upon the second row. The wakes of the second row of turbines exhibit streamwise velocities as low as 0.5 m s^{−1}. The second row of turbine wakes also appear to meander significantly and contain more turbulent motions.

Figure 6*b* shows contours of resolved vertical velocity in the same horizontal plane as figure 6*a*. Because the turbines in case B all rotate clockwise as viewed from upstream, their wakes rotate counter-clockwise. This effect can only be captured using a turbine model that imposes tangential forces on the flow field, like the actuator line or an advanced actuator disc (see Porté-Agel *et al.* [29] for a discussion on the importance of including tangential forces in the turbine model).

To better understand wake evolution with downstream distance, it is useful to look at time-averaged velocity profiles. Figure 7*a* shows the vertical profiles of time-averaged resolved streamwise velocity from the single turbine case A, taken at various downstream locations behind the turbine and aligned with the turbine's centreline. The wake is vertically asymmetric because of the vertical shear in the mean flow. As in figure 6*a*, the velocity peak at the rotor's centre is clearly observable. By about 6*D* downstream, where *D* is rotor diameter, this peak disappears, and the wake more closely resembles that of a bluff body. At 8*D* downstream, the velocity at the centre of the wake is 22 per cent lower than the ambient hub-height average flow speed. In a tidal channel with more ambient turbulent intensity than the 6 per cent simulated in this study, one would expect the wake recovery to be more aggressive than observed here. Figure 7*b* shows horizontal profiles of time-averaged resolved streamwise velocity from case A at various downstream locations and at hub height. The velocity profiles are asymmetric in this plane, as well. Notice that up to 6*D* downstream, there still exist two velocity minima adjacent to the wake centreline. The minima on the left side have a smaller velocity than those on the right. If you imagine the vertical profiles, shown in figure 7*a*, rotated counter-clockwise, as viewed from upstream (the direction that this wake rotates), it is clear that the left minima should have a smaller velocity because it corresponds to the part of the wake that was rotated from the lower-momentum part of the wake nearest the tidal channel lower surface. The profile also is asymmetric at distances past the edge of the wake. The left side of the plot shows a higher than ambient velocity, whereas the right side shows one that is lower. This may have to do with the fact that our domain is not large enough to capture the largest turbulent structures. Figure 3*b*, which presents contours of instantaneous hub-height streamwise resolved velocity from the precursor, shows that the left side of the flow (as viewed from downstream looking upstream) contains a higher-speed structure that is continuously cycled through the flow in the precursor. Likewise, there is a lower-speed structure on the right. When the precursor boundary data are saved, the effect of these structures is saved, and is applied to the turbine array domain. The effect shows up as a consistently higher-speed flow on the left side of figure 7*b* and a lower-speed flow on the right side.

Table 2 outlines the time-averaged power produced, by each turbine and each case collectively, relative to the power produced by a single turbine (case A). Values for turbines T2 and T4 from cases I and J are not shown—they are unreliable because only half the turbine was computed. Turbine T0 in the front row always produces more power than the adjacent turbine T1. The same explanation applies here as that which explains why the left side of figure 7*b* has a time-averaged velocity in excess of the ambient flow, whereas the right has one less than the ambient flow. The high- and low-speed turbulent structures ‘trapped’ on the left and right sides, respectively, of the precursor (as viewed from downstream looking upstream) cause turbine T0 to produce more power and T1 to produce less power.

It is also clear that for the non-staggered, co-rotating cases B–D, increasing the spacing between rows increases the overall array efficiency from 78 to 84 per cent. Staggering the downstream rows, though, even with a downstream spacing equivalent to the non-staggered case B, produces the greatest benefit. The staggered cases I and J perform more efficiently collectively than does the single turbine case A.

Counter-rotating the turbines has little effect on power produced. If consecutive turbines in a cross-stream row rotate counter to one another, but turbines downstream of one another co-rotate, there is no advantage (compare case B with E). If, however, the downstream turbines rotate counter to the next upstream turbines, an efficiency increase of about 3 per cent is observed (compare case B with F). Decreasing the cross-stream spacing of counter-rotating turbines in a row seems to have little effect (compare case B with G and H).

## 5. Conclusions

Our initial work in applying LES to arrays of tidal turbines demonstrates that this method can provide useful unsteady information about wakes and power production. The accuracy of simulations performed within the framework presented here will be assessed as experimental field data become available. Although LES is more computationally expensive than CFD based on the RANS equations, it yields a detailed, time-dependent solution that contains all of the large eddies within the flow. Such simulations could be used to develop improved RANS turbulence models for tidal turbine array calculations.

This work highlights the importance of the application of tangential forces by the turbine model in recreating the wake asymmetry that exists behind actual turbines. However, there is currently a lack of comparison to experimental wake data behind turbines—a deficiency we plan to correct in the future.

Our research also shows that the way in which the turbulent inflow is simulated greatly affects the predicted wake propagation and power production of the array flow. For example, the domain we used here is not long enough to capture the largest scale structures that form in this flow. Therefore, a low-speed structure as long as our domain became ‘trapped’ in the domain and continuously cycled through it during the precursor simulation. This caused the turbines on that side of the domain to consistently produce less power than the turbines on the other side. In the future, we plan to determine the domain size necessary to resolve the largest scales. We also plan to find tidal channel observation data with which to compare our precursor data.

In future studies, we plan to include the physics of salinity and temperature stratification and analyse its effect on the inflow and array performance. It would be interesting to simulate the tidal-cycle variation of the inflow, and its effect on the performance and efficiency of a tidal turbine array.

## Acknowledgements

All computations were performed on the US Department of Energy's Red Mesa high-performance computing system. Financial support for this work was provided through the US Department of Energy Water Power Program. We thank our colleagues: Michael Lawson (NREL) provided useful discussion in performing these simulations, and Tony Martinez (University of Puerto Rico) provided the original development and implementation of our actuator line code.

## 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.