Effects of Support Structures in an LES Actuator Line Model of a Tidal Turbine with Contra-Rotating Rotors

: Computational ﬂuid dynamics is used to study the impact of the support structure of a tidal 1 turbine on performance and the downstream wake characteristics. A high-ﬁdelity computational 2 model of a dual rotor, contra-rotating tidal turbine in a large channel domain is presented, with 3 turbulence modelled using large eddy simulation. Actuator lines represent the turbine blades, 4 permitting the analysis of transient ﬂow features and turbine diagnostics. The following four cases 5 are considered: the ﬂow in an unexploited, empty channel; ﬂow in a channel containing the rotors; 6 ﬂow in a channel containing the support structure; and ﬂow in a channel with both rotors and 7 support structure. The results indicate that the support structure contributes signiﬁcantly to the 8 behaviour of the turbine and to turbulence levels downstream, even when the rotors are upstream. 9 This implies that inclusion of the turbine structure, or some parametrisation thereof, is a prerequisite 10 for the realistic prediction of turbine performance and reliability, particularly for array layouts where 11 wake effects become signiﬁcant. 12


Introduction
The commercial exploitation of tidal energy on a large scale requires the deployment of arrays of full-scale tidal turbines.Given individual turbines of rated power of 1-2 MW, such arrays would have to consist of 50-100 turbines to approach the operating capacities of modern offshore wind farms.
Individual turbines within a farm array will be affected by the wake of any turbines located upstream, and the large-scale environmental flow impact of the farm as a whole must also be understood; thus, modelling tidal arrays becomes a true multiscale problem.The application of computational fluid dynamics (CFD) can shed light in both areas, but this is extremely challenging from a computational perspective.
Wake effects in wind farms have been the subject of many studies.Models range from early empirical linear wake superposition approaches such as the Park model [1], through to Reynolds-Averaged Navier-Stokes (RANS) CFD actuator disc models, large eddy simulation (LES) actuator disc [2][3][4] and actuator line models [5][6][7][8].Detailed reviews of wind turbine and wind farm wake modelling are given by Barthelmie et al. [9], Sanderse et al. [10] and Creech and Früh [11].One striking feature of these models is that, bar a few exceptions [12,13], the turbine support structure is not modelled explicitly, and so only the rotors affect the downwind flow.It is quite likely that, in the mid-to-far wake region, wake effects due to the structure are not important in wind farms; indeed previous, validated studies of single wind turbines [14] and wind farms [4,15] have indicated that the tower and nacelle have negligible impact on the wake and consequently the performance of downwind turbines.The pertinent question here is, then, can the same be said for tidal turbines sited in swiftly flowing water, whose density is over 800 times that of air?At basin scale, it is common to use depth-integrated shallow flow models to assess tidal stream power.In many cases depth-integrated models are used [16][17][18][19], with turbines represented by increased sea bed resistance, and the drag coefficient tuned to include both thrust and structural drag.These representations of turbines enable the thrust to vary with upstream flow speed, but are unable to properly resolve the three-dimensional flow kinematics that occur in the wake of a turbine rotor.Three-dimensional computational fluid dynamics (CFD) is capable of modelling resolved blade motion [20] in good agreement with laboratory experiments [21], but is expensive in terms of computational resources.In such models, simulating wakes over realistic distances downstream (i.e.many multiples of rotor diameter) is extremely challenging, especially with high-fidelity turbulence modelling techniques such as LES, due to the necessity of refining the mesh for the blade boundary layer.Therefore, parameterisation of the blades is required for simulations in larger domains.
An early example of this were LES simulations of a turbine in 800m-long tidal channel using a dynamic actuator disc turbine model [22].This work found that the tidal turbine wake length, when scaled by power output, was on with wind turbines.Others have focussed upon Reynolds Averaged Navier Stokes (RANS) CFD models with actuator disc representations [23][24][25], obtaining good agreement with experimental data.Afgan et al. [26] and Ahmed et al. [27] compared blade-resolved RANS and LES simulations, demonstrating that LES predicts greater fluctuations in blade loads, whilst in similar work McNaughton et al. [28] found that whilst LES produces better agreement with experiments, k − ω RANS models produce acceptable results for far less computational cost.
Churchfield et al. [8] employed actuator line models to produce simulations of four turbines, without support structures, to examine wake effects on downstream turbines.See Section 6 for further discussion of these.
Using LES and the Fluidity CFD software from Imperial College [29], we examine individual and cumulative contributions to the downstream wake of both the rotors and structure in a dual rotor, contra-rotating tidal turbine, located in a large rectilinear channel.The channel is sufficiently large to capture most of the wake, be of representative depth, and contain realistic, fully developed turbulent flow.Whilst computationally demanding, such simulations can provide a wealth of accurate detail, so providing insight into the complex interactions between the rotors and structure.This in turn can inform cheaper, quicker alternative models, such as those used in iterative design and assessment

Initial test cases
For model verification purposes, two preliminary computational tests were conducted to ensure the parameterisations gave realistic results in the absence of the turbine rotors.Simulations were of sheared flow in an empty channel, without, and then with a vertical cylinder surface-piercing present.
The configuration of the rotor and blades is dealt with separately in Section 3.

