Impact of Environmental Turbulence on the Performance and Loadings of a Tidal Stream Turbine

A large-eddy simulation (LES) of a laboratory-scale horizontal axis tidal stream turbine operating over an irregular bathymetry in the form of dunes is performed. The Reynolds number based on the approach velocity and the chord length of the turbine blades is approximately 60,000. The simulated turbine is a 1:30 scale model of a full-scale prototype and both turbines operate at very similar tip-speed ratio of λ ≈ 3. The simulations provide quantitative evidence of the effect of seabed-induced turbulence on the instantaneous performance and structural loadings of the turbine revealing how large-scale, energetic turbulence structures affect turbine performance and bending moments of the rotor blades. The data analysis shows that wake recovery is notably enhanced in comparison to the same turbine operating above a flat-bed and this is due to the higher turbulence levels generated by the dune. The results demonstrate the need for studying in detail the flow and turbulence characteristics at potential tidal turbine deployment sites and to incorporate observed large-scale velocity and pressure fluctuations into the structural design of the turbines.


Introduction
Tidal stream energy has been identified by the European Commission [1,2] as one of the five most promising renewable energy resources. A unique feature of tidal energy is its periodic nature that allows to make reasonably accurate predictions of how much energy can be extracted from a single or array of tidal stream turbines (or marine hydrokinetic turbines) deployed at a chosen location over a given period in time, e.g. for the 20-25 years they are designed for [3,4]. This high predictability in time means that tidal energy is a pseudodeterministic resource and thus ideal to be included in the energy mix with other stochastic renewable energy resources, such as wind and solar energy. Thorsten Stoesser stoesser@cardiff.ac.uk 1 Hydro-environmental Research Centre, School of Engineering, Cardiff University, The Parade, CF24 3AA, Cardiff, UK Despite the advantage of being a predictable source of energy the technology of marine hydrokinetic turbines is considered to be still in their infancy. Their high levelised costof-energy (LCOE) is an issue that needs to be overcome in order to increase their pace of development to become a mature technology. Vazquez et al. [5] acknowledged that a large portion of their LCOE are associated with the turbine's initial one-off manufacturing costs (capital expenditure, CAPEX) and their ongoing maintenance costs (operating expenditures, OPEX). A strategy to reduce the LCOE is to reduce periodic and unscheduled maintenance operations for instance by enhancing the structural response of the device to extreme or fatigue loads and to avoid sudden structural collapse [6]. At present, the knowledge of the long-term effects of hydrodynamic loads on tidal turbines is limited with only a few studies that focused on these matters [7]. Until now, the uncertainty about expected extreme loads has been addressed by turbine manufacturers by using over-estimated safety coefficients which consequently increases the LCOE [8].
The most relevant hydrodynamic features of the marine environment tidal turbines are subjected to are summarised in Fig. 1. Among these, the ones that challenge the turbine's structural integrity the most are energetic turbulent structures of different nature depending on whether they are generated from irregular seabed bathymetry or large-scale turbulence that is convected with the turbulent free-stream flow. For instance, Horse Rock in Ramsey Sound, a submerged obstacle the size of which is of the same order of magnitude as the water depth, is situated in close proximity to a commercial tidal turbine deployment site [9], and turbulent flow structures generated by Horse Rock are expected to cause significant hydrodynamic loads on the turbine's rotor blades. Additionally, waves can cause important velocity fluctuations near the turbine [10]. In addition, operating tidal turbines introduce high levels of turbulence into the flow generating their own energetic large-scale flow structures such as tip vortices due to the shear caused by the blades' tips in the free-steam or the meandering turbine wake that dominates the far-wake featuring large-scale turbulence and substantial velocity fluctuations [11,12]. In the deployment of turbines in arrays, they need to be properly arranged as adverse turbine-to-turbine interaction can trigger fatigue  1 Sketch of the harsh environmental conditions which tidal turbines are subjected to effects and cause a significant reduction of the overall power generation [13,14]. These effects introduced by the turbines alter the marine environment and hence they need to be understood and quantified. Polagye et al. [15] classified the hydrodynamic effects of tidal stream turbines by those relevant in the near-field: water quality (stratification), sediment transport (scour and deposition), and water velocity field (turbulence, acceleration and cavitation); and in the far-field: water quality (dissolved oxygen, nutrients, turbidity), sediment deposition/erosion in intertidal zones, or impact on benthic communities.
During the first stages of the development of tidal turbines a number of experimental campaigns were carried out which focused on understanding their performance [16], e.g. how the performance of the turbine varies depending on its proximity to the free surface, the oncoming velocity or its rotor yaw angle [17][18][19]. The first tidal current turbines were manufactured considering wind turbine rotors as reference point as this is a well-developed technology, although tidal turbine rotors adopted shorter and thicker blades as water forces are much higher than those of wind [6]. The preliminary experiments evidenced that tidal current turbines are a promising novel technology with the tested designs reaching power coefficients above 40%, i.e. similar to their wind counterparts.
A series of studies analysed the effects of turbulence on tidal turbines which were performed in experimental flumes. Mycek et al. [20] compared the performance and wake downstream a Horizontal Axis Tidal Turbine (HATT) considering two turbulence intensity values. They found that for a turbulence intensity (I ) of 3% the wake velocity deficit was not recovered even at 10 diameters downstream whilst for I = 15% the wake velocity at 6 diameters downstream showed no signature of the upstream turbine. Mycek et al. [21] studied two HATTs placed in-line at different distances between them again with turbulence intensities of 3% and 15%. Their results showed that the performance of the downstream turbine was largely affected by the upstream turbine for the low turbulence intensity case regardless of the inter-turbine spacing as a result of a reduced wake velocity recovery. For the high-turbulence case the second turbine was capable of achieving a performance similar to the upstream turbine evidencing that greater turbulence in the ambient flow reduces the negative interaction of multiple turbines. Blackmore et al. [22] investigated several turbulent flow conditions varying the parameters of turbulence intensity and turbulent length scale, and found that the latter had a larger influence on the device's performance and loadings. They also observed that for those cases where the turbulent length scale, which represents the size of the most energetic flow structures, is significantly smaller than the blade length of the turbine the influence of the turbulence is small. Whereas if the turbulent length scale is greater than the turbine's diameter, turbulence structures will cause significant pressure fluctuations acting on the entire device. In view of these experimental campaigns, two main outcomes have been found: (i) turbulence has an effect both on the turbine's performance and its structural loadings and (ii) the turbulence in the ambient flow should be characterised in terms of the turbulence intensity and the turbulent length scale.
At present, a significant challenge for tidal turbine developers is improving the structural design of tidal turbines, e.g. blade thickness or material stiffness, ensuring their long-term integrity avoiding sudden structural failure [6,7,23]. An intrinsic effect of turbulence together with the periodic nature of the tides is the repetitive action of dynamic structural loads which can cause fatigue failure and reduces their survivability. Generally, the industry assumes a lifespan of tidal turbines of 20 years during which they have to stand about 10 7 to 10 8 fatigue cycles [24,25]. Hence, the design of tidal current turbines requires precise knowledge of fatigue loads in order to improve their ability to resist maximum loads (i.e. ultimate limit state).
Results from experimental scale-testing are reliable and these provided invaluable insights into the fluid-structure interaction which enhanced the development of the tidal stream turbine technology. However, main drawbacks of tidal turbines are: (i) relatively high costs that makes extremely difficult to test turbine arrays as large facilities are required [26]; (ii) hydraulic flumes, featuring smooth flat beds which do not represent the conditions found at typical marine deployment locations and the limited width of laboratory flumes which induces a blockage effect that alters the turbine's performance and wake recovery behaviour if compared with boundless environments (e.g. estuaries or sea) [17,27]; and (iii) the experimental data is usually not detailed enough to provide a full understanding of the complex flow-turbine interaction [11,12].
The flow conditions and turbulence characteristics at tidal stream turbine deployment sites are found to be site-specific [28,29], and may be notably different from those in experimental flumes, demanding for detailed site-specific hydrodynamic studies before deploying tidal stream turbines [9]. Milne et al. [30] highlighted the heterogeneity and variability of turbulence properties at different potential deployment locations, with mean turbulence intensities ranging from 5 to 10% whilst instantaneous values can achieve turbulence intensity peaks of up to 20-25% [28]. Reproducing such complex flow conditions in hydraulic flumes is practically impossible, while numerical models can be used as a complementary tool that can overcome some of the experiments' shortcomings.
Blade Element Momentum (BEM) methods are widely used due to their fast calculation of the performance and loadings for a certain turbine design. However, they do not resolve the instantaneous flow field and thus are unable to provide trustworthy information about the effect of seabed turbulence, waves, misaligned flow, or structural behaviour of the turbine under extreme loads. Alternatively, Computational Fluid Dynamics (CFD) is a powerful tool capable of accounting for the previous BEM methods' limitations, and hence, recently, has become an attractive alternative approach to study tidal turbines providing valuable information of the turbine's performance, hydrodynamic loadings and wake. To date, the most commonly used CFD approach to study tidal turbines are those based on the Reynolds Averaged Navier-Stokes (RANS) equations mainly due to its relatively affordable computational cost [31]. RANS-based CFD models compute the time-averaged velocity field and they model all the effects of turbulence. This may be considered suboptimal when computing the turbulent flow around tidal turbines as the flow-turbine interaction appears to be dominated by a wide spectrum of turbulent scales. At the expense of higher computational demand, Large-Eddy Simulation (LES) resolves the large-scale turbulence in the flow, such as those dominating the turbine's wake, while modelling only small-scale turbulence [32,33]. Afgan et al. [34] analysed the results of LES and RANS simulations of a horizontal axis tidal turbine and reported that the method of LES is capable of reproducing quite well the wake downstream of the turbine whilst RANS tends to under-estimate the turbulent kinetic energy in the wake and its velocity recovery.
The computational cost of LES of tidal turbines is affected by the numerical technique adopted to represent the moving geometry. The common strategies used to represent tidal turbines are Geometry Resolved (GR) approaches and actuator techniques, such as actuator disk or actuator line [35]. The most frequent GR method is based on body-fitted meshes [34] in which the turbine geometry is embedded into the fluid domain and it is moved using sliding-mesh techniques [36,37]. Body-fitted GR simulations feature a great accuracy to predict the physics of the flow around tidal turbines although their main inconvenience is the large computational cost due to variable re-allocation processes at each time step. Actuator models are simplified approaches based on the BEM method whose main advantage is their low computational requirements. Actuator line model discretises the turbine blade as a set of discrete points, each representing a section of a blade [38], and assigns tabulated values of lift and drag coefficients to each of these points. This approach has proven to provide a fair estimation of the power produced by one or several turbines in an array and is capable of resolving flow phenomena such as tip vortices generated by the turbine's rotor blades.
Recently, two other GR approaches are used more frequently due to their reduced computational cost compared to body-fitted models: Chimera or overset methods [39,40] and Immersed Boundary (IB) methods [11,12,26]. Both methods use an additional mesh decoupled from the main fluid mesh in order to represent the geometry of the moving turbine, and which is an Eulerian grid in the Chimera approaches and a Lagrangian grid in the IB methods. Kang et al. [35] simulated a HATT adopting a GR approach with the IB method, actuator disk and actuator line models and outlined that only the GR-IB model captured the two main flow structures present in the wake behind tidal current turbines, namely tip vortices and inner vortex or hub-induced wake. The actuator line also accounted for the presence of tip vortices while it failed in predicting the inner wake, while the actuator disk was unable to capture any these flow features. Ouro et al. [11] showed that GR-IB methods can provide accurate results of power performance and relatively accurate estimations of the structural loadings at the root of turbine blades.
This paper presents a LES of a HATT operating in nature-like turbulent flow conditions with the aim of quantifying the impact of bathymetry-induced turbulence on the performance and structural loadings of the turbine. The turbine under investigation is a small-scale version of a tidal turbine turbine prototype that has been deployed at Ramsey sound and was studied numerically in the laboratory under (scaled) geometric conditions similar to the deployment site, however operating above a flat-bed and not above an irregular seabed [11]. This paper will also provide further details of the resulting hydrodynamics and turbulence in the near and far wake of the turbine when operating over irregular bathymetry. The numerical framework and details of the turbine's design and the setup and validation of a precursor simulation performed to generate realistic turbulent inflow conditions are described in Section 2. Results of the instantaneous flow field upstream and downstream the turbine, hydrodynamics and structural loadings, time-averaged wake pattern, and nature of large-scale turbulence are presented and discussed in Section 4. Finally, conclusions are presented in Section 5.

