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

: 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 signiﬁcant e ﬀ ects on their performance. In particular, Oscillating Wave Surge Converters (OWSCs) exploit the horizontal motion of water waves interacting with the ﬂap of the device. In this work, a Boundary Element Method (BEM) is developed, and applied to the investigation of variable bathymetry e ﬀ ects on the performance of a simpliﬁed 2D model of a surge-type wave energy converter excited by harmonic incident waves. Numerical results, illustrating the e ﬀ ects of depth variation in conjunction with other parameters, like inertia and power-take-o ﬀ , on the performance of the device, are presented. Finally, a comparative evaluation of the present simpliﬁed 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.


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. 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: , with g denoting the acceleration of gravity and H the incident wave height. 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.

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 α < x, −h(x) < 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ωt), and the complex wave potential is represented by: where G( x| x 0 ) = (2π) −1 ln x − x 0 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 0 (x,z) indicates the singularity point, and ∂D = ∂D 1 ∪ ∂D 2 ∪ ∂D 3 ∪ ∂D 4 denotes the boundary of the domain decomposed into four sections, as presented in Figure 1: the free surface ∂D 1 , the device part ∂D 2 , the seabed ∂D 3 and the vertical interface ∂D 4 serving as an incidence-radiation boundary. The solution to the potential flow problem is obtained by discretizing the boundary geometry into elements ∂D = ∪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 (σ 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 ∂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: where Q = −igH/(2ω), with g denoting the acceleration of gravity and H the incident wave height. Moreover, R is the reflection coefficient, the function Z (b) 0 (z) = cosh(k 0 (z + h b ))/ cosh(k 0 h b ) and the wavenumber k 0 of the incident-reflected wave is obtained from the dispersion relation formulated at the radiation boundary depth h b : 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: The vertical eigenfunctions {Z with the corresponding parameters ik n , n = 1, 2, . . . , obtained as the imaginary roots of Equation (3). By projecting the two members of Equation (2) in the above orthogonal vertical basis we obtain:  (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 : Substitution of the reflection coefficient's expression in Equation (4) results in the following radiation condition, that applies on the interface (radiation boundary): where ∂ϕ/∂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: where the functions f (z) 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).

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 where M 3 is the number of elements on the seabed, h(x) 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

Discrete BEM Model
By assuming unit source-sink singularity strengths on each element, the induced potential and velocity matrices ϕ k j , U k j , k, j = 1, . . . , 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 k j of the discrete system is obtained, as explained in more detail in the sequel. By denoting with σ 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: Fluids 2020, 5, 99

of 29
Apart from the discrete source singularities' strengths σ j , j = 1, M , the additional unknown parameters are the complex amplitude θ 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: are the data associated with the incident wave on the incidence-radiation interface ∂D 4 , and z k denotes the collocation points, located at the center of the k-element on ∂D 4 . The auxiliary index α, introduced in Equation (9), spans the interval of integers that correspond to the elements distributed across ∂D 4 . Therefore, if δz 4 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: Moreover, the boundary condition on the boundary of the OWSC (see also Appendix A) is: Using the matrix-coefficients for the induced potential and velocity, the above Equation is expressed in discretized form as follows: 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: Fluids 2020, 5, 99 6 of 29 Concerning the additional Equation closing the discrete system, with respect to the additional unknown θ 0 , the dynamic Equation governing the angular harmonic oscillator that models the present 2dOWSC, expressed in the frequency domain, is considered: 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: Using the representation Equation (A5b) of the wave potential in the downwave side of the device to express ϕ(x = 0+, z), 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 δz 2 is the common length of the boundary elements on ∂D 2 , and the result is where the frequency-dependent coefficient Λ(ω) is defined as follows: 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 σ j , j = 1, . . . , M , as well as the complex oscillation amplitude θ 0 of the present 2dOWSC, the reflection coefficient can be calculated from Equation (6) as: 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 and the consideration of energy-flux conservation leads to the following equation: 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.

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:Ĵ P = 1,B P = 1.5,Ĉ 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 1 and ∂D 3 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. The 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.

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: which describes a smooth shoal with maximum bottom slope controlled by the parameter κ, while δh c (x)

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 Fluids 2020, 5, x FOR PEER REVIEW 9 of 31 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 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  In particular, in Figure 3b the calculated real and the imaginary part of the complex amplitude θ 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 ω h a /g ≈ 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. By 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ω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 δh c (x) = 0, κ = 0.09, x b = 200 m, by means of equipotential lines. For the specific period of 8 s (ω = 0.785 rad/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 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. In 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. uniform bathymetry region, whose seabed profile is defined by Equation (20)  ), 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 Fluids 2020, 5, x FOR PEER REVIEW 11 of 31 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.    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) Fluids 2020, 5, 99 11 of 29 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.

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 1 hh  in the region of incidence, to 3 h h  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: x , which is the summation of the radiating fields due to the six modes of oscillation of the WEC: On the wetted surface of the WEC, the following boundary conditions are applied: ,, n n n  n is the normal vector directing inwards the geometry of the floater, and n is the 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: Figure 6. Equipotential lines of the wave field and free surface elevation at different times as calculated by the present BEM in the case of the absorbing-reflecting model in a variable bathymetry region (Parameters :

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: where x = (x, y) denotes the horizontal coordinates and µ = ω 2 /g is the frequency parameter. The wave potential is decomposed on the propagating wave potential ϕ P (x, z), with the absence of the floating body, the diffraction potential ϕ D (x, z) due to the presence of the impermeable and motionless geometry of the WEC and the radiation potential ϕ D (x, z), which is the summation of the radiating fields due to the six modes of oscillation of the WEC: On the wetted surface of the WEC, the following boundary conditions are applied: ∂ϕ D (x, z)/∂n = −∂ϕ P (x, z)/∂n, ∂ϕ (x, z)/∂n = n , = 1, 2, . . . 6, where n = (n 1 , n 2 , n 3 ) is the normal vector directing inwards the geometry of the floater, and n is the 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: In the above expression of heave response amplitude operator-ξ 3 (RAO ξ 3 ), the implemented Froude-Krylov X p and diffraction X D exciting forces are evaluated by integration of the wetted surface of the body ∂D B , as follows: while A(ω) is calculated by: The hydrodynamic coefficients of added mass (α 33 ) and damping (b 33 ) in Equation (26) are obtained by integrating the heave-radiation potential on the submerged boundary surface of the WEC: 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: where η 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: The vertical functions Z n (z; x) are obtained as eigenfunctions of regular Sturm-Liouville problems, formulated at the local depth, while the functions ϕ n (x) are the complex amplitude of the n th -mode, and are found as the solution of the following Coupled Mode System (CMS) 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: 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 = (x, y, z). 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].

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.
where summation refers to all the quadrilateral elements, and Φp and Up are used to denote, respectively, the potential and velocity induced from the p th element of unit singularity distribution to a random field point   ,, x y z  r . 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].

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 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] 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. 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 10(a), 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). 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).

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.

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.

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. 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 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. 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 ˆ1. 19   , justifying the selection of a flap of relative density equal to 0.5 and PTO system parameters 2 , 1. 9, 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.

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 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Ĵ 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:  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  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., T > 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ω ≈ 1.19, justifying the selection of a flap of relative density equal to 0.5 and PTO system parametersĴ P = 2,B P = 1.9,Ĉ 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.

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: ). Figure 13. Distribution of available wave energy at TP5 as a function of the incident wave period.
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  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  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 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 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:

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 m b . 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 ˆhg   , as depicted in Figure 19. 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 ≈ 0.6806 GWh (36)

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 frequencŷ ω = h/g, as depicted in Figure 19.     It can be observed in Figure 19 that the maximum performance achieved is about 70%, forω ≈ 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.
output is the one depicted in Figure 21. The annual energy production (AEP) is 0.581 GW/ℎ. 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 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.

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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Linearized 2dOWSC Model in Constant Water Depth
Assuming, for simplicity, the two-dimensional problem concerning normally incident waves at a constant depth, the linearized model of the considered terminator flap-type model is shown in Figure A1. In particular, except of the standard absorbing-radiating model depicted in Figure A1a, which includes additional effects due to wave transmission in the downwave region behind the flap-board, an absorbing-reflecting model, shown in Figure A1b, is also considered. Although the latter configuration could be represented by appropriately adjusting the PTO parameters to achieve zero transmitted wave, it is still considered here as a special configuration of the 2dOWSC problem, and the solution is used for the validation of the developed BEM. a constant depth, the linearized model of the considered terminator flap-type model is shown in Figure A1. In particular, except of the standard absorbing-radiating model depicted in Figure A1a, which includes additional effects due to wave transmission in the downwave region behind the flapboard, an absorbing-reflecting model, shown in Figure A1b, is also considered. Although the latter configuration could be represented by appropriately adjusting the PTO parameters to achieve zero transmitted wave, it is still considered here as a special configuration of the 2dOWSC problem, and the solution is used for the validation of the developed BEM. Given that the flap motion is generated by incident waves, a fraction of the incoming wave energy (I) received by the device is expected to be reflected (R), while in the case of the more general absorbing-radiating model, another fraction is transmitted (T). Considering that the flap performs harmonic angular oscillations, of maximum width 0  with frequency  , the angle of oscillation at any time is given by: Figure A1. Simplified 2D model of terminator-type OWSC in constant water depth: (a) the standard absorbing-radiating model, and (b) the special case of an absorbing-reflecting model.
Given that the flap motion is generated by incident waves, a fraction of the incoming wave energy (I) received by the device is expected to be reflected (R), while in the case of the more general absorbing-radiating model, another fraction is transmitted (T). Considering that the flap performs harmonic angular oscillations, of maximum width θ 0 with frequency ω, the angle of oscillation at any time is given by: In order to ensure that the system will perform harmonic oscillations, a torsional spring at the hinge point is considered, which restores the flap after its deviation from the vertical position during each period, due to the incident waves. Moreover, the hydraulic generator, which compresses fluid and supplies the hydroelectric power station, uses part of the oscillation's energy, and therefore acts as a damper in the contemplated torsional oscillator. For simplicity, power extracted by the system (PTO) is considered to be converted into electricity at the hydroelectric power station, without any further losses, e.g., due to transmission through the pipeline system, or mechanical losses. The equation for the motion of the angular oscillator is: where J P is the moment of inertia of the flap with respect to the axis of rotation, B P is the damping coefficient associated with the PTO, C P is the spring constant of the installed torsional spring and M(t) is the time dependent torque that the flap receives from the waves. Assuming small angles of deviation from the vertical position (θ 0 ), the horizontal velocity of any point on the flap will be equal to: and L = z + d indicates the lever arm of a given point on the flap, located at depth z, with respect to the center of rotation, and d is the hinge point's depth; see Figure A1. Given the assumption of small angles of deviation from the vertical position (tan(θ 0 ) ≈ θ 0 ), the horizontal excursion of the flap motion at the free surface level is: The general expression of the complex potential ϕ(x, z) describing the wave field in the semi-infinite constant depth strip is as follows: where Q = −igH/2ω and H is the incident wave height. The first term in Equation (A5) corresponds to the left-propagating mode. R represents the complex amplitude of the reflected component and A n , n ≥ 1 represent the corresponding amplitudes of the evanescent modes. In the case of the more general absorbing-radiating model of Figure A1a, the expression of the wave potential in the downwave region behind the flap-board is given by: where T represents the complex amplitude of the transmitted wave and B n , n ≥ 1 are the corresponding evanescent modes. Moreover,Z n (z) = cosh(k n (z + h))/ cosh(k n h) represents the vertical structure of each mode, while k 0 and ik n are obtained as the real and imaginary roots of the dispersion relation: The free surface's elevation is given by: By using the above representations, the free surface elevation in the upwave direction can be expressed as: and in the case of the more general absorbing-radiating model of Figure A1b, the free surface elevation in the downwave direction is: The boundary condition on the surface of the flap specifies that the horizontal body velocity is equal to the fluid's horizontal velocity for z ∈ (−d, 0), and in terms of the complex potential ϕ, it is expressed as follows: The above boundary condition must be satisfied by the wave potential, along with the free surface and bottom boundary conditions, which respectively are: ∂ϕ ∂z − ω 2 g ϕ = 0, z = 0, and ∂ϕ ∂z = 0, z = −h (A10) The x-derivative of the wave potential provides the horizontal component of the fluid's particles' velocity in the flow field. The latter evaluated at the position of the flap x = 0 reads as follows: In the case of the more general absorbing-radiating model of Figure A1a, the x-derivative of the wave potential in the downwave side is obtained from Equation (A5b) as By expressing the function f (z) of the horizontal velocity of the flap in terms of the vertical eigenfunctions Z n (z), and using it in the boundary condition Equations (A9) and (A11), we obtain the following relations for the coefficients R, A n and T, B n , respectively: where M 0 is the moment applied to the flap as a result of the fluid's pressure from the upwave side, and M B is the moment due to pressure forces from the downwave side, which is only included in the open-type model. The latter quantities are expressed in terms of the complex wave potential by integrating the pressure forces evaluated from a linearized Bernoulli's equation (p = iωρϕ), multiplied by the lever arm across the flap's length: which respectively become Using Equation (A12) and substituting in the above the coefficients R and A n , and T and B n , respectively, we finally get: By using the above in Equation (A14) and solving for θ 0 , we obtain for the angular response of the system: O n 2 k n Z n 2 cos 2 (k n h) and the coefficient B = 1 in the absorbing-reflecting model, while B = 1 for the absorbing-radiating model. Finally, by substituting the flap angular response (θ 0 ) in the Equations (A12), the coefficients R and A n , and T and B n , are calculated. The complex wave potential and the free surface elevation are obtained from Equations (A5) and (A8), respectively. The performance of the device is estimated as the fraction of the power input that is converted into beneficial output power. The average wave power (P in ) associated with the incident wave is given by: where c g is the group velocity of water waves in the constant depth strip: The average output power converted by the device is calculated through the damping coefficient (modelling the hydraulic generator) as: Subsequently, the Performance Index (η OWSC ) of the device can be estimated by the quotient of the two above. In addition, for a given incident wave height H, the reflected wave height will be equal to |R|H, and in the case of the open-type model (B = 1), the transmission coefficient is equal to |T|H, and thus, is a relation that can be used to verify conservation of energy in the examined hydromechanical system. In order to illustrate the numerical performance of the above simple model, the following dimensionless parameters are introduced concerning the coefficients of the OWSC Power Take Off model: where ρ m represents the density of the flap's material. It is noted that, since a two-dimensional analogue of the flap type surge converter is being considered, all the terms of the differential equation of motion of the angular oscillator refer to their respective units per unit length in the transverse direction. The extension of the model for the study of three-dimensional phenomena will be examined in future work.
As an example, we consider an OWSC, with non-dimensional PTO parameters equal tô J P = 1,B P = 1.5,Ĉ P = 1 and density ratio ρ m /ρ = 1, operating in water depth equal to h = 10 m, with the device's hinge point located at the depth of 8 m (d/h = 0.8). Numerical results are presented in Figure A2 for incident waveheight H = 1 m, illustrating the calculated performance of the examined system obtained by keeping 10 terms in the modal series Equation (A5), which is found to be enough for the convergence of the results in the whole band of frequencies considered. Both absorbing-radiating and absorbing-reflecting systems 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 A2, respectively. More specifically, the first-row subplots of Figure A2 illustrates the real and the imaginary part of the calculated complex amplitude θ0, as well as its absolute value, as a function of the nondimensional frequency. Moreover, the second-row subplots depict the performance of the studied device (plotted by using thick solid lines) and the reflection coefficient R (plotted by using dashed lines), while in the case of the open-type system the transmission coefficient is also included More specifically, the first-row subplots of Figure A2 illustrates the real and the imaginary part of the calculated complex amplitude θ 0 , as well as its absolute value, as a function of the nondimensional frequency. Moreover, the second-row subplots depict the performance of the studied device (plotted by using thick solid lines) and the reflection coefficient R (plotted by using dashed lines), while in the case of the open-type system the transmission coefficient is also included (by using dotted lines). It is worth noticing that the conservation of energy expressed by Equation (A23) is perfectly satisfied for both examined types and for all frequencies. We also observe in this figure that in the case of the more general absorbing-radiating system, the calculated performance is a few percent lower than in the absorbing-reflecting system, due to additional losses associated with radiating wave energy in the downwave side of the flap-board. On the other hand, the change concerning the angular response of the flap is much smaller.