Basic equations
For all simulations, an incompressible Newtonian fluid was assumed.A control-volume finite element discretisation [30] was used, with first-order, continuous velocity and pressure elements, and a Crank-Nicholson time stepping scheme.Large Eddy Simulation (LES) modelled the effect of unresolved (subgrid) turbulence in fluid flows in the simulations.LES, first developed by Smagorinsky [31], was later adapted to channel flows by Deardroff [32]; the variant applied here within Fluidity takes into account mesh anisotropy [33], as Deardroff's istropic estimate for filter length breaks down as the cell aspect ratio increases [34].Here, the filtered momentum and continuity equations are, respectively, in Einstein notation ∂ ũi where ũi is the filtered (above grid level) i th velocity component, ρ is the fluid density, p is pressure and ν is kinematic viscosity.For the following simulations, ρ = 1027 kg m −3 and ν = 1.831 x10 −6 m 2 s −1 .For application on anisotropic meshes, the subgrid eddy viscosity is represented by a tensor, defined as where C S is the Smagorinsky coefficient, set to 0.1 for all simulations [32], S is the rate-of-strain tensor, and ∆ ij is the element size tensor.More details on the anisotropic LES formulation within Fluidity are found in Bentham [35] and Bull et al. [33].
Before running the 'production run' simulations, a series of test cases were run, and the results compared with published data.The results were used to validate the simulation configurations used, including mesh resolution, turbulence modelling and boundary conditions.For both cases, the inflow boundary conditions were identical.

Specification
Figure 1 shows the idealised channel domain, measuring 1 km x 200 m x 30 m.The chosen depth was close to that of Strangford Narrows where the SeaGen tidal device is situated [36], and the 1 km length allowed the wake behind the turbine to be captured within the model.Furthermore, the domain dimensions would allow large eddies tens of metres across to develop without impingement due to a restrictively small domain.
The surface of the channel was represented as a frictionless, rigid lid, and the lateral walls were also frictionless.Seabed drag was estimated empirically using the quadratic drag law with a bed friction coefficient of C F = 0.005; noting that the quadratic drag law has been found to fit measurements of turbulent tidal flow [37].An open boundary condition was applied at the outflow.The synthetic eddy method (SEM) [38] was applied at the inlet to generate a turbulent inflow.The mean velocity profile was based upon a logarithmic profile, ie.
where u τ is the frictional velocity, K (= 0.41) is the Von Kármán constant and z R is the roughness height of the channel bed, set to 0.05 m.B is a constant, which for turbulent open channels can be taken as B = 8.5 [39].If the flow speed at hub height z H is specified as u H , the frictional velocity can be calculated as where u H is the mean velocity at hub height, set to 2.0 ms −1 .
Building upon previous work [4], both the mean eddy lengthscale and Reynolds stress profiles were specified as a function of height above the seabed for SEM, so that realistic turbulent inflow was generated.Eddy lengthscales were taken from Milne et al. [40], whose measurements from the Sound of Islay agreed with Nezu and Nakagawa [39].This gave the streamwise integral turbulence lengthscales as Cross-stream and vertical components of eddy lengthscale were specified as L v = 0.5 L u and L w = 0.25 L u respectively.The Reynolds stress profiles were taken from Stacey et al. [41], which following Nezu and Nakagawa [39] gives the three diagonal Reynolds stress components for unstratified channel flow as The normalised streamwise comment, R u u, is shown in figure 2.
For the computational mesh, the maximum element dimensions were [2 m, 2 m, 1 m], reduced to [1 m, 1 m, 0.5 m] within a distance of 2 m of the seabed.The overall mesh contained 17.6 million elements, partitioned across 480 computing cores.The time step was fixed at ∆t = 0.3 3 s, with the pressure and velocity fields recorded every 1 second.The model ran initially for 30 minutes of simulation time to 'spin up', followed by another 30 minutes over which flow was to be recorded for analysis.At the time the simulations were carried out, high temporal resolution point probes (detectors) were not functional within Fluidity, which meant that full-domain data outputs were required.This in turn curtailed the sampling period, due to the excessive volume of data produced.Figure 3 shows a typical output of the velocity magnitude distribution throughout the channel at the end of the simulation (at time t = 1800 s).Roll-up of vortical structures can be seen at the bed, consistent with the development of turbulent eddies in open channel flow.