Governing equations and their discretisation
The method of Large-Eddy Simulation (LES) is employed to compute the flow around a horizontal axis tidal stream turbine using the in-house code Hydro3D, which has been successfully applied to the analysis of complex problems such as flow over vegetation [41,42] or flow over dunes [43]. Hydro3D resolves the spatially-filtered three-dimensional incompressible Navier-Stokes equations which read, where u and p are the fluid velocities and pressure, ν and ν t are the viscosity of the fluid and that added by the sub-grid scale model, and f is the source term from the IB method that enforces the no-slip condition at the location of the IB markers [44].
First and second order velocity derivatives are approximated using fourth-order central differences. The fractional-step method is used to advance in time these equations [45] using a low-storage third-order Runge-Kutta algorithm e.g. [46]. The Poisson-type pressure equation is solved using a multigrid approach which provides a scalar field onto which the predicted velocities are projected in order to fulfill the divergence-free condition. The sub-grid scale (sgs) stress tensor is approximated using the WALE model from Nicoud and Ducros [47] which does not explicitly calculate the sgs-viscosities near the wall and hence adapts well to its use in conjunction with moving IB bodies [48]. A direct forcing Immersed Boundary (IB) method originally developed by Uhlmann [44] is implemented in Hydro3D, in which the representation of solid boundaries within the computational domain is accomplished by the discretisation of the geometry in a set of regularly spaced Lagrangian markers constituting the solid domain [11,48]. This Lagrangian grid is detached from the fixed fluid mesh and the interpolation of quantities between fluid and solid domains is accomplished via delta functions [49]. The immersed boundary method is a compromise between numerical accuracy and stability, as moving boundaries in highly turbulent flows pose a challenge to any numerical method. And whilst this method is unable to resolve explicitly the viscous part of the laminar boundary layer developing over the hydrofoil, and hence is unable to resolve the laminar-to-turbulent boundary layer transition, it nevertheless proved to be accurate enough when employing adequate grid spacings and time steps. This LES-IB method has been validated in complex fluid-structure interaction applications such as vortex-induced vibrations [50], dynamic stall in pitching airfoils [51], or flow around vertical axis tidal stream turbines [48]. Most relevant to this study is the validation of the current method for the horizontal axis tidal turbine operating over a flat bed by Ouro et al. [11], which is exactly the same turbine as the one to be studied herein.
Hydro3D uses a domain decomposition technique to divide the computational domain into rectangular sub-domains on which the equations are solved using multiple processors and these communicate via Message Passing Interface (MPI) protocol. Cevheri et al. [52] implemented a local mesh refinement method that allows to refine the fluid mesh in certain areas of interest, e.g. near an immersed solid, while coarser grid resolution can be used far from these thereby reducing the computational load of the LES.

