Interaction of surface waves with an actuated submerged flexible plate: Optimization for wave energy extraction

Journal of Fluids and Structures 81 (2018) 673–692

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Interaction of surface waves with an actuated submerged flexible plate: Optimization for wave energy extraction N. Desmars a,b , J. Tchoufag a , D. Younesian a,c , M.-R. Alam a, * a b c

Department of Mechanical Engineering, University of California, Berkeley, CA 94720, USA LHEEA Lab, UMR CNRS 6598, École Centrale de Nantes, Nantes, 44321, France School of Engineering, Iran University of Science and Technology, Teheran, 16846-13114, Iran

article

info

Article history: Received 8 October 2017 Received in revised form 25 May 2018 Accepted 31 May 2018

Keywords: Wave energy conversion Flexible plate Normal mode decomposition Efficiency optimization Boundary element method

a b s t r a c t We investigate the interaction of linear surface waves with a submerged wave energy converter which consists of a submerged flexible plate actuated by one (or more) units of power take off (PTO), each modeled as a combination of a linear spring and a linear damper. We develop a wave-flexible structure interaction model through modal decomposition of the structure deformations, and use Boundary Element Method to find hydrodynamic coefficients of each deformation mode of the structure. We validate our methodology with existing analytical results of hydro-elasticity. We then perform a case study and obtain the maximum efficiency (or Capture Width Ratio (CWR)) of our wave energy converter when excited by monochromatic waves through a comprehensive parametric study of the device parameters (e.g. rigidity of the plate, and location and characteristics of the PTO units). In the absence of viscous effects, the device can reach an efficiency as high as 80%. We find that the efficiency is more sensitive to the location of the PTOs than to their damping coefficients. For a range of plate length to wavelength ratios (close to unity) and with a single PTO unit, the optimal PTO location is past the middle of the plate (along the incident wave direction). This location is nearly independent of the rigidity of the plate although the resulting CWR depends on the rigidity. When two PTO units are used the optimal configuration of PTOs depends on the plate’s aspect ratio and its placement with respect to the incident wave: the two PTOs may need to be at the same location, lined up along the direction of wave propagation, or placed side-by-side perpendicular to the direction of the wave propagation. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Ocean waves arriving at a shoreline constitute simultaneously a threat of degradation of coastal infrastructures and an opportunity for energy harnessing. The use of flexible membranes as wave barriers for the shore protection has been extensively addressed in previous studies. Early work in attempting to reflect waves mainly dealt with inextensible and vertical flexible membrane, possibly attached to a buoy at its top and hinged to the seafloor (Lee and Chen, 1990; Thompson, 1991; Cho and Kim, 1998; Ul-Hassan et al., 2009; Chakraborty and Mandal, 2014). These studies essentially use the twodimensional linear hydro-elastic theory to formulate the problem. Analytical methods are mainly based on eigenfunction expansion while numerical simulations use Boundary Element Method or Finite Element Method. Supported by experiments, the aforementioned studies showed that a flexible vertical membrane can be an effective passive breakwater due to its

*

Corresponding author. E-mail address: [email protected] (M.-R. Alam).

https://doi.org/10.1016/j.jfluidstructs.2018.05.016 0889-9746/© 2018 Elsevier Ltd. All rights reserved.

674

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

ability to resonate and reflect waves at different frequencies corresponding to its natural modes. Later designs of wave barriers feature porosity on the structures. Therefore, the flexible membrane cannot only scatter and reflect the incoming wave, but also dissipate part of the otherwise transmitted wave energy (Cho and Kim, 2000; Chan and Lee, 2001). While these studies deal with a single structure, the case of multiple vertical flexible membranes moored to the seafloor has also been investigated in the quest for a higher wave blocking efficiency. It was found that dual membrane breakwaters perform significantly (more than twice) better at reflecting waves than a single membrane, for both normal and oblique incident waves (Cho et al., 1998). While this remarkable result was obtained in the asymmetric configuration, i.e when each membrane has a different in-plane tension, the reflection performance is reduced in the symmetric case. Despite their efficiency, such vertical structures in the near shore region of oceans can potentially prevent the circulation of currents in addition to causing animal entanglement, thus influencing the ecosystem. For these reasons, other configurations have been tested and proved to be remarkable wave barriers as well. For instance, through the method of eigenfunction matching, it was shown properly tuned horizontal membranes can generate motion-induced waves that cancel the otherwise transmitted waves (Cho and Kim, 1998). Depending on the membrane parameters, a total reflection configuration can be reached for some incident wave frequencies. Later studies show that adding structural porosity helps reduce waves transmission as well as plate deflection and hydrodynamic loads on the membrane by increasing viscous dissipation (Cho and Kim, 2000; Behera and Sahoo, 2015; Koley and Sahoo, 2016). Therefore, the bulk of studies concerning hydroelastic interaction of deformable structures with waves have regard with shore protection. However, as mentioned earlier, incoming waves also represent a large and reliable source of energy (Cornett, 2008; Mork et al., 2010) that could be harnessed and converted into useful electricity rather than dissipated by perforated structures. To this end, the past decade has witnessed the development of several wave energy converters (WEC) with various operating principles. Despite many WECs made of rigid elements interacting with the incoming waves (Falnes, 2007; Falcao, 2010), the use of flexible bodies for wave energy conversion is not new (Newman, 1979; Farley, 1982). From the start, the motivation for deformable structures was to enhance the performance of WECs by overcoming some intrinsic deficiencies of conventional systems (a single resonance frequency for instance). One of such devices is the ‘‘Wave Carpet’’ (Alam, 2012b), a WEC design inspired by the visco-elastic behavior of muddy seafloors which dissipate the energy of waves arriving a shoreline (Sheremet and Stone, 2003). The ‘‘Wave Carpet’’ device consists in a long submerged flexible and elastic plate equipped with equally spaced power take-off (PTO) units, comprising springs and dampers, to harvest energy. This plate can bend under the action of the wave-induced pressure difference above and below its surface. There are several advantages for this novel design. First, it can be tuned and designed to perform well either for shore protection (wave reflection) or for energy conversion of the transmitted waves. Secondly, in terms of the operational frequency range, flexible systems normally are more appropriate for wider range of frequencies since they contain multiple resonance frequencies. This device can thus be more fitting for ocean waves which are usually polychromatic. Note that other WECs with this particular feature have also been recently developed (Choplain, 2012; Chaplain et al., 2012; Mei, 2014; Babarit et al., 2015). Third, since the Wave Carpet is placed under the water surface, it is more survivable under the action strong waves and during storms. Moreover, the scalability of a concept being another crucial factor in ocean energy devices, the Wave Carpet is modular. Last, in term of visual appearance, it is also unarguably better than the WECs operating at the surface. However, note that in previous analytical studies (Alam, 2012b, a), the carpet was directly mounted on the seafloor and the fluid–solid hydrodynamic interaction was solved for a continuous viscoelastic boundary at the bottom of the sea. In those studies, the carpet flexural stiffness as well as its inertia has been ignored and the WEC is assumed to have an infinite width. Using a finite size Wave Carpet and finite number of localized PTOs in a laboratory wave tank (Börner and Alam, 2015), a real-time hybrid simulation has been performed in order to optimize the power conversion from incoming monochromatic waves. In that work, the design parameters being varied in the optimization procedure were the PTOs damping coefficient and their location along the carpet length. While that study provided optimization parameters for a given setup, it would be very costly if one were to apply such real-time hybrid simulation systematically to various configurations and sea states. The present study aims at providing a reliable wave-structure model that can be systematically used for the optimization of the Wave Carpet in the limit of linear waves. Specifically, we here extend the aforementioned studies on hydro-elastic interaction of the case where the horizontal flexible structure is modified by application of a finite number of localized PTO units. To this end, we work in the frequency domain and make use of the so-called normal mode approach which is widely known in vibration problems and can be extended to study the interaction between waves and deformable bodies (Newman, 1994). Here, we formulate the problem in three-dimension in order to investigate the effects of non-zero angles of the incident waves on the device optimal parameters. The methodology described in this paper can be straightforwardly extended to irregular waves since real seas are known to be well described as superposition of linear monochromatic waves (Falnes, 2004). This paper is organized as follows. We first describe the wave-structure theoretical model and its numerical implementation. In the latter section, we particularly take care of comparing and validating preliminary results with other approaches. Next, we present the results of the modal analysis of the wave-plate system (which acts as a global oscillator at the frequency of the incoming wave) and analyze the influence of the free parameters, being the angle of the incident wave, the number of PTO units, their locations, their damping coefficient, the plate bending modulus, the plate length and aspect ratio, and the in-plane tension on the structure, on the performance of the system.

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