Results
Time-averaged vertical velocity magnitude profiles at a resolution of 0.5 m were taken from the centre of the channel, as shown in Figure 4a.These were calculated at different locations along the centreline of the channel; the distance downstream is plotted in units of D, the rotor diameter of the tidal turbine to be modelled (16 m), with the origin at 250 m downstream of inflow boundary.It can be seen that the time-averaged profile at x = −5 D is very similar to profiles further downstream, with no deviation at any point greater than 0.1 ms −1 at any height or distance downstream, even to x = 20 D.
As a turbulent channel flow with bottom drag, a logarithmic vertical velocity profile should be expected.A logarithmic regression fit was applied to the mean of the velocity profiles in Figure 4a, which gave the following equation u l (z) = 0.26348 ln(z) + 1.29458 (10) Nezu and Nakagawa [39] suggest that the logarithmic law may only be valid in the wall region, and that a power law may be more appropriate.Therefore a power law regression fit was also applied to the mean vertical velocity profile.As the roughness on the channel bottom is not explicitly resolved, the roughness height z R was instead derived from the skin friction coefficient C F [4] for a more appropriate fit.This gave the power law where z R = 0.04852 m.
If we express the exponent as 1/a, then equation ( 11) gives a = 6.48004.Whilst a = 7 is a commonly quoted figure [42], the derived value compares favourably with ADCP profile measurements from Strangford Narrows, from which a = 5 on the flood tide, and a = 7 on the ebb tide [43].Figure 4b superimposes both the log and power law fits on the spatially-averaged velocity magnitude profile.The model profile and the derived log-law match well, apart from a slight overshoot by the log-law near the surface, and a slight undershoot near the channel bottom.This may be due to numerical diffusion arising from insufficient grid resolution.Unfortunately, increasing mesh resolution in this region is presently not an option owing to the prohibitively computational expense; nonetheless, there is good overall agreement, particularly in the mid-region area of interest, near where the turbine rotors will be situated.To quantify the error between the model results and the log plot, the relative 2-norm error was used: x = 0 is where turbine rotors are to be placed, 250 m downstream of the inlet.
where N is the number of sample points in the vertical profile, u m denotes the model results, and u r is the regression fit.
The error norms were determined as l = 0.01918, or under 2% error, and p = 0.03202, or just over 3% error.These were deemed acceptable margins.The turbulent intensity (TI) profiles in Figure 4c show as expected a low TI value (7%) at the surface, which gradually increases towards 15-18% at the channel bed.This compares well with Milne et al. [40], where ADCP measurements in the Sound of Islay gave a TI of 10-11% and a mean flow speed of 1.5 ms −1 , at 5 m above the seabed.The limited sampling frequency of 1 Hz mentioned in Section 2.2.1 means that the higher-frequency turbulence Nezu and Nakagawa [39] found in the lower section of the channel is not detectable.It is likely that if detectors had been available, a more pronounced peak near the bed would have appeared.Even so, the fit of the model data to the log and power law velocity profiles gives confidence that the channel simulation is a reasonable representation of turbulent channel flow.

Channel domain with a cylinder
The purpose of this test was to develop and validate an adequate representation of a structure within the domain, insofar as its effect on the flow is realistic.Flow past a cylinder represents an excellent test case for modelling the flow around a structure, as it is a widely-known problem [44][45][46][47][48] that has been studied extensively using CFD.It is well established that vortex shedding at the cylinder occurs at a predictable frequency for Reynolds numbers within the range 250 < Re < 10 5 ; this behaviour should be observed in the model.

Specification
Previous examples of simulated flow past a cylinder with LES, have involved modelling the boundary layer equations [49,50], or by using Van Driest damping functions to satisfy the zero eddy-viscosity condition at the cylinder surface [51].Neither of these options was practically available due to the size of the domain, and neither was a feature within Fluidity.Instead, an intermediate approach was adopted: to resolve the mesh finely around the cylinder and downstream as far as possible, but to also impose a quadratic drag boundary condition on the cylinder surface.Such a solution would be sensitive to both the mesh resolution near the cylinder and the skin friction coefficient C F chosen.To verify the approach, the Strouhal number St was calculated from the results

St
where f is the frequency of the vortex shedding, D c is the diameter of the cylinder, and u H is the upstream speed of the fluid.By looking at the fluctuations in lift forces acting on the cylinder, the vortex shedding frequency f can be calculated, and so the Strouhal number.
A vertical cylinder of diameter 3 m, similar to the main tower of SeaGen [52], was placed with its centre at [250 m, 100 m, 0 m] on the seabed, extending to the surface 30 m above.As confirmed by the empty channel tests in Section 2.2, this would allow the turbulence sufficient time to develop fully, whilst also avoiding any blockage effects due to narrowing of the passage between the cylinder and the channel walls.Mesh resolution was increased to [0.25m, 0.25m.0.25m] at the surface of the cylinder, as shown in Figure 5.The simulation ran for 1800 s, with a timestep of ∆t = 1 3 s.As with the all the simulations, the velocity and pressure fields were output every 1 s.

Results
To determine the Strouhal number, the frequency of vortex detachment from the cylinder was checked by calculating the lift on a 1-metre thick ring on the cylinder at hub-height, ie.z = 16 m.This was done for two reasons: firstly, as the cylinder is in vertically sheared flow, the Reynolds number can be expected to vary widely from the top to bottom, and secondly, increased turbulence near the seabed would cause large pressure fluctuations not associated with vortex detachment, giving a noisier signal.
The scale of the simulation can be seen in the instantaneous velocity slice in Figure 6, with a close up showing the vortex street caused by shedding in Figure 7.
A fast-Fourier transform (FFT) was applied to the lift force fluctuations.The resulting power spectrum in Figure 8b contains a sharp peak around 0.22-0.225Hz.For u H = 2 ms −1 , equation ( 13), this  Although lower than Achenbach [55] (C D =0.6-0.7), in general there is remarkably good agreement, given the very low blockage ratio presented here, and the differences in oncoming flow profiles.
Therefore, the combination of cylinder surface mesh resolution and quadratic skin drag law was deemed sufficient for simulation of realistic wake effects.

Turbine formulation
The turbine model developed here builds upon previous work, where dynamic torque-controlled actuator discs with active-pitch correction were used to model wind turbines [14] and wind farms [4,15].In the present model, actuator line techniques [5]

