Susan Tully

Ignazio Maria Viola
email: i.m.viola@ed.ac.uk

S Y Mpo

O Tatin

Reducing the wave induced loading of tidal turbine blades through the use of a flexible blade

Keywords:

A study into the loading e ects of wave-current interactions on a blade section of a horizontal axis tidal turbine was performed. Wave-current interactions were calculated based on 2 nd order linear wave theory and a numerical model for estimating loads on both rigid and flexible blade sections is presented, based on a quasi-steady analysis. Results from this model are compared with load measurements on two constant cross-section hydrofoils, one rigid and one flexible, in combined waves and currents in order to assess whether a flexible blade can lead to lower load fluctuations. Particle image velocimetry was used to investigate the flow field surrounding the hydrofoils throughout a wave period in order to better understand the underlying hydrodynamics. The flow experienced by each hydrofoil is found to be highly unsteady with hysteresis e ects resulting in di erent loading profiles than the quasi-steady analysis predicts. The experimental results indicate that the oscillating pressure field, associated with an oscillating free surface, significantly changes the hydrodynamic behaviour of the hydrofoils. The flexible blade was found to reduce the magnitude of load fluctuations in addition to achieving higher lift to drag ratio as compared to its rigid counterpart.

INTRODUCTION

Tidal turbines must operate within an unsteady flow environment. The combination of diurnal currents, sheared current profiles, high levels of turbulence and wave-current interaction leads to highly variable flow across devices (see figure 1), with an associated variation in the loads experienced. Whilst this has ramifications both on the ability of a turbine to generate power e ciently and on the expected lifespan of a device, a lack of understanding of unsteady turbine hydrodynamics has lead to some developers choosing to over-engineer their devices by as much as 30% [START_REF] Marsh | Wave and Tidal Power -an emerging new market for composites[END_REF]. Of particular note is the e ect of cyclical wave-induced loading on devices since their cyclical nature means that blade and rotor fatigue becomes an increasingly important issue. For example, the large fetch provided by the North Atlantic means that UK wave climates are highly energetic, leading to wave-induced velocity oscillations strong enough to propagate through the water column to typical turbine deployment depths. The European Marine Energy Centre (EMEC) reported that wave-induced velocities could be felt as deep as 20 m below the free surface [START_REF] Norris | Update on emec activities, resource description, and characterisation of wave-induced velocities in a tidal flow[END_REF]. Further more, as a turbine blade rotates, and passes through the depthdependent wave-induced flow field, the blade will experience an additional cyclical loading e ect with a period of the blade rotation. These wave-induced loadings can decrease the life of a turbine blade to 30 months, when the expected lifetime is of the order of 30 years.

Recent reviews [START_REF] Lachenal | Review of morphing concepts and materials for wind turbine blade applications[END_REF] and [START_REF] Xiao | A review on flow energy harvesters based on flapping foils[END_REF] indicate that flexible foils have shown promise within the wind energy industry in reducing the increasingly large load fluctuations associated with the ever increasing turbine diameters. Given the high fluid density of water, and the likelihood of tidal turbines experiencing cyclical flow, flexible blades require further investigation in order to decide on their capacity to reduce loads within a tidal environment and thus increase the life time of next generation tidal turbines.

The loading of horizontal axis tidal turbines has been an active research question for a number of years though the vast majority of work has gone into quantifying and understanding steady-state loading. Of the few studies which have looked into unsteady loading behaviour, investigations have focused on torque and thrust measurements of turbines as a whole. Two such investigations by Barltrop et al. [START_REF] Barltrop | Investigation into wave-current interactions in marine current turbines[END_REF] and Galloway et al. [START_REF] Galloway | Studies of a scale tidal turbine in close proximity to waves[END_REF] both found significant variations in instantaneous thrust and torque magnitudes through a wave period, greater than 35% of the mean rotor load; though smaller di erence to average power coe cients as compared to a current only case. Luznik et al. [START_REF] Luznik | The e ect of surface waves on the performance characteristics of a model tidal turbine[END_REF], in corroborating these findings, also noted that the presence of waves raised the lower threshold of the operable tip speed ratio (TSR) range due to the variations in torque produced. Since these studies were investigating rotor-scale loading it is di cult to draw any conclusions on the flow and loading experienced by individual blades.

Milne et al. [START_REF] Milne | Blade loads on tidal turbines in planar oscillatory flow[END_REF] recognised that individual blades were likely to undergo more severe loading cycles than the rotor as a whole, where the torque and thrust behaviour was in e ect an averaging of the forces incident on the three blades. They investigated blade scale loads via instrumented blades on a scaled turbine rotating in planar oscillating flow. While Onset flow experienced by a tidal turbine blade section they found increased loads associated with separated flow, and evidence of approximately 25% greater loads during dynamic stall for certain wave conditions, their use of an unsteady flow which was uniform over the rotor plane meant that their investigation lacked any conclusion on the e ects of each blade's rotating through a depth-dependant wave field.

In 2014 Galloway et al. [START_REF] Galloway | Quantifying wave and yaw e ects on a scale tidal stream turbine[END_REF] instrumented one of the blades of their scaled turbine in order to investigate blade-scale dynamic loading induced by waves and also investigated a depth-dependent wave field. They found that wave frequency, not height, was the driving force behind wave-induced loadings. In addition, they found that dynamic stall had little e ect for the wave conditions they tested, both experimentally and numerically where they used the Boeing-Vertol model developed by Tarzanin [START_REF] Tarzanin | Prediction of control loads due to blade stall[END_REF]. This result is in disagreement with the results of Milne, despite both testing a similar scale device, at similar frequencies and oscillation amplitudes. What di ered between the two analyses is that while Milne tested under purely planar oscillatory flow by towing their turbine at a sinusoidally oscillating velocity in stall water, Galloway towed their device through a wave field. This means that Galloway's tests included the e ects of both a dynamically varying flow field and additionally the oscillations of static pressure from the free surface displacement. Whether this di ering dynamic pressure field was the cause for the di erence in loadings is still an open research question.

