Representing non-linear wave energy converters in coastal wave models

Accepted Manuscript Representing Non-Linear Wave Energy Converters in Coastal Wave Models

Ewelina Luczko, Bryson Robertson, Helen Bailey, Clayton Hiles, Bradley Buckham PII:

S0960-1481(17)31141-2

DOI:

10.1016/j.renene.2017.11.040

Reference:

RENE 9442

To appear in:

Renewable Energy

Received Date:

07 April 2017

Revised Date:

11 October 2017

Accepted Date:

14 November 2017

Please cite this article as: Ewelina Luczko, Bryson Robertson, Helen Bailey, Clayton Hiles, Bradley Buckham, Representing Non-Linear Wave Energy Converters in Coastal Wave Models, Renewable Energy (2017), doi: 10.1016/j.renene.2017.11.040

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Editorial Board Journal of Renewable Energy

April 6, 2017

Re: Journal of Renewable Energy Submission – Representing Non-Linear Wave Energy Converters in Coastal Wave Models To whom it may concern, Please find the manuscript highlights detailed below: Highlights:     

Two novel candidates representations of wave energy converters are applied with SWAN; a coastal wave propagation model. New representations account for both the WEC output power but also the wave energy dissipated through non-recoverable processes; these include hydrodynamic drag and moorings. High fidelity time domain numerical models of a Backwards Bent Ducted Buoy are used to develop new spectral representations of both hydrodynamic and PTO power capture. The resulting power extracted (from the sea state) is significantly increased, while the final power production of the WEC farm is reduced. The far-field environmental impacts of this higher fidelity representation illustrates a great area of disturbance in the lee of the WEC array.

Please do not hesitate to contact me if you have any questions. I look forward to hearing from you. Sincerely yours, Bryson Robertson West Coast Wave Initiative Institute of Integrated Energy Systems (IESVic) University of Victoria www.iesvic.uvic.ca Office: (250) 472-4065 Cell: (250) 891-8117

ACCEPTED MANUSCRIPT

1 2

Representing Non-Linear Wave Energy Converters in Coastal Wave Models

3 4

Ms. Ewelina Luczko†, Dr. Bryson Robertson*†, Dr. Helen Bailey †, Mr. Clayton Hiles††, Dr. Bradley Buckham†

5 6

† Institute for Integrated Energy Systems Victoria, University of Victoria, PO Box 1700 STN CSC, Victoria, BC, Canada.

7

†† Cascadia Coast Research Ltd., Victoria, BC, Canada.

8

*corresponding author: [email protected]

9 10 11 12 13 14 15 16 17 18 19 20 21

Keywords: SWAN, Wave Energy Converter (WEC), coastal modelling, far-field effects, WEC arrays, environmental impacts.

22

1

ACCEPTED MANUSCRIPT

23

Abstract

24 25 26 27 28 29 30 31

To meet the growing global demand for carbon neutral electricity, wave energy is increasingly being identified for its global prevalence and the emerging clean technologies that can harness it. To provide the scale of electricity needed, Wave Energy Converters (WECs) will need to be deployed in large arrays. However, the effects of a WEC array’s layout on the surrounding wave conditions in both the near and far fields is poorly understood. The published literature describes numerous methods to model interactions of incident waves and WEC arrays. Yet, none of the published methods provide a methodology which spectrally resolves the individual WEC’s energy conversion characteristics, and has the flexibility to be applied to any emerging WEC design.

32 33 34 35 36 37 38 39 40 41 42 43

This study develops and implements two novel candidate high-fidelity representations of a Backwards Bent Duct Buoy Oscillating Water Column (BBDB) WEC within a Simulating WAves Nearshore (SWAN) coastal model. Each candidate representation is developed by post-processing time series data generated from a time domain numerical simulation to form a meta-model, referred to as an obstacle case, consistent with the SWAN governing equations and discretization scheme. The study is executed within the framework of SNL-SWAN – a modified version of SWAN that includes a WEC obstacle case that is built on the established concept of using SWAN’s transmission coefficient to emulate the waveWEC interactions. The two new WEC obstacle cases build on SNL-SWAN by considering the WEC’s intercepted power, the rate of energy extraction from the waves, instead of the captured power, the rate of energy conversion by the power-take-off. The two new candidates differ in their level of spectral resolution: one homogenizes the intercepted power across the frequency domain while the other attempts to spectrally resolve the energy transfer mechanisms.

44 45 46 47 48 49 50 51

The two new candidate obstacle cases were compared by considering up to 5 BBDB devices subject to irregular wave conditions. For an individual device, a region of significantly perturbed wave heights, referred to as a wake, stretched 200m in the lee of the WEC, with the maximum impact on wave heights seen 50 m behind the device. The wake behind five devices is largely dependent on the device spacing. The five unit WEC array only presented decreases of wave height between 0.2 and 0.7 metres 200 metres behind the last device in the array. A five unit WEC array model showed that estimates of individual WEC power performance was reduced by 2.3% and 6.0% when implementing the two new obstacle cases.

52 53 54 55 56 57 58 59 60 61 62

Keywords: SWAN, Wave Energy Converter (WEC), Coastal Modelling, Far-field effects, WEC arrays, environmental impacts.

1. Introduction To meet the growing global demand for carbon neutral electricity, wave energy is increasingly being identified for its global prevalence and the emergence of the necessary technologies to extract this energy. However, to provide the scale of electricity needed, Wave Energy Converters (WECs) will need to be deployed in large arrays, or farms. Arrays benefit from increased power production, power smoothing and decreased mooring costs, maintenance costs and electrical connection costs on a per unit basis [1]. Notwithstanding these merits, deployment of a large-scale WEC array faces many obstacles due to unanswered regulatory questions. The energy production of a WEC array and its

2

ACCEPTED MANUSCRIPT

63 64 65

ecological impact must be examined in detail before WEC developers can move forward with array deployments. As such, a predictive WEC array model is a prerequisite to commercial wave energy operations.

66 67 68 69 70 71

Wave Energy Converter array modeling efforts can be divided into two categories. The first are works that develop comprehensive models that fully couple the WEC motion and wave kinematics, but rely on extensive simplifications to the WEC geometry, domain size, bathymetry and wave kinematics [1]. The second sacrifice fidelity in the wave-WEC coupling but consider large numbers of WECs deployed in real world locations subject to realistic irregular wave conditions [2]. Here, we are focussed on the latter category.

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89

