Hydrodynamic Analysis of Surge-Type Wave Energy Devices in Variable Bathymetry by Means of BEM

Abstract

A variety of devices and concepts have been proposed and thoroughly investigated for the exploitation of renewable wave energy. Many of the devices operate in nearshore and coastal regions, and thus, variable bathymetry could have significant effects on their performance. In particular, Oscillating Wave Surge Converters (OWSCs) exploit the horizontal motion of water waves interacting with the flap of the device. In this work, a Boundary Element Method (BEM) is developed, and applied to the investigation of variable bathymetry effects on the performance of a simplified 2D model of a surge-type wave energy converter excited by harmonic incident waves. Numerical results, illustrating the effects of depth variation in conjunction with other parameters, like inertia and power-take-off, on the performance of the device, are presented. Finally, a comparative evaluation of the present simplified surge-type WEC model and point absorbers is presented for a case study in a selected coastal site on the Greek nearshore area, characterized by relatively increased wave energy potential.

1. Introduction

During the last period, marine renewables have been the subject of intensive research worldwide, both by the industry and the academic community. More specifically, water waves provide a sustainable and widely available source of energy, which could be used to enhance the fraction of green energy in the global production. A great variety of Wave Energy Converters (WECs), devices and concepts based on oscillating systems have been proposed (e.g., see [1,2]), many of which have been also experimentally tested and are currently successfully contributing to the energy grid. In particular, Oscillating Wave Surge Converters (OWSCs) exploit the horizontal motion of water waves, interacting with the flap of the device, to generate renewable electricity [3,4,5]. A typical OWSC consists of a frame on which a shaft is mounted and a flap which performs angular oscillations, with the aforementioned shaft being the center of rotation. The frame is fastened at the seabed (or rests on an auxiliary structure at a lower depth) and is fitted with a hydraulic generator. Fluid is compressed into the hydraulic generator due to the angular oscillation, and the pressure is then transferred, by means of an appropriate pipeline system, to a hydroelectric power station to generate electricity. Methods for the analytical and computational modelling of such wave energy systems have been presented by various authors (e.g., see [6]).

As regards the devices that operate in nearshore and coastal regions, effects caused on their response due to finite water depth bathymetric inhomogeneities could be significant, and need to be considered; e.g., see [7]. This has been illustrated in the case of point absorbers operating in variable bathymetry regions, for which models based on the Boundary Element Method (BEM) have been developed, e.g., [8,9], where it is shown that bottom slope and curvature could have an effect on the performance of the systems, especially in shallow water conditions.

In this work, a novel BEM scheme is developed and applied to the investigation of depth inhomogeneity effects on the performance of surge-type wave energy converters excited by harmonic incident waves, studied in the context of linearized water-wave theory. For simplicity, we restrict ourselves to a two-dimensional flap-type configuration and normally incident waves. Although the present simplified 2D BEM model does not accurately treat the 3D effects of an OWSC like the Oyster device, as extensively discussed by Renzi et al. [4], it still provides useful information, especially as concerns the variable bathymetry effects on the performance of the device. The terminator-type model is studied in two dimensions, where the flap of the device is aligned perpendicular to the incident wave direction (e.g., see [4]). As well as the standard absorbing–radiating model depicted in Figure 1a, which includes additional effects due to wave transmission in the downwave region behind the flap-board, an absorbing–reflecting model, shown in Figure 1b, is also considered. The present numerical method is first validated by comparison to an analytical solution, derived in constant water depth and presented in detail in the Appendix A. Although the latter absorbing–reflecting configuration could be represented by the former absorbing–radiating model, with appropriate adjustment of the Power Take Off (PTO) parameters to achieve zero transmitted wave, still it is considered here as a special case, and used for the additional validation of the developed BEM model.

Subsequently, the 2D BEM model is applied to illustrate the effects of depth variation, including sloping and corrugated seabeds, in conjunction with other parameters, like inertia and power-take-off coefficients. Results are used to show and discuss the response of the considered system in inhomogeneous regions. Although nonlinearity associated with large wave heights, wave breaking and viscous fluid effects are significant, a linear model could still provide useful information for the performance and optimization of the devices, as well as estimates of power performance measures that can be derived (e.g., see [10]).

In the last section of this work, a case study is presented providing useful data concerning the available wave energy at a nearshore/coastal site with relatively increased potential in the Aegean Sea. Using the latter data, a comparative evaluation of the present simplified surge-type WEC model and point absorbers is presented, deriving power performance measures that can be further exploited to levelized energy cost analysis. A discussion concerning the validity and shortcomings of the present simplified analysis is also included, especially as concerns the three dimensional and real fluid effects, which will be treated in future work.

2. Hydrodynamic Analysis of A Two-dimensional Flap-Type Model in Variable Bathymetry

In this work, a simplified two-dimensional BEM model is presented for the hydrodynamic analysis of a flap-type device, in the case of normally incident waves in variable bathymetry, shown in Figure 1, in the context of linearized wave theory. In all cases, a zero-thickness flap boundary is considered.

Without loss of generality, a left horizontally semi-infinite domain is considered, $D=\{x_{\alpha }<x,-h\left(x\right)<z<0\}$ (see Figure 1), starting from x_a = 0 where the flap-type device is located. In order to study variable bathymetry effects on the performance of the two-dimensional flap-type model (2dOSWC), a low-order Boundary Element Method (BEM) is developed (e.g., see [11]). The method is first validated by comparison between the numerical and the analytical results in constant water depth, for which an analytical solution is described in Appendix A. Subsequently, the numerical model is used in order to show and discuss the variable bathymetry effects.

The present panel method is based on the approximation of the boundary by linear elements/panels carrying simple source–sink distributions. We assume harmonic time dependence in the form $\exp (-i\omega t)$, and the complex wave potential is represented by:

$$\phi (x)={\displaystyle \underset{\partial D}{\int }\sigma (x_0)G\left(x|x_0\right)}dS(x_0)$$

(1)

where $G\left(x|x_0\right)={\left(2\pi \right)}^{-1}ln\Vert x-x_0\Vert$ denotes the fundamental solution of the Laplace Equation in the two-dimensional space, x(x,z) is a general field point in the domain, while x₀(x,z) indicates the singularity point, and $\partial D=\partial D_1\cup \partial D_2\cup \partial D_3\cup \partial D_4$ denotes the boundary of the domain decomposed into four sections, as presented in Figure 1: the free surface $\partial D_1$, the device part $\partial D_2$, the seabed $\partial D_3$ and the vertical interface $\partial D_4$ serving as an incidence–radiation boundary.