Methodology
The lift and drag force components per unit span acting on a blade are given by where C L (α, Re) and C D (α, Re) are the coefficients of lift and drag respectively, both functions of angle of attack α and Reynolds number Re; ρ is the density of the fluid (for a tidal turbine, seawater) in which the blades move; u rel is the relative speed of the fluid over the blades; and c(r) is the chord thickness as a function of r, the radial distance from the hub centre.In practice, u rel is calculated for each cell point, using the local flow speed.
The Gaussian distribution functions η i at a point in space x for turbine blade i of N blades in the rotor are expressed as where σ is a constant, the standard deviation which controls the filter width, and d i is the ring distance of point x from the actuator line.σ was chosen with care, as too large a value could result in a heavily smeared solution, whereas too small necessitates an extremely fine mesh and very small timesteps.In this case, it was found that one-twentieth of the rotor radius gave an acceptable trade-off between accuracy and computational effort.
Assuming that each blade has identical geometry, aerodynamic characteristics and blade pitch, then We apply this to determine the lift and drag terms as body forces, taking into account tip losses, to give where T is the Prandtl tip-loss factor [14,57].As with previous work [4,14,22], blade-generated turbulence was then added via randomly fluctuating components.These body forces acting on the blades are translated into axial and azimuthal components; from Newton's third law, the consequent body forces on the flow are equal and opposite.Each time-step, these terms are calculated and passed back to the CFD solver to be included in the Navier-Stokes momentum equations.The force terms above are also used to calculate the power output of the turbine.By first calculating the net torque acting on the fluid, we can then calculate the resistive torques turning the generator and blades [4,14], thus turning the drive shaft and the blades.Both the drive chain and the power conversion have associated energy losses, which are written as where P real is the actual power, P ideal is the ideal power without energy losses, E d is the drive train efficiency, and E g the generator and power conversion efficiency.We used E d = 0.94, and E g = 0.96, the values specified in Bedard [58] for the MCT/Siemens' SeaGen device.An active pitching algorithm was used which maximises total lift, matching the behaviour of SeaGen [59].Further details of the numerical model can be found in Creech et al. [4].

Parameterisation
The rotor configuration was based upon that of Marine Current Turbine's SeaGen device [60], ie.
dual rotors, aligned horizontally.As many of SeaGen's technical details are commercially sensitive and not readily available, rotor and performance specifications were sourced from journal papers [52,58,[60][61][62][63].Details often disagreed between papers, so discretion was applied in deciding on the values listed in Table 1.To validate the chosen parameters, candidate models were tested against performance data from SeaGen, as shown in Figure 9.
The aerofoil chosen was a NACA 63-415 type, which has desirable lift characteristics.The lift and drag at limited angles of attack were taken from previous work [4], which based its aerofoil data on that from the Airfoil Catalogue [64].Blade geometry was completely unknown, so as a starting point, the equation for the predicted flow angle at a turbine rotor was taken from Burton et al. [57]: where φ is the predicted inflow angle as a function of radial distance r, λ is the design tip-speed ratio, and µ = r R , where R is the radius of the turbine rotor.If the optimum angle of attack for a given blade α opt , then the blade twist can be given as α opt was calculated from the lift and drag coefficient charts for the chosen aerofoil type, as per Creech et al. [4].This then provided the blade twist angles along the blade.
No information was available on chord thickness, so this was calculated using the equation for an ideal optimised blade derived from blade-element theory [57, Chapter 3], ie.
where σ r = Nc 2πr is the rotor solidity (N is the number of blades, c the local chord thickness), µ = r/R, and C L,opt is the lift coefficient at optimal operation, which is calculated from lift and drag performance data.Rearranging gives  where X(µ) is the right-hand side of (23).The chord length can now be defined as a function of r for a blade with specified characteristics.

Test cases and results
Test cases were devised to check that the specification of the rotor, the blades, and the generator properly represented the real tidal turbine, in a channel flow u H = 2 ms −1 .To lessen computational requirements, a single rotor would be tested in a simulation domain much smaller than the empty channel test case, measuring 250 m x 100 m x 30 m.No support structure was included to further reduce the number of elements required.The turbine rotor was positioned much closer to the inlet, at [50 m, 50 m, 16 m].Such a short domain was acceptable because only the performance of the turbines was of interest, and the wake effects could be ignored.Mesh resolution was increased towards the rotor, so that the actuator volume contained approximately 19 000 points.In the final simulations, the turbine rotors would be operating at below the rated power, so only hub-height flow speeds below 2.5 ms −1 would be tested.Three mean hub-height flow speeds were considered: u H = {1, 1.5, 2} ms −1 .
Figure 9 compares the resulting mean power and thrust from each case with published performance data.The model yields slightly larger time-averaged values for power and thrust at u H = 1 ms −1 when compared to the published data, but it is in good agreement for both power and thrust when u H is at 1.5 and 2 ms −1 , to within a maximum relative error of 8.8%.The cause of the over-performance of the model at the lowest flow speed could be down to minor discrepancies in the rotor, blade or generator specifications, but this cannot be confirmed.Other possible causes may be rotor-rotor wake interaction in the dual rotor configuration; this will be examined in Section 4.
Notwithstanding this, the results provide reasonable confidence in the modelled rotor performance at the target hub height flow speed of u H = 2 ms −1 .

