As ocean current turbines move from the design stage into production and installation, a better understanding of localized loading is required in order to more accurately predict turbine performance and durability. In this study, large eddy simulations (LES) of tidal boundary layers without turbines are used to measure the turbulent bending moments that would be experienced by an ocean current turbine placed in a tidal channel. The LES model captures turbulence due to winds, waves, thermal convection, and tides, thereby providing a high degree of physical realism, and bending moments are calculated for an idealized infinitely thin circular rotor disc. Probability density functions of bending moments are calculated and detailed statistical measures of the turbulent environment are also examined, including vertical profiles of Reynolds stresses, two-point velocity correlations, and velocity structure functions. The simulations show that waves and tidal velocity have the largest impacts on the strength of bending moments, while boundary layer stability and wind speeds have only minimal impacts. It is shown that either transverse velocity structure functions or two-point transverse velocity spatial correlations can be used to predict and understand turbulent bending moments in tidal channels.

Deployments of renewable energy systems—including wind turbines, solar thermal collectors, hydroelectric dams, and many other technologies—have grown at a considerable pace over the past decade, and from 2005 to 2011 renewable energy capacity in the U.S. increased by 300%.1 In order to sustain this rapid increase, however, additional deployments and new sources of renewable energy are needed. To satisfy this need, energy from the ocean has shown considerable promise as a predictable, abundant, and sustainable resource, and it is estimated that ocean energy—comprising wave, tidal, current, and thermal energy—is capable of contributing an additional 400 TWh/yr, or 10% of the current U.S. energy needs, to the U.S. energy portfolio.1 

In order to increase the feasibility of large-scale ocean energy installations, further research is required to understand the turbulent environment in which ocean energy devices operate. Turbulent stresses and bending moments on ocean current turbines, in particular, are not well understood and may contribute to large off-axis loads that can, in turn, lead to unexpectedly short turbine lifetimes due to gearbox and other component failures. Such large off-axis loads arise due to the intermittent nature of turbulence2 and are well known in the wind energy community to often exceed fatigue and extreme loads anticipated during the design process.3,4 Off-axis load cases have not traditionally been included in the design process and such loads can pose problems, in particular, for gears and bearings found in gearboxes.5 

Although ocean current turbines may be placed in open water far from coasts, tidal channels have long been seen as the most viable locations for such turbines due to the strong and highly predictable nature of tidal currents. The flow in tidal channels is, however, complex and chaotic at a wide range of scales due to turbulence production by surface (i.e., wind) and bottom (i.e., no-slip) shears, waves, variable bathymetry (or bottom surface roughness), thermal instabilities, and even turbines themselves (i.e., if a turbine is in the wake of another turbine). Turbulence from all of these sources can create substantial small-scale temporal and spatial variability of the tidal flow field, resulting in potentially large bending moments.

Understanding and predicting bending moments in such turbulent environments require simultaneous two-point spatial measurements across and downstream of the rotor disc. Observational data for such two-point statistics is, however, not yet available in tidal channels over the range of spatial scales (i.e., from meters to hundreds of meters) relevant to tidal turbine arrays. By contrast, high fidelity computer simulations allow comparatively straightforward measurements of two-point statistics over a broad range of scales. Using simulations, the turbulent environment in which tidal current devices operate can be modeled, thereby enabling predictions of long-term turbine performance and turbulent fatigue loading.

Much of the current knowledge of tidal turbine loading is based on experience gained from studies of wind turbines, where attempts have been made to understand turbulent atmospheric environments and to assess the effects of bending moments on turbines. The interactions between the atmospheric boundary layer and utility scale wind turbines have been studied and classified using both experiments and numerical simulations. Using observational data, Mücke et al.,3 Morales et al.,6 and Milan et al.7 have connected extreme and fatigue turbine loads to the intermittent nature of turbulent wind speed fluctuations, and Kelley et al.2 have further examined the impacts of coherent turbulent structures on the dynamic response and bending moments of wind turbines. Wind tunnel boundary layer experiments such as those performed by Chamorro and Porté-Agel8 have revealed additional properties of turbulent wakes formed behind wind turbines, including the creation of non-axisymmetric wakes by inhomogeneous boundary layer inlet flows.

Recently, numerical simulations have begun to provide substantial additional knowledge regarding the connection between turbine loads, atmospheric boundary layer turbulence, and turbine wakes. More specifically, a number of recent studies9–14 have used large eddy simulations (LES) to model arrays of wind turbines in order to determine the extent to which turbine wake interactions reduce power output and increase loading, and to examine how wake recovery occurs in realistic atmospheric boundary layers (refer to Mehta et al.15 for a review of LES methods applied to wind farm aerodynamics). Calaf et al.9 and Meyers and Meneveau,10,14 in particular, have connected wake recovery to momentum entrainment from the atmospheric boundary layer, with the kinetic energy transported vertically into the wake roughly equal to the power extracted by turbines. The connection between properties of the turbulent flow field and turbine loads are often examined using advanced turbine simulators, such as the Fatigue, Aerodynamics, Structures, and Turbulence (FAST) code.16 

Turbine loading and performance in the ocean have also been studied using numerical simulations, often with a focus on accurately modeling boundaries such as channel bathymetry and turbine blade shape. Although there have been many studies of tidal turbines using Reynolds averaged Navier-Stokes (RANS) approaches,17–23 LES has become increasingly common for the study of tidal turbine performance and loading. Churchfield et al.24 and Gebreslassie et al.25 have both used LES to examine arrays of multiple tidal turbines in order to determine the effects of array layout on wake interactions and power output. Other studies have used LES to obtain high-resolution predictions for the near-field flow around single turbines, including the studies of tidal turbine loading by Kang et al.26,27 and Afgan et al.,28 as well as the study of tidal turbine noise generation by Lloyd et al.29 In addition to these computational studies, a number of experimental studies have provided insights into loads, wake properties, power output, and turbulence characteristics for both single tidal turbines30–39 and, more recently, two turbines.40 

By contrast to most prior computational studies of tidal turbine loading, the present paper is specifically focused on understanding the turbulent bending moments that would be experienced by a tidal current turbine placed in a realistic ocean environment, as well as how the moments are connected with fundamental turbulence statistics. Such statistics are critical for predicting moments and, ultimately, turbine performance and fatigue in advance of installation. In order to resolve small-scale turbulent motions and increase the level of physical realism, LES has been performed of an ocean tidal boundary layer for a range of different physical scenarios. Two of the scenarios loosely correspond to prospective tidal energy sites at Admiralty Head and Nodule Point in Puget Sound. LES results are compared with observational measurements of turbulence intensity (TI) profiles at these two locations from Thomson et al.41 The observational measurements have been used previously to examine turbulence intermittency, coherence, and anisotropy42 and are used here to constrain the physical parameters used in the simulations, as well as to gain confidence in the accuracy of the LES.

In the present study, a number of simulations for different physical scenarios have been performed in order to understand the effects of winds, waves, tidal velocity, stability, and tidal channel depth on the moments that would be experienced by a tidal turbine. Changes to the wind and tidal velocities result in changes to the turbulence shear production, waves generate Langmuir turbulence43–45 which increases vertical mixing near the surface, and boundary layer stability strongly affects the creation and properties of convective turbulence. Previously, Li et al.46 have shown using LES that Langmuir turbulence may play a relatively small role in the near-surface turbulence dynamics, and the present study will examine the role of such wave-driven turbulence on tidal turbine loading more specifically. The LES model used in the present study has the capability to represent all of these effects and has been used previously for several high-fidelity process studies of atmospheric and oceanic flows.44,45,47,48 Probability density functions, two-point correlations, velocity structure functions, Reynolds stresses, and other statistics are used to understand the dependence of turbulent bending moments on each of these physical effects, and an analytical connection is made between bending moments and two-point turbulence statistics.

The present study is differentiated from prior research by the extent of ocean physics included in the simulations and by the focus on bending moments that would be experienced by a tidal current turbine. Prior studies of tidal current turbines have often neglected the effects of wind, waves, or tidal motions in order to decrease computational cost or to focus on impacts from other physical processes. Conversely, prior fundamental studies of tidal flow and coastal boundary layers have often not connected physical processes with the loads and turbulent flow fields experienced by tidal turbines. For example, in one of the most comprehensive prior studies of tidal flows, Li et al.49 studied temporally varying tidal currents in an estuarine boundary layer, but did not consider waves, different stability conditions, or varying tidal channel depth, and were not specifically focused on tidal turbine loading. Similarly, the LES studies by Taylor et al.50 and Gayen et al.51 have provided considerable insights into boundary layer structure and turbulence properties for stratified boundary layers in open channel flows, but once again were not focused on bending moments generated in well-mixed tidal boundary layers, such as those examined here.

Finally, it should be noted that the present study does not include turbines in the simulations. Rather, turbines are represented as idealized infinitely thin circular discs during the calculation of bending moments and an attempt is made to understand how the bending moments that would be experienced by a tidal turbine are affected by different physical characteristics of the flow field, independent of the particular choice of turbine design.

Quantification of the expected bending moments on turbines in realistic ocean environments is an important objective of the present study. Prior research has been performed on wind turbine loading for different atmospheric conditions,2,16,52,53 including the specific loading associated with the blade structure. Similar computational and experimental studies have been performed on loads experienced by ocean current turbines,23,24,26–28,32,54–57 but the dependence of such loads on characteristics of the oceanic boundary layer such as wind and wave shear, boundary layer depth, and stability conditions have yet to be examined.

In the present paper, the bending moments (which are sometimes referred to as off-axis loads) that would be experienced by a tidal turbine are considered without restricting the analysis to a specific choice of turbine design. The moments at location x and time t, denoted as M(x, t), are measured for a rigid, infinitely thin circular disc of radius D and with a hub-height above the ocean bottom d, as shown in the schematic in Fig. 1. The direction normal to the disc is the x direction and the resulting moments are calculated from an integral over the area of the disc as