The premise common across these works is that the energy extraction characteristics of an individual WEC can be included in existing near shore wave propagation modeling software through semiempirical models of naturally occurring energy generation, dissipation and transmission. Greenwood et al. used the Danish Hydraulic Institute’s (DHI) MIKE 21 spectral model to simulate WECs as either a source term, an artificial island or a reactive polygon [3]. Folley and Silverthorne used TOMAWAC, a spectral action density model developed at the Electricité de France's Studies and Research Division [4], and also represented WECs as source and sink terms. Their frequency-dependent reflection and absorption terms were based on a WEC’s hydrodynamic coefficients. To determine how large WEC arrays would impact the Portuguese coast, Palha et al. translated knowledge of a Pelamis WEC’s energy conversion performance into a frequency dependent schedule of dissipation factors that were applied inside REFDIF [5] . Millar et al. investigated the effects of WEC arrays on shoreline conditions for the Cornwall WaveHub project in the United Kingdom using the open source wave model SWAN (Simulating Waves Nearshore) by modelling entire arrays as single 4 km long partially transmitting obstacles [6]. The method used a pre-existing function originally intended for modeling linear coastal protection structures, referred to as an obstacle case [7]; a transmission coefficient assigned to the obstacle was used to specify a notional percentage of energy transport to be removed from the modeled wave spectrum between grid points straddling the WEC array. Iglesias and Veigas [8] conducted similar analyses for a wave farm in Tenerife, Canary Islands, Spain.

90 91 92 93 94 95 96 97 98 99 100 101 102 103

The choice of the SWAN framework has dominated very recent works. A modified version of SWAN was developed by Smith et al. in 2012, to allow for user input of frequency and direction dependent transmission coefficients [9]. Carballo and Iglesias used scale model tank tests of the WaveCat WEC, an overtopping type WEC, to establish transmission coefficient values for two wave spectra characteristic of the Iberian Peninsula [10] The experimentally determined transmission coefficients were then applied in SWAN to study the attenuation of wave heights in the lee of two different wave farm layouts under these same two wave conditions. Abanades et al. adopted a single frequency independent value of the WaveCat transmission coefficients in a SWAN study of the impact of a WaveCat array deployed at the Cornwall Wave Hub on beach erosion [11]. In an attempt to account for the dependence of WEC performance on the changing wave conditions, Sandia National Lab has refined the SWAN open source code, referred to as SNL-SWAN, to include a new type of obstacle case in which the transmission coefficient is calculated based on the observed wave spectrum and a user supplied WEC power matrix [12]. This approach allows the transmission coefficient to vary as wave conditions change across relatively large time and spatial scales of a SNL-SWAN analysis.

104 105

However, similar to the strategy of Palha in [5], the WEC specific power matrix that SNL-SWAN depends on provides measures of the power produced by the device, not the total power extracted from the

3

ACCEPTED MANUSCRIPT

106 107 108 109 110 111 112

propagating waves. To differentiate the various stages of the WEC power conversion process, we adopt terminology suggested by Price [13]. The Intercepted Power accounts for the total power removed from the wave field, which is subsequently either converted to a useful output or is lost to the environment. Captured Power refers to the mechanical power captured by the PTO, and Delivered Power refers to the refined and conditioned energy commodity that is delivered to shore. Hydrodynamic Losses comprises the various avenues for power transfer to the environment: eddies due to form drag, turbulence, mooring line drag, etc.

113 114

Figure 1: Power and Loss Terminology. Edited version from Price [13]

115 116 117 118 119 120 121 122 123 124

At present, SNL-SWAN uses an input of Captured Power, which has been measured in a tank or field test and summarized in a power matrix, as a proxy for Intercepted Power, which is not easily measured from a physical trial. Thus, hydrodynamic losses are ignored. The obstacle case framework developed in the current work seeks to remedy this source of error by assimilating knowledge of the wave-WEC hydrodynamic interaction generated using a flexible high-fidelity WEC time domain simulation tool [14] . Two new obstacle cases are developed that use the simulator outputs in place of a standard WEC power matrix. The goal is to provide a flexible process that infuses better accuracy and greater spectral detail in the transmission coefficient calculations, while maintaining the ability of SNL-SWAN to execute on a time scale that allows for it to be used in array siting, layout optimization, a priori assessment of far-field impacts and annual energy production assessments.

125 126 127 128 129

It is important to highlight that, due to the lack of publically available data at the necessary scale, a thorough validation of the sensitivity of the SNL-SWAN calculations to the choice of Captured or Intercepted Power as the driving input is beyond the scope of study. However, the results from this work help inform the prospective methods, measurement metrics and spatial sensitivities required to plan a successful validation study.

130 131 132 133 134 135 136

This paper is structured as follows: Section 2 provides an overview of SWAN, the time domain WEC modelling strategy and the time domain post-processing algorithms that underpin the WEC metamodels, or obstacle cases. Of note, Section 2.2 provides an overview of SNL-SWAN’s existing obstacle case, and the two novel WEC obstacle cases developed in the present work. Section 3 presents the calculated far-field effects and array power production estimates for a test case 5 unit WEC array in SNLSWAN. Finally, Section 4 provides a discussion of the limitations of developed obstacle method and assumptions used in this work, while Section 5 provides a quick concluding synopsis of the differences

4

ACCEPTED MANUSCRIPT

137 138 139

between predictions of WEC array power production and environmental impacts made using the different obstacle cases.

2. SWAN and the representation of WECs

140 141 142 143 144

SWAN is a third generation, implicit, spectral action density model commonly used to model the transformation of offshore waves as they move closer to shore. SWAN’s implicit formulation allows for larger time steps and higher spatial resolution across both structured and unstructured grids. SWAN is capable of accounting for non-linear wave phenomena such as bottom friction, wave breaking, quadruplet and triad wave interactions [15].

145 146 147 148 149 150 151 152

SWAN’s governing spectral action density equation describes how a wave spectrum evolves over the computational domain, as presented in Eq. (1). Note that SWAN uses action density (N), rather than energy density, since action density is conserved in the presence of currents. The first three terms on the left hand side of Eq. (1) denote the propagation of the wave energy spectrum over time, t, and space x and y, respectively. The fourth term represents the depth and current induced refraction. Finally, the fifth term on the left hand side represents the effect of shifting frequencies due to variations in depth and mean currents. The right hand side contains source and sink terms that non-linearly redistribute energy. 2