These investigations into the unsteady loading of tidal turbines focused on investigating the e ects of the largest and longest waves and relied on load measurements alone, with no corresponding flow visualisation techniques used to further elucidate the hydrodynamic cause of the loading. Of course, oscillatory flow over airfoils has been studied extensively, especially for propulsion with most papers taking a Particle Image Velocimetry (PIV) based approach. The majority of studies have either actively pitched or plunged a foil in a steady current or towed it in waves only, as is also the case in each paper listed above. Neither approach is a true representation of the flow field as would be experienced by a tidal turbine blade in real ocean conditions, where higher order wave-current interaction e ects alter the mean vertical velocity profile [START_REF] Olabarrieta | E ects of wave-current interaction on the current profile[END_REF] and the oscillation in free surface height is out of phase with the change in dynamic pressure.

In this paper the hydrodynamics of the wave-blade interaction is studied, for both a rigid and flexible blade section, using a combined wave and current tank which naturally incorporates the higher order e ects expected for real ocean conditions. The aims of this paper are two fold. First is to investigate whether a flexible tidal turbine blade could achieve a more favourable loading profile, as compared to the more standard rigid blade. Second, to answer the question of whether or not the oscillation in static pressure caused by ocean waves has any e ect on the hydrodynamic behaviour of submerged hydrofoils.

To address these aims, this paper investigates the unsteady flow perturbations experienced by individual tidal turbine blade sections when exposed to wave-induced oscillatory flow, utilising both numerically and experimentally derived loadings alongside PIV flow visualisation. In addition a novel flexible blade design is presented which has been designed to reduce the amplitude of load fluctuations experienced through passively flexing the trailing edge section.

METHOD 1.1 Construction of the flow experienced by a blade section

The flow scenario is presented in figure 1, including the coordinate systems used. Whilst the problem is highly three dimensional, insights can be obtained considering the flow incident on a blade section from a rotating reference frame moving with the blade. The parameters of interest are then the relative velocity incident on the hydrofoil U rel , the angle of attack ↵ between the blade chord c and U r el , and the reduced frequency k f of the flow oscillation. U rel varies with time t and it is a function of the current U 0 , the wave orbital velocity components u w and w w in the x and z directions, respectively, the angular velocity of the blade ! and the radial distance r of the blade section from the turbine axis. The reduced frequency k f has been defined as k f = 2⇡ f c/U 0 , where f is the wave frequency. Thus, for k f = ⇡/2 and ⇡, a particle of fluid moving at U 0 would convect x/c ⇡ 4 and x/c ⇡ 2 in one wave period, respectively. Wave velocities, u w and w w are calculated using 2 nd order linear wave theory [START_REF] Dean | Water wave mechanics for engineers and scientists[END_REF], equations ( 1) and ( 2), ensuring that the wavenumber k = 2⇡/ is calculated from the doppler shifted dispersion relation = U 0 k + p gk tanh(kh); where is the wavelength, = 2⇡ f as measured by a static wave buoy, g is the gravity acceleration and h the water depth. If dealing with waves propagating counter to the current direction then u w changes sign, therefore the ± in equation ( 1) stands for opposing and following waves respectively. Only waves which satisfy the Stokes limit of maximum steepness of H s / = 1/7 should be included, while any with k  0 should be assumed blocked. Here, H s refers to the peak-to-trough wave height.

u w = H s 2 gk cosh k (h + z(t)) cosh kh cos (k x ± t) + 3 16
H 2 s k cosh 2k (h + z(t)) sinh 4 (kh) cos 2(k x ± t) (1)
w w = H s 2 gk sinh k (h + z(t)) cosh kh sin (k x ± t) + 3 16
H 2 s k sinh 2k (h + z(t)) sinh 4 (kh) sin 2(k x ± t) (2)
When investigating the flow experienced by a blade cross section as it rotates, x remains constant while z(t) is time dependant to account for the blade rotation through the depth dependant wave field z(t) = z hub +r sin !t. The decay of wave orbital dimensions H (z) and W (Z ) below the free surface are given in equations ( 3) and ( 4) respectively.

H (z) = H s 2 sinh k (h + z) sinh kh (3)
W (z) = H s 2 cosh k (h + z) sinh kh (4)
For stationary, cylindrical coordinates the parameters of interest on the rigid blade section at distance r out from the hub are defined in equations (5 -8).

u axial = U 0 + u w (5) u ✓ = !r + w w cos !t (6)
U rel = q u axial 2 + u ✓ 2 (7) ↵ = 90 arctan u ✓ u axial ( 8
)
where is the static blade pitch, which is the angle between the rotor plane and the blade chord.

The antisymmetric NACA 4415 airfoil profile was chosen for this study given its use in the PerAWAT project [START_REF] Buvat | Identification of test requirements and physical model design[END_REF] for scale model tests designed to replicate the behaviour of Alstom Power's (then Tidal Generation Ltd) turbine. Extensive data on the performance of the NACA 4415, under steady flow conditions and at similar Reynolds number (Re = 8 ⇥ 10 4 ), can be found in [START_REF] Saliveros | The aerodynamic performance of the NACA-4415 aerofoil section at low Reynolds numbers[END_REF].