The solution to the potential flow problem is obtained by discretizing the boundary geometry into elements $\partial D={\displaystyle \cup E_j}$, in conjunction with the approximation of the unknown source–sink function—which is assumed to be piecewise constant on each element—and enforcing the boundary conditions on the control points of each boundary’s part (defined as the elements’ midpoints), from which the source–sink singularity strength on each element $({\sigma }_j)$ is determined. Under the assumption that the variable bathymetry domain, shown in Figure 1, ends in a semi-infinite constant-depth strip, and that the interface $\partial D_4,$ located at $x=x_b$, lies well within the constant-depth region, where the effects of evanescent modes have died out, the complex wave potential consisting of the incident and reflected wave can be expressed by considering the propagating mode only, as follows:

$$\mathsf{\Phi }(x,z)=Q\phi (x,z),\phi (x,z)=\left[\exp (-i{k}_{0}x)+R\exp (i{k}_{0}x)\right]Z_0^{(b)}(z),x\ge x_b,$$

(2)

where $Q=-igH/\left(2\omega \right)$, with $g$ denoting the acceleration of gravity and H the incident wave height. Moreover, R is the reflection coefficient, the function $Z_0^{(b)}(z)=\cosh \left(k_0(z+h_b)\right)/\cosh \left(k_0{h}_b\right)$ and the wavenumber $k_0$ of the incident-reflected wave is obtained from the dispersion relation formulated at the radiation boundary depth $h_b$:

$${\omega }^2=k_{0}g\tanh \left(k_0{h}_b\right)$$

(3)

Using Equation (2), the horizontal wave velocity is calculated by the partial x-derivative of the normalized complex wave potential φ at $x=x_b$, as follows:

$${\frac{\partial \phi }{dx}|}_{x_b}=\left[i{k}_{0}R\exp (i{k}_0{x}_b)-i{k}_0\exp (-i{k}_0{x}_b)\right]Z_0^{\left(b\right)}(z)$$

(4)

The vertical eigenfunctions $\{Z_n^{(b)}(z),n=0,1,2\dots \}$ are derived from the solution of the Sturm–Liouville problem, to which the Laplace equation is reduced in the constant depth strip $-h<z<0,x\ge x_b.$ The latter are $Z_0^{\left(b\right)}\left(z\right),Z_n^{\left(b\right)}\left(z\right)=\cos \left(k_n\left(z+h_b\right)\right)/\cos \left(k_n{h}_b\right),n=1,2,\dots$, with the corresponding parameters $i{k}_n,n=1,2,\dots ,$ obtained as the imaginary roots of Equation (3). By projecting the two members of Equation (2) in the above orthogonal vertical basis we obtain:

$$\langle \phi (x_b,z),Z_0^{\left(b\right)}\rangle =\left(R\exp (i{k}_0{x}_b)+\exp (-i{k}_0{x}_b)\right){\Vert Z_0^{\left(b\right)}\Vert }^2$$

(5)

where $\langle f_1\left(z\right),f_2\left(z\right)\rangle ={\displaystyle \underset{z=-h_b}{\overset{z=0}{\int }}{f}_1\left(z\right)f_2\left(z\right)}dz,\mathrm{and}{\Vert f\Vert }^2=\langle f\left(z\right),f\left(z\right)\rangle$. Consequently, Equation (5) can be solved for R, providing an expression of the reflection coefficient as a function of the complex wave potential in the water column at $x=x_b$:

$$R=\frac{\langle \phi (x_b,z),Z_0^{\left(b\right)}\rangle }{{\Vert Z_0^{\left(b\right)}\Vert }^2\exp (i{k}_0{x}_b)}-\exp (-2i{k}_0{x}_b)$$

(6)

Substitution of the reflection coefficient’s expression in Equation (4) results in the following radiation condition, that applies on the interface (radiation boundary):

$$\frac{\partial \phi }{\partial n}-i{k}_0\frac{\langle \phi (x_b,z),Z_0^{\left(b\right)}\rangle }{{\Vert Z_0^{\left(b\right)}\Vert }^2}{Z}_0^{\left(b\right)}(z)=-2i{k}_0\exp (-i{k}_0{x}_b)Z_0^{\left(b\right)}(z),\mathrm{on}\partial D_4$$

(7a)

where $\partial \phi /\partial n$ denotes the normal derivative of the wave potential on the boundary, with respect to the outward directed unit normal vector $n$. The boundary conditions that apply on the other three sections of the boundary are:

$$\frac{\partial \phi }{\partial n}-\frac{{\omega }^2}{g}\phi =0,\mathrm{on}\partial D_1,\frac{\partial \phi }{\partial n}=f(z),\mathrm{on}\partial D_2\frac{\partial \phi }{\partial n}=0,\mathrm{on}\partial D_3$$

(7b)

where the functions $f(z)=i\omega {\theta }_0\left(z+d\right),-d<z<0$ and $f(z)=0,-h_a<z<-d,$ represent the oscillatory normal component of the velocity due to the flap motion, and are simply zero on the wall below the pivot point of the OWSC, located at z = −d (see Figure 1).

2.1. Mesh Generation

The boundary is discretized by a closed polygonal line, which is formed by distributing nodes across its length. In particular, the free surface nodes are determined by specifying a number of nodes per characteristic wavelength in the domain that is usually larger than 30–40. Next, the free surface nodes are redefined, using a cosine spacing, to allow them to be more densely arranged approaching the corners of the domain. The z-coordinates of the free surface’s nodes are equal to $z_i=0,i=1,M_1+1$, where $M_1$ denotes the number of elements on the free surface. The bottom boundary is discretized accordingly, with the z-coordinates set equal to $z_i=-h\left(x_i\right),i=M_1+M_2+1,M_1+M_2+M_3+1$, where $M_3$ is the number of elements on the seabed, $h\left(x\right)$ is the depth function, and $M_2,M_4$ are the number of elements on the vertical boundaries (OWSC and radiation boundary, respectively). The node spacing is finer on the horizontal boundaries near the corners of the domain, while it is uniform on the vertical sides. The latter is selected to match the fine spacing of the horizontal sides at the corners. Based on the above discretization, the total number of elements is $M=M_1+M_2+M_3+M_4$.

2.2. Discrete BEM Model

By assuming unit source–sink singularity strengths on each element, the induced potential and velocity matrices ${\phi }_{kj},U_{kj},k,j=1,\dots ,M$, from each element to each collocation point, are analytically calculated; see [11]. Using the latter quantities in the boundary conditions, the influence matrix $A_{kj}$ of the discrete system is obtained, as explained in more detail in the sequel. By denoting with ${\sigma }_j$ the strength of the source/sink of the j-element, and by $b_k$ the boundary data on the control point defined as the center of the k-element, the resulting linear system is:

$${\displaystyle \sum _{j=1}^M{A}_k{}_j}{\sigma }_j=b_k j,k=1,2,\dots ,M$$

(8)

Apart from the discrete source singularities’ strengths $\{{\sigma }_j,j=1,M\}$, the additional unknown parameters are the complex amplitude ${\theta }_0$ of the OWSC’s angular motion, calculated as the response of the system, as well as the reflection coefficient R. Given that the reflection coefficient has been expressed as a function of the complex wave potential Equation (6), the extra unknown is the complex oscillation amplitude only. The additional equation required to close the linear system is obtained from the equation of motion of the angular oscillator, expressed in the frequency domain (see also Equation (A14) in the Appendix A).

In order to construct the BEM model, it is necessary to express the boundary conditions in their discrete form. If we denote the unit normal vector on the k-boundary element (directed outwards) by n_k, the boundary condition Equation (7a) on the radiation boundary is expressed as follows:

$$\begin{array}{l}{\displaystyle \sum _{j=1}^M{\sigma }_j}(n_k{U}_k{}_j)-i{k}_0\left\{{\displaystyle \underset{-h_b}{\overset{0}{\int }}{\displaystyle \sum _{j=1}^M{\sigma }_j}{\phi }_{\alpha j}{Z}_0^{\left(b\right)}(z_a)}\right\}{Z}_0^{\left(b\right)}(z_k)=b_k,\\ k=M_1+M_2+M_3+1,\dots M\end{array}$$

(9a)

where

$$b_k=-2i{k}_0\exp (-i{k}_0{x}_b)Z_0^{(b)}\left(z_k\right)$$

(9b)

are the data associated with the incident wave on the incidence–radiation interface $\partial D_4$, and $z_k$ denotes the collocation points, located at the center of the k-element on $\partial D_4$. The auxiliary index α, introduced in Equation (9), spans the interval of integers that correspond to the elements distributed across ∂D₄. Therefore, if δz₄ is the constant length of the boundary elements on the radiation boundary, the integral of the above equation can be replaced by the relation arising from applying a numerical integration rule. Thus, by applying the trapezoidal rule, we get:

$$\begin{array}{l}{\displaystyle \sum _{j=1}^M\left\{(n_k{U}_k{}_j)-i{k}_0\left[{\displaystyle \sum _{a=M_A+1}^M{\sigma }_j{\phi }_{\alpha j}\delta z_4{Z}_0^{(b)}(z_a)}\right]Z_0^{(b)}(z_k)\right\}}=b_k,\\ k=M_A+1,\dots M,M_A={\displaystyle \sum _{i=1}^3{M}_i}\end{array}$$

(10)

Moreover, the boundary condition on the boundary of the OWSC (see also Appendix A) is:

$$\{\begin{array}{c}\frac{\partial \phi }{\partial n}+i\omega {\theta }_0(z+d)=0,-d<z<0\\ \frac{\partial \phi }{\partial n}=0,-h_0<z<d\end{array}$$

(11)

Using the matrix-coefficients for the induced potential and velocity, the above Equation is expressed in discretized form as follows:

$${\displaystyle \sum _{j=1}^M{\sigma }_j}(n_k{U}_k{}_j)+H(z_k+d)\frac{i\omega (z_k+d)}{Q}{\theta }_0=0,k=M_1+1,\dots M_1+M_2$$

(12)

where H denotes the Heaviside step function. Furthermore, the discretized boundary conditions on the free surface and the impermeable bottom, respectively, are expressed as follows:

$${\displaystyle \sum _{j=1}^M{\sigma }_j}(n_k{U}_k{}_j-\frac{{\omega }^2}{g}{\phi }_k{}_j)=0, k=1,\dots M_1$$

(13)

$${\displaystyle \sum _{j=1}^M{\sigma }_j}(n_k{U}_k{}_j)=0, k=1+{{\rm M}}_1+{{\rm M}}_2,\dots {\displaystyle \sum _{i=1}^3{M}_i}$$

(14)

Concerning the additional Equation closing the discrete system, with respect to the additional unknown ${\theta }_0$, the dynamic Equation governing the angular harmonic oscillator that models the present 2dOWSC, expressed in the frequency domain, is considered:

$$\alpha (\omega ){\theta }_0=M_0-B{M}_B,\mathrm{where} \alpha (\omega ){\theta }_0=\left(-{\omega }^2{J}_P-i\omega B_P+C_P\right){\theta }_0$$

(15)

$$\begin{gathered} M_0=-i\omega \rho {\displaystyle {\int }_{-d}^0(z+d)\phi \left(x=0+,z\right)}dz \mathrm{and} \\ {M}_B=-i\omega \rho {\displaystyle {\int }_{-d}^0(z+d)\phi \left(x=0-,z\right)}dz \end{gathered}$$

(16)

and the coefficient B = 1 for the standard absorbing–radiating model, and $B=0$ in the special case of an absorbing-reflecting model (see also Appendix A). In the above equations, $M_0$ is the moment generated by wave pressure forces on the upwave side of the flap, and $M_B$ on the downwave side, respectively. Furthermore, $J_P$, $B_P$, $C_P$ are the coefficients of inertia, damping (representing the power-take-off system) and rotational spring constant, respectively; see the analysis of Appendix A, Equation (A14). Using Equations (15) and (16) we obtain:

$$\alpha (\omega ){\theta }_0+\rho Q{\displaystyle \underset{-d}{\overset{0}{\int }}\left(z+d\right)\left(\phi \left(x=0+,z\right)-B\phi \left(x=0-,z\right)\right)dz}=0$$

(17a)

Using the representation Equation (A5b) of the wave potential in the downwave side of the device to express $\phi \left(x=0+,z\right)$, and setting the water depth $h=h_a$, in conjunction with the expressions of the involved coefficients by Equations (A12b), the above equation can be put in its discrete form, by assuming that $\delta z_2$ is the common length of the boundary elements on $\partial D_2$, and the result is

$$i\omega {\displaystyle \sum _{j=1}^M\left\{{\displaystyle \sum _{a=1+M_1}^{M_1+M_2}H\left(z_a+d\right){\phi }_{aj}\delta z_2(z_a+d)}\right\}}{\sigma }_j+\left[\frac{a(\omega )}{Q\rho }+B\frac{{\omega }^2\mathsf{\Lambda }}{Q}\right]{\theta }_0=0$$

(17b)

where the frequency-dependent coefficient $\mathsf{\Lambda }\left(\omega \right)$ is defined as follows:

$$\mathsf{\Lambda }=\frac{O_0^2}{k_{0}cos{h}^2(k_{0}h){\Vert Z_0\Vert }^2}-{\displaystyle \sum _{n=1}^{\infty }\frac{i{O}_n{}^2}{k_n{[cos(k_{n}h)]}^2{\Vert Z_n\Vert }^2}}$$

(17c)

In conclusion, the discrete BEM model consists of Equations (10), (12)–(14) and (17). After solving the linear system and obtaining the strengths of the unknown source–sink distribution $\{{\sigma }_j,j=1,\dots ,M\}$, as well as the complex oscillation amplitude ${\theta }_0$ of the present 2dOWSC, the reflection coefficient can be calculated from Equation (6) as:

$$R=\frac{{\displaystyle \sum _{a=M_A+1}^M\left({\displaystyle \sum _{j=1}^M{\sigma }_j{\phi }_{aj}\delta z_4{Z}_0(z_a)}\right)}}{{\Vert Z_0\Vert }^2\exp (i{k}_0{x}_b)}-\exp (-2i{k}_0{x}_b)$$

(18)

In the case of the open-type system shown in Figure 1a, the transmission coefficient is calculated by means of Equation (A12b). Finally, the performance of the system is calculated by

$${\eta }_{OWSC}=\frac{P_{out}}{P_{in}}=\frac{4{\omega }^2{B}_p{\left|{\theta }_0\right|}^2}{\rho g{H}^2{c}_g}$$

(19a)

and the consideration of energy-flux conservation leads to the following equation:

$${\eta }_{OWSC}=1-{\left|R\right|}^2-B{\left|T\right|}^2\left(C_{ga}/C_{gb}\right)$$

(19b)

where $C_{ga}$ denotes the group velocity at depth $h_a$ in the far upwave incident region, and $C_{gb}$ the group velocity at depth $h_b$ at the downwave region, which could be used for verifying the conservation of energy of the studied hydromechanical system in variable bathymetry regions.

2.3. Validation of the BEM Model

Applying the present BEM to constant depth regions allows the validation of the numerical method by comparing the results to the analytical solution, described in Appendix A. As an example, we consider the studied system, operating in incident waves of height H = 1 m in constant water depth equal to 10 m, with the hinge point located one meter above the seabed (d = 9m). The non-dimensional PTO parameters (see Appendix A) are set as:${\widehat{J}}_P=1,{\widehat{B}}_P=1.5,{\widehat{C}}_P=1.5$, and the flap’s density is considered as equal to the water’s density.

Figure 2a shows the complex amplitude of the device’s oscillation as a function of the non-dimensional frequency of the incident wave, as calculated both by the analytical and the numerical methods, while Figure 2b illustrates a comparison between the analytical and the numerical approach, regarding the performance of the device as well as the reflection and transmission coefficients. Both the standard absorbing–radiating model and the special case of the absorbing–reflecting model have been considered, and the results concerning the angular response and the reflection–transmission characteristics are plotted in the upper and lower subplots of Figure 2, respectively. The numerical results refer to a 200-m-long domain. The number of wavelengths contained in the domain is a function of the frequency (and the water depth), and thus, is recalculated for each frequency. In order to derive good quality convergent results, the nodes per wavelength, distributed across the free surface and the bottom boundaries, are set to a minimum of 70. In addition, a minimum of 1000 nodes is selected for discretizing the parts ∂D₁ and ∂D₃ of the boundary, in order to avoid coarse meshing at low frequencies (long wavelengths). It is noted that, for uniform bathymetry cases, altering the domain’s length for every iteration/frequency, so that it matches a certain number of wavelengths, does not affect the results. As seen in Figure 2, the differences between the two solutions are negligible. Τhe convergence of the numerical approach can lead to virtually zero differences between the two solutions by using finer meshes of the domain’s boundary. Finally, it is worth noticing that conservation of energy, expressed by Equation (19), is perfectly satisfied for both examined types and for all frequencies.

2.4. Performance of OWSC in Variable Bathymetry

In this section, we present results concerning the calculation of device performance in regions of non-uniform bathymetry, characterized by a depth function h(x). For comparison purposes, all other parameters of the studied 2dOWSC remain the same as in the example previously presented in Figure 2. We first consider a shoaling environment as presented in the top subplots of Figure 3. The variable bathymetry subdomain is defined between two subregions of constant but different depth: the region of wave incidence with depth $h_b=$ 25 m, and the vicinity of the device with depth $h_a=$ 10 m, that remains constant in the downwave region, behind the flap. The water depth in the upwave subregion is defined by the function:

$$z=h(x)=\frac{1}{2}(h_a+h_b)-\frac{1}{2}(h_a-h_b)\tanh \left[\kappa (x-x_m)\right]+\delta h_c\left(x\right), x\in [0,x_b]$$

(20)

which describes a smooth shoal with maximum bottom slope controlled by the parameter $\kappa$, while $\delta h_c\left(x\right)$ is a function of finite support, modelling additional bottom corrugations. The parameter $x_m$ denotes the position of the maximum bottom slope, selected here to coincide with the middle x-coordinate $\left(x_b/2\right)$ of the domain. The seabed profiles shown in the upper left and right subplots of Figure 3 are determined, for $\delta h_c(x)=0$, by setting the parameters as $\kappa =0.03,x_b=400 \mathrm{m}$ and $\kappa =0.09,x_b=200 \mathrm{m}$, respectively.

In particular, in Figure 3b the calculated real and the imaginary part of the complex amplitude ${\theta }_0$, as well as its absolute value, are plotted as a function of the non-dimensional frequency, while the graphs in Figure 3c depict the performance of the device and the reflection coefficient (R) and transmission coefficient (T) as functions of the same parameter. The water depth at the device’s boundary is used for nondimensionalization.

A maximum value of performance about 39% is observed for nondimensional frequency $\omega \sqrt{h_a/g}\approx 1.2$. According to the results shown in the above graphs, the maximum angle achieved during the oscillation is not substantially affected by changes in bathymetry (see Figure 2). On the other hand, the angular responses, especially in the low frequencies, differ significantly compared to the corresponding predictions at constant water depth, due to strong the interaction of the incident waves with the seabed. The observed oscillations in the numerical results are due to variations in the phase of various quantities excited by depth irregularities. By comparing the results presented in Figure 3 to the ones presented in the right column of Figure 2, we observe that as the mean bottom slope reduces, the responses and the performance converge to the results obtained at constant depth. On the other hand, as the bottom slope increases, there is a strong effect on the flap’s oscillation phase, particularly at lower frequencies, due to stronger wave–seabed interaction.

On the contrary, abrupt changes occurring over a short horizontal distance less affect the flap’s oscillation phase, but have a stronger effect on the performance curve, in the bandwidth of frequencies in which there is interaction of the incident waves with the seabed. Βy using an auxiliary grid spanning the entire two-dimensional domain, the distribution of the velocity potential can be calculated by the superposition of the potential induced from each one of the boundary elements to each point of the auxiliary grid. After obtaining the values of the complex wave potential on every point of the 2D grid, the induced field can be represented at any moment (t) by the real part of φ multiplied by $\exp (-i\omega t)$. Figure 4 illustrates the wave field during the OWSC’s operation in a non-uniform bathymetry region, whose seabed profile is defined by Equation (20) for $\delta h_c(x)=0,\kappa =0.09,x_b=200 \mathrm{m}$, by means of equipotential lines. For the specific period of 8 s ($\omega =0.785 \mathrm{rad}/\mathrm{s}$), the resulting reflection coefficient is R = 0.8287, and thus the performance of the device equals P = 31.32%. It can be observed that the wave field resembles an almost stationary wave for fractions of the period, and propagates gradually depending on the flap’s oscillation phase. Furthermore, it can be seen that the wave potential’s field satisfies the Neumann boundary condition Equation (14) on the seabed boundary, which is equivalent to the fact that the equipotential lines intersect the impermeable bottom perpendicularly.

Apart from seabed profiles characterized by depth changes described by hyperbolic tangent functions, the present BEM is applicable to more realistic bathymetric variations which can be determined by considering bottom corrugations (in addition to shoaling effects). Such inhomogeneities can be determined by, among other ways, adding a damped sine or cosine term to the function describing the seabed. Before applying the method, one must ensure that the parameter selection leads the additional term completely dying out, in both the areas of the two vertical boundaries. As an example, we consider in Figure 5a the depth profile [Equation (20)], with $\kappa =0.09,x_b=200 \mathrm{m}$ and bottom corrugations $d{h}_c(x)=3\sin [0.15(x-x_m)]\exp [0.0006{(x-x_m)}^2].$ For comparison purposes, all other parameters remain unchanged. The plot in Figure 5b illustrates the real and the imaginary part of the complex amplitude, as well as its absolute value, as a function of the nondimensional frequency, while Figure 5c of the figure depicts the performance of the device, and the reflection (R) and transmission (T) coefficients, as a function of the same parameter. The water depth at the device’s boundary is used for nondimensionalization. Ιn relation to the corresponding graphs of Figure 3, it can be observed that the current alteration of the seabed profile introduces a significant difference to the device’s performance in the bandwidth of the nondimensional frequency between 0.5 and 1. The wavelengths falling into this bandwidth are 60 m to 170 m long, and thus result in strong interaction with the seabed.

Figure 6 illustrates the calculated wave field during the absorbing–reflecting termination type OWSC operation in the non-uniform bathymetry region, whose seabed profile is defined by Equation (20) for the aforementioned parameter selection, by means of equipotential lines. For the specific period of 8 s, the resulting reflection coefficient is R = 0.8455, and thus the performance of the device equals P = 28.50%. The excellent satisfaction of the bottom boundary condition in the examined irregular seabed topography is clearly observed.

3. WEC Model

In this section, the BEM model developed by the authors of [8,9] is briefly reviewed. This model will be used to obtain numerical results concerning the performance of point absorbers, and present comparative results concerning an application to a nearshore/coastal site in Greece that will be presented in the next section. The point absorber is modelled by an axisymmetric floating body with radius a and draft T, operating in a coastal environment. As before, the bathymetry corresponds to the shoaling seabed profile with depth smoothly ranging from $h=h_1$ in the region of incidence, to $h=h_3$ in the region of transmission. As excitation of the floating device is assumed a harmonic wave of angular frequency ω and incident direction. For a relatively small wave slope, remaining within the limits of the linear theory, the wave potential and free surface elevation are:

$$\mathsf{\Phi }\left(x,z;t\right)=\mathrm{Re}\left\{-\frac{igH}{2\omega }\phi \left(x,z;\mu \right)\cdot \exp \left(-i\omega t\right)\right\},\eta \left(x;t\right)=\mathrm{Re}\left\{\frac{H}{2}\phi \left(x,z;\mu \right)\cdot \exp \left(-i\omega t\right)\right\}$$

(21)

where $x=\left(x,y\right)$ denotes the horizontal coordinates and $\mu ={\omega }^2/g$ is the frequency parameter. The wave potential is decomposed on the propagating wave potential ${\phi }_P\left(x,z\right)$, with the absence of the floating body, the diffraction potential ${\phi }_D\left(x,z\right)$ due to the presence of the impermeable and motionless geometry of the WEC and the radiation potential ${\phi }_D\left(x,z\right)$, which is the summation of the radiating fields due to the six modes of oscillation of the WEC:

$$\phi \left(x,z\right)={\phi }_P\left(x,z\right)+{\phi }_D\left(x,z\right)+\frac{2{\omega }^2}{gH}{\phi }_R\left(x,z\right){\phi }_R\left(x,z\right)={\displaystyle \sum _{\ell =1}^6{\xi }_{\ell }}{\phi }_{\ell }\left(x,z\right)$$

(22)

On the wetted surface of the WEC, the following boundary conditions are applied:

$$\partial {\phi }_D\left(x,z\right)/\partial n=-\partial {\phi }_P\left(x,z\right)/\partial n,\partial {\phi }_{\ell }\left(x,z\right)/\partial n=n_{\ell },\ell =1,2,\dots 6,$$

(23)

where $n=\left(n_1,n_2,n_3\right)$ is the normal vector directing inwards the geometry of the floater, and $n_{\ell }$ is the $\ell$ component of the normal vector in generalized form. Considering that heaving is the only power mode of the WEC, its response in the vertical oscillation can be obtained as:

$${\xi }_3=\left(X_P+X_D\right)/A\left(\omega \right)$$

(24)

In the above expression of heave response amplitude operator-ξ₃ (RAO ξ₃), the implemented Froude–Krylov X_p and diffraction X_D exciting forces are evaluated by integration of the wetted surface of the body $\partial D_B$, as follows:

$$X_P=\frac{\rho gH}{2}{\displaystyle \underset{\partial D_B}{\iint }{\phi }_P{n}_3}dS,X_D=\frac{\rho gH}{2}{\displaystyle \underset{\partial D_B}{\iint }{\phi }_D{n}_3}dS$$

(25)

while $A\left(\omega \right)$ is calculated by:

$$A\left(\omega \right)=-{\omega }^2\left(M+a_{33}\right)-i\omega \left(B_S+b_{33}\right)+\left(C_S+c_{33}\right)$$

(26)

The hydrodynamic coefficients of added mass (α₃₃) and damping (b₃₃) in Equation (26) are obtained by integrating the heave-radiation potential on the submerged boundary surface of the WEC:

$$a_{33}-\frac{1}{i\omega }{b}_{33}=\rho {\displaystyle \underset{\partial D_B}{\iint }{\phi }_3{n}_3}dS$$

(27)

In addition, $c_{33}$ = ρgA_WL is used for the hydrostatic coefficient and A_WL refers to the waterline surface. Finally, Bs and Cs are the damping and spring coefficients of the involving linear PTO system, respectively. The expression for the time-averaged WEC power output of a WEC in heave power mode follows:

$$P\left(\omega ,\theta \right)=\frac{1}{2}{\eta }_{WEC}^{}{\omega }^2\left|B_S{\left({\xi }_3\right)}^2\right|$$

(28)

where ${\eta }_{WEC}^{}$ is the PTO efficiency. The amount of the available wave power that can be harnessed by the device depends on the directional wave spectrum and PTO characteristics; see [8,9].

Evaluation of Wave Field

The total wave field, as already described, can be decomposed into mainly three quantities: the propagation field, the diffracted field and the radiated field. The Coupled Mode Model, developed in [12] and extended to 3D in [13], is used for the evaluation of the propagation field, over bathymetric topographies with even steep depth steps, where analytic solutions are infeasible, ensuring the energy conservation at the same time. The main formula is based on the following local-mode representation, with the enhancement of correction terms in order to globally satisfy the boundary condition of the rigid bottom:

$${\phi }_P\left(x,z\right)={\phi }_{-1}\left(x\right)Z_{-1}\left(z;x\right)+{\displaystyle \sum _{n=0}^{\infty }{\phi }_n\left(x\right)}{Z}_n\left(z;x\right)$$

(29)

The vertical functions Z_n (z; x) are obtained as eigenfunctions of regular Sturm–Liouville problems, formulated at the local depth, while the functions ${\phi }_n\left(x\right)$ are the complex amplitude of the n^th-mode, and are found as the solution of the following Coupled Mode System (CMS)

$${\displaystyle \sum _{n=-1}^{}{A}_{mn}\left(x\right){\nabla }^2{\phi }_n\left(x\right)}+B_{mn}\left(x\right)\nabla {\phi }_n\left(x\right)+C_{mn}\left(x\right){\phi }_n\left(x\right)=0,m=-1,0,1\dots$$

(30)

where the coefficients $A_{mn},B_{mn},C_{mn}$, defined in terms of the vertical eigenfunctions, are described in Table 1 of Reference [13]. The CMS is also supplemented with boundary conditions to treat incidence, reflection and transmission phenomena over regions of general bathymetry.

For the evaluation of the diffraction and radiation potentials, a BEM method, developed in [14] and described in [8,9], is used. More specifically, the induced potential and velocity from the 4-node quadrilateral elements, used for the discretization of all the boundaries (body, free surface and seabed surface), is expressed as:

$$\phi \left(r\right)={\displaystyle \sum _p{F}_p}{\mathsf{\Phi }}_p\left(r\right),\nabla \phi \left(r\right)={\displaystyle \sum _p{F}_p}{U}_p\left(r\right)$$

(31)

where summation refers to all the quadrilateral elements, and Φ_p and U_p are used to denote, respectively, the potential and velocity induced from the p^th element of unit singularity distribution to a random field point $r=\left(x,y,z\right)$. These quantities, induced from each element, can be evaluated by semi-analytical expressions, and the final discrete solution is obtained using the collocation method, aiming to satisfy the boundary conditions at each panel’s centroid.

The aforementioned BEM, used to treat the hydrodynamic problem of a point absorber-type WEC, even in cases of operation over arbitrary bathymetric regions, is developed, validated against analytical solutions whenever possible, and used to obtain useful information regarding the response and power output of the devices, and assist the research towards optimized WECs and efficiency maximization, assuming cases of one-device installation or array deployments; see [8,9].

4. Application to A Specific Coastal Site in Greece

Greece is located at the eastern Mediterranean Basin, and its energy sector is still largely dependent on fossil fuels (coal and oil). On the other hand, the coastline extents to 16,000 km, and therefore marine renewable energy could be a contributor to the power production, especially as concerns the supply of islands. The Aegean Sea region is characterized by high wind potential, which gives rise to relatively intense wave activity. However, the wind fetches are not very long, due to the presence of island complexes. Nonetheless, channeling effects cause certain spots to develop high wave potential [15,16].

Estimation of the mean annual wave power in the Greek Sea region has been presented by various authors, based on calibrated hindcast wave data (e.g., see [15,16]). A nearshore hot spot which is exposed to relatively long wind fetches, corresponding to increased wave potential, is located at the nearshore area of the south-eastern part of Evia Island, in the central Aegean Sea. This region is highlighted in Figure 7.

Wave energy devices are usually located relatively close to the shoreline. Hence, the vast majority of wave energy applications rely on the prediction of wave power availability in nearshore areas. Nonetheless, the data provided by global forecasting models, satellite or buoy measurements are available for off-shore regions. In order to obtain the near-shore wave conditions, the transformation of offshore sea states is required. This is achieved by using wave models for the transformation of waves using offshore wind and wave data, in conjunction with geographical (bathymetric and coastline) information in the examined area (e.g., see [17,18]). The geographical area considered in this study is enclosed by the box [24°10′E, 24°40′E]–[38°05′N, 38°30′N]. The area is shown in the map of Figure 8a, along with the nearest offshore data points (O1, O2 and O3). In a previous work [18], five representative nearshore target points of the specific area were studied in order to transform the offshore sea states and obtain the nearshore wave potential. The detailed analysis showed that the nearshore area around TP5 (North side of Cape Kafireas) is where the available wave potential is maximized. This site is furthermore suitable for the installation of a wave energy park, since it is not very close to villages, and there is the possibility of grid connection. The target points (TP) are shown in Figure 8a using circles. A nautical map illustrating the seabed’s morphology in the area, by means of isobath lines, is presented in Figure 8b.

The wave data used in the study are available from the EUROWAVES database, in the form of time-series covering the 8-year period of 1992–1999, regularly sampled at 6-hour intervals, with a total sample size of 10,189 [18]. The overall mean value of the available wave potential in the offshore area (point O2) is equal to 7.2 kW/m. The calculated nearshore wave climatology is presented in Figure 9, and details concerning the analytical probability models used for univariate and multivariate data can be found in [19,20]. Based on the nearshore wave data, the annual distribution of the wave potential for the point TP5 is shown in Figure 10a, and the mean value is dropped to 3.86 kW/m. The mean monthly distribution is also presented in Figure 10b. The maximum value of available wave power (≈6 kW/m) is observed in winter, while it becomes lowest during spring (≈2 kW/m).

4.1. Estimation of Power Output by Application of OWSC

We next consider an OWSC device positioned along the 20-m isobath in the examined nearshore/coastal area. Due to the fact that the seabed slope in the nearshore area is very mild, we consider as a first approximation that the wave conditions at the OWSC (depth 20 m) remain the same as in TP5 (depth 137 m). Thus, for each pair of mean energy period and significant wave height of the nearshore climatology (Figure 9), the captured wave power per meter of wavefront is calculated by multiplying the resulting wave power flux per meter by the value of the device’s performance curve that corresponds to the given period.

$$P_{out}=P_w{\eta }_{OWSC}(\omega )=\frac{\rho g_{}^2}{64\pi }{H}_s^2{T}_e^{}{\eta }_{OWSC}(\omega ).$$

(32)

Although the present model is restricted to two dimensions, it is approximately used in this section to investigate the effects of the geometrical inhomogeneities of the seabed boundary on the performance of the considered devices. The results concerning power output estimations are considered preliminary, and will be updated by future extensions of the present model to treat realistic three-dimensional and real fluid effects.

Using the analytical model described in the Appendix A, the PTO parameters are first optimized, and at a next step the shoaling effects on the OWSC performance are calculated by means of the present BEM model. The results are obtained under the assumption of a very narrow spectrum, and are based on summing monochromatic components. The device axis of rotation is positioned at a distance of 1 m above the seabed (d/h = 19/20). The site-specific Performance Index (P.I.) of the OWSC is defined as the percentage of the energy captured over the wave power contained in the whole data set, which is calculated by

$$P.I=\frac{{\overline{P}}_{out}}{{\overline{P}}_w}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{P}_{out}(i)}{\left[\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{P}_w(i)}\right]}^{-1}={\displaystyle {\sum }_{i=1}^N{P}_{out}(i)}{\left[{\displaystyle {\sum }_{i=1}^N{P}_w(i)}\right]}^{-1}$$

