Analysis of turbulent bending moments in tidal current boundary layers

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%.¹ 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.¹ 

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 turbulence² 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.⁵ 

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.,³ Morales et al.,⁶ and Milan et al.⁷ have connected extreme and fatigue turbine loads to the intermittent nature of turbulent wind speed fluctuations, and Kelley et al.² 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é-Agel⁸ 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 studies^9–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.¹⁵ for a review of LES methods applied to wind farm aerodynamics). Calaf et al.⁹ 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.¹⁶ 

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.²⁴ and Gebreslassie et al.²⁵ 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.,²⁸ as well as the study of tidal turbine noise generation by Lloyd et al.²⁹ 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 turbines^30–39 and, more recently, two turbines.⁴⁰ 

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.⁴¹ The observational measurements have been used previously to examine turbulence intermittency, coherence, and anisotropy⁴² 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 turbulence^43–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.⁴⁶ 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.⁴⁹ 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.⁵⁰ and Gayen et al.⁵¹ 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

where r = [r_x, r_y, r_z] is the distance from the central point to a location on the disc, u = [u_x, u_y, u_z] is the velocity vector, u_x 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, r_x = 0 and Eq. (1) can be expanded as

The component M_x is about the axis normal to the turbine disc and provides a torque that either accelerates or decelerates the turbine. The components M_y and M_z 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 M_y and M_z 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.⁵ 

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+u′x$⁠, 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)

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

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

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), M_y 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 M_y is termed the mean shear bending moment. By contrast, Eq. (6) indicates that M_z 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 by^58–62 

where u = [u_x, u_y, u_z] is the three-dimensional Eulerian flow velocity, $ω≡∇×u$ is the Eulerian vorticity, u_L = u + u_s is the Lagrangian velocity, u_s is the wave-induced Stokes drift velocity, p is the pressure normalized by the background density ρ₀, b is the buoyancy, f_c 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 u_L and through the forcing term u_L⋅(∇u_s)^⊺ on the right-hand side of Eq. (9), both of which include the wave-induced Stokes-drift velocity u_s. Note that in the present study Coriolis effects are not considered; a similar approach has been used in prior studies of tidal boundary layers⁴⁹ 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,

where β_T = 2 × 10⁻⁴ K⁻¹ is the thermal expansion coefficient and salinity effects have been neglected. The buoyancy is $b=bẑ$ where b = −gρ/ρ₀. Parameter values used to obtain the buoyancy are ρ₀ = 1000 kg/m³, θ₀ = 290.16 K, and g = 9.81 m/s².

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.”⁴³ 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.,⁶³ 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.,⁴⁹ future work may include temporal variations in f_c in order to consider situations in which the tidal stream changes magnitude and direction. Here, however, a time-independent f_c 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 u_s appearing in the Lagrangian velocity u_L. The Stokes drift velocity is modeled as