Quasi-static numerical model

The loads on the rigid hydrofoil are calculated based on the steady-state normal force coe cient, as defined in equation ( 9), for the NACA 4415 profile used in this study, taken from [START_REF] Saliveros | The aerodynamic performance of the NACA-4415 aerofoil section at low Reynolds numbers[END_REF].

C n = F n 1 2 ⇢U 2 rel cb (9)
where C n is set to remain constant, at the minimum and maximum value respectively, to replicate the early stall behaviour. The loads on the rigid hydrofoil are then scaled according to ↵ as calculated in equation ( 8). The load on the flexible blade is computed taking account of the passive deflection of the trailing edge which changes the direction of the chord, which is defined as the line through the leading and trailing edge of the blade section; therefore a displacement of the trailing edge leads to a variation of and thus of ↵. The flexible section of the foil is represented as a cantilever beam, fixed at one end, with a uniformly varying load to mimic the pressure di erence between the pressure and suction sides of the hydrofoil along the flexible length of the chord. A schematic is given in figure 2. The maximum di erence between the pressure coe cient on the suction and pressure sides of the hydrofoil C max P,beam , occurs at the fixed end of the cantilever beam and is related to the maximum coe cient of pressure di erence at the leading edge of the hydrofoil C max P by equation [START_REF] Tarzanin | Prediction of control loads due to blade stall[END_REF], where l = 2 3 c is the length of the flexible part of the foil.

F
C max P,beam = l c C max P ( 10
)
The pressure di erence decays linearly along l to zero at the trailing edge. This triangular C P distribution can be seen in [START_REF] Saliveros | The aerodynamic performance of the NACA-4415 aerofoil section at low Reynolds numbers[END_REF] for the NACA 4415 over a broad range of attached flow conditions with C max

Experimental test

A rigid and a flexible hydrofoil were manufactured on an Objet500 Connex multi-material 3D printer (16 µm layer accuracy) using PolyJet TM VeroWhitePlus RGD835 and Tan-goBlackPlus FLX980 materials. TangoBlackPlus is a flexible, rubber-like material which can be mixed in di erent ratios with VeroWhitePlus to produce di ering grades of flexibility. The flexible blade is developed around a flexible rubber hinge with l = 2 3 c as in the numerical section, allowing the trailing edge to flex in response to pressure changes in the surrounding flow. The first third of the hydrofoil is rigid while the latter two thirds are free to flex by incorporating a thin split on the pressure side such that the trailing edge section can slide back and fourth along the underside of the rigid leading edge section, without major changes to the profile. A schematic of the flexible blade is shown in figure 3 [START_REF] Robinson | The e ect of inlet design on the flow within a combined waves and current flumes, test tank and basins[END_REF]. For this study, only opposing waves are tested experimentally in order to be able to include a series of turbulence reducing meshes near the inlet, reducing the turbulence intensity to approximately 3% at U 0 = 0.5 m/s. Waves are tested at a frequency f = 1 Hz, and wave height H s = 0.09 m leading to a peak-to-trough ↵ oscillation of 24 . The test-section is 5 m downstream of the inlet (4 m downstream of the last mesh screen). The hydrofoil is mounted horizontally at mid-depth, such that the quarter chord position is at z f oil = 0.225 m below the free surface and is sandwiched between two vertical clear acrylic splitter plates, placed to negate tip loss e ects. The foil is attached to the splitter plates via an 8 mm stainless steel rod which threads through both splitter plates and the full span of the foil. There is a 1 mm gap between the hydrofoil and each splitter plate, achieved via placing a washer on the rod between the splitter plate and foil which allows the flexible blade to deflect. There is also a 95 mm gap between the plates and the flume side-walls in order to remain clear of the boundary layer on these walls. Hydrofoils are mounted at = 10 static pitch. The PIV system used consists of a pulsed dual-cavity Nd:YAG laser with 200 mJ at 532 nm and a maximum pulse rate of 15 Hz; with a CCD IMPREX 50MP camera with 2448 x 2050 pixel resolution, mounted outside the flume with the image plane parallel to the free stream flow. The camera has a max frame rate of 50 frames/s but the bu er download speed limits this to 15 frames/s, leading to a maximum PIV image frequency of 7.5 Hz. A Nikkor f /2.8, 50 mm lens is used to achieve a su ciently shallow depth of field, with a 532 nm monochromatic filter installed in order to minimise the noise from ambient light. As a further noise reduction measure the rigid hydrofoil is painted matt black in order to minimise reflections, while the flexible blade's matt black rubber coating means no painting is required. The laser is mounted vertically above the flume, requiring direction over three mirrors and through a beam expander in order to align the laser sheet parallel with the flow direction. The laser sheet is approximately 4 mm thick, centred at the quarter chord depth, with the sheet axis aligned at -5 below the horizontal in order to minimise free surface reflections from the waves. The light sheet is positioned at one third span from the camera facing side of the hydrofoil. Conduct-O-Fil silvered spheres are used to seed the flow with an average diameter of 14 µm and average density of 1.7(g/cc).

