ACTUATOR-LINE CFD MODELLING OF TIDAL-STREAM TURBINES

CFD modelling of tidal turbines is described. For computational efficiency rotors are represented by rotating actuator lines and nacelles by partially-blocked cells. The computational implementation includes an arbitrary Lagrangian-Eulerian (ALE) approach to free-surface movement for wave motion. For a single turbine the model successfully reproduces towing-tank measurements of thrust and power coefficients across a range of TSR. Modelling of two turbines staggered streamwise shows that loads may be reduced or augmented, depending on whether the downstream turbine is in the wake or bypass flow of the upstream turbine. Where the downstream turbine is partially in the wake, individual blades suffer large cyclical load fluctuations. Turbine performance has been simulated under both regular and solitary waves. For long waves, time-varying thrust and power are found to be reasonably well predicted by scaling the non-wave case using the hub-height velocity (wave + current), provided that the reference is adjusted for changing TSR. For short waves this approach over-predicts fluctuating loads.


INTRODUCTION
Tidal resources have long been considered a promising source of renewable energy, offering high energy density and predictable generating periods. Public opposition to the expense and unknown environmental consequences of large barrages have led to tidal-stream devicespredominantly axial-flow turbines -which aim to extract kinetic energy from the tidal current rather than the potential energy built up by impounding. In numerous sites around the world, narrow straits lead to tidal currents in excess of 2.5 m s -1 , where tidal-stream energy becomes commercially viable [1]. Demonstration devices have been tested at the EMEC site in the Orkneys and the FORCE site in Nova Scotia's Bay of Fundy, whilst commercial arrays are under construction in the Pentland Firth off Scotland and Raz Blanchard off Normandy.
Heavily influenced by wind-turbine design, most tidal-stream devices are 3-bladed, horizontal-axis machines. Marine turbines, however, face many additional challenges: the need for a much more substantial nacelle and support; short deployment window; difficulty of access; bio-fouling; marine debris; cavitation; fluctuating loads due to sitespecific turbulence; velocity shear and waves.
Small-scale laboratory studies have been conducted in flumes and towing tanks ( [2], [3]), but these can cover only a limited range of operating conditions and are subject to scale effects. Theoretical models include blade-elementmomentum theory (BEMT), and, increasingly, computational fluid dynamics (CFD). CFD has been used to undertake geometry-resolved simulations of real devices in both low-turbulent flow [4] and realistic site turbulence and velocity shear [5]. However, the computational resources required to model satisfactorily both near and far wake, multiple devices, many operating conditions and interaction with waves is prohibitive. Instead, we have turned to the wind-energy approach ([6], [7]) of replacing a geometrically-resolved turbine ( Figure 1a) by an 'actuator' model -a body-force distribution that provides as closely as possible the same reactive forces as a real turbine. This may be done coarsely with a momentum (and, sometimes, angular momentum) source distributed uniformly over a swept volume to match total thrust and torque (Figure 1b), or, with extra computational expense but much better representation of near wake and no a priori assumption about total loads, by representing the blades as rotating actuator lines (Figure 1c).
Advantages of the rotating-actuator-line approach include a substantial reduction in computer resources, ease of changing rotor design or adding other turbines, and decoupling of turbine rotation from the background grid (allowing, in particular, simulation of waves). Disadvantages include assumptions that the reaction forces are the same as 2-d flow around aerofoils, absence of highfrequency turbulence generated by the blades, and low resolution of the nacelle and support tower. Also, in contrast to momentum sources spread over an actuator disk, rotating actuator lines require a time-dependent calculation.
In the remainder of the paper, Section 2 describes the computational model, including the finitevolume code, representation of rotor and nacelle, and simulation of waves. Section 3 describes validation with a single rotor, whilst Section 4 looks at interactions between multiple turbines. Section 5 simulates turbine behaviour in waves. Finally, Section 6 concludes with key findings and suggested future work.

THE CFD CODE (STREAM)
Calculations here used the in-house finite-volume code STREAM. This solves the Reynoldsaveraged Navier-Stokes (RANS) equations on multi-block, structured, curvilinear meshes. The standard k-ε turbulence model was used throughout. Parallelisation is by blockwise domain decomposition. A 3rd-order flux-limited scheme (UMIST) was used for advection, whilst the second-order Gear's scheme was used for timestepping. In contrast to many commercial codes, STREAM uses a moving-mesh, surfacefitting approach for the free surface (Section 2.4), rather than a volume-of-fluid (VOF) method.