2

𝐾 𝑡 ∂𝑁(𝜎, 𝜃)𝑑𝑦 ∂𝑁(𝜎, 𝜃)𝑑𝜃 ∂𝑁(𝜎, 𝜃)𝑑𝜎 𝑆 ∂𝑁(𝜎, 𝜃) 𝐾 𝑡 ∂𝑁(𝜎, 𝜃)𝑑𝑥 + + + + = ∂𝑡 ∂𝑥 𝑑𝑡 ∂𝑦 𝑑𝑡 ∂𝜃 𝑑𝑡 ∂𝜎 𝑑𝑡 𝜎

(1)

153 154 155 156

where 𝑥 and y correspond to the space and 𝜃 and 𝜎 correspond to directional and frequency space respectively [16]. The source term S is comprised of several contributions - third generation spectral action density models include wave growth by wind, triad wave interactions, quadruplet wave interactions and wave decay due to bottom friction, whitecapping and depth induced wave breaking.

157 158 159 160 161 162

Within the baseline version of SWAN V41.01, an obstacle can be implemented through the transmission coefficient, 0
163

(2)

(𝐻𝑆)𝑖𝑛𝑐 2

Since power transport is proportional to 𝐻𝑆: 2

𝐾𝑡 =

164 165 166 167 168

(𝐻𝑆)𝑙𝑒𝑒

(𝑃)𝑙𝑒𝑒 (𝑃)𝑖𝑛𝑐

=1‒

𝑃𝐼𝐶𝑃 𝐽∙𝑤

= 1 ‒ 𝑅𝐶𝑊

(3)

Where P refers to the total wave energy transport (kW) along the direction linking the two nodes – the x direction in the case shown in Figure 2. Since the obstacle is placed orthogonal to a linkage between grid points, only the component of wave energy transport orthogonal to the obstacle is impacted by the transmission coefficient. The total power intercepted by the WEC is PICP, J is the power transport per unit width of wave crest, w is the physical width of the device and RCW is the relative capture width.

5

ACCEPTED MANUSCRIPT

169 170

Note that here RCW is being calculated using the intercepted rather than the captured power as we are focussed on power being extracted from the waves.

171 172 173 174 175 176 177

It is important to note that the established interpretation of an obstacle as a partial barrier the width of the actual WEC that is placed orthogonal to model grid lines is maintained in this work. As is discussed further in Section 4, this interpretation could be problematic when considering some combinations of directional seas and converter geometries. In this work, the test cases are limited to waves propagating primarily in the x –direction and WECs aligned with that incident wave field. As such, the study focusses on determining the sensitivity of the SWAN outputs to the Kt calculations - not the sensitivity to the underlying physical interpretation of the obstacle.

Figure 2: Representation of an obstacle within SWAN. The device is represented as a line (red, diamond nodes) that interrupts a portion of the action density propagation between nodes of the SWAN model grid (red dots). In the current work the transmission coefficient is applied to the xprojection of the incident energy transport (energy propagating perpendicular to the obstacle) while energy in the y projection is not impacted by the device’s operation. 178 179 180 181 182 183 184 185 186 187 188 189 190

2.1.

WEC Description

To develop the obstacle case representations, the hydrodynamic losses and PTO capture power of a WEC must be determined and charted in each sea state prior to being exogenously applied within each obstacle case. For this study, a reference Backward Bent Duct Buoy (BBDB) Oscillating Water Column (OWC) WEC was chosen to both illustrate the energy balance methodology, and used to provide test case results.

2.1.1. WEC Description The Backward Bent Duct Buoy (BBDB) Oscillating Water Column (OWC) [20], [21], consists of an air chamber, an L-shaped duct, bow and stern buoyance modules, a biradial impulse turbine and a generator. A dimensional drawing is presented in Figure 3 accompanied by a rendering of the WEC in Figure 4. The mooring system was based on the design proposed by Bull and Jacob using the same mooring line material, lengths and positions of floats [22], [14].

191

6

ACCEPTED MANUSCRIPT

Figure 3: Dimensions of the BBDB OWC Reference Model [14]

192 193 194 195 196 197 198 199 200 201 202

Figure 4: OWC and mooring depiction in ProteusDS

2.1.2. WEC Simulation Architecture The time-domain WEC simulation architecture used to generate the new obstacle cases involved the dynamic coupling of three different software packages; WAMIT, ProteusDS and Simulink. The simulator, assumptions and limitations is fully described in [14]. Each software package has a distinct role on developing the necessary high-fidelity model outputs: 1. WAMIT - a linear potential flow model used to establish frequency dependent hydrodynamic coefficients that are determined based on the WEC’s geometry. 2. ProteusDS - a time domain simulator used to calculate a WEC’s dynamics over time using both linear and non-linear force contributions. 3. Simulink – a graphical solver and dynamic simulation tool that is coupled with ProteusDS to characterize the thermodynamics and air turbine dynamics within the OWC’s air chamber [14].

7

ACCEPTED MANUSCRIPT

203 204 205

Figure 5: Summary of the complete WEC modelling architecture and the different software packages used to implement it [14].

206 207 208 209 210 211 212 213 214

The WEC is modeled as two bodies: the main WEC hull and a light piston. The light piston represents the water surface within the air chamber. The model architecture and information passed between the three software packages follows the methods used in [14]. The boundary element method (BEM) code WAMIT is used to calculate the excitation forces and moments, frequency dependent added mass, added damping and the hydrostatic stiffness matrix; which are all used within ProteusDS. ProteusDS calculates the hydrodynamic excitation, radiation, viscous drag, buoyancy forces, mooring forces and the force between the hull and the air chamber pressure. Hydrodynamic and hydrostatic parameters calculated in WAMIT are used as inputs to propagate the WEC’s dynamics over time using an adaptive, variable step Runge-Kutta solver [23], [24].

215 216 217 218 219 220 221 222

