Wave Energy Extraction by Flexible Floaters

Abstract

We present a novel mathematical model to investigate the extraction of wave power by flexible floaters. The model is based on the method of dry modes, coupled with a matched eigenfunction expansion. Our model results compare satisfactorily with preliminary data obtained from a demonstrator device, developed at the University of Groningen. We show that the role of elasticity is to increase the number of resonant frequencies with respect to a rigid body, which has a positive effect on wave power output. The mathematical model is then extended to irregular incident waves, described by a JONSWAP spectrum. Our results show that the peak capture factors decrease in irregular waves, as compared to the monochromatic case. However, the system becomes more efficient at non-resonant frequencies. This work highlights the need to scale-up experimental investigations on flexible wave energy converters, which are still a small minority, compared to those on rigid converters.

1. Introduction

This paper presents novel mathematical and demonstrator models of wave power extraction from flexible floaters. Our groups at Loughborough University and the University of Plymouth, in collaboration with the University of Groningen and Ocean Grazer BV, have recently started an investigation into innovative wave energy converters (WECs), with the goal of decreasing the levelised cost of energy (LCOE) of wave power generation. Indeed, the sheer size and complexity of many of the WEC devices proposed and tested during the past couple of decades has so far hindered their scalability and commercialisation [1]. To overcome such challenges, we consider the possibility of using light and flexible materials (e.g., silicone rubber), instead of bulky metallic components, in the design of the prime mover.

To this aim, Renzi [2] analysed the coupled hydro-electromechanic response of a bimorph plate made by a flexible substrate intertwined with piezoceramic layers. These allow the transformation of the plate elastic motion into useful electricity, by means of the piezoelectric effect. Renzi [2] shows that the piezoelectric plate can extract sufficient energy for low-power devices, like sensors, LEDs, computers and wireless routers. Later, Buriani and Renzi [3] also showed that connecting the flexible piezoelectric device to a vertical wall (e.g., a caisson breakwater) significantly enhances its performance in small-amplitude waves.

The idea of using flexible floaters to extract energy from the ocean has been recently pushed forward by Zheng et al. [4], who investigated the hydrodynamic interaction between water waves and an array of circular porous elastic plates. An important result shown in [4] is that wave power dissipation by the array of elastic plates increases thanks to the constructive interaction between the plates, which suggests an interesting potential for wave power generation.

Further investigations on floating elastic plates include the effects of three dimensional structures on wave energy dissipation [5,6,7], and the interactions between a flexible plate and a bottom ridge [8].

The potential use of arrays of floaters to extract energy from waves has attracted the attention of the wave energy industry as well. For example, the Dutch company Ocean Grazer has recently proposed several versions of its “floater blanket” concept, an array of floater elements each connected to power take-off (PTO) systems [9]. For details on the technology, see [10]. An artist’s sketch of Ocean Grazer’s floater blanket concept is illustrated in Figure 1.

In this paper, we propose a novel mathematical model of wave energy extraction by means of flexible floaters (Section 2). The technology investigated here would correspond to a column of floaters along the incident wave direction, in the case of the Ocean Grazer WEC. The mathematical model investigates a two-dimensional flexible plate floating on the surface of the ocean, and connected to a series of linear PTO devices. By coupling the method of dry modes with a matched eigenfunction expansion, we show that the energy extraction efficiency of the device is enhanced by the bending elastic modes of the plate (Section 3). We also show novel results of a demonstrator model of flexible wave energy device made by silicone sheets (Section 4). The tests were carried out in the Faculty of Science and Engineering, University of Groningen (The Netherlands). To estimate the energy extraction potential of the device, we connected the flexible floater to a 1:35 scale model of the PTO system employed by the Ocean Grazer device. Interestingly, our experimental results show that energy absorption levels are very similar for a continuous floater and a series of single floaters of the same overall length. Finally, the mathematical model is extended to investigate the device performance in random seas. We show that the peak performance decreases in irregular waves, as compared to the monochromatic case. However, the system becomes more efficient at non-resonant frequencies (Section 5). We anticipate that these results will be of interest to wave energy companies working on the development of flexible WECs (Section 6).

2. Mathematical Model

With reference to Figure 2, consider an infinite two-dimensional channel of constant depth h and a rectangular ridge of width $2L$ and height c.

Let us define a Cartesian reference system $\left(x,z\right)$ with the x axis along the undisturbed free surface and the z axis positive upward. At $x\in \left[-L,L\right],z=0$ rests a flexible floater WEC of length $2L$ and thickness $h_p$, allowed to oscillate under the action of incident waves. The WEC is connected to the ridge through a number M of vertical power take-off (PTO) mechanisms, each with damping coefficient $\nu ={\nu }_{PTO}$ and located at $x=x_i$, $i=1,\dots ,M$. We assume $L\gg h_p$, thus the elastic vibration of the floater can be described by the Euler beam equation [11]. We assume also monochromatic incident waves of amplitude A coming from $x\to +\infty$, inviscid fluid and irrotational flow. Hence, the velocity potential $\mathsf{\Phi }\left(x,z,t\right)$ satisfies Laplace’s equation in the fluid domain $\mathsf{\Omega }\left(x,z\right)$. On the free surface, we have the linearised kinematic and mixed boundary conditions

$${\zeta }_t={\mathsf{\Phi }}_z,z=0,x\in [L,+\infty ),$$

(1)

$${\mathsf{\Phi }}_{tt}+g{\mathsf{\Phi }}_z=0,z=0,x\in [L,+\infty ),$$

(2)

where $\zeta$ is the free-surface elevation, g is the acceleration due to gravity and t is time. The subscripts denote differentiation with respect to the relevant variable. We require tangential fluid velocity at the bottom and on the ridge vertical walls, located at $x=\pm L$, i.e.,

$${\mathsf{\Phi }}_n=0,$$

(3)

where n denotes the normal derivative to the relevant surface. The kinematic boundary conditions on the wetted surface of the plate are

$${\mathsf{\Phi }}_z=W_t,z=-d,x\in \left[-L,+L\right],$$

(4)

$${\mathsf{\Phi }}_x=-\left(z-z_g\right)W_{xt},z\in \left[-d,0\right],x=\pm L,$$

(5)

where W is the vertical displacement response of the structure and $z_g$ is the coordinate of the structure’s centre of mass.

Since the system is forced by monochromatic incident waves of frequency $\omega$, we assume harmonic motion

$$\left\{\mathsf{\Phi },\zeta ,W\right\}=\mathrm{Re}\left\{\left(\varphi ,\eta ,w\right)e^{-\mathrm{i}\omega t}\right\},$$

(6)

with $\mathrm{i}$ being the imaginary unit. We now write the governing equations in terms of the spatial variables only

$$\begin{array}{cccc}\hfill {\nabla }^2\varphi & =0,\hfill & \hfill & \mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(7)

$$\begin{array}{cccc}\hfill {\varphi }_z& =-\mathrm{i}\omega \eta ,\hfill & \hfill & z=0,x\in [L,+\infty ),\hfill \end{array}$$

(8)

$$\begin{array}{cccc}\hfill -{\omega }^2\varphi +g{\varphi }_z& =0,\hfill & \hfill & z=0,x\in [L,+\infty ),\hfill \end{array}$$

(9)

$$\begin{array}{cccc}\hfill {\varphi }_z& =-\mathrm{i}\omega w,\hfill & \hfill & z=-dx\in \left[-L,+L\right],\hfill \end{array}$$

(10)

$$\begin{array}{cccc}\hfill {\varphi }_x& =\mathrm{i}\omega \left(z-z_g\right)w_x,\hfill & \hfill & z\in \left[-d,0\right],x=\pm L,\hfill \end{array}$$