(33)

where N is the sample size of the wave climate data at TP5. Figure 11 illustrates the site-specific Performance Indices, as calculated for various combinations of the PTO parameters, as well as the flap material density (defined relatively to the sea water density). Each point of the graphs corresponds to the Performance Index [Equation (33)] of an OWSC positioned at 20 m depth, for the given parameters, as calculated by the analytical model.

In the sequel, in order to present more specific results, we shall restrict the analysis to an OWSC flap of relative density equal to 0.5, and with a nondimensional inertia parameter equal to ${\widehat{J}}_P=2$. Figure 12 illustrates the site-specific performance indices along with resulting performance curves, presented as functions of the non-dimensional frequency. First, it is noted that the value of the nondimensional frequency which corresponds to the most probable period of the incident waves in the area, as obtained from the nearshore/coastal climatology (T ≈ 3.5 s), is equal to:

$$\widehat{\omega }=\frac{2\pi }{T}\sqrt{\frac{h}{g}}=\frac{2\pi }{3.5}\sqrt{\frac{h}{g}}\approx 2.55$$

(34)

Interestingly, we observe in Figure 12 that the optimal (maximum) Performance Index (= 58.14%) does not occur for configurations with eigenperiod coinciding with the above value of the mean period for the incident waves, but for OWSC configurations whose peak performances are achieved for greater oscillation periods of around 6.2 s, which is significantly higher. An explanation for this fact is that, although waves characterized by periods of about T = 3.5 s–4 s occur relatively often due to the climate characteristics of the considered region, these waves are associated with relatively low wave heights, and their energy contents are lower. On the other hand, the occurrence of waves characterized by a reasonably increased value of period of (e.g., Τ > 10 s) is very scarce in the specific geographical area.