Tidal stream turbine implementation
The simulated horizontal axis tidal turbine is the discrete model of a 1/30th-scale prototype comprises three blades and is presented in Fig. 2. The rotor radius R is 0.2 m measured from the centre of the rotor and the radius of the hub r is 0.05 m. A model of the turbine was tested in Cardiff University's hydraulic flume of rectangular cross-section and flat bottom. The Reynolds number based on the chord length, c and freestream velocity, U 0 , is Re c ≈ 60, 00, which is the same as in Ouro et al. [11].
The extracted power of a tidal stream turbine is often obtained for a wide range of rotational speeds or tip speed ratios respectively calculated as λ = R/U 0 , where is the turbine's rotational speed. The optimum rotational speed of a tidal turbine is achieved when power extraction is maximum and in the simulations presented here the tip speed ratio is λ ≈ 3, which is very similar to the expected tip speed ratio of the prototype. The dimensionless performance of the device is computed as the ratio of the power generated by the turbine, P T to that available in the water, P W , represented by the power coefficient C P and defined here as, where Q is the torque generated by all blades, ρ is the fluid density and A = πR 2 is the area swept by the rotor. The flap-wise bending moment is the main contributor to the structural stresses at the blade-hub juncture as its values are generally 3-4 times higher than those from edge-wise bending moments [11]. The forces in the streamwise direction are denoted thrust forces (T ) and are responsible for the flap-wise bending moment M θ T whilst the tangential force Q generated by the blades is responsible for the edge-wide bending moment M θ Q . Following [11], the coefficients of these two structural moments are defined as, Even though the chord Reynolds number of the full-scale prototype is approximately 1-2 magnitudes greater, the lift and drag coefficients (and hence the dimensionless coefficient of power and bending moments) of the simulated turbine are expected to be very similar to the prototype due to the fact that the approach flow is highly turbulent and thus the effects of turbulent boundary layer transition are deemed marginal.