675

Fig. 1. (a) Side and (b) top view of the actuated flexible submerged plate in operating conditions.

2. Problem statement and governing equations We consider the problem of monochromatic waves interacting with a submerged flexible plate connected to one or more PTO units, each composed of a generator modeled as a linear damper (Fig. 1). We consider a Cartesian coordinate system (x, y, z ), with x, y-axes on the mean water surface and z-axis positive upward. A thin flexible plate of dimensions L × W , thickness H, flexural rigidity D and density ρs is placed horizontally at the depth z = −d. Without loss of generality, we assume ( that the ) center of the plate is at x = 0 (see Fig. 1), and that the plate is under the action of a uniform planar tension T⃗ = Tx , Ty , 0 maintained by the simply supported edges at x = ±L/2 and y = ±W /2. The water of constant depth h is assumed to be homogeneous (density ρ ), inviscid, irrotational and incompressible such that potential flow theory applies. Incident waves of amplitude A and wavelength λ = 2π/k arrive at an angle θ with respect to the positive x-axis (Fig. 1b). Each power take off unit ℓ is located at xℓ , yℓ and zℓ = −d and is composed of a generator modeled as a linear dash-pot damper with the damping coefficient cℓ . Later in this work, we shall also consider the case of PTO units composed, in addition to a damper, of a spring of stiffness kℓ . Throughout the paper, we only use dimensionless parameters defined through the following normalization (asterisk ∗ indicates the dimensional variable):

ξ= cℓ =

ξ∗

D=

A cℓ∗

D∗

ρ

kℓ =

√ ρ gh

Tx =

gh4

ρ

gh2

φ∗ φ= √

k∗ℓ

ρg

ρs∗ Ty = ρ = s ρ gh2 √ ρ Ty∗

Tx∗

A gh

ω=ω

∗

h

g

√ t=t

∗

g h

(1)

µ = kh,

with g being the gravity acceleration. Additionally, all parameters having a length dimension (i.e. x, y, z, L, W , H, d, λ) are normalized with the water depth h. Limiting the analysis to the linear free surface waves, governing equations and boundary conditions for the velocity potential φ (x, y, z , t ) read

∆φ = 0