A normalized distribution of the available wave energy, with respect to the incident period associated with the population of waves in the nearshore/coastal area examined, is presented in Figure 13, as obtained from the climatological data of Figure 9. It is clearly seen in this figure that the wave period with the largest energy content for the specific region’s wave climate characteristics is around 6.2 s. The latter result provides useful information concerning the selection of the PTO parameters of the device for the maximization of power capture in the specific area. In particular, the value of the corresponding nondimensional frequency is equal to $\widehat{\omega }\approx 1.19$, justifying the selection of a flap of relative density equal to 0.5 and PTO system parameters ${\widehat{J}}_P=2,{\widehat{B}}_P=1.9,{\widehat{C}}_P=1.25$, which results in a performance curve with a peak of 40%, for a frequency close to the optimal operating point. The latter values will be used in the next part of the present study in order to investigate the effects of variable seabed topography by means of the developed BEM scheme.

4.2. Bottom Topography Effects

As stated earlier, the water depth at the point TP5 is equal to 137 m. Furthermore, as shown in Figure 8b, the distance between point TP5 and the shoreline is approximately equal to 1 NM (≈1850 m). Assuming that the OWSCs array is located 150 m from the shoreline, the upwave length of the domain is 1700 m, and the seabed profile is modelled using Equation (20), as follows:

$$h(x)=\frac{1}{2}(20+137)-\frac{1}{2}(20-137)\tanh \left[0.007(x-x_m)\right],x_m=850m.$$

(35)

As previously discussed, relatively mild bottom slopes do not significantly affect the performance of the device. Applying the BEM, in the above domain and for the specific PTO parameters, results in the performance curve shown in Figure 14a,b. The corresponding results obtained by the analytical model, for a constant depth of 20 m, are presented in Figure 14c,d.

The Performance Index, which results from calculating the percentage of the power that is captured for the given performance curve, was found to be equal to 34.6%, and differs by a small amount from the result at a constant depth (which is 34.8%). This implies that the effects of the seabed’s profile in this case are negligible. As depicted in Figure 14e, significant alterations to the performance curve, due to bathymetry inhomogeneities, occur within the bandwidth of the nondimensional frequency between 0 and 0.5, which corresponds to incident wave periods of 18 s or more. This corresponds to long swells, which result in strong seabed–wave–device interaction. Such data, however, are not present in the climatology of the specific coastal area, and thus, the variable bathymetry causes negligible changes in the performance of the OWSC. Nonetheless, in other cases possibly characterized by different wave statistics and stronger seabed topography variations, the sloping bottom effects are expected to be significant, and can be fully accounted for by the present model.

Based on the above calculated performance of the OWSC, the monthly mean power output, per meter of wavefront, in the examined nearshore/coastal area is presented in Figure 15a. Moreover, the monthly mean percentages of power capture are depicted in Figure 15b. It can be observed that the device’s performance is slightly reduced during the spring season, and appears increased during the winter months. This is mainly due to the variability of the monthly mean incident wave period, which results in a shifting of the device’s operation point on the performance curve, and leads to different fractions of the incident wave energy being captured. Subsequently, we assume an array of three OWSCs positioned along the 20-m isobath in the coastal area, and over the 20-m isodepth in the nearshore area, of the TP5 point, as shown in Figure 16. Each OWSC is considered to be 20 m wide. The assumed device’s spacing is defined in order to ignore, as a first approximation, interaction effects. Furthermore, the device’s width could modify the system’s eigenperiod, and such effects are left to be examined in future work extending the present method to 3D applications.