Generation of Realistic Inflow Turbulence
A precursor simulation is performed initially with the goal to generate a "realistic" fully developed turbulent flow field that is dominated by large-scale, bathymetry-induced turbulence, to be used as inflow condition to the main turbine simulation. As depicted in Fig. 3, the precursor domain comprises one dune and has a streamwise extension of x/k = 20, is  4k) whilst the variable bathymetry of the dune leads to a maximum water depth of 5k found in between 2.5 < x/k < 4. There are two reasons for using this bathymetry: (i) it will lead to significant turbulent structures, similar to those found in real turbine deployment sites [29], and (ii) this flow has been studied extensively in the past [43] and hence the current simulation can be validated with experimental and numerical data, albeit Reynolds numbers are slightly different (Turbine simulation: Re k ≈ 10, 000; Laboratory experiment: Re k ≈ 5, 000, where Re k is based on roughness height and bulk velocity). The flow in the streamwise direction is driven by a pressure gradient which keeps constant the flow rate set to a mean inlet velocity value U 0 = 0.76 m/s. No-slip conditions are set at the side-walls of the domain and the free-surface is treated as shear-free rigid-lid. The boundary conditions are exactly the same as those used in a previous flat-bed channel flow precursor simulations of Ouro et al. [11] so that the wakes of and loadings on a turbine operating over a dune can be compared to the same turbine operating over a flat bed. The dune is represented numerically using the IB method, whose accuracy to predict the flow over a bed of dunes has been proven by Grigoriadis et al. [53]. Cross-sectional planes of instantaneous velocity at the end of the precursor domain are stored and used later as inflow conditions for the main (turbine) domain.
The mesh resolution is kept uniform in the whole domain with x = 4 · 10 −3 m, and consists of 800 × 300 × 200 grid points in the x-, y-and z-directions respectively, yielding to a total of 4.8 · 10 7 fluid elements. The domain of the precursor simulation is divided into 900 sub-domains which are run on 100 Intel Xeon E5-2670 cores using the Supercomputing Wales' high-performance computer. A variable time step is adopted with a fixed CFL value of 0.3 in order to guarantee numerical stability which yielded an average time step of t * = tU 0 /k = 4.2 · 10 −3 . The simulation is initially run for 60 flow-through eddy turn-over times (t e = h/U 0 , where h = 4k), which is approximately 12 flow-through times, until the flow is fully developed. First order statistics are then collected for 250t e before starting to collect second order statistics, and the simulation is continued for another 150t e . The resolution in wall-units is z + = 5.8 which is calculated from the resulting friction velocity u * obtained from the pressure drop along the streamwise direction [54]. The value of z + is lower than 11 justifying the use of no-slip conditions at the side-walls and also along the dune geometry. Here mean (or root mean square) velocity fluctuations are computed from rms(u ) = ((U − u) 2 ) 0.5 , where U stands for the time-averaged velocity and u is the instantaneous streamwise velocity. First-and second-order flow statistics at the profiles L1 to L6, indicated via straight vertical lines in Fig. 4, are presented in Fig. 5 allowing a more quantitative assessment of the precursor simulation and comparisons with experimental data and the LES results of Stoesser et al. [43]. It is worth mentioning that the current simulation feature a larger U 0 and a greater water depth than in [43], which were chosen to match the turbine experiments and simulations of [11]. Figure 5a shows that the present LES nevertheless exhibits good agreement with the experimental and LES data of [43] for normalised time-averaged streamwise velocities (U/U 0 ) up to a water depth of z/k ≈ 3 whilst above this elevation the U -velocity profile predicted by the current LES deviates from the experimental one. These differences are due to the water depth and the presence of side-walls in the current simulation. The aspect ratio (width-to-depth) of the channel is approximately 2 and hence significant secondary currents develop, with the consequence that there is upflow near the sidewalls and downflow in the centre of the channel [55]. The latter pushes the high-momentum fluid at the water surface underneath it and hence there is a velocity dip near the water surface. At profiles L2 to L5, the present LES U -velocity distribution agrees well with the experimental and LES data between 0 < z/k < 1 showing that the phenomenon of flow separation at the recirculating area and boundary layer recovery along the dune stoss side are well captured.
Normalised streamwise (rms(u )/U 0 ) and spanwise (rms(v )/U 0 ) mean velocity fluctuations presented in Fig. 5b show that the results from the present LES match the distribution and values of these turbulent intensities of the reference data at every profile location. In the distribution of vertical mean velocity fluctuations (rms(w )/U 0 ) for the present LES follows the data from [43] at those heights below z/k = 2 coinciding with the areas of flow reversal and boundary layer development (see Fig. 4). Over z/k = 2, the distribution of rms(w ) is slightly different between the compared results as a result of the secondary recirculation due to the use of side-walls and different water depths. Overall, the profiles suggest that the turbulent flow field over the dune(s) is reproducing realistic turbulence typical of such flows and that the present mesh resolution is capable of providing an accurate representation of the complex, highly-turbulent flow over the dune(s). The accurate prediction Also plotted are LES (black line) and experimental data (symbols, for rms(u ) in (b)) from [43] of the velocity fluctuations is crucial for the present case as they are the result of instantaneous turbulence structures. These will be creating "realistic" hydrodynamic loadings on the turbine rotor's blades.

Tidal Turbine Operating Above a Train of Dunes
The HATT operating above a train of dunes is simulated using the full computational domain depicted in Fig. 3, in which two dunes in series comprise the natural bathymetry. The turbine is placed at x hub (x/k = 0, y/k = 3.75, z/k = 3.1) which is close to the crest of the second dune which is where the largest streamwise velocities are expected to occur and which, at this location, matches the water depth of the comparative flat-bed study. The grid is uniform in the three spatial directions and local mesh refinement using x = 0.001 m is employed near the turbine to achieve a fine grid resolution whereas far from it the mesh size is x = 0.004 m, i.e. the same resolution as in the precursor simulation. The domain has a total of 2.41·10 8 fluid elements, and is divided into 2250 sub-domains which are occupying 456 processors on the Supercomputing Wales high-performance computer. The simulation runs for 800,000 iterations for which 760,000 CPU-hours are needed.
A time step of t * = tU 0 /D = 1.71 · 10 −6 is chosen as determined by the time step sensitivity study of Ouro et al. [11], and this is fixed during the simulation to maintain smooth forcing of the immersed boundary method avoiding oscillations and instabilities [48]. The combination of a relatively small time step and fine grid resolution makes the current LES computationally very expensive. Thus only the case with the turbine operating at maximum efficiency is performed, i.e. λ = λ opt = opt R/U 0 , and which is set based on the results of Ouro et al. [11]. Inlet planes from the precursor simulation are fed into the turbine domain as represented in Fig. 3, whilst a convective boundary condition is set at the outlet. A rigid lid is applied at the water surface, which worked well in the previous comparative study [11], and water level increase upstream of the turbine due to flow bloackage are compensated by the presence of a local pressure gradient between upstream and downstream of the turbine. Initially, the simulation is run for 2 flow-throughs before inserting the turbine into the domain. First-order statistics are collected after the device has rotated 4 full revolutions to ensure that the flow around the turbine is fully developed. Averaging of second-order statistics starts after another 4 revolutions, i.e. after the 8th revolution, and is performed for the following 13 revolutions. A total of 21turbine revolutions is simulated. Figure 6 presents iso-surfaces of instantaneous negative pressure fluctuations (p = p − P , where p and P are the instantaneous and time-averaged pressures,respectively) and these are coloured with time-averaged streamwise velocity. Minima of the pressure fluctuation are reliable indicators for vortical motion (e.g. [43,53,56]), and hence allows visualising large-scale turbulence structures of this flow. The flow separates at the crest of the dune and generates a shear layer which, followed by Kelvin-Helmholtz (KH) instabilities, leads to the formation of elongated spanwise vortices denoted as KH rollers. These can grow up to the size of the dune's height and due to the high streamwise velocity are subsequently convected downstream [56]. Along the boundary layer recovery over the stoss side of the dune, Hair-Pin (HP) vortices are similar to those observed by Stoesser et al. [43], and these move upwards interacting with the overlying high-momentum flow and eventually forming Vertical Rollers (VR) that approach the turbine as observed in Fig. 6. Other secondary flow structures, such as streamwise-oriented (S) vortices, are also developed.