Load cell and PIV data are recorded for 20s intervals, capturing 20 wave periods for the 1 Hz wave frequency used. The current is initialised 20 s prior to wave generation and then 6 wave peaks are allowed to pass the hydrofoil before recording begins and the flume is allowed to settle before the next measurement is taken. As the loads measured include a contribution from the splitter plates, loads were measured for the splitter plates alone and the mean, over one wave period, subtracted from the results. The loads are measured at a sampling frequency of 150 Hz in order to allow for direct comparison with the PIV images (7.5 Hz frequency). The average intensity is subtracted from each PIV image set prior post-processing. Vector fields are then calculated using 3pass adaptive correlation with square, 32 pixel, interrogation windows at 50% overlap and any ensemble averages are generated from a minimum of 9 wave periods worth of data.

The conditions tested experimentally are listed in table 1, where ⌫ stands for the kinematic viscosity of the water and the other symbols have been introduced above. H s in table 1 is that set in the numerical model to achieve an ↵ oscillation comparable with the experimental results while U 0 is the mean flow velocity in the x direction measured at z f oil for combined waves and currents. Elastic parameter, rigid F =

E F d 3 F /12 0.5⇢U 2 0 l 3 F 0.13 Elastic parameter, flexible R = E R d 3 R /12 0.5⇢U 2 0 l 3 R 2.3 ⇥ 10 4

RESULTS AND DISCUSSION

Loads on a blade section

Figures 567show the results from the numerical model and experimental test conditions for both the rigid and flexible blades. The markers denote the phase averaged experimental results whilst the predictions from the numerical model are denoted by lines. The results are plotted versus time (t) normalised with the wave period, t w = t f . In addition, the error bars on figures 5 and 6 give the spread of ↵ and U rel respectively at a confidence level of 95%. The turbulence in the flume, combined with the large waves required to generate suitable ↵ oscillations, led to some discrepancy in flow conditions one wave period to the next. This is reflected in the size of the error bars drawn. The computed ↵ oscillation shown in figure 5 for the rigid hydrofoil has been matched to the magnitude of the ↵ oscillation measured in the experiments while the flexible hydrofoil experiences an entirely di erent ↵ oscillation as the deflection of the trailing edge changes ↵. Both hydrofoils were mounted at a static pitch of 10 and the rigid hydrofoil oscillates around this value. However the flexible hydrofoil oscillates around a smaller mean ↵ as the trailing edge initially deflects upwards in response to the force from the current. The rigid blade undergoes a peak-to-peak ↵ oscillation of 25 while the flexible blade only experiences half of this, with a peak-to-trough ↵ oscillation of 12 .

The comparison between measured U rel and that predicted be the numerical model is given in figure 6. The numerical model over predicts the amplitude of the velocity oscillation. Bearing in mind that the numerical wave height was set to match the ↵ oscillation measured, this over-prediction is likely due to the wave velocities being calculated based on linear wave theory where waves are purely sinusoidal. In reality higher H s waves are steeper in the crests and shallower in the troughs. Remembering that the experimental case involved the current and wave field travelling in opposite directions such that the maximum U r el corresponds to a wave trough and the minimum to a wave peak, thus the free surface elevation is minimum at T w = 0 and maximum at T w = 0.5, while U rel is maximum at T w = 0 and minimum at T w = 0.5.

The variation in C n , as defined in equation ( 9), is given in figure 7. The numerical code, based on a quasi-steady analysis predicts that the rigid blade will be in stall for the first half of the wave period, while the experimental results show that C n continues to rise, despite ↵ reaching upwards of 15 , well past the steady state stall angle. This agrees with the current understanding of unsteady loading and dynamic stall in that the normal operational window on ↵ is extended and that separation occurs at higher angles of attack than would be expected from steady state results. In addition for the last quarter of the wave period, the numerical model vastly over predicts the loading on the blade; this is due to flow separation on the upper side of the hydrofoil as will be shown in the following section 2.2.

The variation of the normal force coe cient with the angle of attack is shown in figure 8. The computed trends for the rigid and flexible blades are the same but extend for di erent intervals of ↵, roughly 1 < ↵ R < 22 and 3 < ↵ F < 7 .

C n increases linearly up to a plateau for ↵ > 10 . The flexible blade agrees quite well with the numerical prediction, as also shown in figure 7, and does not present any hysteresis loop. On the contrary, the rigid blade shows a clear hysteresis loop which is traveled anti-clockwise. At T w = 0, under the wave trough, ↵ ⇡ 10 and C n ⇡ 0.5; with the increase of the free surface elevation both ↵ and C n increase and C n overshoots the maximum predicted for quasi-steady conditions. In this initial part of the wave period, C n is much lower than predicted because of the large trailing edge vortex that will be discussed in the next section. After having reached the maximum C n , ↵ begins to decrease and C n decreases as well, but at a lower rate than it increased because trailing edge separation is no longer present on the upper side of the hydrofoil.

While the mean C n for the rigid hydrofoil is higher than that of the flexible hydrofoil, the magnitude of the load fluctuations is also higher. In this respect the flexible hydrofoil has achieved its performance aim, that of reducing the load fluctuations on a tidal turbine blade. The fact that the mean C n for the flexible blade is lower than that of the rigid blade does not imply that there would be a performance loss if a tidal turbine were to use a flexible blade. In fact each blade section should operate at the maximum lift/drag ratio condition in order to maximise the energy extracted, and thus a higher C n is not necessarily desirable. The time average over a wave period of the lift over drag ratio for the rigid blade is L/D R = 8 while for the flexible blade it is L/D F = 10.7. Therefore, while the mean C n is lower for the flexible blade, its lift over drag ratio is higher and so could lead to a higher power output under the right design and operating conditions. Another point to add is that losses due to drag are what limit the maximum possible design TSR [START_REF] Burton | Wind Energy Handbook (Google eBook)[END_REF]. The flexible blade will allow for a higher ceiling on the maximum possible TSR, since its drag losses are lower, allowing for more power output.