The considered OWSC array is illustrated in Figure 16. Given that the results were extracted via a two-dimensional analogue, which refers to unit length in the third dimension, the power output of any unit of the considered OWSC park can be estimated by multiplying the power output per meter of wavefront with the unit’s width. An average of 5% losses was assumed due to the unit’s shut down for maintenance or repair works after failure (η_op = 0.95), and 2% losses due to transmission efficiency (η_tr = 0.98; see [21]). The monthly mean power output values of the OWSC energy park can then be obtained by multiplying the unit’s output power with the number of OWSCs, and the result concerning the estimated corresponding monthly mean power output in [kW] is depicted in Figure 17.

Eventually, the annual energy production (AEP) of the OWSC array is approximately estimated (in kWh) by multiplying the annual averaged electrical power output of the park by the total number of operating hours within a year, as follows:

$$AEP=76.5 kW=8766 h=680,599.00 kWh\approx 0.6806 GWh$$

(36)

4.3. Estimation of Power Output by Point Absorbers and Comparative Evaluation

For comparison purposes, the installation of an energy park, consisting of an array of point absorbers, is examined here, using the model presented in a previous section. The park is located in the exact same geographic area. The array consists of six point absorber buoys, with a diameter of 10 m, so that the total distance covered by wave energy-exploiting devices remains unchanged relative to the previously examined park (60 m), as illustrated in Figure 18. For the purposes of the comparative analysis, an indicative performance curve of these devices was adapted from [8,9].

The performance curve of the WEC point absorber is based on the assumption of a typical PTO damping coefficient, with an indicative value considered as 10 times the mean value of hydrodynamic damping coefficient $b_m$. The expression used in [8,9] for the nondimensionalization of the incident waves’ frequency involves the heaving cylinder’s radius. Therefore, based on the considered diameter, the performance of the buoys can be expressed as a function of the nondimensional frequency $\widehat{\omega }=\sqrt{h/g}$, as depicted in Figure 19.

It can be observed in Figure 19 that the maximum performance achieved is about 70%, for $\widehat{\omega }\approx 1.38$, which implies a natural period of around 6.5 s, compatible with the peak of the distribution of available wave energy at TP5. The monthly mean power output per meter of wavefront, based on the nearshore data calculated for the wave climate in TP5 (Figure 10) and the WEC performance curve of Figure 19, is presented in Figure 20a. Moreover, the monthly mean percentages of power capture are depicted in Figure 20b. It can be observed that the device’s performance differs significantly between months, which is a result of the increased performance in a very limited frequency range, which mainly occurs during winter, based on the available wave climate data.

The power output of any unit of the considered point WEC array can be estimated by multiplying the power output per meter of wavefront with the unit’s diameter. The losses due to device availability and transmission efficiency are considered to be the same as before (η_op = 0.95, η_tr = 0.98). Therefore, for the WEC array of Figure 18, the estimated corresponding monthly mean power output is the one depicted in Figure 21. The annual energy production (AEP) is 0.581 GW/h. The AEP of the point absorbers array is found to be lower than that of the OWSCs park by a factor of around 15%, despite the fact that the performance curve of the heaving buoys reaches a maximum value of 70%, which is greater than that of the OWSC. This is due to the limited frequency range within which the considered WEC performance (without fancy control) is increased, in comparison to the OWSC response and performance, which is more broadband.

Specifically, unlike OWSCs, which are excited by the wave field in a significant fraction of the fluid channel, the point absorber buoys only exploit the power flux concentrated near the free surface’s level. As a result, the performance indices corresponding to frequencies that imply shallow water wave propagation are relatively lower for this type of energy device.

The above analysis is based on a simplified extension of a two-dimensional analogue, and therefore all the effects induced from the three-dimensional nature of the phenomenon, as well as oblique wave incidence, are neglected. The function of individual point absorbers is independent of the propagation direction of the excitation wave field. However, this is not the case for arrays or grids consisting of several such devices, due to the hydrodynamic interactions between the different units. The latter effects are left to be examined in more detail in future works.

In conclusion, the proper design and dimensioning of the wave energy technologies, as the ones studied here, can lead to significant exploitation of marine renewable energy. In addition, with appropriate optimization of device parameters, taking into account the wave climatology of the considered location of installation, significant output can be achieved. As regards Greece’s nearshore/coastal region, the development of wave energy parks would probably lead to increased energy costs, except for sites with increased wave potential. However, as the relative technology and know-how reach new levels of maturity, the costs associated with designing and manufacturing wave energy devices will most likely reduce, and from this perspective, the development of wave energy parks in Greece is likely to take place in the near future. For this purpose, more detailed investigation of various successful devices and concepts is needed in order to identify the most efficient solution.

5. Conclusions

In this work a simplified BEM is developed and applied to the investigation of variable bathymetry effects on the performance of surge-type wave energy converters (OWSC) excited by harmonic incident waves. The terminator-type model is studied in two dimensions, where the flap of the device is aligned perpendicularly to the incident wave direction. On top of the standard absorbing–radiating model, which includes additional effects due to wave transmission in the downwave region behind the flap-board, the special case of absorbing–reflecting model is also considered. Although the present study is restricted for simplicity to two dimensions, the model offers a useful tool that could be systematically applied, with low computational cost, to investigate the effects of the geometrical inhomogeneities of the seabed’s boundary, and could support future extensions into more realistic three-dimensional configurations.

The present method is first validated with the prediction of results against analytical solutions derived at a constant depth. Subsequently, the relation between the various sea states and the device performance is investigated by systematic application of the BEM model, and the effects of several important parameters that characterize the configuration are examined and discussed. Numerical results, illustrating the modification of the wave field and the effects of depth variation in conjunction with other parameters, like inertia and power-take-off, are illustrated.

It is also shown that the effects of bottom slope and curvature, especially on the response and secondly on the performance of the studied device, could be important, especially for low and moderate wave frequencies, and should be taken into account. Finally, the installation of wave energy devices in specific nearshore/coastal regions of Greece, characterized by increased potential, is examined, and the power output from OWSCs is compared with WEC point absorber installations. Future work will be directed towards the extension of the present model into studying the responses and the performance of three-dimensional systems in more realistic, irregular and general-incidence wave conditions, as well as the consideration of non-linear effects.