Overview
The simulations were based upon the tidal channel test case in Section 2.2, as it was sufficiently large (1 km x 200 m x 30 m) to allow realistic turbulence structures to develop, and to capture the extents of the turbine and structure wakes.To isolate and analyse the contributions of the various components of the turbine to the downstream wake, three different scenarios were considered: a) with dual contra-rotating rotors and no structure, b) with the complete support structure only, and c) with both the rotors and the structure present.For each of these scenarios, the mean hub-height flow speed was set to u H = 2.0 ms −1 , with turbulent inflow conditions provided via the synthetic eddy method used in Section 2.2.1.These simulations were run until the turbulent flow was fully developed, and statistical properties of the flow, such as turbulence intensity and time-averaged flow speeds, had become stable; this was a minimum of 2700 s in all runs.Each simulation was then run for a further 900 s, during which full sets of data for the velocity and pressure were saved to disc for analysis each second for post-processing.As detectors were not available, higher frequency sampling of velocity fields for spectral analysis was not possible.
Nevertheless, the 900 s sampling period was sufficient to provide time-averaged velocity profiles that appeared to be statistically stationary.All simulations were run on ARCHER, the UK's national academic supercomputer, using 2400 computing cores each of which typically used 1.5 MAUs of allocation units, with a wall time of 2-3 days.

Configurations
The support structure model is shown in Figure 10, which was based upon SeaGen's design.The model consists of a 30 m high, 3 m diameter monopile that pierces the water surface, with a crossbeam at a height of 15 m above the sea bed.The crossbeam is 27 m broad, and 4 m long, encompassing the monopile, and contains angled sections on either side, that rise 1 m and taper to 3 m long at their furthest extent.Solid nacelle sections measuring 1 m x 3 m x 2 m are located at the ends of the crossbeam.The crossbeam edges have been smoothed to have a curved surface of radius 0.5 m, as have the nacelles.This arrangement closely follows details given by Fraenkel [52], Neill et al. [63] and Fraenkel [60], whilst also taking cues from Keenan et al. [36].The more complex quadropod base was not adopted in the final design, due to the prohibitive mesh refinement and complexity that would have been required.The structure was placed on the seabed in the empty rectilinear channel described in Section 2.2, such that the monopile base was centred at [250 m, 100 m, 0 m].
The contra-rotating rotors case used the design developed in Section 3. in Figure 11.Each rotor was connected to a modelled generator, and like SeaGen, these generators operated asynchronously [65].The mesh resolution was highest near the blades and reduced gradually with distance from the rotors, as can be seen in Figure 12.
The final case combined the dual rotors and structure configurations above.Table 2 lists the mesh resolutions and time-step sizes for each case.

Results
This section examines the time-averaged velocity and turbulence intensity data obtained from the full simulations of turbulent flow in the rectilinear channel with the dual rotors, with the supporting structure on its own, and with both the dual rotors with the support structure.Even with the aforementioned limitations in sampling time and frequency, the results give a good qualitative representation of the persistent flow features.