The hydrodynamic forces acting on the WEC are calculated for each degree of freedom (DOF). The excitation force is the summation of the dynamic pressure across the stationary body from the incoming wave and the resulting diffracted waves [25]. These are obtained directly from WAMIT for different wave frequencies, directions and each DOF. To track the thermodynamic state (temperature, pressure and density) of the air in the chamber, the ProteusDS simulation was dynamically coupled with a Simulink model of the air chamber and turbine. The Simulink model uses knowledge of the platform and water column motions to set the absolute pressure within the air chamber, the differential pressure induced forces on the hull and water column, and the volumetric air flow rate though the air turbine.

223 224 225 226 227 228 229

2.1.3. System External Forces and Power Transmission To refine the calculation of the transmission coefficient Kt, we are concerned with the power transmission from the incident waves to the two WEC bodies (the WEC hull and the light piston), and power dissipated through drag (hydrodynamic losses in Figure 1). Since the WEC simulator internally resolves all of the forces acting between the wave, the WEC hull and the light piston as well as the full 6DOF motion of the two bodies, these power transmissions can be resolved through force-velocity products. Here, a ‘force-velocity’ product refers to a sum of the force-velocity and moment-angular

8

ACCEPTED MANUSCRIPT

230 231 232

velocity products taken across all 6 DOF of the WEC. As an example of the force data that is available from the simulator, consider the equations of motion for the floating, surface piercing, WEC hull in Eq. (4):

(

𝑴𝒙(𝑡) + 𝑨(∞)𝒙(𝑡) +

∫ 𝒌(𝑡 ‒ 𝜏)𝒙(𝜏)𝑑𝜏) + (𝑪𝒙(𝑡)) = (𝑭 (𝑡)) + (𝑭 𝒕

𝑬

𝟎

𝑷𝑻𝑶(𝑡)

(4 ) + ( 𝑭𝒎(𝑡)) + (𝑭𝒗(𝑡)) + (𝑭𝒊𝒓(𝑡) ) )

233 234 235 236

In Eq. (4) 𝑴 is the mass of the object, 𝑨(∞) is the object’s added mass at infinite frequency. The variable 𝒌(𝑡 ‒ 𝜏) is an object’s impulse response kernel, 𝑪 is the object’s hydrodynamic stiffness matrix, and 𝒙 is a vector indicating the WEC’s position in six DOF. The term in parentheses on the left hand side of Eq. (2) is the radiation force on the hull. 𝑭𝑷𝑻𝑶 represents the forces induced on the hull by the pressure in

237

the air chamber; the WECs Power Take Off (PTO). 𝑭𝒎 represents the forces exerted by the moorings. 𝑭𝒗

238

represents the damping forces and moments induced by viscous drag. 𝑭𝒊𝒓 represents the internal

239 240

reaction forces induced by the light piston (the internal water surface) on the WEC hull that can lead to energy being dissipated due to internal drag. Finally, 𝑭𝑬 corresponds to the excitation force and

241

moments induced by the pressure field associated with the scattered wave field.

242

Referring to Eq. (3), the transmission coefficient is characterized by the intercepted power (𝑃𝐼𝐶𝑃)

243

which is defined as the difference between the incident power 𝑃𝐼𝑁𝐶 (‘+’ when flowing from the waves

244

to the two bodies) and the radiated power 𝑃𝑅𝐴𝐷 (‘+’ when flowing from the two bodies to the waves): 𝑃𝐼𝐶𝑃 = 𝑃𝐼𝑁𝐶 ‒ 𝑃𝑅𝐴𝐷

245 246

(5)

The incident power is calculated through knowledge of the excitation and drag forces acting on the hull and piston: 𝑃𝐼𝑁𝐶 = 𝑃𝑖𝑛𝑐,ℎ + 𝑃𝑖𝑛𝑐,𝑝 + 𝑃𝑑𝑟𝑎𝑔,𝑖𝑚

(6)

247 248

All power values are averages calculated over a 20 min duration of steady state operation. In Eq. (6) 𝑃𝑖𝑛𝑐,ℎ and 𝑃𝑖𝑛𝑐,𝑝 are the average values of the products of the excitation forces on the hull and piston,

249

respectively, and the velocities of those bodies, and 𝑃𝑑𝑟𝑎𝑔,𝑖𝑚 is the average power transfer to the WEC

250

hull via the drag force acting on the hull. There is no drag force between the waves and the light piston.

251 252 253 254 255 256 257 258

Figure 6 shows the different manifestations of power transmission due to drag. In the top two scenarios the drag force experienced by the WEC, and the equal and opposite reaction on the waves, result in power being lost from both the WEC and the incident waves in the form of a wake (considered here to be power transfer to ambient fluid that can’t be recaptured); in the bottom two cases, whichever of the WEC or the incident waves that has the larger velocity is the power source, and power flows to the other entity and to the wake. For any instant in the time series, each of the three possible power transmissions can be resolved using knowledge of the drag force and moment existing between the waves and the WEC and the velocities of each entity: 𝑃𝑑𝑟𝑎𝑔,𝑖𝑚 is non zero when the waves and the WEC

259

are moving in the same direction and can be both positive or negative in Eq. (6); 𝑃𝑑𝑟𝑎𝑔,𝑟𝑚 is non-zero

260 261

when the WEC and the waves are moving opposite one another, or when the WEC velocity exceeds the wave velocity; 𝑃𝑑𝑟𝑎𝑔,𝑓 is also non-zero when the WEC and the waves are moving opposite one another,

262

or when the wave velocity exceeds the WEC velocity. Both 𝑃𝑑𝑟𝑎𝑔,𝑟𝑚 and 𝑃𝑑𝑟𝑎𝑔,𝑓 can only exist as power 9

ACCEPTED MANUSCRIPT

263 264 265

transmission to the ambient fluid. As an example on how the different drag induced power transmissions are resolved, consider an instant corresponding to scenario 4 in Figure 6. For a single DOF, if the drag force on the hull is 𝐹𝐷, then: 𝑃𝑑𝑟𝑎𝑔,𝑖𝑚 = 𝐹𝐷 ∙ 𝑣𝐵

(7)

𝑃𝑑𝑟𝑎𝑔,𝑓 = 𝐹𝐷 ∙ 𝑣𝑊 ‒ 𝐹𝐷 ∙ 𝑣𝐵 266

To complete the calculation of 𝑃𝐼𝐶𝑃 in Eq. (5), the total power radiated by the WEC back into the wave

267