Flow fields around a blade section

The time averaged velocity fields around both the rigid and flexible hydrofoils are given in figure 9. The greyed out region on the bottom image covers the region over which the trailing edge deflected throughout the recording. The large degree of trailing edge separation experienced by the rigid hydrofoil can be seen in the velocity deficit behind the hydrofoil which is not present for the flexible blade. This is because the rigid hydrofoil experiences a much larger ↵ oscillation and goes some way towards explaining the lower load fluctuations experienced by the flexible hydrofoil.

To investigate further, figure 10 gives the velocity field around each foil at key points within a single, ensemble averaged, wave period t w . Images are ensemble averages of all flow fields captured in 5 bins at key points through the wave period as indicated in the insert. The angle of attack and velocity is that measured at the leading edge position of the foil. What is not shown in the PIV images is what is happening to the pressure field, along the chord of the hydrofoil, caused by the free surface displacement. The static pressure varies according to the free surface displacement such that the maximum corresponds to a wave crest and the minimum to a wave trough. Therefore the maximum wave-induced pressure gradient (@p w /@ x) is at t w = 0.25 and minimum at t w = 0.75. From the momentum equation it can be observed that @p w /@ x = ⇢Du x /Dt. Therefore, from the derivative of U rel /U 0 in figure 6, it can be estimated that @p w /@ x ⇡ 750 Pam 1 , which is about one q/c, where q = 1/2⇢U 0 2 is the reference dynamic pressure of the tidal stream. Recalling the triangular distribution of C P along the foil chord, it can be estimated that the quasi-static pressure gradient on the upper side of the foil @p 0 /@ x is about 5q/c when ↵ is maximum at t w = 0.25 and negligible when ↵ is minimum at t w = 0.75. Therefore the boundary layer on the upper side of the foil experiences a maximum (adverse) pressure gradient, which promotes separation, of about 6q/c at t w = 0.25 and a minimum (favourable) pressure gradient of about q/c at t w = 0.75. The numerical model accounts only for @p 0 /@ x and neglects @p w /@ x.

The adverse pressure gradient is at its maximum between images B and C and is the cause of the small separation bubble that can be seen in figure C for the rigid blade, roughly half way along the chord. As time progresses, this separation convects downstream, indicated by the point of lowest velocity on images D, E and F. Interestingly, for the rigid blade at image F, the angle of attack is approximately 0 yet there is clearly already a large degree of trailing edge separation, with the velocity at the trailing edge almost zero. This is thought to be due to the combination of the separation that started in image C convecting o the end of the foil, and the trailing edge vortex which begins to develop due to the increasing ↵. The result is that, by the time the flow returns to the beginning of a wave period, there is a large trailing edge vortex generated with counterclockwise circulation which decreases the total foil circulation. This explains the discrepancy between the loading predicted for the final stage of a wave period by the numerical model for the rigid blade and the experimental results, see figure 7.

Comparing with the flexible blade, the same separation, and convection downstream, can be seen in image C, but to a smaller degree and the images are very comparable between points C and E. The di erence in flow behaviour lies in the fact that, given that the trailing edge of the flexible foil is free to deflect, it has deflected upwards in response to the high velocity at point F and the volume of low-speed flow resulting from the separation is significantly lower than for the rigid blade. As a result the experimental loading measured matches well with the attached flow analysis of the numerical model. Further work is required in order to quantify the e ect of the blade rotation, which has been ignored in this work. Furthermore the rough seabed and flood and ebb of the tide results in a constantly varying horizontal shear layer which should be accounted for when carrying out a full analysis of the unsteady loadings. In a similar vein, tower shadow and flow turbulence should also be included. The approach used in this study also relied on superimposing a wave field on a current, with the assumption being that the mean flow field is that of the current alone. Whilst this has given reasonable results here, recent work [START_REF] Olabarrieta | E ects of wave-current interaction on the current profile[END_REF] suggests that surface waves actually alter the mean flow profile which should also be accounted for too.

CONCLUSION

A numerical model for estimating loads on both rigid and flexible blades was presented and compared with load measurements on two hydrofoils, one rigid one flexible, in combined waves and currents with opposite directions. Particle image velocimetry was used to investigate the flow field surrounding the hydrofoils, and the lift and drag on the hydrofoils were measured. The following conclusions are drawn.

• The periodic wave-induced flow fluctuations experienced by a tidal turbine blade, or a generic submerged hydrofoil, can lead to highly unsteady e ects, including chord-normal forces more than 20% higher than in the current only case.

• The quasi-static chordwise pressure gradient on the upper side of the hydrofoil varies due to the waveinduced angle of attack oscillations. The background horizontal pressure field also varies because of the periodic free surface oscillations. These two pressure gradients are in phase, leading to a periodic very large combined pressure gradient which triggers boundary layer separation on the upper side of the hydrofoil at every wave period.

• A passively flexible blade is capable of reducing the amplitude of load fluctuations through decreasing the angle of attack variations, thus reducing separation and drag at high flow angles of attack. This results in a higher lift over drag ratio and so shows promise for improving the design of tidal turbine blades.

Figure 1 .

 1 Figure 1. Onset flow experienced by a tidal turbine blade section

  n is the normal force, ⇢ the density of the fluid, c the chord length and b the span. Under fully steady conditions, with U rel = U 0 , the normal force coe cient C n increases linearly within the attached flow region, from C n = 0.7 at ↵ = 10 to C n = 1.4 at ↵ = 10 . For ↵ < 10 and ↵ > 10 ,