Wake effects
Here, the difference in the wake effects between each of the three cases are considered: rotors only, structure only, and rotors+structure.Two sets of profiles are examined: i) cross-stream (or transect) profiles at rotor hub height z H ; and ii) vertical profiles, on vertical streamwise planes slicing through both rotor hubs at y = {87, 113} m.Both sets of profiles are from x=-1 D upstream to 20 D downstream.
These are augmented by selected instantaneous velocity slices over the full length of the domain (≈47 D downstream).
From the time-averaged velocity profiles for the rotors-only case in Figure 13a, it is clear that by x=20 D, the flow has almost fully recovered to its upstream profile, u H (x = 20 D) having 90% of its upstream value (denoted u 0 ).Immediately downstream of the rotors at 1 D, the deficit matches closely the zone swept out by the blades, beginning at z t =24 m (top of the rotor) and ending at z b =8 m (bottom of the rotor).In the absence of a nacelle, the flow passes through the hub section (z=15.2-16.8m) relatively unimpinged, peaking at 0.85 u 0 .This hub-section flow is still evident at 5 D, but by 10 D has decayed into a quasi-Gaussian deficit.Upstream of the rotors, the vertical profile of turbulence intensity (TI) is similar to that in the undisturbed channel (Figure 4c); downstream it peaks at the hub section and the rotor tips, indicating the presence of tip-vortex shedding.The TI profile maintains this general shape until 10 D, whereupon it starts to begins to decay.By 20 D however, the TI values remain higher than the original upstream profile.The horizontal velocity transect at hub height in   Figure 14a shows a similar pattern, with the same mid-rotor peaks at y={87, 113} m, which eventually become smoothed troughs by 10 D. In the area between the rotors the flow speed increases to 2.2 ms −1 , due to blockage by the rotors and the absence of the support structure.Asymmetry occurs in the rotor deficits from 1-5 D, with the radially-inward tip deficits 0.2 ms −1 higher than the outward tip deficits.
The accelerative effect of the rotor blockage plays an important role here: in the instantaneous velocity field snapshot Figure 15a, there is a jetting phenomenon between the rotors.This has also been noted in previous work modelling offshore wind farms [4].The TI plots also exhibit peaks at the rotor tips and the hub section from 1D onwards; these persist even at 20 D, and do not decay to upstream levels.
For the structure-only case (Figures 13c and 14c), the velocity profiles show a sharp drop immediately behind the crossbeam, with nearly full velocity recovery by 5 D. The horizontal-transect profiles are more complicated, as the transect is at hub-height (z=16 m), and crosses in front of (and behind) the tower, as well as the upper, outer ends of the crossbeam (cf. Figure 10).Furthermore, the nacelles exert a strong influence on turbulence levels, raising the TI to 17.5% at 1 D, 13% at 5 D, before approaching background levels by 10 D. At 20 D, the difference between upstream TI values is negligible.
The velocity profiles for the rotors+structure case (Figures 13c and 14a) are broadly similar to the rotors-only case, but with several important differences.Firstly, the pronounced peak visible in the rotors-only profile (Figure 13a and 14a) has been replaced with flatter troughs.This is due to the influence of the crossbeam and the nacelle sections, causing the flow to accelerate around the solid structure, rather than through the empty hub volume in the rotors as before.This effect can be observed in the instantaneous vertical velocity snapshots in Figures 15c and 15d (zoomed in).Secondly, the flattened velocity peak in the horizontal transects between the rotors no longer occurs at 1D, owing to the presence of the monopile and crossbeam.Instead, two sharp spikes can be seen at approximately u 0 , either side of the central trough at 0.8 ms −1 .By 5 D downstream, these have become one peak at 1.8 ms −1 , somewhat lower than in the rotors case (2.2 ms −1 ).This pattern continues at 10 D and 20 D; the velocity profiles similar but ≈ 0.1 − 0.2 ms −1 lower in the rotors and structure regions.The horizontal velocity snapshot in Figure 15b reflects this, with the pronounced central jet between the rotor wakes in Figure 15a no longer visible.As expected, the vertical TI profile in Figure 13c is broadly similar to that in Figure 13a, but higher turbulence levels occur downstream, particularly behind the structure.At 1 D, the TI is between 20-25% near the nacelle region, much higher than the 12-17% for the rotors only.By 10 D however, TI profiles from both the rotors-only and rotors+structure cases are nearing equivalence.Lastly, of particular note are the TI peaks at the rotor tips, visible in Figures 13c     and 14c, which are 16-40% larger than the rotors case.This suggests that the blades themselves may be subject to, and causing, fluctuations in the flow.Turbulence spectra were calculated at a point 1 D downstream from top of rotor T1 for both cases (cf. Figure 16), which confirm that while both exhibit a peak at 0.4 Hz, this is particularly pronounced when the structure is included.We explore this further in Section 5.2.

Turbine diagnostics
Here we examine how the performance of the turbine is affected by the absence or inclusion of the support structure.Figure 17 presents a selection of diagnostics from the rotors-only test, which show that the power from each rotor fluctuates semi-independently of the other, due to the unsteady, turbulent flow each rotor experiences.We use the term 'semi-independently', because some of the eddies are large enough for each rotor to experience them simultaneously.For both rotors T1 and T2, power output varies approximately ±25 kW from their mean values.The mean power outputs for T1 and T2 differ by 0.25%; for simulations of longer duration these should converge.The weighted angle of attack along the blades, α, shows much less variation, as the pitch control mechanism tries to maximise power output [4,59].The results from the rotors + structure case are broadly similar.
The net time-averaged power output of cases for the single rotor, both rotors only, and then the rotors+structure are compared against published measurements for SeaGen in Table 3.The models are    in close agreement, with a maximum error of 5.1%.It is particularly interesting that the simulations modelling both rotors have more accurately predicted the power output than the single rotor case, suggesting that there is interaction between the two rotors.This may be down to the blockage effects of each individual rotor, which accelerates flow round the edges of the rotors, in theory providing a small performance increase.Indeed, the acceleration effect can clearly discernible in Figure 14a.
Investigating the effect of the support structure on the power output, Table 3 indicates little difference between the rotors and rotors+structure cases, with errors of 4.4% and 4.7% respectively.It is also worth examining the time-series of the power output from each rotor, to see if there are any regular fluctuations as the blades pass in front of the supporting structure (such as the crossbeam).
This was achieved by applying Fast-Fourier transforms (FFT) to the power output.As there was a short sampling period (900 s), and the likelihood of similar spectral characteristics in fluctuations in each rotor, for each simulation, FFTs of the power time-series for each rotor were calculated separately, and then the average of them taken.Figure 18 displays these for the rotors case and the rotors+structure case.Towards the lower end of the frequency range both plots become noisy, as longer wavelengths are not resolved satisfactorily within the sampling period.In Figure 18a, the rotors only case, there is a small peak at 0.4 Hz; however in Figure 18b, when the structure is added, the peak at the same frequency is much larger.The mean rotational frequency of the rotors in one simulation is defined as where ω Tn is a time-averaged rotor angular velocity from the diagnostics data for rotor Tn.
The values for the frequencies in Table 4 demonstrate that both cases are identical to within four significant figures.It should be noted that f C is almost exactly half of 0.4 Hz, shown in Figures 18a     and 18b.This is not surprising, given there are two blades.To determine where in the rotation cycle Figure 19a shows that in the rotors-only case, the power fluctuations have a relatively even spread within a range of ± 2.5 kW, with slightly higher values when either blade points upwards at 0 • .This increase is most likely due to the higher flow speeds and thus the lift which the blades experience at 0 • , as shown in the vertical velocity profile in Figure 13a.Discrepancies between this and the rotor+structure case are quite evident from Figure 19b, where a clear pattern emerges, with fluctuations peaking at a maximum of 7.5 kW at a blade position of 0 • and 180 • , and at a minimum of -7.5 kW at 90 • and 270 • .Table 5 lists these in terms of mean total output per rotor.It is clear that in the rotors+structure simulation, a dip in power output is experienced when the blades are aligned with the cross-beam.By comparing the horizontal velocity profiles at hub height in Figures 14a and 14c, just downstream of the rotors at x=1 D, the velocity deficits are more pronounced when the structure is included.In particular, the absence of the nacelles and tower is noticeable, with peaks flow speed of 1.75 ms −1 in the rotors-only case at y=87 m and y=113 m, and at y=100 m, where the nacelles and tower would be, respectively.This can be attributed to the blockage effect, created by the back thrust of the rotors, causing the flow to accelerate around the rotors and through the centre where the lift is reduced.In contrast, the velocity profiles in Figure 14c show no accelerated flow at the nacelles, and the velocity deficit created by the wake of the tower reaches 0.8 ms −1 , compared to the rotors-only case, which peaks at 2.2 ms −1 , ie. 1.1 u 0 .The culmination of the upstream effect of the supporting structure, in terms of the turbine performance, is a periodic fluctuation in power, responsible for a variation in output of almost 4%.Blade loading has not been analysed, but it can be expected to have more variation than power, due to the system inertia smoothing out fluctuations in power.The effect of inertia on power output can be seen in the slight asymmetries of the polar plot in Figure 19b.