field is the sum of the radiation and cross radiation power contributions: 𝑃𝑅𝐴𝐷 = 𝑃𝑟𝑎𝑑,𝑝 + 𝑃𝑟𝑎𝑑,ℎ + 𝑃𝑟𝑎𝑑,𝑝2ℎ + 𝑃𝑟𝑎𝑑,ℎ2𝑝

(8)

268

where 𝑃𝑟𝑎𝑑,𝑝 is the product of the radiated force of the piston and the velocity of the piston, 𝑃𝑟𝑎𝑑,ℎ is the

269

product of the radiated force of the hull and the velocity of the hull, 𝑃𝑟𝑎𝑑,𝑝2ℎ is the product of the

270 271

radiated force the piston exerts on the hull within the OWC air chamber and the velocity of the piston, 𝑃𝑟𝑎𝑑,ℎ2𝑝 is the product of the radiated force the hull exerts on the piston and the velocity. For the BBDB

272

considered, the dominant cross radiation terms were in heave, surge and pitch degrees of freedom.

Figure 6. The four potential power transmission scenarios between the WEC hull, the Incident waves and the ambient fluid; vB refers to body velocity and vW refers to wave velocity 273 274 275 276 277 278

Equations (5), (6) and (8) can be used to solve for the intercepted power to be applied in Eq. (3). As a means of checking the constituent power calculations, an independent calculation of the intercepted power was also completed following:In Eq. (9), the intercepted power is calculated by summing all of the power leaving the wave-WEC system based on the premise that this exhausted power had to be originally extracted from the waves. In addition to the two drag induced power losses to the ambient fluid 𝑃𝑑𝑟𝑎𝑔,𝑓 and 𝑃𝑑𝑟𝑎𝑔,𝑟𝑚, 𝑃𝑚𝑜𝑜𝑟 is the power extracted from the WEC by the mooring line forces 10

ACCEPTED MANUSCRIPT

279

applied at the mooring termination points on the WEC, and 𝑃𝑖𝑛𝑡,𝑝 is the power lost due to viscous

280

dissipation induced from the relative motion of the light piston and the WEC hull.

281 282

The two techniques for calculating the intercepted power were compared across 66 different sea states and the RMS error in the two estimates of the 𝑃𝐼𝐶𝑃 time series was within 4%. For the remainder of this

283

work 𝑃𝐼𝐶𝑃 was calculated using Eq. (5).

284 285 286 287 288 289 290

2.2.

Generation of obstacles

Within wave propagation models, WECs have previously been modelled as static devices with constant transmission coefficients. Sandia National Laboratories modified the SWAN V41.01 code creating four candidate techniques to calculate the Kt value from user supplied WEC power matrix data – each technique comprising a unique obstacle case [12], [17]–[19]; here two novel obstacle cases (obstacle case five and six) are developed and compared against the baseline obstacle case one of [19].

2.2.1. Obstacle Case One

291 292

For obstacle case one, the calculations within SWAN calculate the transmission coefficient by interpolating the PTO captured power from a power matrix. The captured power, defined as 𝑃𝑃𝑇𝑂 in Eq.

293

(9) and illustrated as the Captured Power in Figure 1 is used in Eq. (3) in place of 𝑃𝐼𝐶𝑃. The same

294 295

coefficient is applied uniformly across all frequencies in the energy density spectrum. This methodology does not account for hydrodynamic losses.

296

2.2.2. Obstacle Case Five

297

In obstacle case five, the transmission coefficient is determined using the actual value of 𝑃𝐼𝐶𝑃 calculated

298 299 300 301 302 303

using Eq. (6). In Eq. (3). As shown in Figure 7 the difference between the captured power and intercepted power is significant. For example, in a sea state of Hs of 1.25 m and Te of 10.5 seconds, obstacle case one would use the captured PTO power of 44 kW (20% of the gross wave power transport) while in obstacle case five an intercepted power of 135 kW (61% of gross) would be applied. Obstacle case five is very similar in structure to that of obstacle case one. A constant transmission coefficient is applied across each frequency and directional bin of the variance density spectrum.

304

(a)

(b)

Figure 7: The captured power (a) and the intercepted power (b) of the BBDB WEC.

11

ACCEPTED MANUSCRIPT

305 306

2.2.3. Obstacle Case Six

307 308 309 310 311 312 313 314 315 316 317

Obstacle case six expands on obstacle case five by accounting for frequency dependent hydrodynamics. Most WECs are tuned to a particular frequency bandwidth and capturing frequency dependence is crucial when characterizing the performance of a WEC. To account for this, Obstacle case six uses a relative capture width (RCW) curve; the RCW captures the frequency dependent performance of the WEC in each incident sea state. Similarly to the transmission coefficient of obstacles one and five, the transmission coefficient is determined by taking the ratio of power absorbed by the device to the power present in the incident wave spectrum, but this calculation is now executed for each frequency bin. The exogenously defined RCW matrix stores RCW values corresponding to a user-defined array of frequencies within each Hs and Tp bivariate bin of the histogram (See Figure 8). For each value of RCW stored in each histogram bin, the Kt value is again arrived at using Eq. (3). Transmission coefficients are then linearly interpolated for each frequency bin within the SWAN simulation.

318 319

Figure 8: Visual representation of the relative capture width (RCW) matrix

320 321 322 323 324 325 326 327 328 329 330 331

To determine the frequency dependent RCW values, the instantaneous force and velocity time series produced using ProteusDS are used to spectrally resolve the intercepted power. For each Hs-Tp histogram bin, a Fast Fourier Transform (FFT) is implemented on all of the simulated forces and velocities produced from a ProteusDS simulation of 20 min of WEC operation (with the initial 200s of data removed to prevent transient start up motions from being considered). The following steps were then used to generate an intercepted power spectrum: 1. For each force and velocity time series considered, sinusoidal force and velocity time series were recreated for each frequency bin using both the phase and amplitude determined from the complex output of the FFT step. 2. For each frequency bin, the product of the phase resolved force and velocity time series is taken to produce an estimate of a phase resolved history of instantaneous power flow. The mean value of that power signal is taken and assigned as the amplitude for the corresponding

12

ACCEPTED MANUSCRIPT

332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347