Instantaneous flow field
The instantaneous flow field behind the turbine is depicted in Fig. 7 in a longitudinal plane at y = y hub , at half width of the domain, with contours of uand v-velocities and instantaneous turbulent kinetic energy, k . Figure 7 allows to identify the main turbulent structures generated by the turbine and the interaction between the wake of the turbine and the near bed flow structures in the trough downstream of the turbine. A recirculation region behind the hub is identified by the streamwise velocity contours of Fig. 7a. Two low-momentum areas are marked in Fig. 7a and these correspond to the recirculation in the dune's trough and in the hub's wake in the region x/D = 1.5 − 4. There is an area of high-momentum flow between the turbine's lower tip and the dune crest. This high-velocity area extends until x/D = 2 − 3 which is where the interaction between the turbine's and dune's wakes takes place and which is marked as an oval in Fig. 7b. Between the turbine's top tip and the free-surface there is another region of high-velocities that interacts with the turbine's wake.   Fig. 7b. These energetic flow structures are generated by the shear induced by the blades' tip into the incident flow. The blades rotate with positive v-velocity when moving above the hub's centre (z > z hub ) while with negative v-velocity for z < z hub . This motion induces the fluid flow to react with an opposite distribution of transversal velocities with negative values over the hub's centre and positive below it. Contours of v-velocity also show the interaction between the wakes generated by the turbine and the second dune. The mixing of the wakes starts approximately after x/D = 2, and this is more noticeable between x/D = 2.5 − 4.0. Noteworthy is the fact that the turbine wake does not exhibit the same characteristics as the wake downstream of a tidal stream turbine placed on a flat bed, where the far wake starts to exhibit a spiraling motion e.g. [12] and [11].
Contours of the instantaneous turbulent kinetic energy (k ) are presented in Fig. 7c and these depict secondary turbulent structures derived from the interaction of the wakes. The area of greatest k is found in the dune's trough as flow reattachment triggers turbulence structures of high turbulence intensity levels (see Fig. 5). The operating turbine introduces significant turbulent kinetic energy by the blades, their tips and the hub.
A three-dimensional view of the turbine-generated turbulence structures is presented in Fig. 8 using iso-surfaces of λ 2 to educe coherent turbulence together with contours of relative pressure, p, in a longitudinal plane in the channel centre. λ 2 is the second eigenvalue of the sum of the symmetric and antisymmetric parts of the velocity gradient tensor and a vortex core is identified by an isosurface of a negative λ 2 . The signature of the tip vortices is clearly visible in the form of low-pressure areas (marked as 1c, 2c etc.) as well as the hub's wake, which starts meandering downstream of the turbine as identified from isosurfaces of a chosen negative pressure fluctuation (very similar to the one shown in [11] and not shown for brevity). The tip vortex denoted as "1a" corresponds to the latest vortex generated by blade 1, or tip vortex denoted "2c" is one that has been shed by blade 2. Shortly downstream of the turbine, bottom tip vortices break down and loose coherence due to their interaction with the ambient flow as just vortices 1a, 2a and 3a are visible from the pressure contours, and vortex 1b is also appreciated but appears to be segregated into 2 cores. Tip vortices near the water surface remain coherent during their advection downstream although the spacing between consecutive vortices is irregular due to their interaction with the turbulent ambient flow, e.g. separation between vortices 2b-3b is larger than between 1c and 2c. The λ 2 isosurfaces show that the characteristic roller shape of the tip vortices [11] is lost in vortex 1b whilst vortex 3a (which is closer to the rotor) is smoother and more compact. Such a loss of coherence is due to their interaction with the turbulent ambient flow. Chamorro et al. [57] visualised experimentally the wake of a HATT using PIV records and also observed that tip vortices feature irregularly spaced patterns in highly turbulent flow.