Discussion and conclusions
Large Eddy Simulation has been used to model a full-scale, dual rotor, contra-rotating turbine, complete with structure, in a realistic-sized channel domain.The results demonstrate that the structure    does have a noticeable effect on performance and the near-wake, causing regular fluctuations in both power output and flow speed.
There are few CFD models of dual rotor configurations in modelling literature.Of these, the tidal turbine array LES simulations of Churchfield et al. [8], show similarities in the downstream wake profiles; however their domain was less than a quarter of the size of the channel domain used here, and so their results may be subject to exaggerated blockage effects due to proximity of the domain walls.
Furthermore, they did not use synthetic eddy methods for realistic inlet turbulence, unlike the present paper.Closer comparisons may be drawn with Afgan et al. [26] and Ahmed et al. [27], who used k − ω SST RANS, and also LES with SEM to model a resolved three-bladed rotor and the support structure.
Their time-averaged LES results of the power coefficient exhibited regular pronounced fluctuations as a function of blade angle with the rotors upstream; in the RANS simulations these were absent.The RANS models of Mason-Jones et al. [66] however, did predict increased fluctuations in torque with blade angle when the stanchion supporting the turbine was included.The lesson here is that care must be used when deploying RANS turbulence schemes to capture transient behaviour in diagnostics.
Moreover, the flow direction experiments of Frost et al. [67] showed that individual blade thrusts varied considerably with blade angle upstream, whereas net rotor thrust and power output did not.It must be noted that their turbine configuration differs somewhat from our model, having 3-blades with one central monopile, versus our dual two-bladed, crossbeam-mounted arrangement.Although there are valuable advances in research using actuator disc approaches [23,25,68,69], particularly in terms of tidal farm modelling, the discretisation of the blades into continuous rings means that such models cannot capture the same fluctuations in power output due to blade-structure interaction.Whether or not this is important for the far-wake of tidal turbines remains an open question, but for mechanical and electrical reliability, these transient features must either be represented or parameterised for accurate simulation.
In terms of wake prediction, measurements from a full-scale turbine would have been useful for model validation.As with the technical specifications however, ADCP measurements of SeaGen's downstream wake were unavailable, so comparison is instead made with the experimental and numerical literature.The water channel tests of Myers and Bahaj [70], despite being of scaled single rotor turbine, support the results in Figures 13b and 14c, that show the structure alone creates substantial turbulence even at 5 D downstream.With the rotors included, our findings agree with those of Batten et al. [24], and Stallard et al. [71,72], in that at 20 D (22 D in Batten), the turbine wake had still not recovered to its original value.They too report higher than background readings for turbulence intensity far downstream.In the near wake, there is good agreement between the cross-stream turbulence profile in Figure 14c and Tedds et al. [73], both showing peaks of 20-25% 1-2 D downstream, as well a peak at the centre near the structure.This is surprising, given the differences in geometry of both rotor and structure; indeed the difference in structure (single stanchion versus monopile+crossbeam for SeaGen) may account for the variance in the central peak.Such conformance across a range of scales and models is encouraging, given that levels of upstream turbulence have been found to influence tidal turbine performance [74]; the turbulence profiles in our rectilinear channel reflect a satisfactory approximation to a generic channel, but they do not reflect any particular tidal site.
The key finding of this paper is that the support structure has a discernible effect on the flow upstream, the downstream wake, and turbine behaviour, for a contra-rotating, dual rotor tidal turbine.
The behavioural changes manifest themselves as regular oscillations in the power output at twice the rotor rotational frequency, and constitute 4% of total output.Furthermore, these fluctuations affect the near-wake of each rotor, where these oscillations are evident in the turbulence spectra, and so undoubtedly causing regular variations in structural loading.Inclusion of the structure in simulations also produces deeper, more persistent wake deficits, as well as substantially higher levels of turbulence downstream, well into the far wake.These results are somewhat dependent on the SeaGen model design, and differences can be expected between it and more conventional single-rotor devices.These matters bear further investigation, particularly for tidal turbine arrays.Whilst array layouts were beyond the scope of the research presented here, it is recommended that future work should investigate the consequences of the resultant complicated wake flows, for the electrical and mechanical systems of tidal turbines.In summary, the present study confirms that if the support structure and individual rotors are not resolved or parameterised appropriately within a numerical model, then the results of such simulations should be held in doubt, from the interrelated perspectives of turbine reliability, performance and fluid dynamics.