frequency bin. This step is based on Parseval’s theorem: only force-velocity products where the force and velocity share a common frequency contribute a non-zero average power transmission. 3. The intercepted power spectrum is calculated by assembling the various power amplitudes according to the right hand side of Eq. (5). 4. In order to match the spectral resolution used in the spectral action density balance of Eq. (1), the intercepted power spectrum is rebinned – multiple frequency bins are combined to form a single bin that matches the frequency bin resolution used in SWAN. 5. The above procedure is conducted for each DOF of the WEC. For each frequency bin, power is an aggregate of contributions from surge, sway, heave, roll, and pitch). An example incident wave spectrum is shown in Figure 9 below in blue. In order to complete the RCW calculations, the incident sea spectrum is represented in terms of omni-directional wave power transport (units of W/m) by scaling the variance density spectrum bins by their corresponding group velocity. The spectrum of intercepted power (normalized by the physical width of the device w for dimensional consistency) determined using the procedure above is overlaid in red. The resulting RCW curve is displayed on the right.

348 349 350 351 352 353 354 355 356 357 358 359 360

Figure 9: Power in a Hs = 1.75m and Te = 10.5s sea state (blue) and intercepted power (red) overlaid on the left. The RCW curve associated with this sea state on the right. To clarify, obstacle cases five and six are similar. Both methods provide a transmission coefficient value for each bin of the Hs-Tp histogram, and both methods account for hydrodynamic losses. However, obstacle case six assigns a different transmission coefficient to each frequency bin used in the SWAN calculations while obstacle case five provides a frequency independent transmission coefficient. The total power removed for any single wave condition is consistent between the two cases.

3. Test Case Description and Results Two test configurations were developed to better quantify how the three WEC representations affect the wave climate and the WEC’s power production. The first configuration examines the impact a single WEC has on the surrounding wave climate. The second configuration is comprised of five WECs in series, positioned directly behind each other, relative to the primary wave heading.

13

ACCEPTED MANUSCRIPT

361 362 363 364 365 366 367

Configurations were run in a 2500 m by 2500 m computational Cartesian domain with a grid resolution of 25 m x 25 m, a 5° directional resolution and a uniform 50 meter depth. Each WEC was assumed to have a 27 m width, represented by a line in SWAN as shown in Figure 2. In the first configuration, a single WEC was positioned at (1215 m, 1215 m). In the second configuration, a 270 m separation distance was chosen between each of the five WEC, with the first WEC positioned at (770 m, 1215 m). Stratigaki [26] recommends the installation of WECs with at least a ten device diameter spacing between devices to limit far-field impacts and reduce intra-array interactions between devices.

368 369 370 371 372

The test case sea state features a 1.75 m significant wave height, 10.5 s wave energy period, an idealized Pierson Moskowitz (PM) spectrum and a cos2 directional spreading. This is most commonly occurring sea state off Amphitrite Bank, Canada in 2006 [23], [27]. To magnify the relative impacts of each obstacle case, the non-linear source terms in SWAN, such as wind growth, triad and quadruplet wave interactions, were disabled.

373 374

The following sections detail the impact of the higher resolution obstacle cases on the surrounding wave field and power production estimates.

375

3.1.

Impact on the Surrounding Wave Field

376 377 378 379 380 381 382

Wave energy transport is proportional to Hs2, so significant wave height is a good metric for gauging the magnitude of the WEC’s far field impacts. A two-dimensional plan view of the Hs in the field surrounding the single WEC is presented in Figure 10. As expected, each of the obstacle case representations (one, five and six) impacted the surrounding wave field. Obstacle cases five and six attenuate Hs more than obstacle case one; indicating that characterizing the WEC interaction with the surrounding wave field using only the captured PTO power underestimates the true impact of the WEC on the environmental conditions.

383 384 385

In obstacle case one, Hs recovers to 1.72 m 100 m behind the WEC. For obstacle case five and six, it takes 330 m and 260 m respectively for this same level of recovery; therefore obstacle case six has less of an impact on the significant wave height than obstacle case five.

14

ACCEPTED MANUSCRIPT

386 387 388 389

Figure 10: Significant wave height comparison demonstrating a WEC's impact on the surrounding wave field when represented with different obstacle cases in the most commonly occurring sea state off of Amphitrite Bank (Hs=1.75 m, Te = 10.5 s).

390 391 392 393 394 395 396

The difference between the obstacle cases becomes more apparent when the wave spectra in the lee of the WEC are studied. Figure 11 presents the omni-directional variance density spectra 5 m, 25 m, 100 m and 200 m behind the WEC. As shown, the variance spectra modified by obstacle cases one and five are scaled forms of the incident variance spectrum due to the use of a frequency independent transmission coefficient. The spectrum corresponding to obstacle case six has been non-uniformly decreased with the most power being removed at the spectrum’s peak (0.09 Hz) and to the right of the peak (0.11 Hz); where the WEC operates most efficiently [14].

397

15

ACCEPTED MANUSCRIPT

398 399 400

Figure 11: Comparison of non-directional variance density spectra at various distances behind the WEC when the WEC is represented with obstacle cases one, five and six.

401 402 403 404 405 406

The recovery of Hs directly in the lee of five WECs is presented in Figure 12. Firstly, all obstacle cases result in a similar reduction pattern in Hs. Obstacle case five consistently reduces the significant wave height more than obstacle case six. Initially, the difference is merely 0.02 m between the two, but by the fifth WEC the difference in significant wave height reduction increases to 0.04 m. Small differences in Hs are observed at the last three WECs. As the incident Hs decreases, a greater transmission coefficient is employed for this particular sea state, thus reducing the power absorbed at each successive WEC.

407 16

ACCEPTED MANUSCRIPT

408 409 410

Figure 12: Comparison of the profile view of the Hs and its recovery in the lee of five WECs when the WECs are represented with different obstacle cases

3.2.

Impact on power production

411 412 413

For the multi-WEC configuration, the captured power and percentage decrease experienced by each WEC is reported in Table 1. The percentage decrease is normalized by the power capture of the first WEC for each obstacle case (88.5 kW).

414 415 416 417 418 419 420 421 422 423

The greatest decrease in captured power is observed in obstacle case five, thus making this obstacle case the most conservative representation for total captured power. For each subsequent WEC, the captured power drop ranges between five and six percent of the first WEC’s power capture. In obstacle case one, the reduction in captured power for each subsequent WEC ranges between 1.5 and 2.5%, while in obstacle case six, the decrease in captured power ranges between 3.5 and 4.5%. While the percentages reported in Table 1 are relatively small, it should be noted that the WEC spacing was set following recommendations on how to mitigate intra-array interactions. Referring to Figure 12, if the 270m spacing was brought down below 200m, more significant performance drops for obstacle cases five and six are expected. Also, differences between the predictions will be exacerbated as the number of devices increases past five; referring to Figure 11 the .