(11)

$$\begin{array}{cccc}\hfill {\varphi }_n& =0,\hfill & \hfill & \mathrm{on}\mathrm{solid}\mathrm{boundaries},\hfill \end{array}$$

(12)

and require that the velocity potential $\varphi$ is outgoing for $x\to -\infty$.

Following Newman [12,13], we now decompose the displacement of the floater into a set of dry modes, i.e., in the absence of fluid or added mass. This allows us to significantly reduce the numerical computations and to obtain a deeper physical insight. A schematic diagram explaining our approach is presented in Figure 3.

Since the plate and the fluid domain are symmetric with respect to the vertical axis $x=0$, we decompose the modal expansion into symmetric and antisymmetric parts, hence

$$w\left(x\right)=\sum _{l=0}^{\infty }\left[{\xi }_l^S{f}_l^S\left(x\right)+{\xi }_l^A{f}_l^A\left(x\right)\right]$$

(13)

where the superscripts S and A denote, respectively, the symmetric and antisymmetric components, while ${\xi }_l^{S,A}$ represents the complex amplitude of the symmetric or antisymmetric lth mode. The plate satisfies the Euler beam equation with free-free end conditions; therefore, the corresponding dry modal shapes are

$$\begin{array}{cc}\hfill f_0^S& =1,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill f_l^S& =\frac{cosh\frac{{\mu }_l^{S}x}{L}}{cosh{\mu }_l^S}+\frac{cos\frac{{\mu }_l^{S}x}{L}}{cos{\mu }_l^S},l=1,\dots ,\infty \hfill \end{array}$$

(15)

and

$$\begin{array}{cc}\hfill f_0^A& =x,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill f_l^A& =\frac{sinh\frac{{\mu }_l^{A}x}{L}}{sinh{\mu }_l^A}+\frac{sin\frac{{\mu }_l^{A}x}{L}}{sin{\mu }_l^A},l=1,\dots ,\infty \hfill \end{array}$$

(17)

where $f_0^S$ and $f_0^A$ correspond to the rigid modes heave and pitch [14], respectively, while the eigenvalues ${\mu }_l^{S,A}$ are the positive real roots of the following eigenvalue conditions

$$tanh{\mu }_l^S+tan{\mu }_l^S=0,$$

(18)

$$tanh{\mu }_l^A-tan{\mu }_l^A=0.$$

(19)

We remark that these modes are orthogonal, i.e., the corresponding shapes satisfy the following properties

$${\int }_{-L}^L\mathrm{d}x{f}_i^S{f}_l^S=\left\{\begin{array}{c}0,l\ne i,\hfill \\ 2L,l=i,\hfill \end{array}\right.,{\int }_{-L}^L\mathrm{d}x{f}_i^A{f}_l^A=\left\{\begin{array}{c}0,l\ne i,\hfill \\ 2L,l=i>0,\hfill \\ \frac{2L^3}{3},l=i=0.\hfill \end{array}\right.$$

(20)

The decomposition into symmetric and antisymmetric parts allows us to analyse the half-problem in the region $x>0$ and simplify significantly the mathematical structure. Let us define the following fluid subdomains

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_1& =\left\{x\in [L,+\infty ),z\in \left[-h,0\right]\right\},\hfill \\ \hfill {\mathsf{\Omega }}_2& =\left\{x\in \left[0,L\right],z\in \left[-h+c,-d\right]\right\},\hfill \end{array}$$

(21)

where ${\mathsf{\Omega }}_1$ represents the domain to the right of the elastic plate, while ${\mathsf{\Omega }}_2$ represents the fluid region between the elastic plate and the ridge. Following the method of [15], we decompose also the velocity potential into diffraction and radiation components, i.e.,

$$\varphi ={\varphi }_D^S+{\varphi }_D^A+\sum _{l=0}^{\infty }{\xi }_l^S{\varphi }_{lR}^S+\sum _{l=0}^{\infty }{\xi }_l^A{\varphi }_{lR}^A,$$

(22)

$${\varphi }_D^{S,A}=\frac{1}{2}{\varphi }_I+{\varphi }_S^{S,A},$$

(23)

where ${\varphi }_S^S$ (${\varphi }_S^A$) is the symmetric (antisymmetric) scattering potential, ${\varphi }_D^S$ (${\varphi }_D^A$) is the symmetric (antisymmetric) diffraction potential satisfying the boundary conditions (8)–(12) with $w=0$, and ${\varphi }_{lR}^S$ (${\varphi }_{lR}^A$) is the lth symmetric (antisymmetric) radiation potential for unit amplitude that satisfies the same conditions with the unknown vertical displacement $w\ne 0$. The incident wave potential is given by

$${\varphi }^I=\frac{-\mathrm{i}Ag}{\omega }\frac{cosh{k}_0\left(h+z\right)}{cosh{k}_{0}h}{e}^{-\mathrm{i}{k}_{0}x},$$

(24)

where the wavenumber $k_0$ is the real root of the dispersion relation

$${\omega }^2=g{k}_{0}tanh{k}_{0}h.$$

(25)

Now omit the superscripts $S,A$ for the sake of brevity, and let ${\varphi }_D^{\left(1\right)}$ (${\varphi }_R^{\left(1\right)}$) be the diffraction (radiation) velocity potential in $\left(x,z\right)\in {\mathsf{\Omega }}_1$ and ${\varphi }_D^{\left(2\right)}$ (${\varphi }_R^{\left(2\right)}$) the diffraction (radiation) velocity potential in $\left(x,z\right)\in {\mathsf{\Omega }}_2$.

The boundary value problem for the subdomain ${\mathsf{\Omega }}_1$ is

$$\begin{array}{cccc}\hfill {\nabla }^2{\varphi }^{\left(1\right)}& =0,\hfill & \hfill & \mathrm{in}{\mathsf{\Omega }}_1,\hfill \end{array}$$

(26)

$$\begin{array}{cccc}\hfill -{\omega }^2{\varphi }^{\left(1\right)}+g\frac{\partial {\varphi }^{\left(1\right)}}{\partial z}& =0,\hfill & \hfill & z=0,x\in [L,+\infty ),\hfill \end{array}$$

(27)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }_{lR}^{S,A\left(1\right)}}{\partial x}& =\mathrm{i}\omega \left(z-z_g\right)\frac{\mathrm{d}{f}_l^{S,A}}{\mathrm{d}x},\hfill & \hfill & z\in \left[-d,0\right],x=\pm L,\hfill \end{array}$$

(28)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }_D^{\left(1\right)}}{\partial x}& =0,\hfill & \hfill & z\in \left[-d,0\right],x=\pm L,\hfill \end{array}$$

(29)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }^{\left(1\right)}}{\partial z}& =0,\hfill & \hfill & z=-h,x\in [L,+\infty ),\hfill \end{array}$$

(30)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }^{\left(1\right)}}{\partial x}& =0,\hfill & \hfill & x=L,z\in \left[-h,-h+c\right],\hfill \end{array}$$

(31)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }^{\left(1\right)}}{\partial x}& =\frac{\partial {\varphi }^{\left(2\right)}}{\partial x},\hfill & \hfill & x=L,z\in \left[-h,0\right],\hfill \end{array}$$

(32)

$$\begin{array}{cccc}\hfill {\varphi }^{\left(1\right)}& ={\varphi }^{\left(2\right)},\hfill & \hfill & x=L,z\in \left[-h+c,-d\right],\hfill \end{array}$$

(33)

where (32)–(33) represent, respectively, continuity of the velocity field and pressure between the fluid domains ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$. The boundary value problem for the subdomain ${\mathsf{\Omega }}_2$ is governed by