Time-averaged flow field
The time-averaged flow field in a longitudinal plane is presented in Fig. 9 with contours and vectors of normalised streamwise velocity (a) and turbulent kinetic energy k (b). Streamwise velocities increase on the stoss side of the upstream dune with the highest velocities found at the dune crest just upstream of the turbine. The turbine extracts energy contained in the approach flow inducing an area of low momentum immediately downstream of it. Outside of the turbine's swept area the flow is accelerated generating high-momentum zones above and below the turbine's blade tips, the locations of which are indicated with dotted lines in Fig. 9. There is a low-pressure area in the dune's trough and the downwards pointing velocity vectors between 1 ≤ z/k ≤ 2 at 1 ≤ x/D ≤ 2 in Fig. 9a indicate that the fast-moving flow found below the turbine's swept area starts to fill this area.
Three wake regions are discerned in Fig. 9: The near-wake extends between 0 < x/D < 1.5 featuring a short recirculation behind the hub. The low-velocity region covering 1.5 < x/D < 4 (blue areas in Fig. 9a) is the characteristic signature of the mid-wake region. Compared to the flat-bed channel results of Ouro et al. [11], see Fig. 10, the presence of an irregular seabed causes the wake to be asymmetric in the vertical direction. The flow decelerates above the dune's trough and the boundary layer deficit downstream of the dune's recirculation zone begins to be filled with high-momentum fluid and hence the turbine's wake is skewed towards the dune trough. Figure 10 plots streamwise velocity contours of both flat-bed (top) and dune-bed (bottom) cases showing that the low-momentum, mid-wake behind the turbine is somewhat wider (due to the flow expanding in the trough) and slightly  Contours and vectors of the mean turbulent kinetic energy (k) are presented in Fig. 9b. Shortly upstream of the turbine at x/D = −0.5, the flow approaching the turbine features larger k values in the lower half of the water-depth, the signature of the turbulence coming from the upstream dune. In the near-wake region, the main regions of high k (red areas in Fig. 9b) coincide with the rotor's swept area as the blades introduce additional turbulence into the flow, in the form of trailing edge and tip vortices. Also, high areas of k are found in the recirculation area of the rotor hub. At 1D downstream of the turbine, the tip vortices and the shear layer above the dune's recirculation zone are responsible for significant levels of k as evidenced by the long k vectors. In the mid-wake the k levels are relatively low and tip vortices diffuse by action of the shear of the free-stream flow. Near the bottom of the channel the zone of high k in the dune's trough merges with the one due to the energetic tip vortices of the near bed blades. At the beginning of the far-wake at x/D ≈ 4, the interaction of the tip vortices with the ambient flow increases again the magnitude of k as observed from the vectors above z/k > 3, while there are also high levels of turbulent kinetic energy at 1 < z/k > 2 due to the interaction between the dune's and turbine's wakes. The two peaks of turbulence kinetic energy at x/D = 4 merge to a large area of high levels of k after x/D ≈ 6. Figure 11 plots contours of the mean turbulent kinetic energy of both flatbed (top) and dune-bed (bottom) cases. The most noteworthy difference is the significantly higher turblence behind the turbine operating over the dune bed and its extension into the far wake. Figure 12a presents profiles of the normalised time-averaged streamwise velocity of the flat-bed and dune-bed cases at selected locations in the centre-line of the channel. This allows more a quantitative comparisons between the two flows. The flow upstream of the turbine, at x/D = −1 is almost identical, however the dune-case features a slightly higher streamwise velocity and a more uniform distribution than the flat-bed case. The first profile downstream of the turbine (x/D = 1) is also quite similar, however slightly higher streamwise velocity is found underneath the rotor in the dune-bed case, the result of the slowly rising channel bed just upstream of the turbine (inducing an acceleration near the bed). The widening of the cross-section in the dune-bed case allows the flow to expand and streamwise velocities in the mid-wake of the dune-case are lower than in the flat-bed case x/D = 4. In the dune bed case the velocity deficit behind the turbine is almost recovered by x/D = 6, whereas the flow over the flat-bed has not fully recovered by x/D = 8. In contrast the streamwise velocity over the dune appears fully recovered by x/D = 8. Figure 12b presents profiles of the normalised streamwise velocity fluctuation of the flat-bed and dune-bed cases at selected locations in the centre-line of the channel. The streamwise turbulence upstream of the turbine, at x/D = −1 is significantly different: the dune-case approach flow having considerably higher streamwise velocity fluctuations than the flow over the flat bed, due to its bedform-induced (or rough-bed) turbulence. In the first profile downstream of the turbine (x/D = 1) the streamwise turbulence is quite similar between the two flows, except for the significant peak in the dune-bed case around z/k = 1, which signifies the location of the shear-layer between the dune's recirculation zone and the fast moving fluid above it. Streamwise turbulence levels remain high near the bed in the dune-bed case and low in the mid-wake x/D = 2, whilst in the flat-bed case the onset of wake meandering signified by high levels of rms(u ) is observed from x/D = 2 onwards. Further downstream of the turbine streamwise turbulence levels are similar (at x/D = 6) before the flow over the dune-bed features aforementioned higher levels of turbulence, and at x/D = 8 this is a remarkable 25% higher streamwise turbulence in the dune-bed flow compared to the flat-bed flow. Figure 13 depicts the distribution of the time-averaged streamwise velocity (left column), mean (k, centre column), and instantaneous (k , right column) turbulent kinetic energy in selected cross-sections located at x/D = 0.5, 1.5, 3.5 and 7.5 downstream of the turbine. In the plane closest to the turbine, the characteristic signature of the turbine near-wake is well appreciated with the circular low-momentum region with negative U/U 0 values behind the hub, and high levels of turbulent kinetic energy behind the hub and along the perimeter Shortly downstream, at x/D = 1.5, the circular shape of the turbine's wake is still visible from the U/U 0 contours that also show the turbulence originating from the dune's trough at z/D < 0.4 in the form of high-and low-speed streaks, previously identified by [43]. The time-averaged k contours suggest that the highest levels of k are found in the dune's trough meanwhile the circular turbulence of the tip vortices is more diffused than at x/D = 0.5 due the loss of coherence of the tip vortices as they interact with (and get sheared by) the highmomentum free-stream flow. This is better illustrated by the instantaneous turbulent kinetic energy contours on the left. The cross-section at x/D = 3.5 is located at the stoss side of the dune where the boundary layer recovery enhances the mixing of the turbine's wake with the highly turbulent ambient flow. The low-momentum region still coincides with the turbine's swept area (dashed lines in Fig. 13) and it has merged with the low-momentum area near the bed. Similar behaviour is observed for the area of highest values of mean turbulent kinetic energy, i.e. merging of turbine-and bed-generated turbulence. The contours of instantaneous k exhibit pockets of high k as a result of the turbine's wake interacting with the turbulent flow generated in the dune's trough. Further downstream of the turbine, at x/D = 7.5, the streamwise velocities appear to be more evenly distributed over the whole cross-section with only minimal signs of the presence of a turbine upstream of this location. The largest timeaveraged k levels are found within the turbine's swept area a signature of wake meandering and this is also visible from the instantaneous turbulent kinetic energy contours. Time-averaged quantities of the normalised streamwise velocity (U/U 0 ) and turbulence intensity (I (%) = (1/3(u + v + w ) 0.5 * 100) along the centreline of the channel at the hub's height, z = z hub are presented in Fig. 14a) and mean velocity fluctuations (rms(u ), rms(v ), rms(w )) and turbulent kinetic energy (k) are presented in Fig. 14b). Upstream of the turbine (x/D < 0), the values of U/U 0 are fairly uniform with U ≈ 1.2U 0 , i.e. approximately 20% higher than the bulk velocity which is owing to the use of sidewalls, hence velocities in the centre of the channel are higher than the cross-sectional average. In the approach flow the turbulence intensity is highest at x/D = −8 with a value of I = 20% and its magnitude diminishes progressively until a minimum of I ≈ 8% is attained at x/D = −2. The fluctuation of vertical velocity (w /U 0 ) is the main contributor to k, especially at −2 < x/D < 0, while velocity fluctuations in x− and y−direction are similar upstream of the turbine. The flow conditions upstream of the turbine at x/D = −0.5 are often adopted as reference values to predict/assess turbine performance. Here, these are U hub ≈ 1.3U 0 and I hub ≈ 14.5%.
The two dashed lines in Fig. 14a and b mark the boundaries of the three characteristic wake regions, as proposed from Fig. 9, and the profiles provide evidence of the chosen extent of each of these regions. The near-wake is characterised by the sudden drop in velocity and spike in turbulence followed by a significant drop-off in turbulence intensity/fluctuations. In the mid-wake, extending between 1.5 < x/D < 4, low streamwise velocities plateau and turbulence levels in this low-momentum zone are relatively constant. In the far-wake, the velocity recovers and the turbulence intensity levels increase again due to high-momentum fluid filling the wake which results in wake meandering and thus high levels of rms(u )/U 0 [35] are observed. The boundary layer recovers along the stoss side of the dune which results in increasing the values of rms(w )/U 0 .