424 425

Table 1: Captured power [kW] produced by five WECs placed ten device widths behind the previous for each obstacle case Captured Power [kW] for each WEC in series (Percentage decrease [%] in captured power) Obstacle case 1 5 6

426

1

2

3

88.50 88.50 88.50

87.13 (1.54%) 83.79 (5.32%) 85.23 (3.69%)

86.69 (2.05%) 83.39 (5.77%) 84.91 (4.06%)

4

5

86.52 (2.24%) 86.45 (2.32%) 83.35 (5.82%) 83.34 (5.83%) 84.82 (4.16%) 84.79 (4.19%)

4. Discussion

427 428 429 430

The test cases presented above quantify the impact a WEC has on the surrounding wave climate while also allowing one to determine how much power was captured by the reference WEC. While these differences are primarily driven by different meta-model representations (see section 2.2), there are nuances in the modelling of WECs as obstacles that must be noted.

431 432 433 434 435 436 437 438

Firstly, the results generated were based on simulations within ProteusDS that used wave spectra with 20 frequency bins and 7 directional bins, while SWAN uses 36 frequency and 72 direction bins. In the generation of the RCW curves for each bin of the histogram, the force and velocity data taken from Proteus were thus relatively discrete – content was concentrated around the center frequencies of the 20 bins used to create the simulated waves. When the intercepted power spectrum are rebinned according the SWAN bin structure, this can lead to somewhat jagged RCW curves. A higher discretization in the ProteusDS spectrum would lead to more robust estimates of the RCW curves, however would incur exponentially increasing computational costs.

439 440

Secondly, the obstacle cases presented here are hindered by the WEC’s representation as a line intersecting the SWAN grid line linking two grid nodes. Using this representation, the obstacle

17

ACCEPTED MANUSCRIPT

441 442 443 444 445 446 447

interrupts energy transport only in the direction of the grid line. However, WECs can interact with waves across all three dimensions (width, length and depth) and it has been shown that some WECs extract more power than is available, when represented solely by their width. The WEC considered in this work never exhibited transmission coefficients in significant excess of unity, but for attenuators this would be expected. Additionally, the use of the current obstacle model limits attempts to resolve the injection of power into the wave field by radiation – as formulated, radiated waves would have to be directed solely along the grid line passing through the obstacle.

448 449 450 451 452 453 454 455 456 457 458

Finally, similar studies have been conducted using different power transfer functions, yet all ultimately converging on similar conclusions. Greenwood et al., found that the wave energy in the lee of a wave surge oscillator modeled in MIKE SW recovered within 100 m when the WEC was represented as a reactive polygon, and 400 - 500 m when represented as an artificial island [3]. Numerical simulations conducted in MILDWAVE for Wave Dragon showed that wave energy reductions of 65% still existed 3000 m leeward in short-crested wave fields, yet were 85% recovered within 1000 metres of the device for long crested waves [28]. The difference in leeward wake is a function of the individual WEC’s geometry, the particular seastate and the means of representation. As these factors change between studies it is difficult to make effective comparisons. A validation exercise is needed in order to help assess the validity, uncertainty and limitations of the newly presented obstacle cases and other previously disseminated techniques.

459

5. Conclusions

460 461 462 463 464 465 466

In order to meet utility scale electricity demands, Wave Energy Converters will need to be deployed in arrays of devices. Array deployments will leverage economies of scale for WEC construction, lower perWEC interconnection costs and result in smoother overall power outputs. However, there currently is a dearth of suitable wave propagation models that allow for arrays of WECs to be accurately modeled, their far field impacts to be studied, and predictions concerning their captured power to be made. This study develops and implements two novel methodologies to represent WEC dynamics within a Simulating WAves Nearshore (SWAN) model.

467 468 469 470

In SWAN, WECs have traditionally been represented by static obstacles with the transfer and extraction of energy defined by a static transmission coefficient. SNL-SWAN, developed by Sandia National Laboratories, is an enhanced version of SWAN and allows for different transmission coefficients to be calculated based on a user supplied power matrix or an RCW curve.

471 472 473 474 475 476 477

In this work, two novel WEC representations are introduced and compared against the baseline SNL transmission coefficient technique. In obstacle case one, the original SNL-SWAN technique uses only the captured PTO power. In obstacle case five and six, the intercepted power is used. Intercepted power accounts for both the captured PTO power but also the power dissipated due to hydrodynamic losses. In obstacle case one and five, a frequency independent transmission co-efficient is used. In obstacle case six, the total intercepted power (captured and dissipated through hydrodynamic losses) is applied via different transmission coefficients for each wave frequency.

478 479 480

As a case study BBDB OWC was used to quantify the effects of various obstacle cases. The BBDB OWC includes many complicating dynamic factors – moorings, surface piercing, and a full 6 DOF range of motion – which make it an ideal test WEC.

18

ACCEPTED MANUSCRIPT

481 482 483 484 485 486 487

For a single WEC, the wave spectrum is most affected within the first 50 m behind the WEC. Once the incident spectrum propagates 200 m behind the WEC, regardless of the representation, the incident spectrum has recovered and the far-field effect limited. The greatest reduction in Hs is observed when the WEC’s intercepted power is described with obstacle case five, a constant transmission coefficient accounting for all intercepted power. Even though obstacle case five creates the greatest decrease in significant wave height in the lee of the incident sea, this representation also causes the surrounding wave field to recover the most quickly.

488 489 490 491 492 493 494 495 496

When an array of five WECs are simulated, the incident Hs to a WEC is reduced with every additional preceding WEC. Obstacle case five reduces the Hs most of the three representations analyzed. Across the three obstacle case representations, the decrease in Hs does not change drastically with each subsequent WEC, however the lateral extent to which the Hs field around the WEC is impacted does increases. When obstacle case one is employed, the power captured by the second WEC is reduced by 1.54% when compared to the first, while the power captured by the fifth WEC is only reduced by 2.30% of the first WEC. When obstacle case five is employed, the power captured by the second and fifth WEC is reduced by 5.56% and 6.00%, respectively. Finally, when obstacle case six is employed, the power captured by the second and fifth WEC is reduced by 2.30% and 2.78% respectively.