$$\begin{array}{cccc}\hfill {\nabla }^2{\varphi }^{\left(2\right)}& =0,\hfill & \hfill & \mathrm{in}{\mathsf{\Omega }}_2,\hfill \end{array}$$

(34)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }_{lR}^{S,A\left(2\right)}}{\partial z}& =-\mathrm{i}\omega f_l^{S,A},\hfill & \hfill & z=-dx\in \left[-L,+L\right],\hfill \end{array}$$

(35)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }_D^{\left(2\right)}}{\partial z}& =0,\hfill & \hfill & z=-dx\in \left[-L,+L\right],\hfill \end{array}$$

(36)

$$\begin{array}{cccc}\hfill \frac{\partial {\varphi }^{\left(2\right)}}{\partial z}& =0,\hfill & \hfill & z=-h+c,x\in \left[-L,+L\right],\hfill \end{array}$$

(37)

and the coupling matching conditions (32)–(33).

In the following sections, we determine the diffraction and radiation potential in ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$, respectively, by matching the potentials at the common boundaries.

2.1. Diffraction Potential Solution

The general solution in ${\mathsf{\Omega }}_1$ is given by

$${\varphi }_D^{S,A\left(1\right)}=\frac{1}{2}{\varphi }_I+\sum _{j=0}^{\infty }{\mathcal{C}}_j^{S,A\left(1\right)}{X}_j^{\left(1\right)}{Z}_j^{\left(1\right)}$$

(38)

where the ${\mathcal{C}}_j^{S,A\left(1\right)}$ are unknown complex constants, $Z_j^{\left(1\right)}$ is the set of orthogonal vertical eigenfunctions

$$Z_j^{\left(1\right)}=cosh{k}_j\left(h+z\right),$$

(39)

$X_j^{\left(1\right)}$ denotes the x dependence

$$X_j^{\left(1\right)}=e^{\mathrm{i}{k}_{j}x},$$

(40)

while the terms $k_j$’s are the roots of the dispersion relation [15]

$$\left.\begin{array}{cc}\hfill & {\omega }^2=g{k}_{0}tanh{k}_{0}h,\hfill \\ \hfill & {\omega }^2=-g{\overline{k}}_{j}tan{\overline{k}}_{j}h,k_j=\mathrm{i}{\overline{k}}_j,j=1,\dots ,\infty \hfill \end{array}\right\}.$$

(41)

The solution in the domain ${\mathsf{\Omega }}_2$ below the elastic plate reads

$${\varphi }_D^{S,A\left(2\right)}=\sum _{j=0}^{\infty }{\mathcal{C}}_j^{S,A\left(2\right)}{X}_j^{\left(2\right)}{Z}_j^{\left(2\right)},$$

(42)

where the ${\mathcal{C}}_j^{S,A\left(2\right)}$ are unknown complex constants and