Figure 2 .

 2 Figure 2. Cantilever beam representation of the flexible blade

P = 3

 3 at ↵ = 10 . The magnitude of the trailing edge deflection Y can then be computed as in equation (11); with inertia I = d 3 /12, beam thickness d and Young's Modulus E. The Young's Modulus of the rigid and flexible blades are 3000 MPa and 0.97 MPa respectively. trailing edge deflection is then used to calculate ↵ for the flexible blade. Since this change in ↵ would lead to a change in hydrodynamic load, an iterative fluid-structure coupling is implemented in order to calculate the load at every time-step. Treating each time step individually and keeping the flow conditions constant, the numerical model loops through calculating the hydrodynamic load, and corresponding deflection, until both are in equilibrium for that time step.

  . Both hydrofoils have c = 0.15 m and b = 0.3 m.

Figure 3 .

 3 Figure 3. Schematic diagram of the cross-section of the flexible blade

Figure 4 .

 4 Figure 4. Schematic of experimental setup: (a) wave absorbing paddle; (b) turbulence reducing meshes x4; (c) hydrofoil; (d) load cells; (e) PIV laser light sheet and image capture window; (f) PIV camera; (g) laser beam expander; (h) laser; (i) wave making paddle.

Figure 5 .

 5 Figure 5. Angle of attack experienced by the hydrofoil during a wave period as measured in front of the hydrofoil (exp rigid); and predicted numerically for a rigid (num rigid) and a flexible foil (num flex).

Figure 6 .

 6 Figure 6. Measured (exp) and predicted (num) flow speed experienced by the hydrofoil during a wave period.

Figure 7 .

 7 Figure 7. Normal force coe cient during a wave period.

Figure 8 .

 8 Figure 8. Normal force coe cient versus angle of attack.

Figure 9 .

 9 Figure 9. Time averaged velocity fields around the rigid and flexible.

Figure 10 .

 10 Figure 10. Velocity field around the rigid and flexible hydrofoil at 6 instants in a wave period. Inserts at the top left of the figure show the corresponding instantaneous angle of attack and flow speed (arrows show the intervals over which the flow fields have been averaged).

ACKNOWLEDGMENTS

The authors would like to acknowledge Dr Ton Van den Bremer for insightful discussions related to this article. This work was supported by the EPSRC, grant number EP/M02038X/1.


Ph Nachtergaele

Daniel J Rixen

A M Steenhoek

Efficient weakly coupled projection basis for the reduction of thermo-mechanical models

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Reducing complex models of multi-physical models is an important step in the design and optimization of high-tech systems. For instance, one is interested in representing the dynamical behaviour of components with low order models in order to integrate these models in a functional analysis of the whole system or to apply optimization strategies. The principle of model reduction has been established long time ago [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF][START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] and substructuring is widely used in structural mechanics. However, novel applications require the reduction methods to be generalized and enhanced in order to efficiently address multidisciplinary problems. In this paper we will discuss an efficient way to build a reduction basis for thermo-mechanical systems. First we shortly summarize the concepts of model reduction and component mode synthesis. Then we outline the thermo-mechanical equations and discuss a special enrichment of the internal modes in order to efficiently represent the dynamic response of a structure to a thermal load. Finally we investigate the efficiency of the proposed basis on simple model of heat actuator.

Concept of model reduction

Let us consider a finite element model characterized by the second order differential matrix equation in which q represents the vector of degrees of freedom (DOF), K , D and M respectively the stiffness, damping and inertia matrices, and f the vector of forces M q + Dq + Kq = f .

(1)

The driving idea of reduction methods is to make an approximation on the model DOF and express them in terms of a linear combination of representative modes. It is clear that in order to obtain a reduced model, the number of modes must be much smaller than the number of DOF in the model. The previous system is then rewritten in terms of ''generalized'' DOF that physically correspond to the amplitude of the considered modes. Mathematically this consists in building an approximation of the DOF vector according to a ''trial'' subspace defined by the modes, assumed to be linearly independent. q = Rq.

(2)

In this approximation, R is a matrix which columns contain the vectors spanning the approximation space. Replacing this expression in the full set of equations leads to the following equation

MR q + DR q + KRq = f + r. (3)
The residual load appears because, due to the approximation on the DOFs, the equation cannot be satisfied exactly in general (unless the approximation space represented by R spans the exact solution, which is never the case in practice). It can be noted that once an approximate solution is computed (as described below) the equation above can always be used to evaluate the associated error in the exact problem. Clearly, there are more equations to satisfy than degrees of freedom in the reduced space. Obviously the idea is to find the set of reduced DOF that minimizes the residue. In the Petrov-Galerkin method one wants to ensure that the residue is null when projected in a given subspace of the problem. Hence in order to get an approximation of the problem, the equation is also projected onto a ''test'' subspace

RT MR q + DR q + KRq = RT {f + r} (4)
and the equation governing the behaviour of the reduced model is obtained by requesting the projection of the residue to be zero

M q + Dq + K q = f (5) M = RT MR D = RT DR K = RT KR f = RT f . (6)
The choice of representative trial and test bases is the key of the method and is directly related to the quality of the approximation. Indeed, if the trial base defines a subspace that contains the exact solution of the problem, then the residual term is null and the reduced model is exact. On the other hand one needs to use a test bases such that the error remaining after solving the reduced problem (the error being in the space orthogonal to R) does not significantly alter the solution.