Large-scale turbulence
The identification of large-scale flow structures can be accomplished via temporal autocorrelations of velocity fluctuations (ρ i ) which allows to determine the integral time and length scales of the turbulent structures present in the flow [28,58], and it is defined as: where i (= 1, 2, 3) refers to the velocity component, t denotes time, and τ is the time lag adopted to measure the phase correlation. The integral time scale, T u i is calculated by integrating the autocorrelation function from τ = 0 until when ρ i = 0, meanwhile the spatial length scale is calculated as L u i = T u i U , where U is the local time-averaged streamwise velocity. Figure 15a shows the autocorrelation function of the three velocity fluctuations obtained at location P1, which is just upstream of the turbine. The largest integral time scale is obtained for T u ≈ 0.17 s indicating that the correlation of the flow in the streamwise direction is more coherent than in the spanwise or transversal directions for which T v ≈ 0.084 s and T w ≈ 0.083 s. The integral length scale value is here L u ≈ 0.15 m and it provides an estimation of the dimension of the large-scale structures in the flow [31]. This value gives L u ≈ 0.75R in terms of rotor radius suggesting that the large-scale flow structures shown in Fig. 6 impact the blades unevenly and hence the instantaneous structural loading on each blade varies significantly (as is demonstrated in Fig. 17 below). Chamorro et al. [59] obtained a similar value of L u ≈ 0.7R indicating that the current case adopting the dune-like bathymetry represents realistic environmental turbulent flow conditions. In order to justify the assumptions made to compute T u i a second location is chosen, P3, to estimate integral time scales. The correlations of Fig. 15b show peaks from which the blade frequency can be inferred, evidence that temporal autocorrelations for these types of (shear) flows are reasonably reliable in estimating large-scale turbulence frequencies and lengthscales. The correlation of Fig. 15a shows a secondary peak at τ ≈ 2.1 s suggesting that a dune generated turbulence structure occurs approximately every two seconds. The anisotropy ratio rms(u ):rms(v ):rms(w ) at location P1, i.e. shortly upstream of the turbine, is found to be 1:0.64:0.88. These values are somewhat dissimilar to the 1:0.75:0.56 found at the Sound of Islay (Scotland) tidal site [28] or to two-dimensional channel flows [60] in which the ratio is 1:0.71:0.55. The present values reflect the influence of the duneinduced turbulence with the vertical velocity fluctuations dominating over those in the transversal direction. According to Mycek et al. [20], ratios of mean turbulence intensities are site-dependent and vary according to the dominating turbulence structures at a given deployment site, e.g. such as the presence of large obstacles in the bathymetry such as Horse Rock in Ramsey Sound [9].
The frequency associated with periodic events due to large-scale flow structures in the velocity time series in Fig. 17a is revealed via Power Spectral Density (PSD) from u− and w-velocities series obtained at four points, denoted P1, P2, P3 and P4 (see Fig. 8), and these are presented in Fig. 16. P1 is shortly upstream of the turbine, P2 and P3 are located at x hub + 0.3D and at x hub + 0.8D, respectively, at a distance of D/2 from the hub's centre, i.e. at the projected turbine swept perimeter that coincide with the path described by the tip vortices, and P4 is aligned with the hub's centre located at a distance of 1.8R downstream. At P1 a low-frequency peak in the u-velocity spectrum is appreciated, labelled as f L and f L = 0.45 s −1 , as a result of the periodic turbulent events generated by the dune upstream of the turbine, such as the vertical rollers depicted in Fig. 6. In the w-velocity spectrum there is also an energy peak at frequency f L although there is another peak with higher associated energy at a frequency of f s = 0.89 Hz, which is linked to bed-induced turbulent events occurring at lower frequency in time than those contributing for f L . Spectra at P1 feature a decay slope of -1 in the production range prior to the inertial sub-range that follows the −5/3 slope expected from homogenous turbulence, which is also commonly found in boundary layer flows [61]. The −5/3-region spans approximately one decade of frequencies and given the relatively high Re this is sufficient for an LES. After the inertial sub-range the slope is significantly steeper and this is due to the sgs-model dissipating energy efficiently.
The turbine introduces turbulent kinetic energy into the flow as result of its operation and the flow structures induced by the turbine blades, e.g. tip vortices, are responsible for increasing the spectral energy associated with frequencies in the energy production range. P2 and P3 are located in the near-wake downstream of the turbine at x = 0.6R and x = Fig. 16 Power spectral density of uand w-velocity signals collected at P1, P2, P3 and P4 for the HATT operating over a train of dunes 1.8R, respectively. The dominating large-scale structures at these locations are identified by the energy peaks occurring at a frequency coincident with the blade passing f b = /2π · N b = 5.75 Hz, and thus these peaks are associated with the tip vortices [62]. Note that the latter energy peaks exhibit higher spectral energy than that found at P1 linked to the low-frequency turbulent structures from the turbulent flow approaching the turbine. The turbine-induced turbulence changes the energy decay rate in the production subrange at P2 and P3 with a flatter slope when compared to the −1 slope found at P1, i.e. in the freestream flow. It is also observed that comparing the u-velocity spectrum at P2 and P3, the slopes at the latter location is steeper compared to that in the former which is closer to the turbine rotor. Jin et al. [58] observed a similar trend from flatter to steeper slopes in the energy cascade when moving from the near-wake to the far-wake of a horizontal axis wind turbine. Energy peaks at blade passing frequency harmonics, i.e. 2f b and 3f b , are observed at P2 especially for the w-velocity spectrum. However, as the tip vortices are convected downstream they interact with the ambient turbulent flow and lose coherence, thus at P3 the spectral energy associated with f b is smaller than that at P2 and energy peaks at f b harmonics are no longer observed. The breakdown of the energetic tip vortices due to the turbulent ambient flow was highlighted in Fig. 8. Turbulence properties in the middle of the near-wake are revealed via spectra at P4, which show an energy peak at f b most noticeably in the w-velocity signal. This is the result of the recurrent shedding of blade root vortices which are generated by the sharp end of the blades near the hub [11].