−h < z < 0, (2a) ∂φ ∂ 2φ =− 2 z = 0, (2b) ∂z ∂t ∂φ =0 z = −h, (2c) ∂z { ∂φ ∂ξ (x, y) ∈ [−L/2, L/2] × [−W /2, W /2] , z = {−d − H , −d} , = (2d) (x, y) = (±L/2, ±W /2) , z ∈ [−d − H , −d] , ∂n ∂t √ ∂φ → iµφ r = x2 + y2 → ∞, (2e) ∂r where ξ (x, y, t ) is the deformation of the plate from its mean location. This deformation is the solution of Kirchhoff’s equation, which reads D∇ 4 ξ −

( Tx

∂2 ∂2 + Ty 2 2 ∂x ∂y

)

ξ + ρs H

∂ 2ξ = P (x, y, t ) + fPTO (x, y, t ) , ∂t2

(3)

with the three boundary conditions ξ (x = ±L/2, y = ±W /2, t ) = 0, ∂x2 ξ (x = ±L/2, y, t ) = 0 and ∂y2 ξ (x, y = ±W /2, t ) = 0. For Eq. (3) to hold, the material damping is neglected and the displacement of the plate is assumed to be small compared to its lateral dimensions. The unsteady component of the net distributed pressure on the plate is given by P (x, y, t ) = −

) ∂ ( φ|z =−d − φ|z =−d−H . ∂t

(4)

676

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

We assume that there are L number of PTO units and express their localized force as fPTO (x, y, t ) = −

L ∑

cℓ

ℓ=1

∂ξ (xℓ , yℓ , t ) δ (x − xℓ , y − yℓ ) , ∂t

(5)

where δ denotes the Dirac delta function. At this point, PTO units are pure dash-pot dampers since structural stiffness is already provided by the plate through its flexural rigidity D and in-plane tension T . Later on, we will see how stiffness from PTO units modifies the efficiency predicted by this pure damping assumption. 3. Normal mode analysis In order to solve this linear wave-structure interaction, we use a normal mode approach and look for time-harmonic (i.e. dependent on time through a factor e−iωt with ω the incident wave frequency) solutions to the governing equations (2). Following this, the plate displacement is of the form

ξ (x, y, t ) =

N M ∑ ∑

Cij Xi (x) Yj (y) e−iωt ,

(6)

i=1 j=1

where the shape functions are Xi (x) = cos (iπ x/L) and Yj (y) = cos (jπ y/W ) for odd mode numbers, and Xi (x) = sin (iπ x/L) and Yj (y) = sin (jπ y/W ) for even mode numbers, (i, j) ∈ [1, N] × [1, M] (we use index i to refer to x-modes and index j to refer to y-modes throughout the paper). Therefore, the odd/even functions correspond to the symmetric/antisymmetric eigenmodes of the unforced (i.e. P = fPTO = 0) thin plate with simply supported edges (Rao, 2007). In the fluid domain, we write the velocity potential φ (x, y, z ) as

φ (x, y, z ) = φI + φS +    φD

N M ∑ ∑

Cij φRij ,

(7)

i=1 j=1

where φI (x, y, z ) is the velocity potential of the incident wave, φS (x, y, z ) the velocity potential of the scattered field due to the plate (at rest), and φRij is the velocity potential of the radiated waves generated by the mode ij of the plate. The sum of the incident and scattered potentials forms the potential of the diffracted wave φD (x, y, z ). Incident wave velocity potential for monochromatic waves satisfying the boundary value problem (2) is

φI (x, y, z ) = −

i cosh [µ (z + 1)]

ω

cosh µ

eiµ(x cos θ+y sin θ ) ,

(8)

where the dimensionless wavenumber µ and frequency ω satisfy the dispersion relation

ω2 = µ tanh µ.

(9)

By substituting the potential decomposition (7) into boundary value problem (2) we obtain one diffraction problem and N × M radiation problems with the following no-penetration conditions on the plate:

∂φD = 0, ∂n

∂φRij = −iωXi Yj . ∂n

(10)

The plate equation of motion, Eq. (3), now turns into N M ∑ ∑

[(

i=1 j=1

D∇ 4 − Tx

) ∂2 ∂2 2 − T − ρ H ω Xi Yj − PRij y s ∂ x2 ∂ y2 −

L ∑

] iωcℓ Xi (xℓ ) Yj (yℓ ) δ (x − xℓ , y − yℓ ) Cij = PD ,

(11)

ℓ=1

where PD = iω φD |z =−d − φD |z =−d−H . A similar expression is obtained for the radiation-problem pressure PRij . Given the orthogonality of the modes, one can recast equation (11) into a matrix form (with diagonal inertia and stiffness matrices) by multiplying it by a given mode shape Xp Yq and integrating over the surface of the plate. The following system of equations is obtained:

(

)

N M ∑ ∑ [ 2( ) ( ) ] −ω aijpq + σijpq − iω bijpq + ψijpq + δijpq + ϵijpq Cij = fpq , i=1 j=1

(12)

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

677

where

δijpq =

∫∫

ϵijpq =

∫∫

σijpq =

∫∫

ω2 aijpq + iωbijpq =

∫∫

ψijpq =

D∇ 4 Xi Yj Xp Yq dS ,

(

( − Tx

L ∑

∫ℓ=∫1 fpq =

)

∂2 ∂2 + Ty 2 2 ∂x ∂y

(13a)

)

Xi Yj Xp Yq dS ,

(

)

(13b)

ρs HXi Yj Xp Yq dS ,

(13c)

PRij Xp Yq dS ,

(13d)

cℓ Xi (xℓ ) Yj (yℓ )

∫∫

δ (x − xℓ , y − yℓ ) Xp Yq dS ,

(PI + PS ) Xp Yq dS .

(13e) (13f)

Eq. (12) can now be solved for coefficients Cij . 4. Solution methodology and validation To solve the boundary value problem (2) for the diffracted and radiated potentials φD and φRij , we use the open source Boundary Element Method solver ‘‘NEMOH’’ (Babarit and Delhommeau, 2015). This linear solver, developed at École Centrale de Nantes, is used to investigate seakeeping problems by calculating the response in frequency domain of a floating or submerged body in specific sea conditions. The three-dimensional (volumic) fluid problem is replaced, through Green theorem, by a two-dimensional (surfacic) problem of a discretized source distribution on a set of panels representing the wetted surface of the body. NEMOH has been extensively used in design optimization of marine energy devices (Babarit et al., 2015; Olondriz et al., 2016; Antonutti et al., 2016, 2014; Andersen et al., 2015; Schmitt et al., 2016) and has undergone numerous validations against experimental (Babarit et al., 2015; Andersen et al., 2015; Schmitt et al., 2016; Crooks et al., 2016) and numerical (Olondriz et al., 2016; Parisella and Gourlay, 2016; Crooks et al., 2016) results. Inputs to NEMOH are wave frequency and angle of incidence, water depth, body geometry, location and number of degrees of freedom. In our case, the degrees of freedom correspond to the normal mode deformations of the plate. The results are given in terms of a pressure field on the body surface and the free surface elevation around the body for each boundary value problem. 4.1. Validation For validation of our approach, which is based on Boundary Element Method and normal mode analysis, we consider the case of a submerged two-dimensional (i.e. W → +∞) membrane (i.e. D = 0) under the action of over-passing surface waves. For this problem, an analytical solution was previously obtained through eigenfunction matching (Cho and Kim, 1998). Here, we compare, in Fig. 2, (i) the total pressure force on the membrane for different tensions in x-direction (Fig. 2a), and (ii) the membrane elevation for different submergence depths (Fig. 2b–e). The total force per unit length is given by Ftot = Fex +

N ∑

Ci FRi ,

(14)

i=1

where Fex = (PI + PS ) dx, is the excitation force and FRi = PRi dx is the radiation force associated with the mode i. Since NEMOH is a three-dimensional solver, we choose a relatively long length (W /L = 16) to minimize the effects of threedimensionality, and also we pick a small thickness (H /L = 10−3 ) to be able to compare our results with theoretical results of Cho and Kim (1998) which are derived under zero-thickness assumption. A convergence study on the number of modes and panels leads to a required minimum number of modes N = 5, and a minimum number of panels of 20, 60 and 3 respectively along the x, y, and z-axes. Fig. 2a presents the normalized total force as a function of the normalized wavenumber for three different tensions. Note that to facilitate the comparison with analytical solution (Cho and Kim, 1998), all parameters are expressed in their dimensional form and normalized explicitly. Clearly theoretical results (lines) and our results (markers) validate each other very well. The peaks of force correspond to resonant wavenumbers of the membrane. As the stiffness of the membrane increases, the dominant resonant wavenumber increases. This can be explained by the fact that on one hand the natural frequencies of the membrane are proportional to the square root of the total stiffness and on the other hand, the excitation frequency of the waves goes as ω ≈ µ1/2 . Therefore, the larger the stiffness, the larger the resonant frequencies, and so the larger wavenumber. Moreover, as shown on Fig. 2b–e where the membrane’s vertical displacement is plotted as a function of the wavenumber and the horizontal coordinate for two different submergence depths, the first resonant peak corresponds to mode 1 (i.e. cos (π x/L)) of the membrane. The vertical displacement from our computations are in good quantitative agreement with analytical results from Cho and Kim (1998) for two different submergence depths.

∫

∫

678

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

Fig. 2. (a) Normalized total force on a two-dimensional membrane plotted as a function of dimensionless wavenumber µ for three different tensions. Analytical results (Cho and Kim, 1998) are shown by lines and results of our numerical simulations (through normal modes analysis) are shown by markers. Other physical parameters are d/h = 0.2, L/h = 1 and ρs = 1. (b–e) Dimensionless amplitude |ξ | /A of the membrane elevation as a function of dimensionless wavenumber µ and dimensionless horizontal coordinate x/L. (b) and (c) correspond to our computations for d/L = 0.2 and d/L = 0.1 respectively (while)(d) and (e) are the analytical results from Cho and Kim (1998) for d/L = 0.2 and d/L = 0.1 respectively. Other physical parameters are L/h = 1, Tx / ρ gh2 = 0.1 and ρs = 1.

4.2. Conservation of energy in computational domain In this work, we are interested in the investigation of the effect of different absorber’s parameters on the performance of a flexible submerged WEC depicted in Fig. 1. As a measure of efficiency of a WEC, Capture Width Ratio (CWR) is widely used in the literature and is defined as CWR =

PPTO P w av e

,

(15)

where Pwav e = W /(4ω) (1 + 2µ/sinh (2µ)) is the dimensionless incident wave energy flux averaged over a period across a wavefront of width W (i.e. the power seen by the device across its width) and PPTO is the power generated during one period

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

679

Fig. 3. Normalized power harvested by the WEC as a function of the normalized damping coefficient of the PTO units cPTO . The blue line is calculated from the motion of the PTO units, whereas the red markers come from the evaluation of the wave energy flux far from the body. The device is mounted on two PTO units located at x1 /L = y1 /W = −0.2 and x2 /L = y2 /W = 0.2, and of damping coefficient c1 = c2 = cPTO . Other physical parameters are λ = 8, d = 0.4, L = λ/2, W /L = 1/2, H /L = 2 × 10−3 , D = 0.03, Tx = Ty = 0.1, θ = 0◦ and ρs = 1.

by the PTOs which is given by PPTO =

L ∑ 1

ℓ=1

2

cℓ ω2 |ξ (xℓ , yℓ )|2 .

(16)

Alternatively, we can calculate energy uptake by our WEC by considering a cylindrical control volume about the device (that extends from the seabed to the free surface), and calculate the net energy flux inward (or outward). With incident waves and in the absence of PTO, this net energy flux is clearly zero (whatever energy gets in from one side, goes out from the other). If the device takes out energy, then the net energy flux inward is positive. This approach was introduced in Newman (1976) and used the same year in Evans (1976) to perform design optimization of oscillating bodies for wave energy extraction. It is now commonly used in WEC performance analysis (Falnes and Hals, 2012; Falnes and Kurniawan, 2015; Rainey, 2012; Farley, 2012) and cloaking phenomenon investigation (Newman, 2013). For each boundary value problem (2), a perturbed potential φp far from the body can be approached (Delhommeau, 1987) by

√ φp (r , β, z ) = 2 where

(

r =

√

2πµ r

H (β) C0

cosh [µ (z + 1)] i(µr + π ) 4 e +O cosh µ

x2 + y2 , β = tan−1 (y/x) , z

)

(

1

) (17)

3

r2

are cylindrical coordinates about the z-axis, C0 = −

µ2 [ 0 ] µ tanh µ µ2 −µ20 +µ

with

µ0 = ω the dimensionless wavenumber in infinite depth, and H (β) is the Kochin function. Eq. (17) shows that H (β) also represents the amplitude of the far-field waves due to the WEC (Newman, 1976). Hence, it contains information about the geometry (in the diffraction problem) and is proportional to the amplitude of the WEC motion (in the radiation problems). In addition to the calculation of the potential of the scattered and radiated waves, the Kochin function is also numerically calculated by NEMOH. The perturbed potential (17) can be used to calculate the wave energy flux through a vertical cylinder surrounding the device (see Babarit, 2015 for numerical details), which corresponds to the absorbed (or generated) power as 2

( PK = J0

8π C0 ωRe (H (θ)) − 8πµC02 ω2

2π

∫

|H (β)|2 dβ

) (18)

0

where J0 = 1/(4ω) [1 + 2µ/sinh (2µ)] is the incident wave energy flux averaged over a period. Let us consider the WEC described previously, located at a submergence depth d = 0.4 and mounted on two PTO units with the same damping coefficients c1 = c2 = cPTO located at x1 /L = y1 /W = −0.2, x2 /L = y2 /W = 0.2. The size of the absorber corresponds to L = λ/2, W /L = 1/2 and H /L = 2 × 10−3 , with a flexural rigidity D = 0.03 and a density ρs = 1. We apply a planar tension Tx = Ty = 0.1 and an incident wave of wavelength λ = 8 is arriving from x = −∞ along the x-axis (θ = 0◦ ). We choose numbers of modes N = M = 6, and numbers of panels 30, 40 and 3 respectively along the x, y and z-axes, for which our simulations are converged. We compare in Fig. 3, PPTO and PK as a function of the PTO damping coefficient cPTO , and excellent agreement is found. By making sure that PPTO and PK are equal, we confirm that the energy absorbed by our WEC does not remain in the fluid domain. In the following, the absorbed power is calculated according to Eq. (16), since it does not require the calculation of the Kochin function, thus saves computational time, and avoids approximation error of the integral term in Eq. (18).

680

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

Fig. 4. Dependence of (a) CWR, (b) optimum location of the PTO along the x-axis, (c) optimum location of the PTO along the y-axis and (d) optimum damping coefficient of the PTO as a function flexural rigidity D and in-plane tension T . At the optimum point CWR = 0.78 for D = 0.018, T = 0, xopt /L ≈ 0.35, yopt = 0 and copt = 92. Other physical parameters are λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , θ = 0◦ and ρs = 1.

5. Results and discussion We proceed in this section to a parametric study of the parameters of the carpet (or flexible plate): the carpet rigidity, the location and the damping coefficient of the PTO units, and the shape of the structure. We first consider a single-PTO configuration, and then add a second unit to the setup. After searching for the optimal design parameters in terms of maximum CWR, we comment on the influence of the carpet length, the incident wave direction and the carpet aspect ratio on the device efficiency. 5.1. Energy harvesting with a single PTO unit Here we consider that one PTO unit is attached to our wave energy device. For our simulations we consider a carpet of density ρs = 1 and with the size L = λ, W /L = 2/3 and thickness H /L = 2 × 10−3 located at the depth d = 0.3. Note that this submergence depth constitutes a trade-off between avoiding large loads existing at the free surface and harvesting wave energy which is more abundant, in a given water column, as one gets close to the free surface. Considering the short distance between the carpet and the water surface and that the carpet is fixed all around its edges, the investigated WEC would be best suitable for intermediate-depth water regimes (relatively shallow) in order to allow for the installation of a supporting structure (to keep the carpet edges fixed) on the seafloor at reasonable manufacturing costs. To maintain the carpet secure at this submergence depth, we also consider the WEC operates under small to moderate amplitudes of sea states. We assume incident monochromatic waves of wavelength λ = 3 arriving normal to the carpet (i.e. θ = 0◦ ). As previously, we choose numbers of modes N = M = 6, and numbers of panels 30, 40 and 3 respectively along the x, y and z-axes. The effect of flexural rigidity D and in-plane tension Tx = Ty = T on the CWR is shown in Fig. 4. For each pair of (D, T ), we find the maximum CWR (Fig. 4a) by optimizing the x-location (Fig. 4b), the y-location of the PTO unit (Fig. 4c) as well as the PTO damping coefficient (Fig. 4d). Note that all optimization processes in this work consist of sweeping PTO parameters over relevant ranges then selecting optimum parameters that yield maximum CWR. We assume that damping coefficient variation leads to a unique maximum CWR, and PTO locations are always swept over their entire meaningful range (x, y) ∈ [−L/2, L/2] × [−W /2, W /2] (i.e. all over the carpet). The optimum y-location has to be at zero as shown

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

681

Fig. 5. (a) CWR, (b) associated PTO amplitude and (c) optimal damping coefficient of the PTO as a function of PTO location on the surface of the carpet. Optimum CWR = 0.71 for (x0 /L, y0 /W ) ≈ (0.35, 0) and copt = 102. Other physical parameters are λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

in Fig. 4c due to the conjunction of the carpet geometry and the incidence angle. On one hand, the geometry of the plate being elongated in the x-direction, it will tend to bend more along this axis. Given that wave power is harvested through the excited natural modes of the structures, one can expect the optimal PTO location to be on a line parallel to the long axis of the plate. On the other hand the carpet being symmetric with respect to the incoming wave direction (here the x-axis), the optimal PTO location will coincidence with the long symmetry axis of the carpet, where y = 0. The optimum x-location, nevertheless, is almost independent of D and T (Fig. 4b), and is always near xopt /L ≈ 0.35, i.e. past the middle section of the flexible board toward the downstream. Optimum damping coefficient copt varies not monotonically, but is usually higher for higher (D, T ) (Fig. 4d). The CWR varies significantly with (D, T ), gaining an impressive maximum of 0.78 for T = 0, D = 0.018 and copt = 92, and decreases as D departs from the optimal value or T increases. The variation with D is intuitive, as larger D means stiffer carpet and at the limit of D → ∞, it is simply too stiff to move at all, hence no energy will be produced. The other extreme limit, when D → 0, the carpet turns into a membrane and does not collect the wave force over the area, hence overall efficiency will be very small. Therefore, there must be an optimum value for D between zero (membrane) and infinity (rigid board), as we obtain in Fig. 4a. Results also suggest that if D is very small, then T takes the control: the optimum CWR is obtained for a finite value of T (cf. the row corresponding to D = 0.005-0.01 in Fig. 4a). But for large D, flexural rigidity dominates the dynamics and no in-line tension is required. In fact, adding a tension decreases the CWR. To gain a better insight into the dynamics of the carpet under the action of waves, we pick a constant flexural rigidity D = 2.25 × 10−2 and an in-plane tension T = 2.7 × 10−2 which give CWR = 0.71. Those values will be kept constant throughout the rest of the paper. We then vary the location of our single PTO unit (whose parameters are denoted by index ‘‘0’’), and for each location (x0 , y0 ) find the optimum damping coefficient copt that gives maximal CWR. We show the variation of the CWR and the amplitude of PTO as a function of its location (x0 , y0 ) in Fig. 5a–b. At any chosen location, we calculate the optimum damping coefficient copt (shown in Fig. 5c) that maximizes the CWR. Clearly, it can be seen that placing a PTO at x0 /L ≈ 0.35 has a very distinct superiority in terms of PTO movement (Fig. 5b) and resulting CWR (Fig. 5a). It is, nevertheless, to be noted that the maximum CWR (obtained at (x0 /L, y0 /W ) ≈ (0.35, 0)) does not occur exactly at the point of maximum PTO elevation (obtained at (x0 /L, y0 /W ) ≈ (0.3, 0)), which is not unexpected as this is the case in almost any vibration system (i.e. one can increase the amplitude by changing system parameters, but energy damped or generated may decrease). Assuming that the PTO is at y0 = 0 (i.e. on the centerline of the carpet), we plot in Fig. 6 the variation of the CWR as a function of x0 and c0 (the latter being the damping coefficient of the PTO). Clearly and as before, whatever the value of c0 , x0 has to stay toward the downstream part of the carpet to get the optimal CWR. Moreover, when c0 increases, the x0 -bandwidth of high CWR gets smaller and the PTO has to get closer to the downstream edge of the carpet (x0 /L → 0.5).

682

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

Fig. 6. CWR as a function of damping coefficient c0 and location x0 of the PTO. Optimum CWR = 0.71 is obtained for x0 /L ≈ 0.35 and c0 = 97. Other physical parameters are y0 = 0, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

Fig. 7. (a, b) CWR and (c, d) optimal damping coefficient as a function of PTO location on the surface of the carpet for a submergence depth of (a, c) d = 0.2 and (b, d) d = 0.4. Other physical parameters are λ = 3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

Finally, we show in Fig. 7 the maximal CWR and optimal damping coefficient of a single PTO as a function of its location (equivalent to Fig. 5a, c) for two submergence depths d = d∗ /h = 0.2 and 0.4, d∗ being the dimensional submergence depth. Clearly the two configurations lead to different carpet behaviors, with a symmetric response with respect to x0 /L ≈ 0 for d = 0.2 and a very similar response to the case of d = 0.3 for d = 0.4. We find that the maximum CWR is 0.27 (with copt = 125) and 0.81 (with copt = 90) for the submergence depth d = 0.2 and 0.4, respectively. However, it is not possible to conclude from this data that the optimal CWR of a WEC increases monotonically with the submergence depth since these results have been obtained with for fixed values of structural parameters (D, T ). Given that these stiffness parameters modify the resonance frequencies of the carpet, the optimization of the CWR with respect to (D, T ) at each submergence depth could

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

683

Fig. 8. CWR as a function of PTO locations x1 , x2 on the surface of the carpet. Two hotspots corresponding to optimum CWR ≈ 0.8 are located on each side of the symmetry line x1 = x2 , e.g. (x1 /L, x2 /L) ≈ (−0.35, 0.35) and (0, 0.35). Other physical parameters are y1 = y2 = 0, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

lead to a different conclusion. In fact, the value of d has a non-trivial effect upon the performance of the carpet since on one hand, the wave excitation force decreases with the submergence depth, thus decreasing the motion of the structure. But on the other hand, the hydrodynamic added mass and damping also tend to decrease with the submergence depth (Mihelcic, 1989; Zheng et al., 2007), thereby favorising the response of the carpet. Therefore, one can expect the value of d to influence the performance of the carpet in a non-monotonic manner. While the methodology developed in this paper work is relevant for large to moderate values of the submergence depth (smaller values leading quickly to nonlinear effects), the results from the three values of d (0.2, 0.3, 0.4) investigated here need to be supplemented with more computations to conclude on the effect of the submergence depth, a task we leave for future work. 5.2. Energy harvesting with two PTO units Now let us consider two PTO units, denoted PTO1 and PTO2 . If waves arrive normal to the carpet (i.e. θ = 0◦ ), then due to symmetry, PTOs at their optimum location are either on the centerline of the carpet (y1 = y2 = 0), or side by side (y1 = −y2 , x1 = x2 ). Here, since the carpet is elongated along the incoming wave direction (i.e. L > W ) and considering previous results about optimum y-location with one PTO (Fig. 4c), both PTOs are assumed to be on the centerline. For each location of the two PTOs, i.e. (x1 , x2 ), we calculate the optimum c1 , c2 that yield the maximum CWR. We obtain c1 = 297 and c2 = 85. This maximum CWR as a function of x1 and x2 is presented in Fig. 8. Note that as expected the plot is symmetric about the line x1 = x2 as areas about this line basically represent the same configuration with index of PTO1,2 switched. Two hotspots on each side of the symmetry line can be identified corresponding to (x1 /L, x2 /L) ≈ (−0.35, 0.35) and (0, 0.35) for which CWR ≈ 0.8. Similarly to the case of a single PTO, the common position of these hotspots, namely when x2 ≈ 0.35, is close to the location of maximum displacement of the plate, as we later show in Fig. 11a–b. In many practical applications, the limiting factor is the damping coefficient of the PTO system, which has a limited freedom in tuning. The location of the two PTO units, assuming that the damping coefficient of PTOs is given and fixed, influence the overall performance. We show the overall CWR as a function of locations of the two PTOs (x1 , x2 ), for four different damping coefficients c1 = c2 = 20, 51, 200 and 1000 respectively in Fig. 9a–d. For a relatively small damping coefficient c1 = c2 = 20, the optimum location is at x1 /L = x2 /L = 0.32 for which CWR = 0.59 (Fig. 9a) suggesting that at least one of the PTOs must be nearly halfway on the downstream side, but the optimum obtains for both at the same location. That is, the optimum case is a single-PTO configuration with a damping coefficient twice as high. For a higher damping of c1 = c2 = 51, the CWR distribution is close to that of in Fig. 8, with five highlighted subareas with highest CWR. As the damping ratio further increases (e.g. c1 = c2 = 200, Fig. 9c), the optimum locations of the two PTOs move toward the two ends: one PTO must be very close to the upstream end and the other must be close to the downstream end. The overall performance decreases under a very high damping coefficients of c1 = c2 = 1000 as is shown in Fig. 9d where only at corners CWR may be large. In fact, for the case of Fig. 9d, the maximum CWR = 0.69 for x1 /L = 0.46 and x2 /L = −0.38. 5.2.1. Influence of the length Here we consider the previous configuration, and we show the evolution of the CWR when the PTOs are set to their optimal damping coefficient which are no longer assumed to be equal (i.e. c1 = 297, c2 = 85) and their positions (x1 , x2 )

684

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

Fig. 9. CWR as a function of PTO locations x1 , x2 for (a) c1 = c2 = 20, (b) c1 = c2 = 51, (c) c1 = c2 = 200, and (d) c1 = c2 = 1000. Other physical parameters are y1 = y2 = 0, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

are swept (Fig. 10a), and when the PTOs are at their optimal location (i.e. (x1 /L, x2 /L) = (−0.36, 0.34) from Fig. 8) and their damping coefficients c1 , c2 are swept (Fig. 10b). As Fig. 10a shows, having different damping coefficients breaks the symmetry found previously (Fig. 9), although the highest CWR (= 0.83) is obtained for a similar configuration, that is having one PTO on each side of the carpet, with (x1 /L, x2 /L) = (−0.36, 0.34), (c1 , c2 ) = (297, 85). Note that PTO2 (located downstream) has the lowest damping coefficient, and c2 has a significant influence on the CWR, unlike c1 which influence is very limited (Fig. 10b). This indicates that PTO1 does not extract power, and only has a role of conditioning the carpet to maximize the power extracted by PTO2 . To get a better appreciation of the dynamics of the PTOs, we consider the same case as Fig. 10a and plot the amplitude of the displacement of PTO1 (Fig. 11a). Clearly, the optimal locations (x1 /L, x2 /L) = (−0.36, 0.34) correspond to the minimum PTO1 displacement. We show in Fig. 11b that, at those optimal PTO locations, a significant part of the carpet (i.e. x/L < −0.25) is almost static. This suggests not only that one PTO would be enough to harvest essentially all the power, but also that there exists an effective length of the depicted flexible carpet material beyond which the WEC would not produce a much higher power. To investigate this issue, we evaluate the power extracted by a moving PTO, denoted with index ‘‘m’’ and having a damping coefficient cm = 85, as a function of its position along the x-axis in two cases, that are with and without a second fixed PTO with damping coefficient c0 = 297 located at x0 /L = −0.36. Both PTOs are on the centerline of the carpet. As we can see in Fig. 12, having a fixed PTO located upstream (continuous curve) leads to a higher maximum power extracted by the moving one, than not constraining the upstream part of the carpet (dashed curve). This shows that, more than having PTOs, it is the tuning of the upstream motion of the carpet that leads to an overall better WEC efficiency. In addition, the ‘‘rigid board’’ tuning appears to be the most efficient one. Physically, that is the effect of a brutal apparent water depth reduction, leading waves amplitude to increase as they are going over the carpet (i.e. shoaling). Therefore, the pressure differential between the top and the bottom gets larger with distance from the upstream edge, and likewise for the carpet excitation. This effect being directly related to the constraint applied on the fluid motion, a rigid carpet yields larger downstream excitation. In

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

685

Fig. 10. CWR as a function of (a) PTO locations x1 , x2 for c1 = 297, c2 = 85, and (b) damping coefficients c1 , c2 for (x1 /L, x2 /L) = (−0.36, 0.34). The global optimum CWR = 0.83 is obtained for (x1 /L, x2 /L) = (−0.36, 0.34), c1 = 297 and c2 = 85. Other physical parameters are y1 = y2 = 0, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

Fig. 11. (a) Amplitude of the upstream PTO (PTO1 ) as a function of PTO locations x1 , x2 , and (b) elevation of the carpet for (x1 /L, x2 /L) = (−0.36, 0.34). Other physical parameters are c1 = 297, c2 = 85, y1 = y2 = 0, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

real working conditions, the stability of the device is reduced when the board gets longer as it becomes more easily subject to potential buckling and other structural instabilities. It accordingly requires higher manufacturing costs, resulting in the existence of an optimal length providing highest energy production–manufacturing costs ratio, which determination is out of the scope of this paper. So far, we limited our study to a carpet length equal to the wavelength of incoming waves. However, in operational conditions, the sea state is a spectrum that would lead to different values of L/λ, the carpet length to wavelength ratio. To investigate the performance of our WEC when varying this ratio, we consider the optimal configuration found previously, that is two PTOs on the centerline of the carpet with a fixed one located at (x0 /L, y0 /W ) = (−0.36, 0) with a damping coefficient c0 = 297. We plot in Fig. 13 the power harvested by the moving PTO of damping coefficient cm = 85 as a function of its position along the x-axis and the ratio L/λ ∈ [0.5, 2]. It shows that the extracted power is maximal for relative lengths close to L/λ = 1. This could be expected because it is for this ratio that the CWR has been optimized with respect to the structural rigidity (D, T ) of the carpet. For other values of the relative length, the pair (D, T ) has been held fixed and the extracted power optimized only with respect to the second PTO position. Fig. 13 also shows that the optimal PTO location varies considerably with L/λ. When this relative length is different from unity, the optimal PTO locations are no longer concentrated around x/L ≈ 0.35 but are rather distributed over a wider range of the carpet length. Even though for a large part of the relative length (L/λ ∈ [0.5, 1.5]), the spot of optimal PTO location is downstream, there is a transition for longer carpets L/λ > 1.5 where the spots of largest CWR are found not only to broaden around the center position x/L ≈ 0, but also to show a few oscillations as a function of xm /L. These oscillations for longer carpets and the lack thereof for shorter ones can be related to the excited modes of the WEC. Clearly, every other parameter being held fixed, longer carpets experience the incoming wave with an apparent higher frequency. The response of such carpets would be thus dominated by the higher

686

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

Fig. 12. Value of the dimensionless power extracted by a moving PTO with damping cm = 85 as a function of its position along the x-axis in two different cases: with (continuous curve) and without (dashed curve) a fixed PTO located at x0 /L = −0.36 with damping c0 = 297. Other physical parameters are y0 = ym = 0, λ = 3, d = 0.3, L = λ, W /λ = 2/3, H /λ = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

Fig. 13. Value of the dimensionless power extracted by a moving PTO with damping cm = 85 as a function of its position along the x-axis for different carpet lengths. Another fixed PTO is at (x0 /L, y0 /W ) = (−0.36, 0). Other physical parameters are ym = 0, λ = 3, d = 0.3, W /λ = 2/3, H /λ = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

deformation modes, while shorter carpets are dominated by the first modes. In conjunction with this, the presence of the upstream PTO inducing a ‘‘rigid board’’ zone and the shoaling effect leads to the optimal PTO locations to be concentrated (first vibration modes) downstream for short carpets, while they spread over longer carpets (higher vibration modes). 5.2.2. Influence of the angle of incidence In actual operation conditions, ocean waves arrive in the form of a directional frequency spectrum. It is thus essential to be able to predict the behavior of the structure for the corresponding range of incident wave angles. In this section, we investigate the influence of the incident angle on the optimal PTO locations by considering two PTOs no longer assumed to be on the centerline of the carpet. However, we fixed their damping coefficient to half of the optimal damping found with one single PTO and fixed (D, T ) (Fig. 5c), that is c1 = c2 = 102/2 = 51. Four incident angles are tested: 0◦ , 30◦ , 60◦ and 90◦ . In the zero degree case, optimal locations are on the centerline of the carpet (i.e. y1 = y2 = 0), one PTO being upstream and the second one downstream (Fig. 14a–b). For positive angle, the optimal locations move toward the positive values of the y-axis (Fig. 14c–d). This is in agreement with previous observation that PTOs perform better downstream, despite the fact that now the centerline symmetry is broken (i.e. y1 ̸ = 0 and y2 ̸ = 0). For all angles, the optimal PTO locations are found to be on lines almost parallel to the long axis of the carpet, in accordance with the explanation we provided earlier on the bending tendency of an elongated plate. Yet, the incidence direction of the wave being allowed to vary from the long axis, the optimal PTO locations will differ from the y = 0 axis except for the zero angle case. For 90◦ angle, optimal locations are now symmetric with respect to the y-axis (i.e. x1 = −x2 and y1 = y2 ). Since the incident wave is unidirectional and the in-plane tension identical in both x and y-directions, 90◦ angle case only differs from 0◦ angle case in the carpet aspect ratio W /L seen by the wave. That means that there exists a critical aspect ratio corresponding to a transition of optimal locations of the PTOs from the centerline to the sides of the carpet. We investigate this phenomenon in next section. Note that for the 90◦ angle, the highest CWR is greater than 1 due to the definition of the CWR (15) using the initial width of the carpet W . This definition assumed that the width W provided the largest exposure of the carpet to the incident waves, which is strictly

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

687

Fig. 14. CWR as a function of one PTO location while the other one stays at its optimal location found after the simultaneous optimization of both PTO locations. Left row (a, c, e, g): PTO1 location varies while PTO2 remains fixed; right row (b, d, f, h): PTO2 location varies while PTO1 remains fixed. Four incident wave directions from θ = 0◦ to 90◦ are taken. (a) and (b) correspond to a 0◦ incident angle, with (x2 /L, y2 /W ) = (0.32, 0) for (a) and (x1 /L, y1 /W ) = (−0.33, 0) for (b). (c) and (d) correspond to a 30◦ incident angle, with (x2 /L, y2 /W ) = (0.28, 0.18) for (c) and (x1 /L, y1 /W ) = (−0.02, 0.15) for (d). (e) and (f) correspond to a 60◦ incident angle, with (x2 /L, y2 /W ) = (0.17, 0.2) for (e) and (x1 /L, y1 /W ) = (−0.13, 0.2) for (f). (g) and (h) correspond to a 90◦ incident angle, with (x2 /L, y2 /W ) = (0.13, 0.2) for (g) and (x1 /L, y1 /W ) = (−0.13, 0.2) for (h). Other physical parameters are c1 = c2 = 51, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 and ρs = 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

true for the zero angle case only. For 90◦ incident waves however, the carpet exposure length is L, which is larger than W . Therefore, the maximum CWR is likely to exceed 1. Additionally, color contrast between Fig. 14a, c, e and 14b, d, f indicates that location of the downstream PTO is much more influent on the CWR value than the upstream one, which concurs with previous observation that the upstream PTO is used to constrain the beginning of the carpet to maximize power extraction by the downstream one.

688

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

Fig. 15. CWR as a function of incident wave angle, with and without PTO location optimization. Without optimization, the PTO locations correspond to the optimal ones for the 0◦ incident angle, i.e. (x1 /L, y1 /W ) = (−0.33, 0) and (x2 /L, y2 /W ) = (0.32, 0). Other physical parameters are c1 = c2 = 51, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 and ρs = 1.

Hereafter PTO locations are fixed, and thus could not be optimized for each incoming wave angle. In order to quantify the effect of such an optimization on the WEC performance, we plotted in Fig. 15 the CWR as a function of incident angle in two configurations: with and without PTO locations optimization. When locations are not optimized, we take the optimal ones for the 0◦ case, i.e. (x1 /L, y1 /W ) = (−0.33, 0) and (x2 /L, y2 /W ) = (0.32, 0). As we could anticipate from the displacement of optimal PTO locations, the difference between optimized and non-optimized CWR increases as the angle gets larger. Hence, we expect the CWR to significantly decrease when the system is interacting with directionally distributed waves. 5.2.3. Influence of the aspect ratio As for the angle of incidence, we investigate the optimal PTO positions for different carpet aspect ratios. We keep the same setup as previously, i.e. we have two PTOs with c1 = c2 = 51. Five aspect ratios are taken, from a square (W /L = 1) to W /L = 1.45, while the initial carpet area is kept constant. The definition of the capture width ratio is adapted for each case by considering in Eq. (15) the actual width W of the carpet. For low aspect ratios (i.e. W /L ≤ 1.4) the optimal PTO positions are on the centerline of the carpet (Fig. 16a–h) and the patterns are only slightly contrasted (CWR always between 0.65 and 0.5). That means that if one PTO is well positioned (i.e. downstream), the other will not have a significant influence on the overall CWR. However and as mentioned previously, by increasing its aspect ratio, the carpet will tend to bend more along the y-axis (direction of its elongated dimension), and the resulting optimal PTO positions will eventually both be downstream on an axis parallel to the y-axis (Fig. 16i–j). In this latter configuration, both PTOs act in an exact same way on each side of the carpet, thus they equally participate to the power extraction. An inappropriate location of one of the PTOs can make the CWR drop to 0.42, while in optimal configuration it reaches 0.7. This indicates that for a carpet shape elongated in a direction normal to the incident wave one, having multiple PTOs becomes essential for an efficient power extraction. 6. Added PTO stiffness In order to be able to tune the WEC response to the excitation frequency, PTO units usually include a system whose stiffness can be adjusted. In this section, we investigate the effect of an added linear spring in parallel to our dash-pot damper. Taking the spring restoring force into account, the localized PTO force becomes fPTO (x, y, t ) =

L ∑

ℓ=1

( ) ∂ − cℓ + kℓ ξ (xℓ , yℓ , t ) δ (x − xℓ , y − yℓ ) , ∂t

(19)

where kℓ is the spring stiffness. The equation of motion (12) turns into N M ∑ ∑ [ 2( ) ( ) ] −ω aijpq + σijpq − iω bijpq + ψijpq + δijpq + ϵijpq + κijpq Cij = fpq ,

(20)

i=1 j=1

where κijpq is the PTO stiffness matrix, and is of the form:

κijpq =

L ∑

kℓ Xi (xℓ ) Yj (yℓ )

∫∫

δ (x − xℓ , y − yℓ ) Xp Yq dS .

(21)

ℓ=1

Similarly, as in Section 5.1, maximum CWR is found by optimizing the stiffness and damping coefficient of one single and fixed PTO (x0 /L = 0.35, y0 = 0) in a specific range of (D, T ). We notice in Fig. 17a that adding a localized PTO stiffness does

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

689

Fig. 16. CWR as a function of one PTO location while the other one stays at its optimal location found after the simultaneous optimization of both PTO locations. Left row (a, c, e, g, i): PTO1 location varies while PTO2 remains fixed; right row (b, d, f, h, j): PTO2 location varies while PTO1 remains fixed. (a) and (b) correspond to the aspect ratio W /L = 1, with (x2 /L, y2 /W ) = (0.32, 0) for (a) and (x1 /L, y1 /W ) = (0.2, 0) for (b). (c) and (d) correspond to the aspect ratio W /L = 5/4, with (x2 /L, y2 /W ) = (0.27, 0) for (c) and (x1 /L, y1 /W ) = (0.13, 0) for (d). (e) and (f) correspond to the aspect ratio W /L = 4/3, with (x2 /L, y2 /W ) = (0.24, 0) for (e) and (x1 /L, y1 /W ) = (0.12, 0) for (f). (g) and (h) correspond to the aspect ratio W /L = 1.4, with (x2 /L, y2 /W ) = (0.23, 0) for (g) and (x1 /L, y1 /W ) = (0.15, 0) for (h). (i) and (j) correspond to the aspect ratio W /L = 1.45, with (x2 /L, y2 /W ) = (0.2, 0.12) for (i) and (x1 /L, y1 /W ) = (0.2, −0.12) for (j). Other physical parameters are c1 = c2 = 51, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , D = 2.25 × 10−2 , Tx = Ty = 2.7 × 10−2 , θ = 0◦ and ρs = 1.

690

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

Fig. 16. (continued)

not change the range of optimal (D, T ) (optimum found for D ≈ 0.02 and T ≈ 0 likewise Fig. 4a). That means that for this range of flexural rigidity and tension, the structural rigidity of the carpet dominates its dynamics, and adding a localized PTO stiffness does not have any effect on the optimal structural properties. Nevertheless, it increases the CWR for each pair (D, T ). Even if the gain is small for the maximum CWR (from 0.78 for a simple damper to 0.82 with a PTO stiffness), it leads to a significant increase when (D, T ) get larger (from 0.16 to 0.32 for D = 0.05, T = 0.5). The optimal stiffness is negative for each pair (D, T ) and gets more negative as (D, T ) get larger (Fig. 17b). That denotes that the carpet is over-stiff, and the PTO response needs to be detuned from the initial resonant frequency determined by the structural rigidity and damping coefficient to a lower one. The optimal damping coefficient is higher for higher (D, T ) (Fig. 17c), as previously noticed in Section 5.1. However, it now varies monotonically due to the presence of the spring. 7. Concluding remarks In this work, we showed that the use of a properly tuned flexible carpet for wave energy conversion can reach a capture width ratio close to 80%. A three-dimensional linear wave-flexible structure interaction model was developed by using a BEM solver. Considering a horizontal submerged flexible plate, we computed the diffraction and radiation solutions in the case of monochromatic wave fields. The validity of the results was then verified by comparing with a two-dimensional analytical formulation based on the linear hydro-elastic theory, and excellent agreements were obtained. Another validation was done through the assessment of energy balance in the domain. The effects of the modification of the main Wave Carpet parameters were explored for a typical implementation case of the device in the ocean, and optimization procedures allowed us to underline the optimal setup to harvest a maximal amount of energy by a comprehensive parametric study. First, we showed that the flexural rigidity of the carpet plays an important role in the efficiency, and there is an optimal range in which the power extracted from the incident wave is maximal. Also, the location of the PTOs has a more significant effect on the efficiency than their damping coefficient. Then, we pointed up

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

691

Fig. 17. Dependence of (a) CWR, (b) optimum damping coefficient of the PTO and (c) optimum stiffness of the PTO as a function flexural rigidity D and in-plane tension T . At the optimum point CWR = 0.82 for D = 0.021, T = 0, kopt = −60 and copt = 85. Physical parameters are x0 /L = 0.35, y0 = 0, λ = 3, d = 0.3, L = λ, W /L = 2/3, H /L = 2 × 10−3 , θ = 0◦ and ρs = 1.

that, for a range of carpet lengths around the incoming wavelength, the optimal location of the PTOs is downstream (past the middle of the carpet along the incident wave direction), while the almost stationary upstream part of the carpet has a role of conditioning the incoming wave. Furthermore, a carpet with a rigid upstream region leads to a higher CWR than that of a fully flexible plate because the former amplifies the shoaling effect to the incoming waves, due a more brutal change of water depth. In addition, it was shown that the optimal PTO locations are sensible to the incident wave direction as well as the carpet aspect ratio, and the power extracted by the device can be significantly affected. Finally, we highlighted that adding a negative PTO stiffness results in higher efficiency by detuning PTO response from initial resonance frequency. As the model developed works in frequency domain, it can be easily extended to a polychromatic and multidirectional excitation wave field with the exact same formulation by using a directional wave spectrum as an input. Using a numerical solver (as Finite Element Method) to determine structure shape functions could extend our approach to any carpet geometry and boundary conditions. The main limitation being the linearity assumption on the fluid and structural motion, next logical steps would be to take into account relevant nonlinear effects. Viscous effects, for instance, could easily be approached in frequency domain by considering a linear or quadratic viscous damping term derived from an estimated drag coefficient. Later, a time domain resolution would allow to increase order of nonlinearities, either in the fluid and solid domains. Acknowledgments JT and ND would like to thank the Cyclotron Road program at the Lawrence Berkeley National Lab for the partial support through DE-AC02-05CH11231 awarded by the Department of Energy. MRA would like to acknowledge a partial support provided under ABS Endowed Chair in Ocean Engineering. ND would like to thank the Agence National de la Recherche for their financial support for international mobility through the program ‘‘Investissements d’avenir’’ ANR-10-LABX-19-01. Authors would like to thank M. Lehmann, and T. Boerner for interesting and useful discussions throughout the project, and A. Babarit for technical support pertaining to the NEMOH package.

692

N. Desmars et al. / Journal of Fluids and Structures 81 (2018) 673–692

References Alam, M.-R., 2012a. A flexible seafloor carpet for high-performance wave energy extraction. In: Proc. of the 31st International Conference on Ocean, Offshore and Artic Engineering (OMAE). Alam, M.R., 2012b. Nonlinear analysis of an actuated seafloor-mounted carpet for a high-performance wave energy extraction. Proc. R. Soc. A 468, 3153– 3171. Andersen, M.T., Hindhede, D., Lauridsen, J., 2015. Influence of model simplifications excitation force in surge for a floating foundation for offshore wind turbines. Energies 8, 3212–3224. Antonutti, R., Peyrard, C., Johanning, L., Incecik, A., Ingram, D., 2014. An investigation of the effects of wind-induced inclination on floating wind turbine dynamics: heave plate excursion. Ocean Eng. 91, 208–217. Antonutti, R., Peyrard, C., Johanning, L., Incecik, A., Ingram, D., 2016. The effects of wind-induced inclination on the dynamics of semi-submersible floating wind turbines in the time domain. Renew. Energy 88, 83–94. Babarit, A., 2015. Mémoire d’habilitation à diriger des recherches, Tech. rep., École Centrale de Nantes. Babarit, A., Delhommeau, G., 2015. Theoretical and numerical aspects of the open source BEM solver NEMOH. In: Proc. of the 11th European Wave and Tidal Energy Conference (EWTEC2015). Babarit, A., Gendron, B., Singh, J., Mélis, C., Jean, P., 2015. Hydro-elastic modeling of an electro-active wave energy converter In: Proceedings of the ASME 2013 32nd International Conference on Ocean, Offshore and Artic Engineering. Behera, H., Sahoo, T., 2015. Hydroelastic analysis of gravity wave interaction with submerged horizontal flexible porous plate. J. Fluids Struct. 54, 643–660. Börner, T., Alam, M.R., 2015. Real time hybrid modeling for ocean wave energy converters. Ren. Sus. En. Rev. 84, 207–217. Chakraborty, R., Mandal, B.N., 2014. Scattering of water waves by a submerged thin vertical elastic plate. Arch. Appl. Math. 84, 207–217. Chan, A.T., Lee, S.W.C., 2001. Wave characteristics past a flexible fishnet. Ocean Eng. 28, 1517–1529. Chaplain, J.R., Heller, V., Farley, F.J.M., Hearn, G.E., Rainey, R.C.T., 2012. Laboratory testing of the anaconda. Phil. Trans. Roy. Soc. Lon. 370, 403–424. Cho, I.H., Kee, S.T., Kim, M.H., 1998. Performance of dual flexible membrane wave barriers in oblique waves. J. Waterway 124, 21–30. Cho, I.H., Kim, M.H., 1998. Interactions of a horizontal flexible membrane with oblique incident waves. J. Fluid Mech. 367, 139–161. Cho, I.H., Kim, M.H., 2000. Interactions of a horizontal porous flexible membrane with waves. J. Waterway, Port, Coastal, Ocean Eng. 126, 245–253. Choplain, N., 2012. Interactions of a Submerged Membrane with Water Waves and its use in Harnessing Nearshore Wave Power (Ph.D. thesis), University of Southampton. Cornett, A.M., 2008. A global wave energy resource assessment. In: Proc. of the 18th International Society of Offshore and Polar Engineers Conference. Crooks, D., Hoff, J.V., Folley, M., Elsaesser, B., 2016. Oscillating wave surge converter forced oscillation tests. In: ASME 2016 35th International Conference on Ocean, Offshore and Arctic Engineering. Delhommeau, G., 1987. Les problèmes de diffraction-radiation et de résistance de vagues : étude théorique et résolution numérique par la méthode des singularités (Ph.D. thesis), École Nationale Supérieure de Mécanique. Evans, D.V., 1976. A theory for wave power absorption by oscillating bodies. In: Proc. 11th Symp. on Naval Hydrodynamics, pp. 503–515. Falcao, A.F. de O., 2010. Wave energy utilization: A review of the technologies. Ren. Sus. En. Rev. 14, 899–918. Falnes, J., 2004. Ocean Waves and Oscillating Systems: Linear Interactions Including Wave-energy Extraction. Cambridge University Press. Falnes, J., 2007. A review of wave-energy extraction. Mar. Struct. 20, 185–201. Falnes, J., Hals, J., 2012. Heaving buoys, point absorbers and arrays. Phil. Trans. R. Soc. A 370, 246–277. Falnes, J., Kurniawan, A., 2015. Fundamental formulae for wave-energy conversion. Roy. Soc. Open Sci. 2, 140305. Farley, F.J.M., 1982. Wave energy conversion by flexible resonant rafts. Appl. Ocean Res. 4 (1), 57–63. Farley, F.J.M., 2012. Far-field theory of wave power capture by oscillating systems. Phil. Trans. R. Soc. A 278–287. Koley, S., Sahoo, T., 2016. Oblique wave scattering by horizontal floating flexible porous membrane. Meccanica 1–14. Lee, J.F., Chen, C.J., 1990. Wave interaction with hinged flexible breakwater. J. Hyd. Res. 28, 283–295. Mei, C.C., 2014. Nonlinear resonance in anaconda. J. Fluid Mech. 750, 507–517. Mihelcic, C.S., 1989. Hydrodynamic Force Coefficients of a Vertical Circular Cylinder (Ph.D. thesis), The University of British Columbia. Mork, G., Barstow, S., Pontes, M.T., Kabuth, A., 2010. Assessing the global wave energy potential. In: Proc. of the 29th International Conference on Ocean Offshore Mechanics and Artic Engineering. Newman, J., 1976. The interaction of stationary vessels with regular waves. In: Proc. 11th Symp. on Naval Hydrodynamics, pp. 491–501. Newman, J.N., 1979. Absorption of wave energy by elongated bodies. Appl. Ocean Res. 1 (4), 189–196. Newman, J.N., 1994. Wave effects on deformable bodies. Appl. Ocean Res. 16, 47–59. Newman, J.N., 2013. Cloaking a circular cylinder in water waves. Eur. J. Mech. B. Fluids (Enok Palm Memorial Volume) 47, 145–150. Olondriz, J., Elorza, I., Trojaola, I., Pujana, A., Landaluze, J., 2016. On the effects of basic platform design characteristics on floating offshore wind turbine control and their mitigation. J. Phys. Conf. Ser. 753, 052008. Parisella, G., Gourlay, T., 2016. Comparison of open-source code NEMOH with WAMIT for cargo ship motions in shallow water, Tech. rep., Center for Marine Science and Technology - Curtin University. Rainey, R.C.T., 2012. Key features of wave energy. Phil. Trans. R. Soc. A 425–438. Rao, S.S., 2007. Vibration of Continuous Systems. John Wiley & Sons. Schmitt, P., Windt, C., Nicholson, J., Elsässer, B., 2016. Development and validation of a procedure for numerical vibration analysis of an oscillating wave surge converter. Eur. J. Mech. B/Fluids 58, 9–19. Sheremet, A., Stone, G., 2003. Observations of nearshore wave dissipation over muddy sea beds. J. Geophys. Res. 108, 1–11. Thompson, G.O., 1991. Flexible Membrane Wave Barrier (Masters thesis), Oregon State University. Ul-Hassan, M., Meylan, M.H., Peter, M.A., 2009. Water-wave scattering by submerged elastic plates. Quart. J. Mech. Appl. Math. 62, 321–344. Zheng, Y.H., Shen, Y.M., Tang, J., 2007. Radiation and diffraction of linear water waves by an infinitely long submerged rectangular structure parallel to a vertical wall. Ocean Eng. 34, 69–82.