Often (such as in symmetric problems) the test and trial bases are taken as identical. This case corresponds to a Galerkin technique (equivalent to a Rayleigh-Ritz method for dynamic problems) and can be interpreted as a virtual work principle (see [START_REF] Géradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]).

Static and dynamic condensation

The well-known Guyan-Irons static condensation algorithm [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF], can be seen as a projection method using static deformation shapes as trial and test bases. Nevertheless, in the case of dynamic analysis, the dynamic behaviour of the internal domain of the model, that means the condensed DOF, is not taken into account. The improvement proposed in [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] consists in adding to the trial and test bases a given number of modes directly related to the dynamic behaviour of the internal part of the model. Craig proposed to consider the internal eigen modes φ of the model, computed by constraining the retained DOF to zero. These modes indeed represent the dynamic neglected in the quasi-static modes. In a general way, the modes added to the trial and test bases, respectively Φ and Φ, can be different. Nevertheless, when using internal eigen modes, the common practice consists in considering the same modes in the two bases. In the Craig-Bampton method the following trial and test bases are considered

q = q r q c = R q r γ = Rq (7) R = I 0 -K -1 cc K cr Φ RT = I -K rc K -1 cc 0 ΦT (8)
where we denote with a subscript c the DOF that are condensed (i.e. the ones internal to a substructure) and by r the remaining ones (the DOF on an interface).

Reduction of thermo-mechanical models

In the case of thermo-mechanical problems, the DOF representing the model contains structural and thermal DOF, respectively noted ü and θ , and reordering the equations according to the physical nature of the DOF leads to the classical formulation of thermo-mechanical problems

M uu 0 0 0 ü θ + D uu 0 D θu D θθ u θ + K uu K uθ 0 K θθ u θ = f u f θ . (9)
In the previous expression, M uu , D uu and K uu correspond to the structural inertia, damping and stiffness matrices, that are symmetric. The D θθ and K θθ matrices are symmetric, and, in the context of second order equations in time, can be seen as damping and stiffness terms. Physically however they describe respectively the heat capacitance and conduction effects. The two remaining matrices D θu and K uθ characterize the thermo-mechanical coupling. Due to these terms, the complete system is unsymmetric. Partitioning again the DOF vector in terms of the retained part q r and condensed part q c , each of them containing both structural and thermal DOF, the equations of the complete system can thus be reordered and rewritten

   M uu rr 0 0 0 0 0 0 0 0 0 M uu cc 0 0 0 0 0        qu r qθ r qu c qθ c     +     D uu rr 0 D uu rc 0 D θu rr D θθ rr D θu rc D θ θ rc D uu cr 0 D uu cc 0 D θu cr D θθ cr D θu cc D θ θ cc         qu r qθ r qu c qθ c     +     K uu rr K uθ rr K uu rc K uθ rc 0 K θ θ rr 0 K θ θ rc K uu cr K uθ cr K uu cc K uθ cc 0 K θ θ cr 0 K θ θ cc         q u r q θ r q u c q θ c     =    f u r f θ r f u c f θ c    . (10)

Static projection modes

We want to stress that, unlike for purely structural models, the static modes in the trial and test matrices are not identical when handling thermo-mechanical models. Indeed, in (8) the block -K -1 cc K cr is not simply the transpose of K rc K -1 cc since K cc is not symmetric and since K rc is not the transpose of K cr . Also a particular attention must be paid to the structure of the matrices in order to implement an efficient reduction algorithm. Indeed, it is preferable to factorize both K uu cc and K θ θ cc rather than factorize the entire K cc matrix. It must also be noted that despite the symmetry of the mass matrix in the complete model, the reduced mass matrix of a thermo-mechanical system is non-symmetric. This comes from the fact that the test and trial spaces are not identical.

Internal projection modes

The choice for the internal modes could be based on the state-space form of the internal problem, similarly to what was proposed for highly damped structures [START_REF] Craig | Generalized substructure coupling procedure for damping systems[END_REF]. However, our main concern was to deal with systems composed of a structural and a thermal part that have well separate time constants, namely the period of the lowest pure structural frequencies is much smaller than the time constant of the thermal diffusion in the system. In this case, a quasi-static behaviour of the structure can be assumed as a good approximation. Therefore, we conjectured that the dynamic coupling between thermal and structural problem can be neglected in the process of internal mode computation. This last one is thus performed by solving two decoupled eigen problems. The structural symmetric eigen problem characterizes the dynamic response of the system to structural time dependent loads.

K uu cc Φ u u = M uu cc Φ u u Ω u . (11)
The ''thermal-driven'' internal modes are computed in two steps. First, the decoupled thermal eigen problem, also symmetric, that is representative of the thermal response to thermal loading is solved

K θθ cc Φ θ θ = D θθ cc Φ θ θ Ω θ . (12)
By doing this, the coupling between the two physics is not considered, leading to poor quality approximations. Therefore, in accordance with the assumption of quasi-static behaviour of the structure with respect of the thermal excitations, a correction is proposed. This new approach, which is innovative regarding the current practice, consists in adding a coupling term Φ u θ to the thermal internal modes that corresponds to the static structural response to the thermal eigen modes.

Φ u θ = -K uu cc -1 K uθ cc Φ θ θ . (13)
According to the quasi-static approximation, the trial and test basis are defined by the following projection operators in which the last column of the matrices will be called the ''thermal-driven modes''