where u_s(z) is the Stokes drift vertical profile and $ϑs$ is the wave direction. Using the Donelan empirical spectrum^64,65 for the ocean wave field, the resulting vertical profile u_s(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 u_s(z) are determined by setting the 10 m wind speed, denoted U₁₀, in the Donelan spectrum. Profiles of u_s(z) used in the present simulations are shown in Fig. 2. These profiles correspond to U₁₀ = 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.⁶⁶ 

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., u_x and u_y) 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 u_z 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, Q₀, 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 z₀ = 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.⁶³). 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.⁴¹ 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 L_x = 600 m and L_y = 300 m with depth H = 39 m, a tidal velocity of U_t = 2 m/s, upper wind forcing resulting in u_τ = 6.3 mm/s, waves created by a 10 m wind of U₁₀ = 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 N_x = 512 and N_y = 256, resulting in a horizontal grid resolution of Δ_x,y = 1.2 m. The vertical grid size is N_z = 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,⁴⁴ although it may be larger than the values found in typical tidal channels.

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 L_x, L_y, N_x, and N_y 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.⁴¹ 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 U_t result in large changes in flow properties, and so the T₋ and T₊ simulations correspond to values of U_t that are only 1 m/s smaller and larger, respectively, than the baseline value of U_t = 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 u_s(z) used for the B₊, B₋, and Wa₊ simulations are shown in Fig. 2. Finally, ocean surface cooling of Q₀ = −5.0 W/m² has been found previously⁴⁵ 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 u_x 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 u_x 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 [L_x, L_y] = [600m, 300 m].

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.⁴¹ Prior work by Sullivan and Patton in the atmospheric convective boundary layer⁶⁷ indicates that mesh sizes of 256³ in a domain of size (L_x, L_y, L_z) = (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, $ux′2¯>uy′2¯>uz′2¯$ and $u′xu′z¯<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 atmosphere^47,67 and tidal channels.⁶⁸ 

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 as^63,68

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 studies^63,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

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 $u′xu′z¯$ 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.⁴¹ are made using the TI, which is defined as

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).⁴¹ 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 $ux′2¯$ increases and $u¯x$ decreases (note that immediately at the bottom boundary TI = 0 since $ux′2¯=0$⁠). While the roughness length z₀ used in the simulations (see Table I) was chosen to give a close match to the observed TI profiles from Thomson et al.,⁴¹ 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.

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.⁴⁵ and explained by McWilliams and Fox-Kemper.⁶⁹ 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.

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 $ux′2¯$ compared to the base case, although little change is seen in either $uz′2¯$ or $u′xu′z¯$ (note also that the increase in $ux′2¯$ is relatively small for In₊). An increase in $ux′2¯$ 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 $ux′2¯$ Reynolds stress, while halving the mean velocity results in roughly a factor of four decrease in $ux′2¯$⁠.

The mean temperature profiles remain close to the background value θ₀ 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₋), $u′xθ′¯<0$⁠, corresponding to an anti-correlation between $u′x$ and $θ′$⁠; while for the stable case (In₊), $u′xθ′¯>0$⁠. In Fig. 10(c), $u′zθ′¯>0$ for the unstable case, corresponding to a downward flux of low temperature fluid as the surface is cooled; while $u′zθ′¯<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 $u′xθ′¯$ and $u′yθ′¯$ 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.

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, f₁₁, and transverse velocity correlation, f₁₂, as functions of height, where the correlation functions are given by

and the corresponding integral scales Λ₁₁ and Λ₁₂ are defined here as

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.

For the base case simulation B, Fig. 11 shows that over the full tidal channel depth the longitudinal integral scale Λ₁₁ is between roughly 20 and 50 m and the transverse integral scale Λ₁₂ 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 Λ₁₁(z) with increasing z/h. For the transverse correlation f₁₂ in Fig. 11(b), Λ₁₂(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 f₁₂, which was shown in Section II to be closely connected to the bending moment M_z, 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 f₁₂ occur as the wave forcing varies, with approximately a factor of 2 increase in Λ₁₂ 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%.

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 S_1j defined as

where x is a position in space and r_j 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 r^1/3 scaling, in accordance with classical Kolmogorov theory,⁷⁰ although the departures from this scaling are greatest for S₁₁ and S₁₂ near the bottom boundary. In both cases, the structure functions have reduced slope, indicating greater prominence of small scale motions relative to large scales.

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 M_z, 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.

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 M_y and eddy moments M_z. 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(10 m)$ real-world turbines, such as the 11 m diameter Seaflow turbine,⁷¹ the 16 m diameter OpenHydro turbine,⁷² and the 20 m diameter Seagen turbine.⁷³ 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 M_y 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.

Changes to the M_y 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 M_y 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 M_y, 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.

The eddy bending moment M_z, shown for the baseline case in Fig. 15(b), is driven by horizontal imbalances in u_x 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 M_z is similar to the magnitude of M_y. Much like M_y, the magnitude of the 90th percentile bending moments is approximately double the mean. For M_z, 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 M_z.

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.⁴¹ 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,⁵ and also reflect intermittent variations in the speed of the incoming flow, which can affect blade, tower, and other structural loads.³ Consistent with prior studies^2,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.