M(x,t)=Aρ[r×u(x+r,t)]ux(x+r,t)drydrz,
(1)

where r = [rx, ry, rz] is the distance from the central point to a location on the disc, u = [ux, uy, uz] is the velocity vector, ux is the velocity normal to the disc, and A is the cross-sectional area of the disc. Since the disc is assumed to be infinitely thin, rx = 0 and Eq. (1) can be expanded as

Mx(x,t)=Aρ[ryuz(x+r,t)rzuy(x+r,t)]ux(x+r,t)dA,
(2)
My(x,t)=Aρrzux2(x+r,t)dA,
(3)
Mz(x,t)=Aρryux2(x+r,t)dA.
(4)

The component Mx is about the axis normal to the turbine disc and provides a torque that either accelerates or decelerates the turbine. The components My and Mz are bending moments (or off-axis loads) about the horizontal and vertical bisects of the turbine disc, respectively (see Fig. 1). The expressions in Eqs. (3) and (4) are used to calculate My and Mz in Section IV. Bending moments are of particular concern in the design of turbine gearboxes, since such off-axis loads can place unexpectedly large stresses on gears and bearings that are not typically addressed during the design process.5 

FIG. 1.

Schematic of an idealized circular rotor disc of diameter D and hub-height d, showing the direction of bending moments Mx, My, and Mz defined in Eqs. (2)–(4). The hub-height d is measured from the bottom boundary of the tidal channel.

FIG. 1.

Schematic of an idealized circular rotor disc of diameter D and hub-height d, showing the direction of bending moments Mx, My, and Mz defined in Eqs. (2)–(4). The hub-height d is measured from the bottom boundary of the tidal channel.

Close modal

Bending moments are related to the characteristics of the turbulent flow field, and in particular to transverse velocity differences over different spatial separations across the rotor disc. This can be seen by assuming the rotor disc to be a circle with diameter D that is symmetric about the y and z axes and by decomposing the velocity into mean and fluctuating parts as ux=u¯x+ux, where (·)¯ indicates an average over time and horizontal directions. The resulting bending moments from Eqs. (3) and (4) are then written as (the x coordinate is suppressed for brevity)

My(x,t)=ryrz0ρα1rz[Δrzux(y+ry,z,t)]+ρα2rz[Δrzu¯x(z)]drzdry,
(5)
Mz(x,t)=rzry0ρα3ry[Δryux(y,z+rz,t)]drydrz,
(6)

where the αi variables have been introduced to simplify the notation and are given by

α1ux(y+ry,z+rz,t)+ux(y+ry,zrz,t)+2u¯x(z+rz),α2u¯x(z+rz)+u¯x(zrz)2ux(y+ry,zrz,t),α3ux(yry,z+rz,t)+ux(y+ry,z+rz,t)+u¯x(z+rz).
(7)

The velocity difference terms introduced in Eqs. (5) and (6) are defined as

Δrzux(y+ry,z,t)=ux(y+ry,z+rz,t)ux(y+ry,zrz,t),Δryux(y,z+rz,t)=ux(y+ry,z+rz,t)ux(yry,z+rz,t),Δrzu¯x(z)=u¯x(z+rz)u¯x(zrz),
(8)

where the last expression involving u¯x is connected to the mean vertical shear du¯x/dz. Note that since the tidal channel is assumed here to be horizontally homogeneous, mean flow statistics depend only on z and other shear components such as du¯x/dy are zero. Moments of the fluctuating velocity differences in Eq. (8) are, in turn, connected to turbulent structure functions. Properties of structure functions in tidal channels and their connection to bending moments will be discussed in more detail in Section IV.

From Eq. (5), My is seen to depend on two primary effects: a contribution due to variations in turbulent fluctuations over the disc and a contribution due to variations in the mean velocity over the disc. The latter effect is otherwise known as the mean shear, and thus My is termed the mean shear bending moment. By contrast, Eq. (6) indicates that Mz depends only on the differences in turbulent fluctuations over the disc; this moment is thus termed the turbulent eddy bending moment. In Section IV, Eqs. (5) and (6) will be used to understand bending moments in tidal channels for a range of physical scenarios.

The National Center for Atmospheric Research (NCAR) LES model is used to perform the simulations.44,47 In order to represent the effects of winds, waves, and tides, the simulations solve the forced Wave-Averaged Boussinesq (WAB) equations given by58–62 

ut+uL·u=p+buL·(us)+fc+sgs,
(9)
θt+uL·θ=sgs,
(10)
·u=0,
(11)

where u = [ux, uy, uz] is the three-dimensional Eulerian flow velocity, ω×u is the Eulerian vorticity, uL = u + us is the Lagrangian velocity, us is the wave-induced Stokes drift velocity, p is the pressure normalized by the background density ρ0, b is the buoyancy, fc is a driving term used to create the tidal current, θ is the potential temperature, and sgs denotes subgrid-scale (SGS) terms introduced by the LES modeling. The WAB equations account for the effects of wave forcing through advection by the Lagrangian velocity uL and through the forcing term uL⋅(∇us) on the right-hand side of Eq. (9), both of which include the wave-induced Stokes-drift velocity us. Note that in the present study Coriolis effects are not considered; a similar approach has been used in prior studies of tidal boundary layers49 and rotational effects are not expected to make a significant contribution to the dynamics at the range of scales examined here. An equation of state is used for the density ρ, namely,

ρ=ρ0[1+βT(θ0θ)],
(12)

where βT = 2 × 10−4 K−1 is the thermal expansion coefficient and salinity effects have been neglected. The buoyancy is b=bẑ where b = −/ρ0. Parameter values used to obtain the buoyancy are ρ0 = 1000 kg/m3, θ0 = 290.16 K, and g = 9.81 m/s2.

The WAB equations in Eqs. (9)–(12) have been used previously in a number of computational studies of wave-driven Langmuir turbulence.44,45,48 Langmuir turbulence consists of disordered collections of counter-rotating vortical cells (typically called “Langmuir” cells). These cells are typically 10 m deep in the vertical and can be up to 1 km long in the horizontal direction (although typically they are much shorter, as in the present simulations). They create surface convergence zones where foam, plankton, and other debris collect, resulting in characteristic “windrows.”43 The primary effect of Langmuir turbulence on the flow field is to create more intense vertical mixing, which may substantially alter the loads experienced by a turbine in a tidal channel.

The NCAR LES model solves Eqs. (9)–(12) on a structured, rectilinear grid. A spectral method is used in horizontal (x, y) directions and a second-order finite difference method is used in the vertical (z) direction for the velocities. A third-order finite difference method is used in the vertical direction for θ. A three-step explicit Runge-Kutta method is used to advance the solution in time, and the Poisson pressure equation associated with Eq. (11) is solved iteratively at each time step. A variable time step is used with a constant Courant-Friedrichs-Lewy (CFL) number of 0.5. The SGS model used in the simulations is described by Sullivan et al.,63 and takes into account not only the subgrid-scale influences proportional to the resolved-scale strain rate but also the influences proportional to the horizontal average strain rate. At the bottom boundary, the SGS model is calibrated to match Monin-Obukhov similarity theory. This SGS model has been used in previous studies of Langmuir turbulence.44,45,48

The forcing term, fc=fcx̂, in Eq. (9) is used to create a tidal current in the x direction and is constant in both space and time. Following Li et al.,49 future work may include temporal variations in fc in order to consider situations in which the tidal stream changes magnitude and direction. Here, however, a time-independent fc is used in the simulations to create uniform tidal velocities of 1.0 m/s, 2.0 m/s, and 3.0 m/s.

The effects of wave forcing are parameterized in Eqs. (9) and (10) via the Stokes drift velocity us appearing in the Lagrangian velocity uL. The Stokes drift velocity is modeled as

us=us(z)[cos(ϑs)x̂+sin(ϑs)ŷ],
(13)

where us(z) is the Stokes drift vertical profile and ϑs is the wave direction. Using the Donelan empirical spectrum64,65 for the ocean wave field, the resulting vertical profile us(z) decreases rapidly with depth, as shown in Fig. 2. In the simulations, the Stokes drift velocity is assumed to be the same at all horizontal (x, y) locations with an angle of ϑs=0°. That is, the wave field is assumed to be perfectly aligned with the tidal direction. The magnitude of the Stokes drift and the corresponding profile us(z) are determined by setting the 10 m wind speed, denoted U10, in the Donelan spectrum. Profiles of us(z) used in the present simulations are shown in Fig. 2. These profiles correspond to U10 = 5.0 m/s and 10.0 m/s, which give spectral significant wave heights of 0.7 m and 2.7 m, respectively. These wave heights have been chosen based on the typical (i.e., less than 1 m) and maximum (i.e., 2.3 m) wave heights observed at Admiralty Inlet in Puget Sound.66 

FIG. 2.

Profiles of the Stokes drift magnitude us(z) for the B (solid black lines), B (red dashed lines), B+ (blue dashed-dotted lines), and Wa+ (heavy black lines) simulations. The profiles for the B, B, and B+ simulations correspond to U10 = 5.0 m/s, and the profile for the Wa+ simulation corresponds to U10 = 10.0 m/s, as summarized in Table I. The inset shows the same profiles with a logarithmic axis for us(z).

FIG. 2.

Profiles of the Stokes drift magnitude us(z) for the B (solid black lines), B (red dashed lines), B+ (blue dashed-dotted lines), and Wa+ (heavy black lines) simulations. The profiles for the B, B, and B+ simulations correspond to U10 = 5.0 m/s, and the profile for the Wa+ simulation corresponds to U10 = 10.0 m/s, as summarized in Table I. The inset shows the same profiles with a logarithmic axis for us(z).

Close modal