Figure 1 .
Figure 1.The idealised tidal channel, showing boundary conditions and dimensions.

Figure 2 .
Figure 2. Vertical profile of normalised streamwise Reynolds stress at the inlet, as a function of height.

Figure 3 .
Figure 3. Velocity magnitude distribution of the empty channel at t=1800s.

Figure 4 .
Figure 4. Turbulent flow in an empty rectilinear channel, showing vertical profiles for a) time-averaged velocity magnitude, b) spatially-averaged velocity magnitude versus ideal log and power law profiles, and c) calculated turbulent intensity.Units for x are in D, the diameter of the turbine rotors (16m).x = 0 is where turbine rotors are to be placed, 250 m downstream of the inlet.

Figure 5 .
Figure 5. Close up of a horizontal slice through the mesh at hub-height, showing the mesh resolution around and downstream of the cylinder.The resolution of the mesh at the cylinder's surface is 0.25 m; this increases to approximately [2 m, 2 m, 1 m] over a distance of 10 m upstream, 20 m cross-stream, and 250 m downstream.

Figure 6 .
Figure 6.Horizontal slice through the velocity field at z=16 m and t=900s, showing the full extent of the wake behind the cylinder at the scale of the channel.

Figure 7 .
Figure 7. Horizontal slice through the vorticity field at z=16 m and t=900 s, showing the vorticity generated by flow past the cylinder.
have been used to represent the rotor, whereby the blades themselves are not resolved, but the forces exerted by them on the fluid are still present.In actuator line theory, the blade is replaced by the actuator line, and the forces are spread spatially via a normalised Gaussian distribution function to become body forces.In this implementation, we use a two-dimensional Gaussian distribution function, which is described below.The code for the model has been released as open-source, under the Lesser GNU Public License, version 2.1 [56].

Figure 8 .
Figure 8.The time-series of the lift force on the cylinder is shown in a), with the FFT of lift in b) showing a pronounced peak at about 0.22-0.225Hz.

Figure 9 .
Figure 9. Time-averaged power (blue) and thrust (red) for a single rotor of SeaGen, as a function of mean hub-height flow speed u H .The solid lines represent data sourced from Douglas et al. [62] and Fraenkel et al. [52,60]; the squares and triangles represent simulation results.

2 .
They were positioned in front of the solid nacelle structures, with the first rotor T1 at [247 m, 87 m, 16 m], and the second rotor T2 at [247 m, 113 m, 16 m].The blades were oriented on T1 and T2 so that the lift-induced torque would cause clockwise and anti-clockwise rotations respectively, thus forming the contra-rotating pair shown

Figure 10 .
Figure 10.The turbine structure: a) a perspective view illustrating the main components, and b) a head-on view indicating important dimensions.

Figure 12 .
Figure 12.Horizontal slice at hub height through the mesh for the dual rotors only case.Mesh resolution increases sharply toward the rotors, and reduces gradually downstream to the resolution used in the empty channel case.

Figure 13 .
Figure 13.Vertical profiles of time-averaged velocity magnitude and turbulence intensity for rotors, structure, and rotors with structure.The units for x are D, ie. the rotor diameter (=16 m).

Figure 14 .
Figure 14.Horizontal transects of time-averaged velocity and turbulence intensity at z=16 m for each full-scale case.Units for x are D.

Figure 15 .Figure 16 .
Figure 15.Instantaneous snapshots of the velocity field for rotor and rotor+structure cases, at the end of the simulation.

Figure 17 .
Figure 17.Diagnostic output from both rotors, T1 and T2, showing power, and α, the weighted mean angle of attack (see Creech et al. [4] for details) for each.The dashed lines indicate mean power output.

Figure 18 .Figure 19 .
Figure 18.Rotor-averaged FFT plots of the power time-series, for (a) rotors only, and (b) rotors+structure.There is a pronounced spike at 0.4 Hz for the simulation that includes the support structure.

Table 1 .
General specifications for modelled turbine.

Table 3 .
Comparison of power output values

Table 4 .
Mean rotational frequencies of the rotors in each simulation case.

Table 5 .
High frequency power fluctuation ranges for rotor and rotor+structure cases.CaseFluctuation range ( P max − P min ) (kW) % mean rotor power