$$Z_j^{\left(2\right)}=\left\{\begin{array}{cc}1,\hfill & \hfill j=0\\ cos\frac{j\pi \left(d+z\right)}{c+d-h},\hfill & \hfill j=1,\dots ,\infty \end{array}\right.,$$

(43)

$$X_j^{S\left(2\right)}=\left\{\begin{array}{cc}1,\hfill & \hfill j=0\\ cosh\frac{j\pi x}{c+d-h},\hfill & \hfill j=1,\dots ,\infty \end{array}\right.,$$

(44)

$$X_j^{A\left(2\right)}=\left\{\begin{array}{cc}x,\hfill & \hfill j=0\\ sinh\frac{j\pi x}{c+d-h},\hfill & \hfill j=1,\dots ,\infty \end{array}\right..$$

(45)

Usage of the boundary conditions (32)–(33) gives

$$\sum _{j=0}^{\infty }{\mathcal{C}}_j^{S,A\left(1\right)}{\left.X_j^{S,A{\left(1\right)}^{\prime }}\right|}_{x=L}{Z}_j^{\left(1\right)}-{\mathcal{C}}_j^{S,A\left(2\right)}{\left.X_j^{S,A{\left(2\right)}^{\prime }}\right|}_{x=L}{Z}_j^{\left(2\right)}={\left.\frac{\mathrm{i}{k}_0}{2}{\varphi }_I\right|}_{x=L},$$

(46)

$$\sum _{j=0}^{\infty }{\mathcal{C}}_j^{S,A\left(1\right)}{X}_j^{S,A\left(1\right)}\left(L\right)Z_j^{\left(1\right)}-{\mathcal{C}}_j^{S,A\left(2\right)}{X}_j^{S,A\left(2\right)}\left(L\right)Z_j^{\left(2\right)}={\left.-\frac{1}{2}{\varphi }_I\right|}_{x=L}.$$

(47)

Multiplying each side of (46) and (47) by $Z_i^{\left(1\right)}$ and $Z_i^{\left(2\right)}$, respectively, and integrating over the relevant intervals, $z\in \left[-h,0\right]$ and $z\in \left[-h+c,-d\right]$, yields the following inhomogeneous linear systems in the unknown coefficients ${\mathcal{C}}_j^{S,A\left(1\right)}$ and ${\mathcal{C}}_j^{S,A\left(2\right)}$

$$\left\{\begin{array}{c}\sum _{j=0}^{\infty }{\mathcal{C}}_j^{S\left(1\right)}\mathrm{i}{k}_j{e}^{\mathrm{i}{k}_{j}L}{b}_{ij}{\delta }_{ij}-\frac{{\mathcal{C}}_j^{S\left(2\right)}j\pi }{c+d-h}sinh\left(\frac{j\pi x}{c+d-h}\right)a_{ij}=\frac{Ag{k}_0{b}_{ii}{\delta }_{i0}{e}^{-\mathrm{i}{k}_{0}L}}{2\omega cosh{k}_{0}h}\hfill \\ \sum _{j=0}^{\infty }{\mathcal{C}}_j^{S\left(1\right)}{e}^{\mathrm{i}{k}_{j}L}{a}_{ji}-\sum _{j=1}^{\infty }{\mathcal{C}}_j^{S\left(2\right)}\frac{h-d-c}{2}cosh\left(\frac{j\pi x}{c+d-h}\right){\delta }_{ij}-{\mathcal{C}}_0^{S\left(2\right)}\left(h-d-c\right)=\frac{\mathrm{i}Ag{a}_{0i}}{2\omega cosh{k}_{0}h}\hfill \\ i=0,\dots ,\infty \hfill \end{array}\right.,$$

(48)

$$\left\{\begin{array}{c}\sum _{j=0}^{\infty }{\mathcal{C}}_j^{A\left(1\right)}\mathrm{i}{k}_j{e}^{\mathrm{i}{k}_{j}L}{b}_{ij}{\delta }_{ij}-\sum _{j=1}^{\infty }\frac{{\mathcal{C}}_j^{A\left(2\right)}j\pi }{c+d-h}cosh\left(\frac{j\pi x}{c+d-h}\right)a_{ij}-{\mathcal{C}}_0^{A\left(2\right)}{a}_{i0}=\frac{Ag{k}_0{b}_{ii}{\delta }_{i0}{e}^{-\mathrm{i}{k}_{0}L}}{2\omega cosh{k}_{0}h}\hfill \\ \sum _{j=0}^{\infty }{\mathcal{C}}_j^{A\left(1\right)}{e}^{\mathrm{i}{k}_{j}L}{a}_{ji}-\sum _{j=1}^{\infty }{\mathcal{C}}_j^{A\left(2\right)}\frac{h-d-c}{2}sinh\left(\frac{j\pi x}{c+d-h}\right){\delta }_{ij}-{\mathcal{C}}_0^{A\left(2\right)}L\left(h-d-c\right){\delta }_{i0}=\frac{\mathrm{i}Ag{a}_{0i}}{2\omega cosh{k}_{0}h}\hfill \\ i=0,\dots ,\infty \hfill \end{array}\right.,$$

(49)

where ${\delta }_{ij}$ is the Kronecker delta, while

$$a_{ij}=\frac{k_i\left(c+d-h\right)\left[{\left(-1\right)}^{1+j}sinh{k}_{i}c+sinh{k}_i\left(h-d\right)\right]}{{\left(c+d-h\right)}^2{k}_i^2+j^2{\pi }^2},$$

(50)

$$b_{ij}=\left\{\begin{array}{c}\frac{1}{2}\left(h+\frac{sinh2k_{i}h}{2k_i}\right),i=j\hfill \\ 0,i\ne j\hfill \end{array}\right..$$

(51)

Systems (48) and (49) are solved by truncating and numerically solving an $N\times N$ system of equations. The singularity at the bottom edges of the WEC is weaker than that of objects characterised by sharp corners, thus the numerical convergence is fast [16,17,18]. In Section 2.3, we check the numerical computations through theoretical integral relations, such as the Haskind–Hanaoka formula.

2.2. Radiation Potential Solution

Since the problem is linear, the solution in ${\mathsf{\Omega }}_1$ can be written as

$${\varphi }_{lR}^{S,A\left(1\right)}=\sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{S,A\left(1\right)}{X}_j^{\left(1\right)}{Z}_j^{\left(1\right)}$$

(52)

where the ${\mathcal{D}}_{lj}^{S,A\left(1\right)}$ are unknown complex constants, while the eigenfunctions $X_j^{\left(1\right)}$ and $Z_j^{\left(1\right)}$ are given by (39) and (40).

The solution in ${\mathsf{\Omega }}_2$ is given by the homogeneous part and a particular solution that accounts for the plate vibration in $z=-d$

$${\varphi }_{lR}^{S,A\left(2\right)}=\sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{S,A\left(2\right)}{X}_j^{\left(2\right)}{Z}_j^{\left(2\right)}+{\tilde{\varphi }}_{lR}^{S,A},$$

(53)

where ${\mathcal{D}}_{lj}^{S,A\left(2\right)}$ are unknown complex constants, while $X_j^{\left(2\right)}$ and $Z_j^{\left(2\right)}$ are expressed by (43) and (45). The particular solution for the rigid heave mode reads

$${\tilde{\varphi }}_{0R}^S=\frac{-\mathrm{i}\omega }{2\left(h-c-d\right)}\left[z^2-x^2+2z\left(h-c\right)\right],$$

(54)

while the particular solution for the pitching mode is given by [14]

$${\tilde{\varphi }}_{0R}^A=\frac{\mathrm{i}\omega x}{6\left(h-c-d\right)}\left[x^2-3z^2-6z\left(h-c\right)\right].$$

(55)

For the symmetric and antisymmetric lth bending mode (15)–(17), the particular solutions are

$${\tilde{\varphi }}_{lR}^S=\frac{-\mathrm{i}\omega L}{{\mu }_l^S}\left[\frac{cos\frac{{\mu }_l^{S}x}{L}cosh\frac{{\mu }_l^S\left(z+h-c\right)}{L}}{cos{\mu }_l^{S}sinh\frac{{\mu }_l^S\left(h-c-d\right)}{L}}-\frac{cosh\frac{{\mu }_l^{S}x}{L}cos\frac{{\mu }_l^S\left(z+h-c\right)}{L}}{cosh{\mu }_l^{S}sin\frac{{\mu }_l^S\left(h-c-d\right)}{L}}\right],$$

(56)

$${\tilde{\varphi }}_{lR}^A=\frac{-\mathrm{i}\omega L}{{\mu }_l^A}\left[\frac{sin\frac{{\mu }_l^{A}x}{L}cosh\frac{{\mu }_l^A\left(z+h-c\right)}{L}}{sin{\mu }_l^{A}sinh\frac{{\mu }_l^A\left(h-c-d\right)}{L}}-\frac{sinh\frac{{\mu }_l^{A}x}{L}cos\frac{{\mu }_l^A\left(z+h-c\right)}{L}}{sinh{\mu }_l^{A}sin\frac{{\mu }_l^A\left(h-c-d\right)}{L}}\right].$$

(57)

As in the previous subsection, continuity of pressure and velocity (32) and (33) gives

$$\sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{S,A\left(1\right)}{\left.X_j^{S,A{\left(1\right)}^{\prime }}\right|}_{x=L}{Z}_j^{\left(1\right)}-{\mathcal{D}}_{lj}^{S,A\left(2\right)}{\left.X_j^{S,A{\left(2\right)}^{\prime }}\right|}_{x=L}{Z}_j^{\left(2\right)}={\left.\frac{\partial {\tilde{\varphi }}_{lR}^{S,A}}{\partial x}\right|}_{x=L},$$

(58)

$$\sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{S,A\left(1\right)}{X}_j^{S,A\left(1\right)}\left(L\right)Z_j^{\left(1\right)}-{\mathcal{D}}_{lj}^{S,A\left(2\right)}{X}_j^{S,A\left(2\right)}\left(L\right)Z_j^{\left(2\right)}={\left.{\tilde{\varphi }}_{lR}^{S,A}\right|}_{x=L}.$$

(59)

As in the case of the diffraction potential, by multiplying each side of (58) and (59) by $Z_i^{\left(1\right)}$ and $Z_i^{\left(2\right)}$, respectively, and integrating over the relevant intervals, $z\in \left[-h,0\right]$ and $z\in \left[-h+c,-d\right]$, yields the following inhomogeneous linear systems in $D_{lj}^{S,A\left(1\right)}$ and $D_{lj}^{S,A\left(2\right)}$

$$\left\{\begin{array}{c}\sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{S\left(1\right)}\mathrm{i}{k}_j{e}^{\mathrm{i}{k}_{j}L}{b}_{ij}{\delta }_{ij}-\frac{{\mathcal{D}}_{lj}^{S\left(2\right)}j\pi }{c+d-h}sinh\left(\frac{j\pi x}{c+d-h}\right)a_{ij}=c_{il}^{S\left(1\right)}\hfill \\ \sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{S\left(1\right)}{e}^{\mathrm{i}{k}_{j}L}{a}_{ji}-\sum _{j=1}^{\infty }{\mathcal{D}}_{lj}^{S\left(2\right)}\frac{h-d-c}{2}cosh\left(\frac{j\pi x}{c+d-h}\right){\delta }_{ij}-{\mathcal{D}}_{l0}^{S\left(2\right)}\left(h-d-c\right)=c_{il}^{S\left(2\right)}\hfill \\ i=0,\dots ,\infty \hfill \end{array}\right.,$$

(60)

$$\left\{\begin{array}{c}\sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{A\left(1\right)}\mathrm{i}{k}_j{e}^{\mathrm{i}{k}_{j}L}{b}_{ij}{\delta }_{ij}-\sum _{j=1}^{\infty }\frac{{\mathcal{D}}_{lj}^{A\left(2\right)}j\pi }{c+d-h}cosh\left(\frac{j\pi x}{c+d-h}\right)a_{ij}-{\mathcal{D}}_{lj}^{A\left(2\right)}{a}_{i0}=c_{il}^{A\left(1\right)}\hfill \\ \sum _{j=0}^{\infty }{\mathcal{D}}_{lj}^{A\left(1\right)}{e}^{\mathrm{i}{k}_{j}L}{a}_{ji}-\sum _{j=1}^{\infty }{\mathcal{D}}_{lj}^{A\left(2\right)}\frac{h-d-c}{2}sinh\left(\frac{j\pi x}{c+d-h}\right){\delta }_{ij}-{\mathcal{D}}_{l0}^{A\left(2\right)}L\left(h-d-c\right){\delta }_{i0}=c_{il}^{A\left(2\right)}\hfill \\ i=0,\dots ,\infty \hfill \end{array}\right..$$

(61)

In the latter, $a_{ij}$ and $b_{ij}$ are given by (50) and (51), while the constant terms on the right hand side read

$$c_{il}^{S,A\left(1\right)}={\int }_{-h+c}^{-d}\mathrm{d}z{Z}_i^{\left(1\right)}{\left.\frac{\partial {\tilde{\varphi }}_{lR}^{S,A}}{\partial x}\right|}_{x=L}+{\left.\mathrm{i}\omega \frac{\mathrm{d}{f}_l^{S,A}}{\mathrm{d}x}\right|}_{x=L}{\int }_{-d}^0\mathrm{d}z{Z}_i^{\left(1\right)}\left(z-z_g\right),$$

(62)

$$c_{il}^{S,A\left(2\right)}={\int }_{-h+c}^{-d}\mathrm{d}z{Z}_i^{\left(2\right)}{\left.{\tilde{\varphi }}_{lR}^{S,A}\right|}_{x=L}$$

(63)

As in the previous section, the linear systems (60) and (61) are solved numerically by truncating the series at $n=N$ and $i=N$.

2.3. Structural Response and the Haskind–Hanaoka Formula

The vibration of the floating elastic plate is governed by the following Euler dynamic equation

$$-EI{W}_{xxxx}-\rho {\mathsf{\Phi }}_t-\sum _{m=1}^M\delta \left(x-x_m\right){\nu }_{PTO}{W}_t-\rho gW-{\rho }_p{h}_p{W}_{tt}=0,$$

(64)

where E is the elastic modulus of the plate, I is the area moment of inertia of the plate and $\delta$ denotes the Dirac delta function. The first term in the equation above represents the flexural rigidity, the second term is the dynamic pressure exerted by the diffracted and radiated wave fields, the third term represents the effect of localised forces due to the PTO system, the fourth term is the hydrostatic contribution, while the last term represents the inertia of the plate.

By expanding W through the dry mode decomposition (13) and recalling that for a free-free beam in the absence of applied loads

$$f_{l_{xxxx}}^{S,A}=f_l^{S,A}{\left(\frac{{\mu }_l^{S,A}}{L}\right)}^4,$$

(65)

we obtain

$$\begin{array}{cc}\hfill & \sum _{l=0}^{\infty }\left\{\left[EI{ϵ}_l{\left(\frac{{\mu }_l^S}{L}\right)}^4+\rho g-{\omega }^2{\rho }_p{h}_p-\sum _{m=1}^M\delta \left(x-x_m\right){\nu }_{PTO}\mathrm{i}\omega \right]f_l^S{\xi }_l^S-\mathrm{i}\omega \rho {\left.{\varphi }_{lR}^{S\left(2\right)}\right|}_{z=-d}{\xi }_l^S\right.\hfill \\ \hfill & +\left.\left[EI{ϵ}_l{\left(\frac{{\mu }_l^A}{L}\right)}^4+\rho g-{\omega }^2{\rho }_p{h}_p-\sum _{m=1}^M\delta \left(x-x_m\right){\nu }_{PTO}\mathrm{i}\omega \right]f_l^A{\xi }_l^A-\mathrm{i}\omega \rho {\left.{\varphi }_{lR}^{A\left(2\right)}\right|}_{z=-d}{\xi }_l^A\right\}\hfill \\ \hfill & =\mathrm{i}\omega \rho {\left.\left({\varphi }_D^{S\left(2\right)}+{\varphi }_D^{A\left(2\right)}\right)\right|}_{z=-d},\hfill \end{array}$$

(66)

where, ${ϵ}_l$ is defined by

$${ϵ}_0=0,{ϵ}_l=1,l=1,\dots ,\infty .$$

(67)

The complex coefficients ${\xi }_l^S$ and ${\xi }_l^A$ can be found by multiplying the latter equation by $f_i^S$ and $f_i^S$, respectively, and then integrating along the total length of the plate $x\in \left[-L,L\right]$. Truncating the series at $l=J$, we obtain a $\left(2J+1\right)\times \left(2J+1\right)$ non-homogeneous system in ${\xi }_l^S$ and ${\xi }_l^A$ that can be written in compact form

$$\begin{array}{cc}\hfill & \sum _{l=0}^{\infty }\left\{\left[K_{il}^S-{\omega }^2\left(I_{il}^S+M_{il}^S\right)-\mathrm{i}\omega \left(C_{il}^S-{\nu }_{PTO}{D}_{il}^S\right)\right]{\xi }_l^S\right.\hfill \\ \hfill & \left.+\left[K_{il}^A-{\omega }^2\left(I_{il}^A+M_{il}^A\right)-\mathrm{i}\omega \left(C_{il}^A-{\nu }_{PTO}{D}_{il}^A\right)\right]{\xi }_l^A\right\}=F_i,i=0,\dots ,\infty ,\hfill \end{array}$$

(68)

where $K_{il}$, $I_{il}$, $M_{il}$, $C_{il}$ and $D_{il}$ are the generalised stiffness matrix, mass matrix, added mas matrix, radiation damping matrix and PTO damping matrix, while the term at the right-hand side represents the exciting force vector. Their respective expressions are given in Appendix A.

The structure of the latter equation suggests that the floating plate behaves as a linear forced harmonic oscillator. The natural modes of the WEC are then evaluated from the following homogeneous unforced and undamped system

$$\begin{array}{cc}\hfill & \sum _{l=0}^{\infty }\left\{\left[K_{il}^S-{\omega }^2\left(I_{il}^S+M_{il}^S\right)\right]{\xi }_l^S+\left[K_{il}^A-{\omega }^2\left(I_{il}^A+M_{il}^A\right)\right]{\xi }_l^A\right\}=0,i=0,\dots ,\infty .\hfill \end{array}$$

(69)

By equating to zero the determinant of the coefficient matrix, we obtain the eigenfrequencies ${\omega }_i$ and the respective modal forms.

Now, we apply the Haskind–Hanaoka formula valid for two-dimensional domains to check the numerical computations of the diffraction and radiation velocity potentials. Its expression reads [15]

$$F_i^{S,A}=-2\mathrm{i}\rho gA{\mathcal{A}}_i^{S,A^{-}}{C}_g{\omega }^{-1},$$

(70)

where the term on the left hand side represents the exciting force given by (A21) and (A22), i.e.,

$$F_i^{S,A}=\mathrm{i}\omega \rho {\int }_{-L}^L\mathrm{d}x{f}_i^{S,A}{\left.\left({\varphi }_D^{S\left(2\right)}+{\varphi }_D^{A\left(2\right)}\right)\right|}_{z=-d},$$

(71)

$C_g$ is the group velocity,

$$C_g=\frac{\omega }{2k_0}\left(1+\frac{2k_{0}h}{sinh2k_{0}h}\right),$$

(72)

while ${\mathcal{A}}_i^{S,A^{-}}$ is the amplitude of the radiated waves at large distance from the plate, for unit modal amplitude and in the direction opposite the incident waves

$${\mathcal{A}}_i^{S,A^{-}}=\frac{\mathrm{i}\omega cosh{k}_{0}h{D}_{i0}^{S,A\left(1\right)}}{g},$$

(73)

in which $D_{i0}^{S,A\left(1\right)}$ is the first complex coefficient of the radiation potential in ${\mathsf{\Omega }}_1$ (52). Substitution of (73) in (70) gives

$$F_i^{S,A}=2\rho A{C}_g{D}_{i0}^{S,A\left(1\right)}cosh{k}_{0}h.$$

(74)

The latter relates the diffraction to the radiation potential, and it is used to perform numerical check evaluations.

2.4. Wave Power Extraction and Theoretical Maximum Efficiency

Once the complex coefficients ${\xi }_l^S$ and ${\xi }_l^A$ are determined, the average power absorbed over a wave period $T=2\pi /\omega$ by the plate in monochromatic waves can be calculated as

$$P=\frac{1}{T}{\int }_0^T{\nu }_{PTO}\sum _{m=1}^M{\left(\frac{\mathrm{d}W\left(x_m\right)}{\mathrm{d}t}\right)}^2\mathrm{d}t=\frac{1}{2}{\nu }_{PTO}{\omega }^2\sum _{m=1}^M{\left|w\left(x_m\right)\right|}^2,$$

(75)

which, after the substitution of (13), becomes

$$P=\frac{1}{2}{\nu }_{PTO}{\omega }^2\sum _{m=1}^M{\left|\sum _{l=0}^{\infty }{\xi }_l^S{f}_l^S\left(x_m\right)+{\xi }_l^A{f}_l^A\left(x_m\right)\right|}^2.$$

(76)

Then, we define the capture factor as the ratio between the power output P and the incident wave energy flux per unit width

$$C_F=\frac{P}{{CE}_g},$$

(77)

where

$$E=\frac{1}{2}\rho g{A}^2,$$

(78)

is the total energy and $C_g$ is the group speed.

We now turn to the evaluation of the theoretical maximum capture factor. Using the radiated wave amplitudes (73), the most general expression of the capture factor (77) for two-dimensional flexible floaters becomes

$$C_F=-\mathrm{Re}\left\{2\sum _l{\mathcal{A}}_l^{-}{\xi }_l^{*}+\sum _l\sum _i{\xi }_l{\xi }_i^{*}\left({\mathcal{A}}_l^{-}{\mathcal{A}}_i^{-*}+{\mathcal{A}}_l^{+}{\mathcal{A}}_i^{+*}\right)\right\},$$

(79)

where ${\left(·\right)}^{*}$ denotes the complex conjugate of $\left(·\right)$, while the superscripts $A,S$ denoting symmetric and antisymmetric components are omitted for brevity. If there is one degree of freedom, the latter becomes

$$C_F=-2\mathrm{Re}\left\{{\mathcal{A}}_l^{-}{\xi }_l^{*}\right\}-2{\left|{\xi }_l{\mathcal{A}}_l^{-}\right|}^2,$$

(80)

which is maximised when $\mathrm{Re}\left\{{\mathcal{A}}_l^{-}{\xi }_l^{*}\right\}=-\left|{\mathcal{A}}_l^{-}{\xi }_l^{*}\right|$. The maximum value of the capture factor occurs when $\mathrm{d}{C}_F/\mathrm{d}{\xi }_l=0$, i.e., for

$$\left|{\xi }_l\right|=\frac{1}{2\left|{\mathcal{A}}_l^{-}\right|},$$

(81)

and is equal to $C_F=0.5$. This result can be derived in a different way directly from the equation of motion of a two-dimensional rigid absorber properly constrained [15,19].

We remark that the theoretical maximum of the capture factor for a two-dimensional WEC cannot be larger than 1 because of conservation of energy. The plate considered in this work is elastic and characterised by two rigid modes (heave and pitch) and an infinite set of bending modes, thus $C_F$ can be maximised several times within the range of frequencies of interest.

For example, let us consider two modes that dominate the dynamics with respect to the others, one symmetric and the other one antisymmetric. Recalling that wave energy extraction is optimised when the total scattered and radiated waves are maximised in the direction opposite to the incident wave field, we assume ${\mathcal{A}}_l^{S-}={\mathcal{A}}_l^{A-}$, ${\mathcal{A}}_l^{S+}=-{\mathcal{A}}_l^{A+}$, i.e., the radiated wave amplitude of each mode in the direction opposite to the incident waves is the same. The corresponding capture factor becomes

$$C_F=-2\mathrm{Re}\left\{{\mathcal{A}}_l^{S-}\left({\xi }_l^{S*}+{\xi }_l^{A*}\right)\right\}-2{\left|{\xi }_l^S{\mathcal{A}}_l^{S-}\right|}^2-2{\left|{\xi }_l^A{\mathcal{A}}_l^{S-}\right|}^2.$$

(82)

The capture factor is maximised when the first term on the right-hand side is real and negative and when the modal coefficients ${\xi }_l^{S,A}$ satisfy the following condition

$$\left|{\xi }_l^{S,A}\right|=\frac{1}{2\left|{\mathcal{A}}_l^{S-}\right|}.$$

(83)

Substitution of (83) into (82) gives $C_F=1$, a value independent of the WEC size. This result has been obtained from the simplified assumption of WEC motion dominated by two modes, while the flexible floater considered in this work includes rigid and bending elastic modes as well. This aspect potentially implies multiple optimisation and consequent larger efficiency with respect to standard rigid devices.

3. Results and Discussion

In this section, we investigate the effects of the plate geometry, ridge height c and PTO distribution on the hydrodynamic behaviour and efficiency of the system. We choose the following parameters: $A=1$ m, $h=5$ m, $\rho =1000$ kg m${}^{-3}$, $EI=6.9\times {10}^4$ kg m${}^3$s${}^{-2}$, and $L=10$ m. Since in the expressions for the velocity potentials and bending modes there are infinite terms, the summations must be truncated up to a limiting value in the computations. Here, we use $j=10$ and consider the first 5 dry modes to reach good accuracy [13]. Calculations were carried out using the MATLAB software.

3.1. Effects of the PTO

Here, we investigate the effects of the PTO coefficient ${\nu }_{PTO}$ and PTO distribution on the power extraction efficiency. For the sake of brevity, let us assume $c=0$ m.

Figure 4 shows the surface plot of the capture factor versus the incident wave frequency and the PTO coefficient for different PTO distributions. Figure 4a refers to the case of two PTO systems located at the ends of the plate $x_i=\pm L$, while Figure 4b shows the case of 5 PTOs equally spaced located at $x_i=\pm L,\pm L/2,0$. When the number of the PTOs increases, the bandwidth of the capture factor increases and the system becomes more efficient. Several peaks with value $C_F\sim 1$ are shown, therefore the theoretical maximum of a two-dimensional WEC in a channel can be almost reached. This result confirms the theoretical predictions obtained in Section 2.4. Note also that the maxima are located in correspondence to the first eigenfrequencies of the system, i.e., ${\omega }_i=0.99,1.55,2.24,2.77,3.66,5.06$ Rad s${}^{-1}$. As in the case of oscillating wave surge converters and oscillating water columns [17,18,20,21,22,23,24], resonance of natural modes is beneficial in terms of power extraction efficiency.

3.2. Effects of the Ridge Height

Now we analyse the effects of the ridge height c on the capture factor. Let us fix the PTO distribution to be equally spaced every 5 m and consider the same floater geometry analysed in the previous section. By comparing Figure 4b (case without ridge) with Figure 5a,b, we note that when the ridge height increases, the overall efficiency at small frequencies decreases, while a narrow peak appears around $\omega =2.2$ Rad s${}^{-1}$. This behaviour is mainly due to the reflection of the incident waves with increasing c and the shifting of the eigenfrequencies towards smaller values. In any case, this analysis shows that if a bottom structure with a large height c is needed, the floater geometry and the PTO distribution can be still optimised to maximise $C_F$ up to values around one.

3.3. Effects of the Plate Stiffness

In order to evaluate the effects of the plate stiffness on the generated power, a parametric analysis is performed for a softened plate characterised by a smaller value of the stiffness factor $EI=6.9\times {10}^3$ kg m${}^3$s${}^{-2}$ and a rigid plate. Figure 6a shows that when the flexural rigidity of the plate decreases, the efficiency of the system can increase. This is due to the shifting of the bending mode eigenfrequencies towards smaller frequencies and the fact that these frequencies tend to be much close to each other. Differently from the case shown in Figure 4b, there are now four peaks in which the capture factor $C_F$ is close to one.

When the plate is rigid, or characterised by very large stiffness, there are no contributions from the bending modes and the dynamics is governed by pitching and heaving only. Figure 6b shows that there is one maximum around $\omega =2.2$ Rad s${}^{-1}$ with value $C_f\sim 0.8$. This frequency falls within the eigenfrequencies of the heave mode ${\omega }_i=1.07$ Rad s${}^{-1}$ and pitching mode ${\omega }_i=1.67$ Rad s${}^{-1}$, respectively. The overall efficiency is clearly smaller with respect to the cases shown so far because we reduced the number of eigenfrequencies, the resonances of the natural modes and the possible modal optimisations. This highlights the beneficial effects of the bending elastic modes on the power extraction efficiency.

4. Comparison with Preliminary Demonstrator Data

The flexible floater demonstrator is made by two layers of Polymax SILO-CELL silicone sponge sheet, each of dimensions 2 m (length) × $0.2$ m (width) × $0.01$ m (thickness), see Figure 7.

The floater is attached to a laboratory model (scale 1:35) of the multi-pump, multi-piston power take off system (MP${}^2$PTO) designed by the University of Groningen and Ocean Grazer, shown in Figure 8.

The MP${}^2$PTO system is installed inside a wave tank, which is $1.20$ m high, $0.77$ m wide and 10 m long. Two transparent lateral walls are also installed inside the wave tank, along the direction of the incident waves, restricting the width to $0.2$ m. This effectively creates a channel inside the tank, and the flexible floater is then installed inside this channel. The channel width matches the floater width ($0.2$ m), hence the dynamics inside the channel are two-dimensional. The water level in the tank is $h=0.9$ m. The waves are generated by a flap paddle driven by a rotating-arm engine, located at one side of the tank. The engine frequency and the length of the rotating arm can be set to a maximum of 60 Hz and $0.25$ m, respectively. At the other side of the tank, a plane beach induces wave breaking, thus dissipating wave energy and reducing reflection. To calibrate the wave maker, digital particle image velocimetry (DPIV) measurements were taken using a high-definition camera. The camera recorded the motion of tracer particles (polyamide particles) seeded in the tank, which was illuminated by means of a laser sheet. For details on the measurement procedure and associated errors, we refer the interested reader to Refs. [25,26] and references cited therein. The system of pistons and cables was calibrated using load cells, as described in Ref. [27].

4.1. Description of the Demonstrator Device

The flexible floater is connected to the MP${}^2$PTO system in the tank via high-performance polyethylene cables, which are in turn connected to pistons. Each piston is located inside a $0.057$ m wide cylinder. As the floater deforms under the action of the incident waves, it transmits its motion to the pistons via the cables. In turn, the pistons pump water inside the cylinders, as shown in Figure 9.

The power extracted by each PTO element can be approximated by calculating the work done by the piston against the force of gravity to lift the water column of weight $\rho gAH$, where A is the cross-sectional area of the cylinder and H is the maximum hydraulic head over a cycle. Thus, the extracted power is given by

$$P=\rho gA{H}^2/\tau ,$$

where $\tau$ is the duration of the upstroke motion of the piston.

Four sets of demonstrations are conducted. In two sets (S1 and S2), the floater is a continuous flexible silicone sheet of 2 m length, whereas in two sets (S3 and S4), the floater is cut into 10 pieces, each being $0.2$ m long and connected to a single piston. For each setup, the maximum hydraulic head in each cylinder is measured eight times, and the average of such measurements is then taken. The engine and wave parameters for the four sets of tests are reported in

Table 1. Parameters for the four sets of experiments on the flexible floater device.

Set	Floater	Wave Height	Wave Period
S1	Continuous	$0.06$ m	$1.36$ s
S2	Continuous	$0.08$ m	$1.62$ s
S3	Discontinuous	$0.06$ m	$1.36$ s
S4	Discontinuous	$0.08$ m	$1.62$ s

. We use Froude’s law to scale between test and model results. The Froude number is

$$Fr=\frac{U^2}{gD},$$

(84)

where U is the fluid speed and D any linear dimension of the system. Using Equation (84) and a 1:35 scale ratio, the wave periods $1.36$ s and $1.62$ s scale up to 8 and $9.5$ s, respectively, whereas the wave heights $0.06$ m and $0.08$ m scale up to $2.1$ m and $2.8$ m, respectively. These values are consistent with an oceanic wave climate, such as the North Atlantic Ocean [28].

4.2. Demonstrator Results

Figure 10 shows the average power at each of the 10 pistons in the MP${}^2$PTO laboratory model, for the four sets of tests detailed in Table 1. Piston 1 is the closest to the generator, whereas piston 10 is the closest to the absorbing beach, on the other side of the tank. Note that, in all configurations, the power absorbed by the device first increases moving towards the centre of the floater, and then decreases moving further on. Therefore, in all cases, the maximum power output is achieved towards the middle of the floater. Note also that the power has a second maximum at the tail of the floater (piston 10), which results from end reflections by the inclined beach. This is especially visible in the more energetic sea state (configurations S2 and S4), whereas it is almost negligible in the milder sea state (configurations S1 and S3).

As expected, the configurations relevant to the more energetic seas, S2 and S4, are also those with the largest power. Figure 10 also shows that the maximum power attained by the continuous floater (configurations S1 and S2) is similar to that achieved by the discontinuous floater (configurations S3 and S4). Interestingly, in the less energetic state, the continuous floater (configuration S1) performs slightly better than the discontinuous one (configuration S3) at almost all pistons. In the more energetic sea state, the continuous floater (configuration S2) performs better than the discontinuous one (configuration S4) towards the front (pistons 1,2 and 3), and worse towards the middle (pistons 4–9). Since there is no substantial advantage in choosing the continuous configuration over the discontinuous one (at least in the sea states analysed here), our results suggest that in practice one can opt for either configuration based on other criteria, such as installation and maintenance costs. In this sense, the discontinuous configuration appears more advantageous to reduce maintenance costs, for example due to the possibility of replacing only one element, instead of the whole floater, in case of localised damage. The positive practical implication of our preliminary findings prompts the need to scale-up the demonstrator device to a full-fledged experimental model. This is the object of our ongoing research effort.

4.3. Comparison with Mathematical Model

In this section, we show a preliminary comparison between the results of the experimental and mathematical models. We remark that the PTO system is modelled as a linear damper in the mathematical model, whereas it is nonlinear in the experimental model. For the sake of comparison, we selected a PTO coefficient for the model which generates the same total power output as that in the demonstrator device. This allows some quantitative comparisons and discussions.

Figure 11 shows the behaviour of the power output for each piston, in both the mathematical model and the tests, for the same configuration as Figure 10. The overall behaviour is captured well and, in general, the comparison is satisfactory. Both models show that the maximum power output is achieved by those pistons located towards the front of the device. Some differences (especially in the maximum power output) still remain, and these are likely due to the use of a linear PTO in the mathematical model.

We remark that, to date, very few studies have investigated the non-linear dynamics of wave-plate systems analytically, see for example Refs. [29,30]. However, these dealt with the interaction of waves with ice sheets. On the contrary, to the best of our knowledge, no application of nonlinear theories to wave power extraction from flexible plates has been made so far. This highlights the need for developing higher-order mathematical models to achieve a more accurate description of the power extraction dynamics.

5. Power Extraction in Irregular Waves

In this section, we investigate the effect of irregular sea waves on the floating plate dynamics and power extraction efficiency. Let us assume the following JONSWAP spectrum function [31]

$$S_{\zeta }(\omega ,{\omega }_p)=\frac{\alpha H_s^2}{\omega }{\left(\frac{{\omega }_p}{\omega }\right)}^4\exp \left[-1.25{\left(\frac{{\omega }_p}{\omega }\right)}^4\right]{\gamma }^{\exp \left[-{(\omega /{\omega }_p-1)}^2/\left(2{\sigma }^2\right)\right]},$$

(85)

in which $H_s$ is the significant wave height, ${\omega }_p$ denotes the peak frequency and

$$\alpha =\frac{0.0624(1.094-0.01915ln\gamma )}{0.23+0.0336\gamma -0.185{(1.9+\gamma )}^{-1}},\sigma =\left\{\begin{array}{c}\hfill 0.07:\omega \le {\omega }_p\\ \hfill 0.09:\omega >{\omega }_p\end{array}\right.,\gamma =3.3.$$

(86)

Since the problem is linear, the oscillation of the flexible plate can be written as

$$W\left(t,{\omega }_p,x,{\nu }_{PTO}\right)=\sum _{n=1}^{\infty }\sqrt{2S_{\zeta }\left({\omega }_n\right)\mathsf{\Delta }\omega }\mathrm{RAO}\left({\omega }_n,x,{\nu }_{PTO}\right)cos\left({\omega }_{n}t+{\delta }_n\right),$$

(87)

where ${\omega }_n$ is the nth component of the discretised spectrum, $\mathsf{\Delta }\omega$ is the frequency interval, ${\delta }_n$ is a random phase related to ${\omega }_n$, whereas RAO is the response amplitude operator for the plate, i.e.,

$$\mathrm{RAO}\left({\omega }_n,x,{\nu }_{PTO}\right)=\left|w\left({\omega }_n,x,{\nu }_{PTO}\right)\right|.$$

(88)

Then, the instantaneous generated power by the entire system is

$$\begin{array}{cc}\hfill P_s\left(t,{\omega }_p,{\nu }_{PTO}\right)& ={\nu }_{PTO}{\left[\sum _{m=1}^M\sum _{n=1}^{\infty }\sqrt{2S_{\zeta }\left({\omega }_n,{\omega }_p\right)\mathsf{\Delta }\omega }\mathrm{RAO}\left({\omega }_n,x_m,{\nu }_{PTO}\right){\omega }_{n}sin\left({\omega }_{n}t+{\delta }_n\right)\right]}^2.\hfill \end{array}$$

(89)

From the previous expression, we obtain the averaged generated power [17,20,22]

$${\overline{P}}_s\left({\omega }_p,{\nu }_{PTO}\right)=\underset{\tau \to \infty }{lim}\frac{1}{\tau }{\int }_0^{\tau }{P}_s\mathrm{d}t={\nu }_{PTO}\sum _{m=1}^M\sum _{n=1}^{\infty }{S}_{\zeta }\left({\omega }_n,{\omega }_p\right)\mathsf{\Delta }\omega {\mathrm{RAO}}^2\left({\omega }_n,x_m,{\nu }_{PTO}\right){\omega }_n^2,$$

(90)

whose expression in the limit $\mathsf{\Delta }\omega \to 0$ becomes

$${\overline{P}}_s\left({\omega }_p,{\nu }_{PTO}\right)={\nu }_{PTO}\sum _{m=1}^M{\int }_0^{\infty }{S}_{\zeta }\left(\omega ,{\omega }_p\right){\mathrm{RAO}}^2\left(\omega ,x_m,{\nu }_{PTO}\right){\omega }^2\mathrm{d}\omega .$$

(91)

Defining $P_{\zeta }$ as the total incident wave power per unit crest width

$$P_{\zeta }\left({\omega }_p\right)={\int }_0^{\infty }\rho g{C}_g\left(\omega \right)S_{\zeta }\left(\omega ,{\omega }_p\right)\mathrm{d}\omega ,$$

(92)

the capture factor in irregular seas $C_{F\zeta }$ can then be written as

$$C_{F\zeta }\left({\omega }_p,{\nu }_{PTO}\right)=\frac{{\overline{P}}_s}{P_{\zeta }}.$$

(93)

The latter expression gives the capture factor for any sea state characterised by significant wave $H_s$, peak frequency ${\omega }_p$ and PTO coefficient ${\nu }_{PTO}$.

Let us investigate the capture factor in irregular waves of the same plates analysed in Section 3.3. Figure 12 shows the behaviour of $C_{F\zeta }$ versus peak frequency ${\omega }_p$ and PTO-coefficient ${\nu }_{PTO}$ for the softened flexible plate and the rigid plate. The softened plate has stiffness factor $EI=6.9\times {10}^3$ kg m${}^3$s${}^{-2}$, while the rigid plate has $EI\to \infty$.

As in the case of monochromatic incident waves, the flexible plate results in being more efficient than the rigid plate and can be optimised for several values of ${\omega }_p$ and ${\nu }_{PTO}$. Indeed, the rigid plate shows a single peak, while the flexible plate shows three maxima and a much larger bandwidth.

In addition, by comparing Figure 6 and Figure 12, we note that the maxima are reduced with respect to the case of monochromatic waves, whereas the system can be more efficient outside the resonant frequencies. This is mainly due to the coupling between the broadband incident waves and the eigenfrequencies of the system. Similar results were already obtained in the context of oscillating wave surge converters and oscillating water columns [17,20,22].

6. Conclusions

This paper has presented novel results on analytical and demonstrator models of wave energy converters made by flexible floating plates. The mathematical model is based on a linearised potential-flow theory, whereby the method of dry modes is combined to matched eigenfunction expansions, in order to solve the hydrodynamics of the converter. The main results of the analytical model are:

• The effect of the plate elasticity is to increase the number of the resonant frequencies with respect to a rigid plate, while wave power extraction and the bandwidth of the capture factor become larger. The same result has been obtained both in monochromatic and irregular waves.
• The PTO distribution plays a significant role, and it is seen that, by increasing the number of PTO devices, modal optimisation occurs and the overall efficiency of the system improves.
• We also investigated the effect of the ridge height below the plate. Analytical results showed that if a bottom structure is needed, the floater can be properly designed to maximise power extraction, despite reduced incident wave transmission. This aspect has potentially strong implications for the design of nearshore structures for coastal protection.
• We analysed the plate response to irregular waves described by a JONSWAP spectrum. We showed that the presence of a broad wave frequency range reduces the maximum resonant peaks of the system. However, away from resonance, the efficiency can be larger than that of the monochromatic case and the benefit of irregular waves is significant.

A preliminary demonstrator model was also realised, by connecting a flexible Polymax SILO-CELL silicone sheet to the multi-piston MP${}^2$PTO system developed at the Faculty of Science and Engineering, University of Groningen. A comparison between the mathematical model results and the demonstrator data was encouraging. It is to be noted that the demonstrator results were characterised by abnormal power output values for the last piston in the PTO system. This was due to the inclined beach at the end of the tank not being able to suppress end reflections. Further experiments will need to be undertaken in a larger wave tank, with a more sophisticated absorption system, to confirm the trend identified in this work. Given the promising efficiency levels shown by the mathematical model, our results highlight the need to scale-up experimental investigations on flexible wave energy converters, which are still a small minority, compared to those on rigid converters.