All of the simulations performed in the present study are fully periodic in the horizontal (x, y) directions for all simulation variables. With the tidal forcing fc=fcx̂, this results in an inflow along the x direction that is directly tied to the outflow. On the top boundary, the horizontal velocities (i.e., ux and uy) are driven by a wind shear stress of magnitude τ in the ϑw direction. In the simulations, ϑw is always chosen to be in the positive x direction, resulting in a wind shear friction velocity uτ in the x direction. The vertical velocity uz is constrained to be zero at the top boundary in all simulations. Depending on the intended thermal stability of the simulations, either a cooling, warming, or adiabatic temperature flux, Q0, is applied at the top boundary. At the bottom boundary, a no-slip condition is used for velocities in all three directions and the temperature flux is zero. The surface roughness for the bottom no-slip boundary is z0 = 0.001 m, and the SGS model used in the LES matches the log-law for the mean velocity at the lower boundary (as described in Sullivan et al.63). Note that the bottom boundary is assumed in this study to be flat and absent of any bathymetric features.

The physical setup and conditions used in the numerical simulations have been loosely guided by the wind, wave, and tidal conditions found at the Admiralty Head and Nodule Point locations in Puget Sound. These sites have been characterized in Thomson et al.41 and recently collected data allows some limited comparisons to be made between the simulation and experimental results. Since the present simulations are, however, more broadly focused on the effects of winds, waves, boundary layer depth, instabilities, and tidal velocity, a series of 13 simulations corresponding to 11 different physical scenarios have been performed, as summarized in Table I. The baseline simulation (denoted “B” in Table I) has horizontal lengths Lx = 600 m and Ly = 300 m with depth H = 39 m, a tidal velocity of Ut = 2 m/s, upper wind forcing resulting in uτ = 6.3 mm/s, waves created by a 10 m wind of U10 = 5.0 m/s, and no heat flux at the surface (corresponding to a neutral boundary layer). Figure 2 shows the Stokes drift profile used in the baseline case. The horizontal dimensions of the computational grid in the baseline case are Nx = 512 and Ny = 256, resulting in a horizontal grid resolution of Δx,y = 1.2 m. The vertical grid size is Nz = 128 giving a vertical resolution of Δz = 0.30 m. The turbulent Langmuir number in the baseline simulation is Lat=uτ/us(0)=0.25, which is a standard value observed in realistic ocean environments,44 although it may be larger than the values found in typical tidal channels.

TABLE I.

Summary of the physical setup used in the simulations. Horizontal and vertical domain sizes are denoted as Lx, Ly, and h, computational domain sizes are given by Nx, Ny, and Nz, Ut is the tidal velocity, uτ is the top surface friction velocity, τ is the wind stress, ϑw is the wind direction with respect to the x-axis, U10 is the wind speed 10 m above the ocean surface, us(0) is the Stokes drift velocity at the surface, Lat=uτ/us(0) is the turbulent Langmuir number, ϑs is the Stokes drift direction with respect to the x-axis, Q0 is the surface heat flux, and z0 is the roughness height of the bottom boundary. Parameters varied in each of the simulations are shown in bold.

