CFD modelling of a tidal stream turbine subjected to profiled flow and surface gravity waves

This study used computational ﬂuid dynamics to investigate the effect of waves and a velocity proﬁle on the performance of a tidal stream turbine (TST). A full scale TST was transiently modelled operating near its maximum power point, and then subjected to waves both in and out of phase with its period of rotation. A proﬁle was then added to one of the wave models. For this set of conditions it was found that the longer period and in-phase wave had a signiﬁcant effect on the power range ﬂuctuations, with more modest variations for thrust and the average values, although this is dependent on the tur- bine tip speed ratio. The addition of the proﬁle had a strong effect on the bending moment. It has been concluded that a naturally varying sea state may yield a smoothing effect in this turbine response, but that with further structural investigation it may be that some mea-suring and mitigation techniques are required in the event of a predominantly single long period, in-phase wave.


Introduction
Tidal current and wave energy resources around the UK could provide 70 TWh/yr of electricity, contributing towards 20% of the total demand [1]. With increasing pressure to implement renewable energy measures across the EU towards a goal of 20% by 2020 [2], or 15% by 2020 for the UK [3], ocean resources will play a vital role if they are correctly exploited. Tidal power may offer specific advantages over other renewable energy technologies due to its predictable nature [4] and thereby offer some balancing contribution against more variable resources. Two main technology groups are in development; barrage or impoundment schemes such as the EDF Energy La Rance power plant [5], and tidal stream turbines, which generate electricity from the kinetic energy of the tidal stream. Advantages of this are more localised and less intrusive deployment, allowing an installation to be tailored to the local surroundings.
According to the UK's Department for Energy and Climate Change, tidal stream energy has considerable potential to play a part in combating the rising global energy crisis [6], although tidal stream turbine (TST) technology is still in the prototype stage, and has been shown to exhibit unprecedented operating and maintenance costs [7]. Tidal stream velocities exhibit considerable fluctuations [4], and may present operating conditions far in excess of those for which a device may be designed. Excessive loading may cause failure in the mechanical power transmission or electrical generating equipment, or even catastrophic failure of the turbine blades and components [8]. This effect is magnified for offshore applications, where access to devices is not arbitrary and requires expertise and expensive access equipment [9].
The conditions in which a horizontal axis TST may be expected to operate vary depending on the location. Other than the predictable diurnal tidal flow, they may be subjected to local turbulence, excessive storm loads and depth penetrating surface wave effects [10]. As opposed to the expected gravitational changes in tidal velocities, these other variations are harder to predict [11]. In addition, frictional effects with the seabed and water surface mean that a TST is actually subjected to a velocity profile rather than a uniform plug flow [9], which means that it is subjected to varying loads on each blade as it rotates. O'Doherty et al. [7] examined this effect with ADCP data from the Severn Estuary, and confirmed its existence, with local turbulence and oscillations also having an effect; subsequently Mason Jones et al. [12] verified that it has an effect on the turbine's loading characteristic and power output.
The UK has a very strong tidal energy resource and many potential deployment areas [1], with optimal sites having a free stream velocity of 2-3 m/s and a depth of 20-30 m. The technology is still in early development, with many prototype and demonstration projects being implemented in the UK, such as the DeltaStream project in Pembrokeshire [13], and the array, of SeaGen S turbines, in Anglesey [14]. Design decisions and economics can drastically benefit from information obtained before the prototype stage, which is where comparatively inexpensive numerical techniques such as computational fluid dynamics become crucial.
Much previous CFD work has assumed a uniform plug flow in a non-free surface body of fluid [15]. The effects of high shear profiles and turbulence has been carried out by Mason-Jones et al. [16] and Morris [17], and the next stage of development is to start examining the effect of surface waves on turbine performance. Surface gravity waves exhibit significant penetration into the water column of $50% of their wavelength [18], meaning their influence becomes greater with longer periods and as the turbine moves towards the surface. It is also possible that the wave may not be incident on the turbine perpendicular to the current flow direction, causing further modelling complications [19]. The aim of this study was to examine the effect of a wave and velocity profile on a tidal stream turbine, considering the case where the wave was both in and out of phase with the period of turbine rotation.