R =     I 0 0 I 0 0 0 0 - K uu cc K uθ cc 0 K θθ cc -1 K uu cr K uθ cr 0 K θθ cr Φ u u 0 Φ u θ Φ θ θ     (14) RT =     I 0 0 I - K uu rc K uθ rc 0 K θθ rc K uu cc K uθ cc 0 K θθ cc -1 0 0 0 0 Φ u u T 0 Φ u θ T Φ θ θ T     .
(15)

Applications

The efficiency of this approach has been evaluated using a simplified model of heat actuator, illustrated on Fig. 1. The model is made of thermo-mechanical volume elements and a symmetry condition has been applied on the lower face so that the motion of the actuator is in the plane of the structure. The two anchors are clamped and an iso-temperature constraint is applied on each of them.

The whole structure is thermally loaded using volume heat source to simulate a Joule heating effect. To simulate the physics, a different value of the heat generation rate is considered in the two arms of the structure. In terms of model reduction, two mechanical DOF and two thermal DOF are retained. The mechanical retained DOF correspond to the in-plane displacements of the actuation point located at the corner of the structure tip, in the symmetry plane. The retained thermal DOF, named T 1 and T 2, are the anchor's temperatures. Note that the thermal DOF retained for the reduction will no longer be unknowns when the reduced model is used for the analysis since the anchor temperatures are then imposed.

Projection modes

The projection modes corresponding to our reduction approach are presented on Fig. 2. The first two modes correspond to the retained mechanical DOF modes (''static modes''). As these modes represent a static mechanical behaviour, they do not involve any thermal component.

The internal thermal-driven projection modes correspond to the eigen modes of the thermal condensed DOF to which the corresponding quasi-static mechanical response is added. The reduced model involves ten internal thermal-driven modes even if only the four first modes are illustrated.

Transient response

In order to evaluate the efficiency of the proposed reduction basis, transient analyses has been performed for different configuration of the model. A geometrical scaling of the model and two values of the heat capacitance of the material have been considered. Table 1 summarizes the mechanical and thermal time constants of the model in the different configurations of analyses.

For each studied configuration, the lateral displacement of the actuation point is compared for the initial 3D model, the reduced model generated using thermal-driven internal modes and a reduced model generated using pure thermal internal modes (no structural modes are used here). The simulation of the models in meter demonstrates that, even if the ratio of time constants is huge, the thermal-driven approach provides significantly more accurate results only when a low value of heat capacitance is considered, as illustrated on Fig. 3.

The simulation of the models in micron for real value of heat capacitance, presented on Fig. 4, show that the two reduction methods lead to accurate results, even if the ratio of time constants is small. However, in the same geometrical configuration, if the heat capacitance is significantly reduced the thermal-driven approach provides again more accurate results, as illustrated on Fig. 5. Firstly, these results demonstrate that, even if the ratio between mechanical and thermal time constants correspond to a relevant property in terms of thermoelastic damping [START_REF] Lifshitz | Thermoelastic damping in micro-and nano-mechanical systems[END_REF], this value is not a significant indicator in terms of model reduction strategy. Secondly, for usual values of heat capacitance, the two reduction methods provide similar results very close to the reference values.

However, for low values of heat capacitance, the thermal-driven approach reveal significantly more accurate than the reduction strategy based on uncoupled pure thermal internal modes. Therefore, taking into account the low additional computation cost required by the thermal-driven approach, we suggest using thermal-driven modes when building reduced models of thermo-mechanical systems.

Conclusion

In this paper we have proposed an extension of the classical Craig-Bampton reduction method to thermo-mechanical systems. This new algorithm introduces an enrichment of the reduction basis allowing a more accurate approximation of the thermo-mechanical response. The idea consists in adding to the pure thermal modes their corresponding quasi-static structural response. It has been shown using a simple model of heat actuator that the quasi-static structural enrichment of the thermal modes can significantly improves the fidelity of the reduced model dynamic response. However, it also appears that the efficiency of this improvement depends on the system configuration that is not only characterized by the ratio between structural and thermal time constants. Further investigations are required on this subject to understand in which configuration the thermal-driven mode reduction method is particularly well suited.

Fig. 1 .

 1 Fig. 1. Schematic view of the thermal actuator model.

  (a) Mechanical static mode (U x ). (b) Mechanical static mode (U y ). (c) Thermal static mode (T 1 ). (d) Thermal static mode (T 2 ).(e) Thermal-driven mode 1.(f) Thermal-driven mode 2.(g) Thermal-driven mode 3.(h) Thermal-driven mode 4.

Fig. 2 .

 2 Fig. 2. Projection modes (temperature on deformed shape).

5

 5

  model using 3 thermal-driven modes Reduced model using 3 pure thermal modes Thermo-mechanical model in meters, low heat capacitance (a) Realistic heat capacitance. (b) Low heat capacitance.

Fig. 3 .

 3 Fig. 3. Transient analysis of 200 m model.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Transient analysis of 200 µm model with real heat capacitance.

Table 1

 1 Mechanical and thermal time constants of the model.

	Model length	Heat capacitance	Mechanical time constant	Thermal time constant	Ratio
	200 m	1.64 × 10 6	1.22929 × 10 1	5.5022 × 10 8	10 7
	200 m	1.64 × 10 3	1.22929 × 10 1	5.5022 × 10 5	10 4
	200 µm	1.64 × 10 6	1.22929 × 10 -5	5.5022 × 10 -4	10 1
	200 µm	1.64 × 10 3	1.22929 × 10 -5	5.5022 × 10 -7	10 -2