TURBINE REPRESENTATION
Drawing on the actuator-line model of [6], with recommendations of [7], the turbine rotor is replaced by rotating lines of actuator points ( Figure 1c). The body forces associated with these are derived from 2-d aerofoil theory, knowing blade geometry (shape, chord and blade-angle) at each radius. At each inner iteration the local flow velocity U a at a point on a blade is compounded with the local blade velocity (Ωr tangentially) and projected onto the plane perpendicular to the blade to determine the relative velocity U rel . Given blade angle β(r) this determines angle of attack α and thence, using chord c(r) and blade-section lift and drag coefficients C L (α,r) and C D (α,r), an actuator force for radial length Δr of ) )( where D ê and L ê are unit vectors along and normal to the local relative velocity U rel (Figure 2). Each isolated point force is distributed over a wider volume of scale σ as volume force density: where x is cell centre and x p an actuator point.
Here, the search for influence is cut off at 2.5σ. [7] provides guidance on σ. In this paper a (fixed) value for σ of roughly 2 times local cell size was used, with 30 actuator points along each blade.
Look-up tables are created for chord c(r), blade angle β(r) and lift and drag coefficients C L (α,r) and C D (α,r). Manufacturer's data prescribes the first two; where aerofoil data is not available we have also developed vortex-panel and 2-d viscous solvers to determine the force coefficients.
Timestep-dependence tests yielded a suitable Δt = T rot /250, where T rot is the turbine rotation period, which varies with tip-speed ratio (TSR). With this value of Δt, the distance swept by a rotor tip in one timestep (ΩRΔt) was approximately the same as the grid spacing.

NACELLE REPRESENTATION
The nacelle for a tidal-stream turbine is far more substantial than a wind turbine. Ignoring it will fail to incorporate its blocking effect and provide a low-resistance route for streamlines to bypass the actuators representing the blades. However, creating a grid that resolves and aligns with the nacelle runs counter to the desire to decouple the turbine from the background flow grid.
In these simulations we adopt a blocking-out-cells technique, calculating for any cell the fraction f which is covered by a nacelle. Here, the latter are represented as cylindrical tubes of appropriate diameter and length. The overlapping fraction in any cross-stream plane is approximated by that for a polygon ( Figure 3).
To match dimensions and force each velocity component  to a small value in the required regions we introduce an implicit force density for a wholly or partially-covered cell of the form: where u* is a velocity scale significantly larger than any in the domain (say, 10 times inflow velocity) and L is a typical cell scale (e.g., V 1/3 ). Simulations showed that, for a bluff body, this approach generated a velocity and pressure field (including recirculating flow) very similar to those on a geometry-resolving multiblock mesh.

FREE-SURFACE MODELLING
In contrast to many commercial codes using the volume-of-fluid (VOF) method on a fixed grid, STREAM employs the arbitrary Lagrangian-Eulerian (ALE) approach to solve the finitevolume equations of motion on a free-surfaceconforming moving mesh [8].
Free-surface vertices are interpolated from intermediate control points (Figure 4a), which are iteratively adjusted over a timestep to satisfy nonet-volume-flux through each surface cell face ( Figure 4b):  Wave motions are input by specifying timevarying inlet surface elevation η(t) and velocity (u(t),v(t),0). On the outlet plane, wave reflection is minimised by solving a 1-d wave equation for all field variables and surface elevation. This is designed to allow forward-travelling waves with absolute speed U 0 + c to pass, but not backwardtravelling waves.

SINGLE TURBINE
Here, we simulate the experiments of [2] who made load measurements on a three-bladed axialflow turbine in a towing tank. The cross-section of the computational domain matched that of the tank. The upper surface was treated as a symmetry plane as surface displacements due to the presence of the turbine are small at low Froude numbers.
Blade twist β(r) and chord c(r) distributions were as specified in Bahaj et al.'s paper (for a 20º blade pitch angle), whilst aerofoil section coefficients c D (α,r) and c L (α,r) were generated by XFoil. The diameter and streamwise length of the nacelle were the same as those in the experiments, but no attempt was made to include the vertical support.
In design conditions (TSR  6), numerical tests showed that satisfactory grid independence for loads was obtained with a mesh extending 3D upstream and 7D downstream of the rotor (where D is rotor diameter), with a background mesh of 690000 cells, clustered at the centre of the rotor and allowed to expand uniformly in all directions. With a low-turbulence, uniform approach flow, the total number of cells was found to be less important than the size of cells in the vicinity of the rotor (here, 0.01D at rotor centre). For parallelisation on a multi-core desktop PC the domain was divided into 8 streamwise blocks. Figure 5 illustrates wake and vortex structures obtained at tip-speed ratio (TSR) = 6. Vortex structure is revealed by an isosurface of the vorticity indicator Q, where and S ij and Ω ij are strain and vorticity tensors.  Figure 6 compares the turbine thrust and power coefficients: across a range of TSR with the measurements of [2]. (Since power is torque  rotation rate, both C T and C P are referred to in this paper as "load" coefficients). Thrust and torque are directly related to momentum deficit and angular momentum in the near wake. In Figure 6 (but not our subsequent sections) corrections for the 7.5% blockage have been performed in the same way as [2], because their original uncorrected load coefficients were not provided. Computed thrust is slightly low, but power is well-predicted across most of the range tested. In particular, actuator-line computations demonstrate the significant fall-off in power at offdesign TSR that simple BEMT models [9] are not able to replicate.

MULTIPLE TURBINES
An actuator approach makes it straightforward to accommodate and redistribute additional turbines. Figure 7 shows predicted wakes for two turbines separated streamwise by 5D and laterally by 0, ½ and 1D (between centrelines). Figure 8 shows the resulting power and thrust coefficients. Turbine geometry and channel dimensions are the same as Section 3, but the computational domain has been extended downstream. Fixed rotation speeds with TSR of 6 based on onset flow far upstream is used for both turbines, although the downstream turbine will experience a very different onset velocity, dependent on its position relative to the upstream turbine. In the low-turbulence flow here, wakes develop slowly, so that turbine interactions are exaggerated; in real marine flows, turbulence and velocity shear promote more rapid wake recovery.  In all cases the upstream turbine experiences essentially the same load coefficients as the individual turbine of Section 3. However, the lateral displacement of the second turbine leads to significant differences in load.
In Figure 7a the second turbine is directly downstream of the first turbine and its load coefficients are substantially reduced, albeit slightly unrealistically: for a more turbulent flow, wakes would spread and recover more rapidly.
In Figure 7b the downstream turbine is half "shaded" by the upstream turbine. Whole-rotor power coefficients in Figure 8 show the expected reduction in thrust and mean power. However, Figure 9, which plots power coefficients for the whole rotor and an individual blade on each turbine (scaled by factor 3), indicates how individual blades on the downstream turbine rotate in and out of the upstream turbine's wake and experience huge cyclical changes in loading, with implications for fatigue damage.
At lateral spacing 1D ( Figure 7c) the downstream turbine actually experiences greater thrust and produces more power than the upstream turbine ( Figure 8). Here, the downstream turbine is effectively in the bypass flow of the upstream turbine. Because of the relatively large channel blockage ratio (0.075 for each turbine) this results in a higher onset velocity to the downstream turbine. This illustrates the possibility of enhancing power production in a narrow channel by judicious use of blockage effects [10]. Figure 9. Power coefficients for turbines staggered streamwise by 5D and laterally by ½D; also included, the power coefficient (3) for a particular blade.

TURBINES IN WAVES
To date, two wave types have been simulated: regular waves and a solitary wave.
For common rotor sizes in realistic water depths and current speeds, turbine rotation periods are smaller than typical mean wave periods and an initial approach to predicting changes in load might be to adopt a "quasi-static" approach, whereby time-varying thrust and power are scaled from the current-only case (adjusted for the change in TSR) according to the square or cube of the fluctuating onset velocity (either at nacelle height, or averaged over the rotor). However, this neglects important effects due to streamwise acceleration, wave-induced pressure differences, wave-induced vertical velocity shear and a nonzero vertical velocity component. A key objective of these simulations is to establish how accurate the quasi-static approach is.
The simulations presented below are nominally for a turbine of diameter D = 8 m in water depth 18 m and a current U 0 of 2 m s -1 . Geometry and relative sizes and positions are the same as those in Section 3, in order to have background non-wave data for comparison. Results are presented nondimensionally, using velocity scale U 0 , length scale D and density of fluid ρ. Assuming a TSR of 6 in a steady current, the rotation period is 2.09 s (non-dimensional period T rot U 0 /D = 0.524) and the inertia of the rotor is assumed to be such that it maintains a constant rotation rate (but not constant TSR) as a wave passes.

REGULAR WAVES
Linear wave theory admits surface displacements ) ω cos( ω is the wave frequency, and t U x x 0    the displacement, relative to the current U 0 . The frequency in an absolute frame (such as the computational domain) is ω a = ω +kU 0 . For all calculations here the wave celerity ( k c / ω  ) far exceeds the current speed. Here also, all simulations are for waves travelling in the same direction as the current.
Realistic mean wave periods are of order 8 s in a frame moving with a current in the same direction (or 6.75 s in a stationary frame, which is what is set at inlet). This is more than three times the turbine rotation period. To investigate the effect of the time-varying wave motions we consider both this wave period and also a shorter period of 4 s relative to current (3.03 s in a stationary frame). A summary of the relevant time periods and wavelengths/wavenumbers is given in Table 1 below, where λ is wavelength, k is wavenumber and d is water depth. "Long waves / shallow water" typically implies kd < π/10, whilst "short waves / deep water" usually refers to kd > π [11] Both wave sets fall between between these bounds; however, for convenience, we shall refer to them as "long" and "short" waves, the adjectives applying both to period and wavelength.   (b) short waves. Figure 11 shows power coefficients for the two regular wave cases. Superposed are power coefficients obtained by scaling a current-only case according to the time-varying hub-height velocity (current + wave), assuming that power scales as velocity cubed. Two scalings are shown: one with a fixed current-only power coefficient (C P = 0.493) and the other adjusting the currentonly power coefficient, using curve fits to the computed characteristics in Figure 6, to reflect the change in TSR. For the longer waves, scaling reproduces the computed C P variation satisfactorily, provided that the second scaling is used. For the shorter waves, however, power coefficients from the simulations produce significantly smaller excursions about mean power than predicted from the quasi-static assumption with either scaling. Possible reasons include wave-induced streamwise acceleration and pressure differences across the rotor. For these steeper waves there are also small differences in phase, due to linear wave theory under-predicting the phase velocity of the waves, which therefore reach the turbine slightly earlier than anticipated.

SOLITARY WAVE
Solitary wave theory [11] admits waves of fixed shape: Here, we consider a solitary wave of amplitude A = 2.4 m (= 0.3 D) passing the turbine. The wave celerity is 7.07 times the current velocity, the length of wave (taken between points where rise is 1% of maximum) is 5.99D, and this length passes over the turbine in 1.42 turbine rotations. Figure 12 shows the wave during its passage over the turbine. Turbine effects on the wave are negligible, but wave effects on the turbine are substantial, disrupting the vortex structures in the core of the wake and increasing wake spread. Figure 13 shows the transient thrust and power coefficients. The quasi-static approach provides acceptable agreement with both, but only if changes to TSR are reflected in the constantcurrent load coefficients before scaling by powers of velocity.

CONCLUSIONS
Rotating actuator lines are a computationallypractical means of combining real rotor geometries for one or more devices with complex fluid flows, notably waves. Applications include optimising the layout of multiple turbines, predicting the impact of complex flows (bathymetry, turbulence or waves) and devising control strategies.
Turbine characteristics have been successfully computed and validated against experimental data across a range of TSR. Two-turbine calculations have shown the (positive and negative) effects of interacting turbines. Wave calculations demonstrate significant fluctuations in load that must be absorbed in the power train for realistic wave periods. For wave periods significantly longer than one turbine rotation a quasi-static approach based on scaling load coefficients by the relative power of the onset velocity works well, provided that the constant-current coefficients from which they are scaled are continuously adjusted to reflect changing TSR. However, for short-period waves the quasi-static approach overpredicts loading, possibly as a result of streamwise flow acceleration or wave-induced pressure differences across the rotor.
Further work may include: • effects of onset turbulence and velocity shear; • realistic wave series, based on measured spectra; • fluctuating speeds and operational control of devices, with drive-train models.
Modelling improvements should include: • better understanding of 3-d (and transient) effects on the flow about blade sections; • body-source terms for turbulent kinetic energy as well as momentum.