CaseBaselineWindWaveDepthStabilityTide
SimulationBBB+WiWi+WaWa+HH+InIn+TT+
Lx (m) 600 600 600 600 600 600 600 600 600 600 600 600 600 
Ly (m) 300 300 300 300 300 300 300 300 300 300 300 300 300 
h (m) 39 39 39 39 39 39 39 22 56 39 39 39 39 
Nx 512 256 1024 512 512 512 512 512 512 512 512 512 512 
Ny 256 128 512 256 256 256 256 256 256 256 256 256 256 
Nz 128 64 256 128 128 128 128 72 184 128 128 128 128 
Ut (m/s) 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 3.0 
uτ (mm/s) 6.3 6.3 6.3 0.0 12.6 6.3 6.3 6.3 6.3 6.3 6.3 6.3 6.3 
τ (N/m20.04 0.04 0.04 0.0 0.08 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 
ϑw (deg) 
U10 (m/s) 5.0 5.0 5.0 5.0 5.0 0.0 10.0 5.0 5.0 5.0 5.0 5.0 5.0 
us(0) (m/s) 0.10 0.10 0.10 0.10 0.10 0 0.20 0.10 0.10 0.10 0.10 0.10 0.10 
Lat 0.25 0.25 0.25 0.0 0.35  0.18 0.25 0.25 0.25 0.25 0.25 0.25 
ϑs (deg) 
Q0 (W/m25.0 5.0 
z0 (m) 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 
CaseBaselineWindWaveDepthStabilityTide
SimulationBBB+WiWi+WaWa+HH+InIn+TT+
Lx (m) 600 600 600 600 600 600 600 600 600 600 600 600 600 
Ly (m) 300 300 300 300 300 300 300 300 300 300 300 300 300 
h (m) 39 39 39 39 39 39 39 22 56 39 39 39 39 
Nx 512 256 1024 512 512 512 512 512 512 512 512 512 512 
Ny 256 128 512 256 256 256 256 256 256 256 256 256 256 
Nz 128 64 256 128 128 128 128 72 184 128 128 128 128 
Ut (m/s) 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 3.0 
uτ (mm/s) 6.3 6.3 6.3 0.0 12.6 6.3 6.3 6.3 6.3 6.3 6.3 6.3 6.3 
τ (N/m20.04 0.04 0.04 0.0 0.08 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 
ϑw (deg) 
U10 (m/s) 5.0 5.0 5.0 5.0 5.0 0.0 10.0 5.0 5.0 5.0 5.0 5.0 5.0 
us(0) (m/s) 0.10 0.10 0.10 0.10 0.10 0 0.20 0.10 0.10 0.10 0.10 0.10 0.10 
Lat 0.25 0.25 0.25 0.0 0.35  0.18 0.25 0.25 0.25 0.25 0.25 0.25 
ϑs (deg) 
Q0 (W/m25.0 5.0 
z0 (m) 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3 

A pair of simulations with identical physical parameters to the baseline simulation but with doubled (B+) and halved (B) resolution have been performed in order to ascertain grid convergence. All other simulations have identical values of Lx, Ly, Nx, and Ny to the baseline case (i.e., simulation B), but h is varied to 22 m and 56 m in two of the simulations (denoted H and H+ in Table I). These two values of h correspond to the locations at Nodule Point and Admiralty Head, respectively, allowing qualitative comparisons of the simulations with observations from Thomson et al.41 Note, however, that the comparisons with observations are limited since the present LES only matches the tidal channel depths at the Nodule Point and Admiralty Head locations. More rigorous validation of the simulations requires additional information from the observational locations about realistic wind and wave forcing, bathymetry, and inlet flow properties. As such, the focus of the present study is on understanding the effects of different physical scenarios on current turbine bending moments, as opposed to rigorously modeling specific real-world tidal boundary layer locations.

As summarized in Table I, the tidal velocity (T and T+), wind strength (Wi and Wi+), wave strength (Wa and Wa+), and surface heat flux (In and In+) have each been varied in the simulations to be successively smaller and larger than the baseline values. Section IV will show that relatively small changes in the tidal velocity Ut result in large changes in flow properties, and so the T and T+ simulations correspond to values of Ut that are only 1 m/s smaller and larger, respectively, than the baseline value of Ut = 2.0 m/s. In order to determine wind and wave effects, the surface wind was turned off for the Wi simulation and the Stokes drift forcing was turned off in the Wa simulation. The Wi+ and Wa+ simulations correspond to surface wind and wave forcings, respectively, which are twice as large as the baseline values. Stronger forcings produce even greater changes in flow properties, but become increasingly less realistic in comparison to the standard values used for the baseline case. Profiles of the Stokes drift us(z) used for the B+, B, and Wa+ simulations are shown in Fig. 2. Finally, ocean surface cooling of Q0 = −5.0 W/m2 has been found previously45 to create an unstable mixed layer in the presence of Stokes drift forcing, and the same cooling is used here for the In simulation. In order to perform a more meaningful examination of stability effects, heating with the same magnitude but opposite sign is used for the In+ simulation.

After an initial spin-up period during which the boundary layer turbulence is allowed to develop, each of the simulations summarized in Table I were run for a total of 6 h of virtual time, with an average timestep of 0.15 s for simulations with the baseline resolution (i.e., the resolution of case B). The higher resolution simulation B+ required a smaller average timestep and thus more overall timesteps were necessary to reach 6 h of virtual time. The total simulated time in each case was sufficient to allow turbulence to fully develop, and the ux velocity and Reynolds stress profiles were found to be stable and unchanging with additional simulation time.

The basic flow state in each of the simulations at the end of 6 h is represented by the snapshot of the ux velocity field for case B shown in Fig. 3. As expected for a tidal boundary layer, the velocity is generally smallest at the bottom boundary and greatest near the surface, and there is substantial spatial variability resulting from the turbulent nature of the flow field. The domain shown in Figure 3 is also sufficiently large in the horizontal (x, y) directions to contain many statistically independent turbulent structures; this will be confirmed quantitatively in Section IV C where it will be shown that the turbulent integral scales are roughly 50 m in the x direction and 10 m in the y direction, both of which are small relative to the domain size [Lx, Ly] = [600m, 300 m].

FIG. 3.

Three-dimensional snapshot of the horizontal velocity field ux after 6 h of virtual time for the B simulation.

FIG. 3.

Three-dimensional snapshot of the horizontal velocity field ux after 6 h of virtual time for the B simulation.

Close modal

In the following, a validation of the simulation results is outlined, including a study of grid convergence. Single point statistics of the velocity and temperature for the different simulations are then presented, followed by an analysis of two point turbulence statistics. Finally, off-axis bending moments for each of the different physical scenarios summarized in Table I are outlined.

Validation of the physical parameters and numerical methods used in the simulations is critical for grounding the simulations in reality and for identifying any potential numerical issues. Validation of the simulations is performed here through comparison with theory, prior studies, and observational data from Thomson et al.41 Prior work by Sullivan and Patton in the atmospheric convective boundary layer67 indicates that mesh sizes of 2563 in a domain of size (Lx, Ly, Lz) = (5120, 5120, 2048)m are sufficient for achieving scale separation between the large scale energy containing eddies and the filter cutoff scale, and the present simulations have a similar dynamic range in a smaller domain with finer grid resolution. This computational grid size has been used as a baseline for the present simulations (roughly corresponding to the size of case B), and in the following validation results are presented for three different grid resolutions in order to determine grid convergence.

Figure 4 shows vertical profiles of u¯x along with diagonal and off-diagonal Reynolds stresses for the B, B−, and B+ simulations. Here, (·)¯ denotes an average taken over horizontal (x, y) directions as well as time t. The profile of u¯x shows good agreement with the theoretically predicted log-layer profile and the qualitative behavior of the diagonal and off-diagonal Reynolds stresses is consistent with prior studies of wall bounded shear flows. That is, ux2¯>uy2¯>uz2¯ and uxuz¯<0, with the greatest magnitudes of all stresses occurring close to the bottom boundary. The Reynolds stress profiles also show good agreement with previous studies of the atmosphere47,67 and tidal channels.68 

FIG. 4.

Vertical profiles of (a) u¯x, (b) diagonal Reynolds stresses, ux2¯ (solid), uy2¯ (dashed), and uz2¯ (dashed–dotted), and (c) off-diagonal Reynolds stresses, uxuz¯ (solid), uxuy¯ (dashed), and uyuz¯ (dashed–dotted). Each panel shows values for three different grid resolutions, corresponding to simulations B (black lines), B (blue lines), and B+ (red lines), as a function of z/h. Additionally, a theoretical log layer velocity profile, u=(u*/κ)ln(z/z0), is shown in grey in (a) (this profile is obscured by the overlapping lines for the B, B, and B+ simulations).

FIG. 4.

Vertical profiles of (a) u¯x, (b) diagonal Reynolds stresses, ux2¯ (solid), uy2¯ (dashed), and uz2¯ (dashed–dotted), and (c) off-diagonal Reynolds stresses, uxuz¯ (solid), uxuy¯ (dashed), and uyuz¯ (dashed–dotted). Each panel shows values for three different grid resolutions, corresponding to simulations B (black lines), B (blue lines), and B+ (red lines), as a function of z/h. Additionally, a theoretical log layer velocity profile, u=(u*/κ)ln(z/z0), is shown in grey in (a) (this profile is obscured by the overlapping lines for the B, B, and B+ simulations).

Close modal

From a more quantitative perspective, deviations of the simulation results from the theoretical log-layer profile are characterized by the quantity ϕm, which is defined as63,68

ϕm(z)κzu*Uz,
(14)

where U(u¯x+u¯y)1/2,κ=0.41, and u* is the friction velocity at the bottom no-slip boundary. While the ϕm profile for case B, shown in Fig. 5, is not precisely unity, as would be expected in the case of a perfect log-layer profile, the small deviations are consistent with those found in previous studies63,68 of atmospheric and oceanic boundary layers. Figure 5 shows that the deviations are largest near the channel bottom, approaching a high value of 1.2 and a low value of 0.7. In the simulations, grid anisotropy (i.e., the ratio of Δx,y to Δz) was found to have a large effect on the ϕm profile, with higher values of the ratio (greater than 10) resulting in more substantial deviations from unity. A grid anisotropy ratio of approximately four, as used in the present simulations, results in ϕm profiles that are in line with previous studies.63,68

FIG. 5.

Vertical profiles of ϕm, defined in Eq. (14), for three different grid resolutions, corresponding to the simulations B (solid black lines), B (blue dashed lines), and B+ (red dashed-dotted lines), as a function of z/h.

FIG. 5.

Vertical profiles of ϕm, defined in Eq. (14), for three different grid resolutions, corresponding to the simulations B (solid black lines), B (blue dashed lines), and B+ (red dashed-dotted lines), as a function of z/h.

Close modal

Figures 4 and 5 also show the effects of varying grid resolution on the velocity and Reynolds stress profiles for low (B), base (B), and high (B+) resolutions. The vertical profiles of u¯x in Fig. 4(a) indicate that there is good agreement between the base and high resolution runs, with slight deviations near the top of the boundary layer for the low resolution run. The diagonal Reynolds stresses in Fig. 4(b) show moderate variations for the three resolutions along the entire depth of the domain. The uxuz¯ Reynolds stress in Fig. 4(c) shows good agreement for z/h between 0.2 and 1, with more substantial deviations in the lower resolution run. The ϕm profiles in Fig. 5, show similar deviations from unity for all three resolutions, but the deviations occur at heights proportional to the grid size. In particular, the ϕm profile for the highest resolution run does not show substantially less deviation from unity, just a deviation at a height closer to z/h = 0. Overall, the simulations with different grid resolutions indicate that the resolution of the base run (B) captures most statistics with acceptable accuracy relative to a simulation with grid resolution doubled in each direction. In the interest of exploring as large a parameter space as possible, the resolution of the base simulations thus appears to be appropriate for developing insights into the variation of off-axis loads for varying winds, waves, stability, and other physical parameters.

Finally, comparisons with the observational data of Thomson et al.41 are made using the TI, which is defined as

TIux2¯1/2u¯x.
(15)

The resulting vertical profiles of TI are shown in Fig. 6 for the simulations B (corresponding to a 39 m depth), H (corresponding to a 22 m depth at the Nodule Point location), and H+ (corresponding to a 56 m depth at the Admiralty Head location).41 Figure 6 shows that the simulations correctly capture the basic shape of the observational TI profiles, but slightly underestimate the magnitude of TI (note that observational data were only available up to roughly 18 m above the ocean bottom). Increasing the grid resolution (i.e., simulation B+) resulted in a slight (roughly 5%) increase in the magnitude of the TI profile, but the LES results for all values of h do display the same qualitative behavior seen in the observations; namely, a large increase in TI near the bottom of the ocean boundary layer as ux2¯ increases and u¯x decreases (note that immediately at the bottom boundary TI = 0 since ux2¯=0). While the roughness length z0 used in the simulations (see Table I) was chosen to give a close match to the observed TI profiles from Thomson et al.,41 more complete agreement with the observations most likely requires more realistic modeling of bathymetry and inlet flow conditions, both of which are beyond the scope and intent of the present study. Such improved agreement would also be necessary to obtain more accurate predictions of wake expansion rate in simulations that include models for tidal turbines. Nevertheless, the good qualitative agreement and approximate quantitative agreement between the LES results and observations shown in Fig. 6 provides confidence in the validity of the simulations.

FIG. 6.

Turbulence intensity (TI) as a function of z/h for h = 22 m, 39 m, and 56 m, corresponding to simulations H, B, and H+, respectively. Observational results from Thomson et al.41 are also shown for the 22 m Nodule Point location (red dashed line) and the 56 m Admiralty Head location (red dashed-dotted line).

FIG. 6.

Turbulence intensity (TI) as a function of z/h for h = 22 m, 39 m, and 56 m, corresponding to simulations H, B, and H+, respectively. Observational results from Thomson et al.41 are also shown for the 22 m Nodule Point location (red dashed line) and the 56 m Admiralty Head location (red dashed-dotted line).

Close modal

Single point statistics are commonly used for understanding and characterizing tidal current environments and provide valuable information to developers, designers, and modelers regarding mean flow magnitudes, turbulence intensities, and transient effects. However, the ability of single point statistics to predict bending moments and off-axis loads is limited. In this section, the impacts of various physical parameters on single point statistics are analyzed, and in Sections IV D and V the connections between both single and two point statistics and bending moments are examined in more detail.

Figures 7(a) and 7(c) show that profiles of the mean velocity u¯x do not change substantially as either the wind or the stability conditions change. Changes in the wave strength do, however, affect u¯x near the surface, where u¯x for Wa+ decreases below the corresponding values for B and Wa, as shown in Fig. 7(b). The increased wave strength for the Wa+ simulation, modeled as a stronger Stokes drift velocity (see Table I) results in an anti-Stokes effect, as previously observed by Hamlington et al.45 and explained by McWilliams and Fox-Kemper.69 Profiles of u¯x also change for T and T+ as compared to the base case due to the large changes in mean tidal velocity. Overall, the general uniformity of the mean velocity profile with respect to changes in the wind, wave, and stability conditions is significant for the understanding of turbine loads, since a drastic change in vertical velocity profile for different physical parameters would result in added complexity towards addressing loads associated with the mean vertical shear.

FIG. 7.

Vertical profiles of u¯x for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. In each panel, solid black lines are baseline results for case B, dashed blue lines show Wi, Wa, In, and T, and red dashed-dotted lines show Wi+, Wa+, In+, and T+.

FIG. 7.

Vertical profiles of u¯x for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. In each panel, solid black lines are baseline results for case B, dashed blue lines show Wi, Wa, In, and T, and red dashed-dotted lines show Wi+, Wa+, In+, and T+.

Close modal

Vertical profiles of Reynolds stress are shown in Fig. 8 and provide insights into the level of turbulent activity within the flow. Both the Wa+ and In+ simulations show increases in ux2¯ compared to the base case, although little change is seen in either uz2¯ or uxuz¯ (note also that the increase in ux2¯ is relatively small for In+). An increase in ux2¯ represents an increase in the largest component of turbulence kinetic energy, and will impact not only the bending moments discussed in this study but also transient torque loads. Large deviations from the base Reynolds stress profiles are present in the T and T+ profiles, and the magnitudes of the changes are worth noting: a 50% increase of mean velocity approximately doubles the ux2¯ Reynolds stress, while halving the mean velocity results in roughly a factor of four decrease in ux2¯.

FIG. 8.

Vertical profiles of Reynolds stresses ux2¯ (solid), uxuz¯ (dashed), and uz2¯ (dashed-dotted) for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. In each panel, black lines are baseline results for case B, blue lines show Wi, Wa, In, and T, and red lines show Wi+, Wa+, In+, and T+. Note that the horizontal axis of panel (d) is roughly twice as large as the horizontal axes for panels (a)–(c), all of which are equal.

FIG. 8.

Vertical profiles of Reynolds stresses ux2¯ (solid), uxuz¯ (dashed), and uz2¯ (dashed-dotted) for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. In each panel, black lines are baseline results for case B, blue lines show Wi, Wa, In, and T, and red lines show Wi+, Wa+, In+, and T+. Note that the horizontal axis of panel (d) is roughly twice as large as the horizontal axes for panels (a)–(c), all of which are equal.

Close modal

The mean temperature profiles remain close to the background value θ0 for all but the unstable (In) and stable (In+) cases. Figure 9 shows that there is a sharp increase in the mean potential temperature θθ0¯ for In+ due to heating at the surface while there is a sharp decrease for In due to surface cooling. These large changes in the temperature profiles are accompanied by corresponding changes in the temperature fluxes, as shown in Fig. 10. For the unstable case (In), uxθ¯<0, corresponding to an anti-correlation between ux and θ; while for the stable case (In+), uxθ¯>0. In Fig. 10(c), uzθ¯>0 for the unstable case, corresponding to a downward flux of low temperature fluid as the surface is cooled; while uzθ¯<0 for the stable case, corresponding to a downward flux of high temperature fluid as the surface is heated. The horizontal fluxes in Fig. 10(b) for the stable and unstable cases remain close to zero. It should be noted that there are large variations in uxθ¯ and uyθ¯ for the base case (B). These variations are due primarily to large θ2¯ in the base simulation resulting from no imposed constraint on the correlation between the horizontal turbulence velocities and θ, other than the condition that the vertical flux of θ be zero at the surface.

FIG. 9.

Vertical profiles of θθ0¯ for unstable In (blue dashed line), baseline B (black solid line), and stable In+ (red dashed-dotted line) simulations.

FIG. 9.

Vertical profiles of θθ0¯ for unstable In (blue dashed line), baseline B (black solid line), and stable In+ (red dashed-dotted line) simulations.

Close modal
FIG. 10.

Vertical profiles of (a) uxθ¯, (b) uyθ¯, and (c) uzθ¯ for In (blue dashed lines), B (black solid lines), and In+ (red dashed-dotted lines) simulations.

FIG. 10.

Vertical profiles of (a) uxθ¯, (b) uyθ¯, and (c) uzθ¯ for In (blue dashed lines), B (black solid lines), and In+ (red dashed-dotted lines) simulations.

Close modal

Two-point velocity correlation functions and integral length scales have been shown in prior studies to be accurate predictors of turbine loading.2,3,6 In particular, these studies have shown that when the integral length scale of the turbulence is similar to the diameter of the turbine rotor, the loads on the turbine are greatest; turbulence with this characteristic length scale results in large velocity gradients across the rotor swept area and correspondingly large bending moments. In Section II, it was explicitly shown that two-point statistics such as turbulence structure functions across the rotor disc are closely connected to bending moments.

Figure 11 shows the longitudinal velocity correlation, f11, and transverse velocity correlation, f12, as functions of height, where the correlation functions are given by

f11(z,Δx)=ux(x,y,z)ux(x+Δx,y,z)¯ux2¯,f12(z,Δy)=ux(x,y,z)ux(x,y+Δy,z)¯ux2¯,
(16)

and the corresponding integral scales Λ11 and Λ12 are defined here as

Λ11(z)=0Lx/2f11(z,x)dx,Λ12(z)=0Ly/2f12(z,y)dy.
(17)

Note that the upper bounds of the integrals in Eq. (17) are taken to be half the full horizontal domain in the x and y directions in order to avoid spuriously high integral scales due to the periodicity of the simulations in horizontal directions.

FIG. 11.

Vertical profiles of the (a) longitudinal, f11(z, Δx), and (b) transverse, f12(z, Δy), correlation functions for the baseline simulation (B). Colors show values of f11 and f12 at each z and Δx or Δy, and solid black lines show the integral scales Λ11 and Λ12 defined in Eq. (17).

FIG. 11.

Vertical profiles of the (a) longitudinal, f11(z, Δx), and (b) transverse, f12(z, Δy), correlation functions for the baseline simulation (B). Colors show values of f11 and f12 at each z and Δx or Δy, and solid black lines show the integral scales Λ11 and Λ12 defined in Eq. (17).

Close modal

For the base case simulation B, Fig. 11 shows that over the full tidal channel depth the longitudinal integral scale Λ11 is between roughly 20 and 50 m and the transverse integral scale Λ12 is between 0 and 15 m. Figure 11(a) shows that, generally, the turbulence remains longitudinally correlated for larger Δx as z/h increases, although the correlation approaches zero at the upper surface; these variations are reflected in the increase of Λ11(z) with increasing z/h. For the transverse correlation f12 in Fig. 11(b), Λ12(z) increases monotonically with increasing z/h, reflecting increased transverse correlations for larger Δy from the bottom to the top of the tidal channel.

As shown in Fig. 12, the transverse correlation f12, which was shown in Section II to be closely connected to the bending moment Mz, is qualitatively similar for all of the simulations performed in the present study. As in Figs. 7 and 8 for the mean velocities and Reynolds stresses, respectively, some of the largest changes in f12 occur as the wave forcing varies, with approximately a factor of 2 increase in Λ12 at all z/h for the Wa+ simulation. This indicates that the greatest changes in bending moments may occur as the wave forcing is varied. The transverse length scale also varies for the other simulations, but magnitude changes are more modest, on the order of 10%.

FIG. 12.

Vertical profiles of the transverse correlation f12(z, Δy) for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. The top row shows Wi, Wa, In, and T, and the bottom row shows Wi+, Wa+, In+, and T+. Colors show values of f11 and f12 at each z and Δy, solid black lines are Λ12 for each simulation, and dashed black lines are Λ12 for the base case (B).

FIG. 12.

Vertical profiles of the transverse correlation f12(z, Δy) for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. The top row shows Wi, Wa, In, and T, and the bottom row shows Wi+, Wa+, In+, and T+. Colors show values of f11 and f12 at each z and Δy, solid black lines are Λ12 for each simulation, and dashed black lines are Λ12 for the base case (B).

Close modal

In addition to understanding the scale of the turbulent fluxes present in the ocean boundary, it is also important to understand the magnitude and distribution of velocity variations—these two properties feed directly into the distribution of bending moments experienced by the turbine, as discussed in Section II. To understand these spatial velocity differences, Fig. 13 shows structure functions S1j defined as

S1j(z,rj)=|ux(x+rj)ux(xrj)|¯ux2¯1/2,
(18)

where x is a position in space and rj is a vector of length r in the j (i.e., x, y, or z) direction. Each of the structure functions shown in Fig. 13 follows an approximate r1/3 scaling, in accordance with classical Kolmogorov theory,70 although the departures from this scaling are greatest for S11 and S12 near the bottom boundary. In both cases, the structure functions have reduced slope, indicating greater prominence of small scale motions relative to large scales.

FIG. 13.

Structure functions (a) S11, (b) S12, and (c) S13 for the baseline simulation (B), where the structure functions are defined in Eq. (18). Results are shown for eight different heights, for z/h between the bottom boundary at z/h = 0 (blue) and the top surface at z/h = 1 (red). Black dashed lines show r1/3 scaling relations, corresponding to the prediction from Kolmogorov.70 

FIG. 13.

Structure functions (a) S11, (b) S12, and (c) S13 for the baseline simulation (B), where the structure functions are defined in Eq. (18). Results are shown for eight different heights, for z/h between the bottom boundary at z/h = 0 (blue) and the top surface at z/h = 1 (red). Black dashed lines show r1/3 scaling relations, corresponding to the prediction from Kolmogorov.70 

Close modal

Finally, Fig. 14 shows probability density functions (pdfs) of the transverse velocity increments Δryux for different values of Δry, where Δryux is defined in Eq. (8). These increments play a central role in determining the magnitude of Mz, as shown by the integral in Eq. (6). In general, the pdfs become increasingly non-Gaussian as Δry decreases, and the pdfs for Δry at typical rotor radii (∼5 m) show indications of non-Gaussianity and intermittency. The implications of this are significant for models that assume a normal distribution of velocity differences, since Fig. 14 indicates that extreme events, leading to greater bending moments, will occur more frequently than expected from a normal distribution.

FIG. 14.

Probability density functions (pdfs) P(Δryux) for (a) Δry=1,5, and 25 m at z = 15 m and (b) Δry=5m at z = 5, 15, and 25 m. Each of the pdfs is shifted by two orders of magnitude for clarity and solid black lines show Gaussian distributions.

FIG. 14.

Probability density functions (pdfs) P(Δryux) for (a) Δry=1,5, and 25 m at z = 15 m and (b) Δry=5m at z = 5, 15, and 25 m. Each of the pdfs is shifted by two orders of magnitude for clarity and solid black lines show Gaussian distributions.

Close modal

The final component of the analysis regards bending moments, which are outlined analytically in Section II. Bending moments are primary drivers of turbine design, and understanding the magnitude and physical properties of such moments will aid in the design and performance prediction of ocean current turbines. Two primary types of bending moments are considered here: shear moments My and eddy moments Mz. The eddy moments are so named since they have zero means for the x–y homogeneous tidal channels examined in the present study. These moments are calculated using Eqs. (3) and (4) for a circular infinitely thin rotor disc of diameter D = 10 m. This diameter was chosen to be approximately representative of O(10m) real-world turbines, such as the 11 m diameter Seaflow turbine,71 the 16 m diameter OpenHydro turbine,72 and the 20 m diameter Seagen turbine.73 The moments are calculated at a range of different hub heights d (see Fig. 1) between 5 m and 34 m.

Figure 15(a) shows that My is strongly dominated by mean shear in the flow and increases substantially as z/h approaches zero where the mean shear is largest. The magnitude of the median bending moment at the channel bottom is three times larger than the median bending moment at the surface, with half of the change occurring in the bottom 10% of the boundary layer. The 90th percentile bending moment is approximately double the median bending moment, and the two vertical profiles are qualitatively very similar.

FIG. 15.

Probability density functions (pdfs) of bending moments (a) My(z) and (b) Mz(z) as functions of height z/h for the baseline simulation (B). The magnitude of the pdf is represented by the color, the median bending moments are shown as solid black lines, and 10th and 90th percentiles are shown as black dashed-dotted lines. The logarithm of pdf magnitudes is shown for clarity.

FIG. 15.

Probability density functions (pdfs) of bending moments (a) My(z) and (b) Mz(z) as functions of height z/h for the baseline simulation (B). The magnitude of the pdf is represented by the color, the median bending moments are shown as solid black lines, and 10th and 90th percentiles are shown as black dashed-dotted lines. The logarithm of pdf magnitudes is shown for clarity.

Close modal

Changes to the My bending moments for the different simulation cases, as shown in Fig. 16, are subtle but important. In Fig. 16(b), increasing the strength of the wave forcing results in an increase in My near the surface and a small decrease in the middle half of the boundary layer. This change is attributable, in part, to a similar change in the mean velocity profile, shown in Fig. 7, where the increased wave forcing case was the only case to substantially alter the mean velocity profile. Changes in the surface wind forcing and boundary layer stability have little impact on My, as shown in Figs. 16(a) and 16(c), respectively. The bending moments are very sensitive to changes in mean tidal velocity, as shown in Fig. 16(d), with significantly weaker moments in the T case and larger moments in the T+ case.

FIG. 16.

Probability density functions (pdfs) of shear bending moments My for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. The top row shows Wi, Wa, In, and T and the bottom row shows Wi+, Wa+, In+, and T+. The median bending moments for each simulation case are shown as solid black lines, and 10th and 90th percentiles are shown as black dashed-dotted lines. Baseline results for the median moments are shown as black dashed lines. The logarithm of pdf magnitudes is shown for clarity.

FIG. 16.

Probability density functions (pdfs) of shear bending moments My for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. The top row shows Wi, Wa, In, and T and the bottom row shows Wi+, Wa+, In+, and T+. The median bending moments for each simulation case are shown as solid black lines, and 10th and 90th percentiles are shown as black dashed-dotted lines. Baseline results for the median moments are shown as black dashed lines. The logarithm of pdf magnitudes is shown for clarity.

Close modal

The eddy bending moment Mz, shown for the baseline case in Fig. 15(b), is driven by horizontal imbalances in ux velocities across the rotor disc. The median moments at the channel bottom are nearly double the surface value, and the increase between the surface and the bottom is nearly linear. At mid-depths, the magnitude of Mz is similar to the magnitude of My. Much like My, the magnitude of the 90th percentile bending moments is approximately double the mean. For Mz, waves again have the most noticeable impact on the magnitude and vertical profile of the bending moment. In particular, increasing the wave strength (bottom panel of Fig. 17(b)) decreases the bending moments by 10%–15% through the entire domain. Changes in wind shear (Fig. 17(a)) and stability (Fig. 17(c)) result in minimal changes to the bending moment, while changes to the mean velocity (Fig. 17(d)) again result in large changes to the magnitude of Mz.

FIG. 17.

Probability density functions (pdfs) of eddy bending moments Mz for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. The top row shows Wi, Wa, In, and T and the bottom row shows Wi+, Wa+, In+, and T+. The median bending moments for each simulation case are shown as solid black lines, and 10th and 90th percentiles are shown as black dashed-dotted lines. Baseline results for the median moments are shown as black dashed lines. The logarithm of pdf magnitudes is shown for clarity.

FIG. 17.

Probability density functions (pdfs) of eddy bending moments Mz for (a) Wi±, (b) Wa±, (c) In±, and (d) T±. The top row shows Wi, Wa, In, and T and the bottom row shows Wi+, Wa+, In+, and T+. The median bending moments for each simulation case are shown as solid black lines, and 10th and 90th percentiles are shown as black dashed-dotted lines. Baseline results for the median moments are shown as black dashed lines. The logarithm of pdf magnitudes is shown for clarity.

Close modal

In Sections II–IV, methodologies for simulating flows in ocean tidal channels have been outlined and results pertaining to loading of ocean current turbines in these environments have been presented. These methods and results have implications for turbine designers, observational campaign managers, and modelers; implications of the present results for each of these three groups are outlined in more detail below.

Based upon the present results, the most important consideration for a turbine designer is the choice of hub height. Turbines mounted near the ocean floor not only see much larger bending moments (Fig. 15) but also see much lower velocities (Fig. 7(a)), while turbines higher in the ocean boundary layer harness power from much higher velocity flows with lower magnitude bending moments. One of the main drivers of turbulence near the ocean surface—strong waves—have been shown to actually slightly decrease the bending moment at moderate depths (Section IV D); it should be noted, however, that this analysis only considers bending loads and not transient torque loads, which may depend differently on waves and other physical characteristics.

Effects of length scale and structure function magnitude have also been shown to combine in interesting ways that warrant additional investigation by turbine designers. The inclusion of stronger waves in the simulations results in a large increase in the transverse length scale, which is a likely contributor to overall decreased bending moments at many depths for the strong wave simulation. Design flexibility in rotor diameter could allow for an opportunity to decrease overall loads by ensuring rotor diameter is of a different order of magnitude than the characteristic scale of turbulent eddies.

Finally, turbine designers should note that the underlying drivers of bending moments have been shown to deviate from a normal distribution (Fig. 14). Extreme loads will occur more frequently than a Gaussian distribution would imply. Designers should account for these higher probabilities in assessing turbine loads.

Observational campaigns are critical for understanding the specific environment in which a proposed turbine will be installed. These campaigns allow for very accurate understanding of mean velocities and the primary component of turbulence kinetic energy. Of all physical properties investigated in this paper, the mean tidal velocity has shown the largest impact on overall bending moment magnitude.

Through the analytical derivation (Section II), loads have been shown to tie closely to two point correlations and structure functions. For these properties to be understood in a specific installation environment, two point measurements are needed—a single point measurement device cannot capture transverse length scales and structure functions. While capturing a variety of measurement widths and physical conditions will help in understanding bending moments, the present study indicates that priority should be given to capturing information on transverse lengths of the same scale as the turbine diameter.

Finally, as waves have been shown to have a significant impact on bending moments (Section IV D), campaign managers should consider capturing wave specific data, such as amplitude and period, to accompany data on mean velocities and turbulence kinetic energy. These wave properties will help designers and modelers better understand the turbulent ocean environment. Stability was found to be less important in understanding loads for this idealized analysis, but could be of substantial importance when understanding wake recovery in arrays of real-world turbines.

Tidal channel simulations have the potential to supplement and enhance measurement campaigns in the design process. However, for maximum effect, simulations should include the physics most relevant to ocean current boundary conditions.

In the present analysis of bending moments, it was found that the most important variables to consider in the simulation process were mean tidal velocity and waves. Simulations should assess a variety of tidal velocities that are anticipated at the installation site, not just the mean velocity. Waves were found to be important, particularly in the top 70% of the boundary layer, and should be included in simulation physics.

One parameter that was found to have little effect on the resulting bending moments was stability. The present set of simulations use the Boussinesq hypothesis (as described in Section III), but for estimating mean and 90th percentile loads, stability effects were negligible. However, this study has only investigated moderate levels of surface heating and cooling; for areas with more substantial heating and cooling effects, additional analysis is required. Stability could also be of primary importance when considering wake recovery in turbine arrays.

Comparisons of the simulation results to the measurement campaign by Thomson et al.41 show that the simulations slightly underestimate the turbulence intensity profile in the tidal channel. Possible explanations for this include the influences of bathymetry and coastal features on the inlet and bottom boundaries. In the case of a complex surrounding environment, free stream turbulence in the simulation may need to be supplemented to account for these non-local effects.

Finally, it was found in the simulations that high levels of grid anisotropy (the ratio of Δx,y to Δz) resulted in substantial deviations from a log-layer profile. While the ocean channel environment lends itself well to anisotropic grids (since the horizontal domain length is typically much larger than the vertical depth), efforts should be made to either keep the grid anisotropy below four or to adapt the SGS model to account for the grid anisotropy.

A series of large eddy simulations have been performed in order to better understand the tidal channel environment and turbulent bending moments experienced by ocean current turbine rotors, with a focus on the physical parameters that impact these moments. The simulations were validated against theoretical predictions as well as compared to observational data from the Admiralty Head and Nodule Point inlets off the coast of Washington. In addition to simulation results and statistics related to the understanding of bending moments, an analytical derivation was presented that explicitly links two-point turbulence statistics to bending moments.

The simulations were analyzed to provide information on mean velocity, Reynolds stresses, longitudinal and transverse correlations, length scales, structure functions, and bending moment distributions. Bending moments were calculated for an idealized infinitely thin circular rotor disc and vertical and horizontal bending moments were found to be of the same order of magnitude at moderate depths, although horizontal bending moments increased substantially near the ocean floor. The mean tidal velocity profile and wave strength were found to be the most important physical factors in determining the bending moment magnitude and distribution.

In future work, extending the analysis to include statistics and distributions on transient torque loads is an important task. The analysis could also be expanded to run the resulting simulation velocity snapshots through an advanced turbine simulator to provide statistics on all forms of turbine loads. As with off-axis loads, torque loads can have negative effects on gearbox and drivetrain components,5 and also reflect intermittent variations in the speed of the incoming flow, which can affect blade, tower, and other structural loads.3 Consistent with prior studies2,3,6,7 for wind turbines, the present study indicates that intermittency can also play an important role in determining bending moments in tidal boundary layers. As opposed to the idealized analysis of turbine loading presented here, using velocity snapshots from the current study as input to an advanced turbine simulator will therefore provide a more complete picture of the blade, drivetrain, and tower loads experienced by realistic ocean current turbines. Additionally, comparisons against other observational campaigns could further validate the simulation parameters while also providing guidance on improved modeling of inflow turbulence conditions for different types of environment. Finally, future work will include realistic models for tidal turbines in order to examine wake effects in turbine arrays. The present study has specifically focused on the effects of different physical sources of turbulence on bending moments, for which it was convenient to remove turbine-turbine interactions, but waking of one turbine by another is a primary real-world contribution to bending moments that deserves further study.

Helpful discussions with Dr. Peter Sullivan, Dr. Baylor Fox-Kemper, Dr. Jim Thomson, and Dr. Katherine McCaffrey are gratefully acknowledged. The work of Katherine Smith in performing final test simulations is also greatly appreciated. The authors were partially supported by NSF 1258995. This research utilized the Janus supercomputer, which was supported by NSF (Award No. CNS-0821794) and the University of Colorado at Boulder. The Janus supercomputer is a joint effort of the University of Colorado at Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research.

1.
K.
Burman
and
A.
Walker
, “
Ocean energy technology overview
,”
Technical Report No. DOE/GO-102009-2823, U.S. Department of Energy, Federal Energy Management Program
,
2009
.
2.
N. D.
Kelley
,
B. J.
Jonkman
,
G. N.
Scott
,
J. T.
Bialasiewicz
, and
L. S.
Redmond
, “
The impact of coherent turbulence on wind turbine aeroelastic response and its simulation
,” in
Windpower 2005 Conference Proceedings
(
2005
).
3.
T.
Mücke
,
D.
Kleinhaus
, and
J.
Peinke
, “
Atmospheric turbulence and its influence on the alternating loads on wind turbines
,”
Wind Energy
14
,
301
316
(
2011
).
4.
A.
Sathe
,
J.
Mann
,
T.
Barlas
,
W. A. A. M.
Bierbooms
, and
G. J. W.
van Bussel
, “
Influence of atmospheric stability on wind turbine loads
,”
Wind Energy
16
,
1013
(
2013
).
5.
F.
Oyague
, “
Gearbox modeling and load simulation of a baseline 750-kW wind turbine using state-of-the-art simulation codes
,”
NREL Technical Report No. NREL/TP-500-41160
,
2009
.
6.
A.
Morales
,
M.
Wächter
, and
J.
Peinke
, “
Characterization of wind turbulence by higher-order statistics
,”
Wind Energy
15
,
391
406
(
2012
).
7.
P.
Milan
,
M.
Wächter
, and
J.
Peinke
, “
Turbulent character of wind energy
,”
Phys. Rev. Lett.
110
,
138701
(
2013
).
8.
L. P.
Chamorro
and
F.
Porté-Agel
, “
A wind-tunnel investigation of wind-turbine wakes: Boundary-layer turbulence effects
,”
Boundary-Layer Meteorol.
132
(
1
),
129
149
(
2009
).
9.
M.
Calaf
,
C.
Meneveau
, and
J.
Meyers
, “
Large eddy simulation study of fully developed wind-turbine array boundary layers
,”
Phys. Fluids
22
,
015110
(
2010
).
10.
J.
Meyers
and
C.
Meneveau
, “
Large eddy simulations of large wind-turbine arrays in the atmospheric boundary layer
,” AIAA Paper 2010-827,
2010
.
11.
F.
Porté-Agel
,
Y.-T.
Wu
,
H.
Lu
, and
R. J.
Conzemius
, “
Large-eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms
,”
J. Wind Eng. Ind. Aerodyn.
99
,
154
168
(
2011
).
12.
M. J.
Churchfield
,
S.
Lee
,
P. J.
Moriarty
,
L. A.
Martinez
,
S.
Leonardi
,
G.
Vijayakumar
, and
J. G.
Brasseur
, “
A large-eddy simulation of wind-plant aerodynamics
,” AIAA Paper 2012-537,
2012
.
13.
A. W.
Lavely
,
G.
Vijayakumar
,
J. G.
Brasseur
,
E. G.
Paterson
, and
M. P.
Kinzel
, “
Comparing unsteady loadings on wind turbines using TurbSim and LES flow fields
,” AIAA Paper 2012-0818,
2012
.
14.
J.
Meyers
and
C.
Meneveau
, “
Flow visualization using momentum and energy transport tubes and applications to turbulent flow in wind farms
,”
J. Fluid Mech.
715
,
335
358
(
2013
).
15.
D.
Mehta
,
A. H.
van Zuijlen
,
B.
Koren
,
J. G.
Holierhoek
, and
H.
Bijl
, “
Large eddy simulation of wind farm aerodynamics: A review
,”
J. Wind Eng. Ind. Aerodyn.
133
,
1
17
(
2014
).
16.
J. M.
Jonkman
and
M. L.
Buhl
, Jr.
, “
FAST user's guide
,”
Technical Report No. NREL/EL-500-38230, National Renewable Energy Laboratory
,
2005
.
17.
S.
Gant
and
T.
Stallard
, “
Modelling a tidal turbine in unsteady flow
,” in
Proceedings of the Eighteenth (2008) International Offshore and Polar Engineering Conference
(
2008
), pp.
473
479
.
18.
W. M. J.
Batten
,
M. E.
Harrison
, and
A. S.
Bahaj
, “
Accuracy of the actuator disc-RANS approach for predicting the performance and wake of tidal turbines
,”
Philos. Trans. R. Soc. A
371
,
20120293
(
2013
).
19.
M. E.
Harrison
,
W. M. J.
Batten
,
L. E.
Myers
, and
A. S.
Bahaj
, “
Comparison between CFD simulations and experiments for predicting the far wake of horizontal axis tidal turbines
,”
IET Renewable Power Gener.
4
,
613
627
(
2010
).
20.
J. H.
Lee
,
S.
Park
,
D. H.
Kim
,
S. H.
Rhee
, and
M.-C.
Kim
, “
Computational methods for performance analysis of horizontal axis tidal stream turbines
,”
Appl. Energy
98
,
512
523
(
2012
).
21.
X.
Sun
,
J. P.
Chick
, and
I. G.
Bryden
, “
Laboratory-scale simulation of energy extraction from tidal currents
,”
Renewable Energy
33
,
1267
1274
(
2008
).
22.
S. R.
Turnock
,
A. B.
Phillips
,
J.
Banks
, and
R.
Nicholls-Lee
, “
Modelling tidal current turbine wakes using a coupled RANS-BEMT approach as a tool for analysing power capture of arrays of turbines
,”
Ocean Eng.
38
,
1300
1307
(
2011
).
23.
R.
Malki
,
A. J.
Williams
,
T. N.
Croft
,
M.
Togneri
, and
I.
Masters
, “
A coupled blade element momentum-computational fluid dynamics model for evaluating tidal stream turbine performance
,”
Appl. Math. Modell.
37
,
3006
3020
(
2013
).
24.
M. J.
Churchfield
,
Y.
Li
, and
P. J.
Moriarty
, “
Large-eddy simulation study of wake propagation and power production in an array of tidal-current turbines
,”
Philos. Trans. R. Soc. A
371
,
20120421
(
2013
).
25.
M. G.
Gebreslassie
,
G. R.
Tabor
, and
M. R.
Belmont
, “
Investigation of the performance of a staggered configuration of tidal turbines using CFD
,”
Renewable Energy
80
,
690
698
(
2015
).
26.
S.
Kang
,
I.
Borazjani
,
J. A.
Colby
, and
F.
Sotiropoulos
, “
Numerical simulation of 3D flow past a real-life marine hydrokinetic turbine
,”
Adv. Water Resour.
39
,
33
43
(
2012
).
27.
S.
Kang
,
A.
Lightbody
,
C.
Hill
, and
F.
Sotiropoulos
, “
High-resolution numerical simulation of turbulence in natural waterways
,”
Adv. Water Resour.
34
(
1
),
98
113
(
2011
).
28.
I.
Afgan
,
J.
McNaughton
,
S.
Rolfo
,
D. D.
Apsley
,
T.
Stallard
, and
P.
Stansby
, “
Turbulent flow and loading on a tidal stream turbine by LES and RANS
,”
Int. J. Heat Fluid Flow
43
,
96
108
(
2013
).
29.
T. P.
Lloyd
,
S. R.
Turnock
, and
V. F.
Humphrey
, “
Assessing the influence of inflow turbulence on noise and performance of a tidal turbine using large eddy simulations
,”
Renewable Energy
71
,
742
754
(
2014
).
30.
A. S.
Bahaj
,
W. M. J.
Batten
, and
G.
McCann
, “
Experimental verifications of numerical predictions for the hydrodynamic performance of horizontal axis marine current turbines
,”
Renewable Energy
32
,
2479
2490
(
2007
).
31.
A. S.
Bahaj
,
A. F.
Molland
,
J. R.
Chaplin
, and
W. M. J.
Batten
, “
Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank
,”
Renewable Energy
32
,
407
426
(
2007
).
32.
L.
Myers
and
A. S.
Bahaj
, “
Near wake properties of horizontal axis marine current turbines
,” in
Proceedings of the 8th European Wave and Tidal Energy Conference
(
2009
), pp.
558
565
.
33.
B.
Gaurier
,
P.
Davies
,
A.
Deuff
, and
G.
Germain
, “
Flume tank characterization of marine current turbine blade behaviour under current and wave loading
,”
Renewable Energy
59
,
1
12
(
2013
).
34.
L. P.
Chamorro
,
C.
Hill
,
S.
Morton
,
C.
Ellis
,
R. E. A.
Arndt
, and
F.
Sotiropoulos
, “
On the interaction between a turbulence open channel flow and an axial-flow turbine
,”
J. Fluid Mech.
716
,
658
670
(
2013
).
35.
L. P.
Chamorro
,
D. R.
Troolin
,
S.-J.
Lee
,
R. E. A.
Arndt
, and
F.
Sotiropoulos
, “
Three-dimensional flow visualization in the wake of a miniature axial-flow hydrokinetic turbine
,”
Exp. Fluids
54
,
1459
(
2013
).
36.
E.
Fernandez-Rodriguez
,
T. J.
Stallard
, and
P. K.
Stansby
, “
Experimental study of extreme thrust on a tidal stream rotor due to turbulent flow and with opposing waves
,”
J. Fluids Struct.
51
,
354
361
(
2014
).
37.
P.
Mycek
,
B.
Gaurier
,
G.
Germain
,
G.
Pinon
, and
E.
Rivoalen
, “
Experimental study of the turbulence intensity effects on marine current turbines behaviour. Part I: One single turbine
,”
Renewable Energy
66
,
729
746
(
2014
).
38.
I. A.
Milne
,
A. H.
Day
,
R. N.
Sharma
, and
R. G. J.
Flay
, “
Blade loading on tidal turbines for uniform unsteady flow
,”
Renewable Energy
77
,
338
350
(
2015
).
39.
T.
Stallard
,
T.
Feng
, and
P. K.
Stansby
, “
Experimental study of the mean wake of a tidal stream rotor in a shallow turbulent flow
,”
J. Fluids Struct.
54
,
235
246
(
2015
).
40.
P.
Mycek
,
B.
Gaurier
,
G.
Germain
,
G.
Pinon
, and
E.
Rivoalen
, “
Experimental study of the turbulence intensity effects on marine current turbines behaviour. Part II: Two interacting turbines
,”
Renewable Energy
68
,
876
892
(
2014
).
41.
J.
Thomson
,
B.
Polagye
,
V.
Durgesh
, and
M. C.
Richmond
, “
Measurements of turbulence at two tidal energy sites in Puget Sound, WA
,”
IEEE J. Oceanic Eng.
37
(
3
),
363
374
(
2012
).
42.
K.
McCaffrey
,
B.
Fox-Kemper
,
P. E.
Hamlington
, and
J.
Thomson
, “
Characterization of turbulence anisotropy, coherence, and intermittency at a prospective tidal energy site: Observational data analysis
,”
Renewable Energy
76
,
441
453
(
2015
).
43.
I.
Langmuir
, “
Surface motion of water induced by wind
,”
Science
87
,
119
123
(
1938
).
44.
J. C.
McWilliams
,
P. P.
Sullivan
, and
C. H.
Moeng
, “
Langmuir turbulence in the ocean
,”
J. Fluid Mech.
334
(
1
),
1
30
(
1997
).
45.
P. E.
Hamlington
,
L. P.
Van Roekel
,
B.
Fox-Kemper
,
K.
Julien
, and
G. P.
Chini
, “
Langmuir-submesoscale interactions: Descriptive analysis of multiscale frontal spin-down simulations
,”
J. Phys. Oceanogr.
44
,
2249
2272
(
2014
).
46.
S.
Li
,
M.
Li
,
G. P.
Gerbi
, and
J.-B.
Song
, “
Roles of breaking waves and Langmuir circulation in the surface boundary layer of a coastal ocean
,”
J. Geophys. Res.: Oceans
118
,
5173
5187
(
2013
).
47.
C.-H.
Moeng
, “
A large-eddy-simulation model for the study of planetary boundary-layer turbulence
,”
J. Atmos. Sci.
41
(
13
),
2052
2062
(
1984
).
48.
L. P.
Van Roekel
,
B.
Fox-Kemper
,
P. P.
Sullivan
,
P. E.
Hamlington
, and
S. R.
Haney
, “
The form and orientation of Langmuir cells for misaligned winds and waves
,”
J. Geophys. Res.-Oceans
117
,
C05001
(
2012
).
49.
M.
Li
,
S.
Radhakrishnan
,
U.
Piomelli
, and
W. R.
Geyer
, “
Large eddy simulation of the tidal cycle variations of an estuarine boundary layer
,”
J. Geophys. Res.
115
,
C08003
, doi: (
2010
).
50.
J. R.
Taylor
,
S.
Sarkar
, and
V.
Armenio
, “
Large eddy simulation of stably stratified open channel flow
,”
Phys. Fluids
17
,
116602
(
2005
).
51.
B.
Gayen
,
S.
Sarkar
, and
J. R.
Taylor
, “
Large eddy simulation of a stratified boundary layer under an oscillatory current
,”
J. Fluid Mech.
643
,
233
266
(
2010
).
52.
L. Y.
Pao
and
K. E.
Johnson
, “
A tutorial on the dynamics and control of wind turbines and wind farms
,” in
American Control Conference, 2009. ACC'09
(
2009
), pp.
2076
2089
.
53.
J. H.
Laks
,
L. Y.
Pao
, and
A. D.
Wright
, “
Control of wind turbines: Past, present, and future
,” in
American Control Conference, 2009. ACC'09
(
2009
), pp.
2096
2103
.
54.
A. S.
Bahaj
,
W. M. J.
Batten
,
A. F.
Molland
, and
J. R.
Chaplin
, “
Experimental investigation into the hydrodynamic performance of marine current turbines
,”
Sustainable Energy Series Report No. 3
,
2005
.
55.
L.
Myers
and
A. S.
Bahaj
, “
Simulated electrical power potential harnessed by marine current turbine arrays in the Alderney Race
,”
Renewable Energy
30
(
11
),
1713
1731
(
2005
).
56.
A. S.
Bahaj
and
L.
Myers
, “
Analytical estimates of the energy yield potential from the Alderney Race (channel islands) using marine current energy converters
,”
Renewable Energy
29
(
12
),
1931
1945
(
2004
).
57.
R.
McSherry
,
J.
Grimwade
,
I.
Jones
,
S.
Mathias
,
A.
Wells
, and
A.
Mateus
, “
3D CFD modelling of tidal turbine performance with validation against laboratory experiments
,” in
Proceedings of the 9th European Wave and Tidal Energy Conference
, University of Southampton, UK (
2011
).
58.
A. D. D.
Craik
and
S.
Leibovich
, “
A rational model for Langmuir circulations
,”
J. Fluid Mech.
73
,
401
426
(
1976
).
59.
I.
Gjaja
and
D. D.
Holm
, “
Self-consistent Hamiltonian dynamics of wave mean-flow interaction for a rotating stratified incompressible fluid
,”
Physica D
98
,
343
378
(
1996
).
60.
D. D.
Holm
, “
The ideal Craik-Leibovich equations
,”
Physica D
98
(
2–4
),
415
441
(
1996
).
61.
J. C.
McWilliams
,
J. M.
Restrepo
, and
E. M.
Lane
, “
An asymptotic theory for the interaction of waves and currents in coastal waters
,”
J. Fluid Mech.
511
,
135
178
(
2004
).
62.
N.
Suzuki
and
B.
Fox-Kemper
, private communication (2015).
63.
P. P.
Sullivan
,
J. C.
McWilliams
, and
C.-H.
Moeng
, “
A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows
,”
Boundary-Layer Meteorol.
71
(
3
),
247
276
(
1994
).
64.
M. A.
Donelan
,
J.
Hamilton
, and
W. H.
Hui
, “
Directional spectra of wind-generated waves
,”
Philos. Trans. R. Soc. London Ser. A
315
(
1534
),
509
562
(
1985
).
65.
A.
Webb
and
B.
Fox-Kemper
, “
Wave spectral moments and Stokes drift estimation
,”
Ocean Modell.
40
(
3–4
),
273
288
(
2011
).
66.
B.
Polagye
and
J.
Thomson
, “
Tidal energy resource characterization: Methodology and field study in Admiralty Inlet, Puget Sound, WA (USA)
,”
Proc. IMechE Part A: J. Power Energy
227
,
352
(
2013
).
67.
P. P.
Sullivan
and
E. G.
Patton
, “
The effect of mesh resolution on convective boundary layer statistics and structures generated by large-eddy simulation
,”
J. Atmos. Sci.
68
(
10
),
2395
(
2011
).
68.
M. J.
Churchfield
,
Y.
Li
, and
P. J.
Moriarty
, “
A large-eddy simulation study of wake propagation and power production in an array of tidal current turbines
,”
NREL, Technical Report No. CP-5000-51765
,
2011
.
69.
J. C.
McWilliams
and
B.
Fox-Kemper
, “
Oceanic wave-balanced surface fronts and filaments
,”
J. Fluid Mech.
730
,
464
490
(
2013
).
70.
A. N.
Kolmogorov
, “
The local structure of turbulence in incompressible fluid for very large Reynolds number
,”
Dokl. Akad. Nauk SSSR
30
,
299
303
(
1941
).
71.
European Commission, Directorate General for Research and European Commission, Directorate General for Research, Energy, New and Renewable Energy Sources, and IT Power (Organization), “SEAFLOW: World's first pilot project for the exploitation of marine currents at a commercial scale: Contact JOR3-CT98-0202: Final Publishable Report,” Office for Official Publications of the European Communities, 2005.
72.
OpenHydro, http://fundyforce.ca/technology/openhydro-nova-scotia-power/ for Fundy Ocean Research Center for Energy,
2015
.
73.
SeaGen-S 2 MW, MCT Product Brochure, Marine Current Turbines,
2013
.