Impact of bathymetry-induced turbulence on hydrodynamic loadings
Survivability of tidal turbines relies on the ability of their material to sustain large load fluctuations due to sudden velocity and pressure variations in the flow field. Figure 17 presents the time history of the streamwise velocity upstream of the turbine obtained at P1 (located  Fig. 6. The instantaneous power and thrust coefficients drop markedly when these events pass through as depicted in Fig. 17b, as a result of the sudden drop in speed of the approaching flow. After each event, the power generation of the device is quickly recovered as a result of the following increase of streamwise velocity, suggesting significant load fluctuations on all turbine components. Noteworthy is the fact that a drop in streamwise velocities is associated with a sudden increase in vertical velocity (not shown for brevity), the result of turbulent ejection events which travelling on a upwards path towars the water surface [43]. Figure 17c and d present the bending moment coefficients for each of the blades comprising the rotor and highlights that each blade is affected differently during the turbulent Table 1 Peak, mean and range of the bending moments of the three blades of the turbine operating over a bed of dunes, and those from the flat-bed channel simulations with artificial turbulent inflows of Ouro et al. [11] Bending events. For instance, at time t/T ≈ 4.8, blade 1 is almost in vertical position and thus its tip is near to the water surface, whilst blades 2 and 3 are below half water depth and hence their tips are closer to the channel bed. The bathymetry-induced turbulence structures are more energetic near the bed and thus induce on blades 2 and 3 a larger impact than on blade 1. The results confirm that the size of the turbulent structures is significantly smaller than the turbine diameter as it was verified from the integral spatial length scale. Table 1 presents averaged/maximum values of structural bending moments for the three blades together with the results of Ouro et al. [11] for the flat-bed channel simulations. The "range" quantity, here calculated as the difference between maximum and minimum values, is a key parameter in the prediction of fatigue loads [7].
The computed time-averaged structural coefficients of the present simulation are very similar to those found by [11] for the cases of the turbine operating above a flat-bed channel subjected to different levels of turbulence intensity. The mean flapwise bending moment coefficient C T M θ = 0.1012 is only approx. 0.101 which is similar to the 0.0964 and 0.0975 values obtained for the I = 10% and 20% cases, respectively. Such small variation under different turbulence conditions follow the findings of Blackmore et al. [22] who stated that even though turbulence intensities in the flow vary considerably, mean bending moments stay within certain bounds, i.e. no great differences are expected irrespective of the turbulence intensity in the flow. An analogous trend is found for C Q M θ with a value of C Q M θ = 0.0259 for the present case with a bed of dunes whereas it is C Q M θ = 0.0230 or C Q M θ = 0.0235, respectively for two artificial turbulent inflow cases [11].

Conclusions
A horizontal axis tidal turbine operating over an irregular bathymetry has been simulated with the goal to study the impact of realistic turbulent flow conditions on the performance and loadings of the device. A precursor simulation featuring a single dune with periodic boundary conditions was adopted to validate first-and second-order statistics predicted by the large-eddy simulation which used immersed boundaries to represent the dune's geometry. The velocity field in a cross-sections of the precursor simulation were then fed into a second domain featuring two consecutive dunes with the operating turbine placed at the crest of the second dune.
The wake developed downstream of the turbine features three distinct regions that were identified from the flow characteristics of the mean velocity field and were found to be similar to the one observed over a turbine operating over a flat-bed channel. The analysis of the instantaneous velocity field revealed that dune-induced turbulence enhances wake deficit recovery. The mean streamwise velocity in the wake attained values close to the free-stream condition at a distance of 7-8D downstream of the turbine whilst in absence of bathymetry-induced turbulence the wake recovery is achieved further downstream. The fact that momentum is recovered faster in the presence of irregular bathymetry is encouraging in light of future design of turbine arrays as irregular seabed bathymetry is the norm at potential deployment sites. What follows is that it allows to deploy tidal turbines in closer proximity of each other and this will save time and money or will allow placing more turbines within a given area.
It is shown that the instantaneous flow contains energetic large-scale turbulence structures that are generated by the bathymetry, mainly in the form of roller and hairpin vortices and that these have a profound impact on the tidal turbine. In particular, significant fluctuations in the velocity and pressure fields were observed upstream of the turbine. These induce sudden drops in the turbine's instantaneous performance as well as large fluctuations in the hydrodynamic loadings on the blades, the latter being a potential risk of fatigue failure of the blades. Integral turbulent length scales were computed and revealed that the size of the most energetic large-scale structures appears to be smaller than the turbine radius which results in a variation of structural bending moments in each blade to every turbulent event.
This work gives further insights into the challenges tidal turbines face when operating in highly turbulent flows and in the presence of bathymetric seabed features such as dunes or rocks underlines the importance of considering velocity and turbulence data at future deployment sites of tidal turbines in their structural design.