497

Bibliography

498 499

[1]

B. F. M. Child and V. Venugopal, “Optimal configurations of wave energy device arrays,” Ocean Eng., vol. 37, no. 16, pp. 1402–1417, 2010.

500 501 502 503

[2]

M. Folley, A. Babarit, B. F. M. Child, D. Forehand, L. O’Boyle, K. Silverthorne, J. Spinneken, V. Stratigaki, and P. Troch, “A Review of Numerical Modelling of Wave Energy Converter Arrays,” in ASME 2012 31st International Conference on Ocean, Offshore and Arctic Engineering, 2012, pp. 1–12.

504 505 506

[3]

C. E. Greenwood, D. Christie, and V. Venugopal, “The Simulation of Nearshore Wave Energy Converters and their Associated Impacts around the Outer Hebrides,” 10th Eur. Wave Tidal Energy Conf. (EWTEC 2013), 2013.

507 508

[4]

K. E. Silverthorne and M. Folley, “A New Numerical Representation of Wave Energy Converters in a Spectral Wave Model.”

509 510 511

[5]

A. Palha, L. Mendes, C. J. Fortes, A. Brito-Melo, and A. Sarmento, “The impact of wave energy farms in the shoreline wave climate: Portuguese pilot zone case study using Pelamis energy wave devices,” Renew. Energy, vol. 35, no. 1, pp. 62–77, 2010.

512 513

[6]

D. L. L. Millar, H. C. M. C. M. Smith, and D. E. E. Reeve, “Modelling analysis of the sensitivity of shoreline change to a wave farm,” Ocean Eng., vol. 34, no. 5–6, pp. 884–901, 2007.

514

[7]

Delft University, “USER MANUAL SWAN - Cycle III version 41.01A,” 2009.

515 516

[8]

G. Iglesias and M. Veigas, “Wave farm impact and wave-WEC interaction : a sensitivity analysis,” Proc. 10th Eur. Wave Tidal Energy Conf., Aalborg, Denmark, 2013.

517 518 519

[9]

H. C. M. Smith, C. Pearce, and D. L. Millar, “Further analysis of change in nearshore wave climate due to an offshore wave farm: An enhanced case study for the Wave Hub site,” Renew. Energy, vol. 40, no. 1, pp. 51–64, 2012.

520 521

[10]

R. Carballo and G. Iglesias, “Wave farm impact based on realistic wave-WEC interaction,” Energy, vol. 51, pp. 216–229, 2013.

19

ACCEPTED MANUSCRIPT

522 523

[11]

J. Abanades, D. Greaves, and G. Iglesias, “Wave farm impact on the beach profile: A case study,” Coast. Eng., vol. 86, pp. 36–44, 2014.

524

[12]

K. Ruehl, C. Chartrand, and A. K. Porter, “SNL-SWAN User’s Manual,” 2014.

525

[13]

A. A. E. Price, “New perspectives on wave energy converter control,” 2009.

526 527

[14]

H. Bailey, B. Robertson, and B. J. Buckham, “Wave-to-Wire Simulation of a Floating Oscillating Water Column Wave Energy Converter.” 2016.

528 529

[15]

SWAN, “SWAN Scientific and Technical Documentation,” Delft University of Technology, The Netherlands, 41.01A, 2006.

530

[16]

SWAN, “SWAN Technical Documentation,” Delft Univ. Technol. Netherlands, 2007.

531 532

[17]

A. Porter, K. Ruehl, C. Chartrand, and H. Smith, “Development and release of the open-source wave climate environment assessment tool SNL-SWAN,” pp. 1–5, 2015.

533 534

[18]

A. Porter, K. Ruehl, and C. Chartrand, “Further Development of Snl - Swan , a Validated Wave Energy Converter Array Modeling Tool,” 2nd Mar. Energy Technol. Symp., pp. 1–9, 2014.

535 536 537

[19]

K. Ruehl, A. Porter, A. Posner, and J. Roberts, “Development of SNL-SWAN , a Validated Wave Energy Converter Array Modeling Tool,” Proc. 10th Eur. Wave Tidal Energy Conf., Aalborg, Denmark, 2013.

538 539 540

[20]

D. Bull, C. Smith, D. S. Jenne, P. Jacob, A. Copping, S. Willits, A. Fontaine, D. Brefort, G. Copeland, M. Gordon, and R. Jepsen, “Reference Model 6 ( RM6 ): Oscillating Wave Energy Converter,” vol. 6, no. September, 2014.

541

[21]

U.S DOE, “Reference Model Project (RMP),” 2014. .

542 543

[22]

D. Bull and P. Jacob, “Methodology for creating nonaxisymmetric WECs to screen mooring designs using a Morison Equation approach,” in Oceans. IEEE., 2012.

544 545 546

[23]

B. Robertson, H. Bailey, D. Clancy, J. Ortiz, and B. Buckham, “Influence of wave resource assessment methodology on wave energy production estimates,” Renew. Energy, vol. 86, pp. 1145–1160, 2016.

547

[24]

Dynamic Systems Analysis, “ProteusDS 2015 Manual,” Victoria, BC, 2015.

548 549

[25]

J. Falnes, Ocean Waves and Oscillating Systems: Linear Interactions Including Wave Energy Extraction. Cambridge: Cambridge University Press, 2004.

550 551 552 553

[26]

V. Stratigaki, P. Troch, T. Stallard, D. Forehand, J. Kofoed, M. Folley, M. Benoit, A. Babarit, and J. Kirkegaard, “Wave Basin Experiments with Large Wave Energy Converter Arrays to Study Interactions between the Converters and Effects on Other Users in the Sea and the Coastal Area,” Energies, vol. 7, no. 2, pp. 701–734, 2014.

554 555

[27]

E. Luczko, “Assimilating a Higher Fidelity Representation of Wave Energy Converters in a Spectral Model,” M. A. Sc. Thesis University of Victoria, 2016.

556 557 558

[28]

C. Beels, P. Troch, K. De Visch, J. P. Kofoed, and G. De Backer, “Application of the time-dependent mild-slope equations for the simulation of wake effects in the lee of a farm of Wave Dragon wave energy converters,” Renew. Energy, vol. 35, no. 8, pp. 1644–1661, 2010.

559

20