Modelling approach
A suitable representation of a TST deployed in 35 m of water was modelled with a domain of rectangular cross section of 50 m by 50 m, and 350 m in length. A free surface was defined at 35 m, with the turbine situated 100 m downstream of the inlet. The distances from the walls, inlet and outlet were located to preclude any interference or blockage effects [15]. This tidal stream turbine also ensures total flow development before interaction with the turbine and full wake development before the fluid exits the model [17]. A rotating cylindrical domain incorporating the turbine geometry was used to represent the rotation of a turbine due to the action of the tidal current. This model set up has been extensively applied and experimentally validated using a recirculating water flume and non-free surface (''closed") CFD models [12,20,21], as well as some free surface work [17].

Turbine modelling
The turbine model was developed from an experimental device [22] with a diameter of 0.5 m and based on a Wortmann FX 63-137 profile, with a 33o twist from the blade root to tip [23]. A support structure was added based on an existing design [24], resulting in the geometry in Fig. 1. The turbine was meshed with a tetrahedral scheme employing a higher element concentration around the blade tips. The rotational cylindrical domain comprised 2.8 million elements, as shown in Fig. 2 [24].

Volume of fluid technique
The majority of previous work [15,17,23,24] has been based on a closed model approach, neglecting gravity. In order to model surface effects, a volume of fluid technique was used, with activated buoyancy. To achieve mesh refinement in the region of the free surface, the 50 m by 50 m by 350 m domain was created with two separate boxes, the intersection of which was defined at the water depth of 35 m. The turbine model was incorporated into this domain as shown in Fig. 3. The surface of intersection between the two boxes was meshed with constant thickness inflation layers of total 4 m in each direction, to incorporate the maximum amplitude of the applied waves to ensure surface resolution. The remainder of the body section was meshed to yield just over 1 million elements, as shown in Fig. 4. The z direction is the rotational axis of the turbine and the y direction is the gravitational direction perpendicular to the free surface.

Wave and profile modelling
The waves and profiles were added to the model mathematically. The velocity profile was specified as a simple power law (Eq. (1)) that was assumed to be uniform across the width of the domain, and only varying in the direction of the water depth The waves were applied using a modified linear Airy wave theory to apply surface elevation, sub-surface orbital fluid motion, and to extrapolate the kinematics at the mean water height throughout the water column. This was because when the wave heights and water depth used in this study were non-dimensionalised and compared to Fig.5 [25] or similar representations [26], it was found that they lay outside of the linear regime, which assumes that the amplitude of the wave is very small. Rather than specifying a complex and time consuming Stokes model [27], a modification known as Wheeler stretching [28] was applied, which was found to yield an acceptable approximation.

Boundary conditions, viscous model and solution
The Reynolds Averaged Navier-Stokes (RANS) equations were used to relate the Reynolds Stresses to the mean velocity gradients, with closure achieved with the Shear Stress Transport (SST) viscous model. Previous studies have shown the Reynolds Stress Model (RSM) to be particularly suitable for TST modelling [15]; however more recent and complex modelling has used the SST model due to added complexity and computational expense due to transient [17], FSI [24] and free surface modelling [17], as well as a more representative turbulence dissipation rate [17].
The fluid entered the model perpendicular to the front face at a uniform velocity of 3.086 m/s, upon which the wave and profile were superimposed. The fluid exited the model through the rear face, which was set as having a corresponding downstream pressure profile due to the action of the fluid under the influence of gravity (Fig. 6). Both of these boundaries were specified as openings, which are able to tolerate potential reverse flow in the model; a likely occurrence due to the forward and back motion of the subsurface fluid. The seabed was set as a 'No Slip' boundary, and the side walls were specified to be 'Free Slip' to simulate a body of fluid in open sea. The top was also set as an opening to simulate the free surface.
The turbine domain was set to rotate about its axis at an angular velocity of 2.25 rad/s, which is located very close to peak power on the power curve for this particular design [15].

Velocity profile
The tenth power law velocity profile was based on site data [29] and has been used extensively in previous work [23]. To allow comparison of results between the profiled model and the plug flow model, the volumetric flow rate through the turbine was matched in both cases, which was achieved with an iterative methodology whereby the turbine swept area was split into segments, the flow rate through each calculated and a total volumetric flow rate calculated such that the average flow velocity across the turbine was 3.086 m/s [15,23]; resulting in Eq. (2).

Wave characterisation
Sea wave characterisation was carried out with real data, an example of which is shown in Fig. 7, from the British Oceanographic Database [19] in order to determine the parameters for a realistic but extreme case sea wave. This was achieved by analysing the significant periods and wave heights of various waves measured at multiple locations around the Welsh and West coasts.
Cases with a wave both in phase and out of phase with the rotational period of the turbine were investigated. With the turbine rotation specified as 2.25 rad/s, this equates to a multiple of the time taken for one rotation, which is $2.79 s. Due to added complexities in modelling very short period waves, in that the wave amplitude becomes to large compared to the wavelength to maintain approximate linear wave theory, it was decided that a wave period corresponding to two turbine rotations would be selected as the in phase wave would be generated, which gives a wave period of $5.58 s. For the out of phase wave, a wavelength of 30 m was selected that would allow sufficient depth penetration [18]. Calculation of the Once the period was selected or calculated from Eq. (4), an inverse tanh function based on a set of data pairs for (kh) and tanh(kh) was used to determine the wave number from Eq. (5), given that the right hand terms are known or can be calculated. The resulting conditions are shown in Table 1. kdtanhðkdÞ ¼

Results
The models were solved using Ansys CFX, and post processed in Ansys Post to obtain curves for thrust, power and bending moments. The following graphs are plotted against turbine angular position, with 0°defined as the position where the turbine blade 1 is in the top ''legs down" configuration as shown in Fig. 8.
This means that the sequence of blades passing the top in the legs down position is Blade 1 -Blade 3 -Blade 2. The results are presented in graphs split into 120o intervals, where each position is specified above the graph with use of a colour-coded turbine to demonstrate which blade is in the top location at that position. A dimensionless representation of the wave is depicted in black in the graphs to show the wave time period and the relative locations of peaks and troughs.

Result stabilisation
Figs. 9-12 show that there is a period of settling after which the results become consistent and can be analysed. This time is required for the flow field and waves to stabilise and for the results to become suitably converged. With this mesh and model set up, the residual error came down to 10-6 for each time step; an acceptable level [32]. Sections from the converged and stabilised portions of the graphs were analysed, with the x axis showing the actual angular position at the time. The exception is for the combined wave and profile model (OP + P), as it was found that the combination caused considerable instability and required much longer to settle into a consistent pattern. Therefore, the first 19 rotations of the turbine are omitted from the results and the rotational angle reset to 0 at this point as shown in Fig. 11.

Steady state power curves
A series of steady state models were run at varying tip speed ratios (TSR) to produce a total power curve for the turbine design, which matched previous studies very well [15,23]. This plot is shown in Fig. 13, with the average power from the stabilised portion of the transient curves added at the TSR corresponding to a turbine angular velocity of 2.25 rad/s. Fig. 14 shows the averages of the transient models along with the range from the minimum to the maximum.    An out of phase wave of the given period therefore causes a decrease in average power of $5.3%; the same wave combined with a tenth power law velocity profile results in a decrease in average power of $8.5%, and the wave in phase with two turbine rotations gives a decrease of $3.4%. The power range in the flat case is 17.5 kW, which increases by a huge 1194% and 1317% for the OP and OP + P cases respectively. This range then almost doubles again in the 2IP case. The effect of adding the profile shifts the range of the OP model down by $30-50 kW as seen in Fig. 14, whereas the 2IP model gives an increased and relatively even spread.

Transient performance curves: thrust
The Flat thrust curves show a repeating, steady pattern for each blade, with minimal fluctuation in the total thrust value. The peak thrust for each blade occurs at $100°after the blade passes the top position, and the minimum $200°after it passes the top position. The blade values fluctuate from $90 kN to 103kN, and the repeating fluctuations that are due to the effect of the blades passing the stanchion and their position in the water column damp out the oscillations in the total thrust value, which range from 300 kN to 306 kN (Fig. 15). However, when the out of phase wave is added in Fig. 16, the repeating pattern is lost as the effect of the wave far outweighs the natural fluctuations due to blade rotation. The blade values now vary on an inconsistent basis from 80 kN to 115 kN, and the total from 260 kN to 340 kN. There are still local peaks due to blade position, but the main peaks occur at locations equivalent to the wave peaks. In Fig. 17 the profile has been added, which has very minimal effect on the total thrust values, with a slight overall decrease.
A much greater effect is seen in the 2IP model, with the re-emergence of a repeating blade and total thrust pattern as the wave is in phase with the turbine rotation. This means there is an identically repeating pattern every second rotation. Again the local blade maxima occur in the same locations, but are almost not visible due the overriding effect of the in phase wave. The ranges have widened for the blades from 69 kW to 122 kN, and for the total from 215 kN to 364 kN. The wave peak is

Transient performance curves: power
The Flat power curves again show a repeating, steady pattern for each blade, with ranges as shown in Tables 2 and 3. The application of the out of phase wave again destroys this pattern, with a further reduction in power on applying the profile as previously discussed. The difference in this case is that the addition of the 2IP wave causes the power on blades 1 and 3 to peak at the wave peak, as opposed to the thrust case where the peak was on blades 1 and 2. Another notable difference is that the power increase in terms of averages and ranges is much greater than the effect on thrust.

Transient performance curves: bending moments
The presence of the waves has a significant effect on bending moment, but the addition of the profile clearly has the strongest effect (Fig. 23). The in phase wave yields a greater increase in resultant bending moment than the out of phase case, and there is a very significant negative spike in the bending moment angle for the 2IP wave.  The presence of the waves has a significant effect on bending moment, but the addition of the profile clearly has the strongest effect. The in phase wave yields a greater increase in resultant bending moment than the out of phase case, and there is a very significant negative spike in the bending moment angle for the 2IP wave.

Frequency analysis
Frequency analysis via the Fast Fourier Transform (FFT) has been undertaken to highlight the interaction between load fluctuations due to the rotational frequency of the turbine and those resulting from the presence of the ocean wave. In each case FFTs (Figs. [24][25][26][27][28][29] were used to calculate the thrust and power spectrums for each individual blade and on the total rotor.  For each condition the rotational frequency of the turbine and the frequency of the wave loading are readily detectable. The rotational frequency of the turbine in general is accompanied by thrust peaks at integer multiples (harmonics) of the rotational frequency of the turbine, which decay rapidly until the third or fourth harmonic after which point they are no longer detectable.

No wave case
Using thrust as an example for analysis, Fig. 24 shows the thrust spectrum under the Flat condition for each individual turbine blade. The degree of load fluctuation on each blade is approximately 17%. There is a high degree of symmetry highlighted by the similarity in the observed thrust spectrum in each case. Fig. 24 also shows the resultant total thrust spectrum for the Flat case. The phase relationship between the individual blade spectrums leads to a construction and deconstruction of the harmonics shown in the total thrust spectrum. Specifically the 1st and 2nd harmonics observed in the individual blade spectrums are cancelled out due to the 120°phase relationship,  whereas the 3rd harmonic is attenuated and superimposed in the total spectrum. In the total thrust spectrum the degree of fluctuation relative to the mean total thrust is 1.68% or 5 kN.
The analysis of the power spectrum for the Flat case is shown in Fig. 25. This again shows a level of symmetry between the blades as would be expected. Each blade generates $155.5 kW of power, but with a fluctuation $28%. Due to the absence of any waves or velocity profile through the water column this fluctuation is due to shadowing effects as a result of the presence of the support structure.

Out of phase wave case
The thrust spectrum calculated with the OP wave for each individual blade is shown in Fig. 26. The frequencies relating to the rotational velocity of the turbine are again clearly observable. In the OP case the symmetry between each individual blade has been reduced due to the presence of the wave. The wave loading frequency for each blade is clearly observable. The degree of fluctuation due to the shadowing effect, characterised by the relationship between the spectrum frequency and the rotational frequency of the turbine ranges between 14% and 22% of the mean thrust on each blade. The fluctuation due to the wave loading on an individual blade basis is less than the shadowing fluctuation, ranging between 9.5% and 10.5% of the mean thrust loading per blade. Fig. 26 shows the total resultant thrust spectrum for the OP case. Again, the construction at three times the rotational frequency due to the phase relationship between the individual blade spectrums is observable. The effect of the wave has Table 2 Averages and ranges of power from the transient models.

Model designation
Average summed relative to the wave loading on each individual blade and appears as the source of highest load fluctuation. In the total thrust case the fluctuation due to the shadowing effect is 7.15%, whereas the fluctuation due to the wave loading is 9.92%. Fig. 27 shows the power spectrum for the OP case. The frequencies follow a similar trend to those of the thrust. The power generated is reduced from the Flat case from $115.5 kW to $147 kW but also have a large level of fluctuations over 47%. Overall the turbine now generates $441 kW with a fluctuation of $11.5% due to the shadowing effects. The introduction of waves provides the largest effect on the power generation at over 17%. Fig. 28 shows the thrust spectrum for each individual blade under the 2IP wave. Again the harmonics due to the shadowing effect at the blade pass frequency and its harmonics are visible. However, the spectrum is dominated by the effect of the in phase wave. In this case the shadowing effect is between 14% and 16% of the mean thrust on each individual blade. The depth of fluctuation due to the wave is significantly higher, ranging from 32% and 36% of the mean thrust on each blade.

In phase wave case
The effect of the wave has had a compounded effect on the total thrust spectrum for the 2IP wave case as shown in Fig. 28 The phase relationship due to the shadowing effect observed in the data can still be noted by the presence of the thrust peak  at three times the rotational frequency of the turbine, which $3.21% of the total thrust on the turbine. The wave effect in this case is a source of extreme load fluctuation and is 34.05% of the total thrust on the turbine rotor. Fig. 29 shows the power spectrum for each individual blade under the 2IP wave. As with the thrust data the harmonics due to the shadowing effect at the blade pass frequency and its harmonics are visible. The spectrum is again dominated by the effect of the in phase wave. The fluctuations due to the shadowing effect are however much higher for the power at     between 34% and 55% on each individual blade. The fluctuations due to the wave are significantly higher again, ranging from 53% and 69% on each blade.
The effect of the wave has also had a compounded effect on the total power spectrum for the 2IP wave case as shown in Fig. 29. Again the phase relationship due to the shadowing effect observed in the data can still be noted by the presence of the power peak at three times the rotational frequency of the turbine, which $4.8% of the total power on the turbine. The wave effect in this case is a source of extreme fluctuation and is 61.3% of the total power generated by the turbine rotor.

Total value comparisons
The effect of adding an out of phase wave of the period in this study increases the amplitude of the power fluctuations from the turbine by over 200%, whilst decreasing the average power by $5.3%. There is a similar trend for thrust, but the increases and decreases are on a much smaller scale. This is possibly due to the selection of the turbine rotational velocity to be just off the maximum power point. As the wave increases and decreases the local velocity, parts of the turbine that coincide with the increased velocity due to the wave peaks are moving through the maximum power point and exhibiting a very large increase in power; whereas the maximum thrust point is at a different location. If the turbine were operating at a different tip speed ratio it is therefore possible that the effect of a wave may be to drastically increase the thrust and reduce power.
Once the wave becomes in phase with the turbine rotation and the period increases, the effect on the power vastly increases, with the range doubling from the out of phase case. In the model run, the wave peak corresponded to a particular blade being in the top position, meaning that the combined effect of this additional increase and the longer period led to much greater power fluctuations. Again, the effect on thrust, while a large increase in range, was less dramatic. It is expected that an increase in time period would have a greater effect on the turbine, as the wave has a greater depth penetration effect and causes greater changes in subsurface velocities. It is expected that as the period increases further, more significant power fluctuations will be observed, although as discussed this will depend on the tip speed ratio of the turbine, and could potentially lad to vast thrust variations instead.
If the peak of the wave were lined up with where the power of each blade peaks; approximately 100°after it passes the top position, a significant increase in fluctuations on a single blade could be observed. This uneven distribution of power and thrust spikes and decreases could lead to non-uniform loading and fatigue, and could result in much greater bending moments. However, considering the OP case, the natural variation in power and thrust across the blades as a result of the blade positions and wave peaks coming at different times has a smoothing effect, meaning that no single point or blade is subjected to such a repeating pattern. In a real sea state, the incident waves will be composed of many simultaneous long and short period waveforms, which will prevent sharp repetitive spikes.
There is the risk though of a single dominant, long period wave in phase with the turbine rotation, and techniques to predict this and mitigate the risk, such as altering the tip rotational velocity of the turbine, load shedding or pitch alteration would be required. The total thrust and power values are shown in Figs. 30 and 31 for the same time period, and the dramatic effect the longer period and in phase wave has on the turbine is apparent. There is some smoothing effect from the out of phase nature of the OP case. On a modelling note, it is clear that the generated wave is shifted to the left with the addition of the profile, and is thought to be due to the velocity dragging effect of the profile and its different flow velocities. This must be considered when designing future models, but does not affect translation of results to a real sea state.

Bending moment comparisons
The bending moment graphs in Fig. 23 show that clearly the greatest effect on resultant bending moment occurs with the presence of the velocity profile. The profile used was fairly modest compared to some site data [29], and was approximated Fig. 30. Total thrust for the three wave models. as a smooth power law, which again is not representative of reality. This effect in a real sea state will therefore be amplified considerably, and the combination of a higher period, in phase wave with a high shear, irregular velocity profile could result in very considerable fluctuations and bending moments.

Conclusions
In conclusion, the long period and in phase wave would seem to constitute the greatest risk to turbine survivability, with modest changes in average power and thrust, but a startling increase of over 1000% in the power range. This was confirmed in the frequency analysis of the results, showing the wave effect become dominant over the blade shadowing. This effect may depend on the tip speed ratio of the turbine, as its operating location on the thrust and power curve will determine whether effects from local wave velocities will move the performance towards or away from the maximum operating point. Therefore in this case, the predominant effect was on power range, but in other cases this may be different.
A real sea state and out of phase wave will have a smoothing effect on this phenomenon, especially when considering irregular and composite waves. However, the risk of predominant long period and in phase ways may be significant, and will possibly require techniques for measurement and mitigation, such as changing the turbine tip speed ratio.
Further study will be conducted to determine if this increasing trend continues with increased period and is exacerbated by in phase waves, with models to be run incorporating composite waves, which will possibly smooth out the fluctuations. Work is required to determine the effect that these waves have on the turbine loading, and will be carried out by way of fluid-structural interaction modelling. The immediate requirement is to investigate the potentially severe effects of combining a long period wave, such as the 2IP configuration, with a high shear velocity profile.