A Novel Method for Estimating Wave Energy Converter Performance in Variable Bathymetry Regions and Applications

Abstract

A numerical model is presented for the estimation of Wave Energy Converter (WEC) performance in variable bathymetry regions, taking into account the interaction of the floating units with the bottom topography. The proposed method is based on a coupled-mode model for the propagation of the water waves over the general bottom topography, in combination with a Boundary Element Method for the treatment of the diffraction/radiation problems and the evaluation of the flow details on the local scale of the energy absorbers. An important feature of the present method is that it is free of mild bottom slope assumptions and restrictions and it is able to resolve the 3D wave field all over the water column, in variable bathymetry regions including the interactions of floating bodies of general shape. Numerical results are presented concerning the wave field and the power output of a single device in inhomogeneous environment, focusing on the effect of the shape of the floater. Extensions of the method to treat the WEC arrays in variable bathymetry regions are also presented and discussed.

1. Introduction

Interaction of the free-surface gravity waves with floating bodies, in water of intermediate depth and in variable bathymetry regions, is an interesting problem with important applications. Specific examples concern the design and evaluation of the performances of special-type ships and structures operating in nearshore and coastal waters; see, e.g., [1,2]. Also, pontoon-type floating bodies of relatively small dimensions find applications as coastal protection devices (floating breakwaters) and they are also frequently used as small boat marinas; see, e.g., [3,4,5,6,7]. In all these cases, the estimation of the wave-induced loads and motions of the floating structures can be based on the solution of the classical wave-body-seabed hydrodynamic-interaction problems; see, e.g., [8,9]. In particular, the performance of the Wave Energy Converters (WECs) operating in nearshore and coastal areas, characterized by variable bottom topography, is important for the estimation of the wave power absorption, determination of the operational characteristics of the system and could significantly contribute to the efficient design and layout of the WEC farms. In this case, wave-seabed interactions may have a significant effect; see [10,11].

In the above studies the details of the wave field propagating and scattered over the variable bathymetry region could be important in order to consistently calculate the responses of the floating bodies. For rapidly varying seabed topographies, including steep bottom parts, local or evanescent modes may have a significant impact on the wave phase evolution during propagation. Such a fact was demonstrated through the interference process in one-directional wave propagation as observed for either varying topographies (see e.g., [12,13]) or abrupt bathymetries including coastal structures (see e.g., [14,15,16]). For such problems, the consistent coupled-mode theory has been developed in [17], for the water waves propagation in variable bathymetry regions. Furthermore, it was subsequently extended for 3D bathymetry in [18], and applied successfully to treat the wave transformation over nearshore/coastal sites with steep 3D bottom features, like underwater canyons; see, e.g., [19,20].

In recent works [21,22] the coupled-mode model is further extended to treat the wave-current-seabed interaction problem, with application to the wave scattering by non-homogeneous current over general bottom topography. The problem of the directional spectrum transformation of an incident wave system over a region of strongly varying three-dimensional bottom topography is further studied in [22]. The accuracy and efficiency of the coupled-mode method is tested, comparing numerical predictions against experimental data by [23] and calculations by the phase-averaged model SWAN [24,25]. Results are shown in various representative test cases demonstrating the importance of the first evanescent modes and the additional sloping-bottom mode when the bottom slope is not negligible.

In this work, a methodology is presented to treat the propagation-diffraction-radiation problem locally around each WEC, supporting the calculation of the interaction effects of the floating units with variable bottom topography at a local scale. The method is based on the coupled-mode model developed by [17], and extended to 3D by [18], for water wave propagation over general bottom topography, in conjunction with the Boundary Element Method (BEM) for the hydrodynamic analysis of floating bodies over general bottom topography [15] and the corresponding 3D Green’s function [26]. An important feature of the present method is that it is free of mild-slope assumptions and restrictions and it is able to resolve the 3D wave field all over the water column, in variable bathymetry regions including the interactions of floating bodies of general shape. Numerical results are presented and discussed concerning simple bodies (heaving vertical cylinders) illustrating the applicability of the present method.

2. Formulation

We consider here the hydrodynamic problem concerning the behavior of a number N of identical cylindrical-shaped WECs, $D_B^{\left(k\right)},k=1,N$, of characteristic radius a and draft d, operating in the nearshore environment, as shown in Figure 1.

The variable bathymetry region is considered between two infinite sub-regions of constant, but possibly different depths $h=h_1$ (region of incidence) and $h=h_3$ (region of transmission). In the middle sub-region, it is assumed that the depth $h$ exhibits an arbitrary variation. The wave field is excited by a harmonic incident wave of angular frequency ω, propagating with direction θ; see Figure 1. Under the assumptions that the free-surface elevation and the wave velocities are small, the wave potential is expressed as follows:

$$\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\},$$

(1)

where $x=\left(x_1,x_2\right),$ and satisfies the linearized water wave equations; see [27]. In the above equation $H$ is the incident wave height, g is the acceleration due to gravity, $\mu ={\omega }^2/g$ is the frequency parameter, and $i=\sqrt{-1}$. The free surface elevation is then obtained in terms of the wave potential as follows:

$$\eta \left(x;t\right)=-\frac{1}{g}\frac{\partial \mathsf{\Phi }\left(x,z=0\right)}{\partial t}.$$

(2)

Using standard floating-body hydrodynamic theory [8], the complex potential can be decomposed as follows:

$$\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 _{k=1}^N{\displaystyle \sum _{\ell =1}^6{\xi }_{k\ell }}}{\phi }_{k\ell }\left(x,z\right)$$

(3)

where ${\phi }_P\left(x,z\right)$ is the normalized propagation wave potential in the variable bathymetry region in the absence of the WECs, ${\phi }_D\left(x,z\right)$ is the diffraction potential due to the presence fixed (motionless) bodies $D_B^{\left(k\right)},k=1,N$, that satisfies the boundary condition $\partial {\phi }_D\left(x,z\right)/\partial n_k=-\partial {\phi }_P\left(x,z\right)/\partial n_k$ on k-WEC, where $n_k={\left(n_1,n_2,n_3\right)}_k$ the normal vector on the wetted surface of the k-body, directed outwards the fluid domain (inwards the body). Furthermore, ${\phi }_{k\ell }\left(x,z\right),k=1,N$, denotes the radiation potential in the non-uniform domain associated with the $\ell$-oscillatory motion of the k-body with complex amplitude ${\xi }_{k\ell }$, satisfying $\partial {\phi }_{k\ell }\left(x,z\right)/\partial n=n_{k\ell }$, equal to the $\ell$-component of generalized normal vector on the wetted surface of the k-WEC ($n_{k\ell }={\left(r\times n_k\right)}_{\ell -3}$ for $\ell =4,5,6$).

In the case of simple heaving WECs, only the vertical oscillation of each body is considered ${\xi }_k={\xi }_{k3}$, which is one of the most powerful intensive modes concerning this type of wave energy systems. In the present work we will concentrate on this simpler configuration, leaving the analysis of the more complex case to be examined in future work. For an array of heaving WECs the hydrodynamic response is obtained by:

$${\xi }_{k3}={\left(A_{km}^{}\right)}^{-1}\left(X_{Pm}+X_{Dm}\right),k,m=1,\dots N,$$

(4)

where $X_{Pm}+X_{Dm}$ denote the exciting vertical force on each WEC due to propagating and diffraction field, respectively, and the matrix coefficient $A_{km}^{}$ is given by:

$$A_{km}^{}=-{\omega }^2\left(M+a_{km}\right)-i\omega \left(B_m{\delta }_{km}^{}+b_{km}\right)+\left(C_m+c\right){\delta }_{km}^{},$$

(5)

where ${\delta }_{km}^{}$ denotes Kronecker’s delta and M is the body mass (assumed the same for all WECs). The hydrodynamic coefficients (added mass and damping) are calculated by the following integrals:

$$a_{km}-\frac{1}{i\omega }{b}_{km}=\rho {\displaystyle \underset{\partial D_{Bm}}{\iint }{\phi }_{k3}{n}_{m3}}dS,k,m=1,\dots N,$$

(6)

of the heaving-radiation potential of the k-WEC on the wetted surface $\partial D_{Bm}$ of the m-WEC. Moreover, $c=\rho g{A}_{WL}$ is the hydrostatic coefficient in heaving motion with $A_{WL}$ the waterline surface, and $B_m,C_m$ are characteristic constants of the Power Take Off (PTO) system associated with the m-th degree of freedom of the floater. The components of the excitation (Froude-Krylov and diffraction) forces are calculated by the following integrals of the corresponding potentials:

$$X_{Pm}=\frac{\rho gH}{2}{\displaystyle \underset{\partial D_{Bm}}{\iint }{\phi }_P{n}_m}dS,X_{Dm}=\frac{\rho gH}{2}{\displaystyle \underset{\partial D_{Bm}}{\iint }{\phi }_D{n}_m}dS,m=1,\dots N$$

(7)

on the wetted surface $\partial D_{Bm}$ of the m-WEC. The total power extracted by the array is obtained as:

$$P\left(N;\omega ,\theta \right)=\frac{1}{2}{\omega }^2\left|{\displaystyle \sum _{k=1}^N{\eta }_{eff}^m{B}_m{\left({\xi }_k\right)}^2}\right|,$$

(8)

where ${\eta }_k$ indicates the efficiency of the PTO associated with the k-th degree of freedom (that could be a function of the frequency ω). Finally, the q-index can be estimated by:

$$q_N\left(\omega ,\theta \right)=N^{-1}P\left(N;\omega ,\theta ,H\right)/P_1\left(\omega ,\theta ,H\right),$$

(9)

where $P_1\left(\omega ,\theta ,H\right)$ indicates the output of a single device operating in the same environment and wave conditions. Obviously, the calculated performance depends on the frequency, direction and height of the incident wave, as well as on the physical environment and the positioning of the WECs in the array (farm layout). Finally, the operational characteristics of the farm, in general multi chromatic wave conditions, characterized by directional wave spectrum, could be obtained by appropriate spectral synthesis; see, e.g., [20,22].

3. Propagating Wave Field

The wave potential ${\phi }_P\left(x,z\right)$ associated with the propagation of water waves in the variable bathymetry region, without the presence of the scatterer (floating body), can be conveniently calculated by means of the consistent coupled-mode model developed [17], as extended to three-dimensional environments by [18]. This model is based on the following enhanced local-mode representation:

$${\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).$$

(10)

In the above expansion, the term ${\phi }_0\left(x\right)Z_0\left(z;x\right)$ denotes the propagating mode of the generalized incident field. The remaining terms ${\phi }_n\left(x\right)Z_n\left(z;x\right),$ $n=1,2,\dots ,$ are the corresponding evanescent modes, and the additional term ${\phi }_{-1}\left(x\right)Z_{-1}\left(z;x\right)$ is a correction term, called the sloping-bottom mode, which properly accounts for the satisfaction of the Neumann bottom boundary condition of the non-horizontal parts of the bottom. The function $Z_n\left(z;x\right)$ represents the vertical structure of the n-th mode. The function ${\phi }_n\left(x\right)$ describes the horizontal pattern of the n-th mode and is called the complex amplitude of the n-th mode. The functions $Z_n\left(z;x\right)$, $n=0,1,2\dots$, are obtained as the eigenfunctions of local vertical Sturm-Liouville problems formulated in the local vertical intervals $-h\left(x\right)\le z\le 0$, and are given by:

$$Z_0^{}\left(z;x\right)=\frac{\cosh \left[k_0^{}\left(x\right)\left(z+h\left(x\right)\right)\right]}{\cosh \left(k_0^{}\left(x\right)h\left(x\right)\right)},Z_n^{}\left(z;x\right)=\frac{\cos \left[k_n^{}\left(x\right)\left(z+h\left(x\right)\right)\right]}{\cos \left(k_n^{}\left(x\right)h\left(x\right)\right)},n=1,2,\dots$$

(11)

In the above equations the eigenvalues $\left\{i{k}_0^{}\left(x\right),k_n^{}\left(x\right)\right\}$ are obtained as the roots of the local dispersion relations:

$$\mu h\left(x\right)=k_0\left(x\right)h\left(x\right)\tanh \left[k_0\left(x\right)h\left(x\right)\right],\mu h\left(x\right)=-k_n\left(x\right)h\left(x\right)\tan \left[k_n\left(x\right)h\left(x\right)\right]n=1,2,\dots$$

(12)

The function $Z_{-1}\left(z;x\right)$ is defined as the vertical structure of the sloping-bottom mode. This term is introduced in the series in order to consistently satisfy the Neumann boundary condition on the non-horizontal parts of the seabed. It becomes significant in the case of seabottom topographies with non-mildly sloped parts and has the effect of significant acceleration of the convergence of the local mode series Equation (10); see [17]. In fact, truncation of the series (10) keeping only a small number 4–6 totally terms have been proved enough for calculating the propagating wave field in variable bathymetry regions with bottom slopes up to and exceeding 100%. For specific convenient forms of $Z_{-1}\left(z;x\right)$ see the discussion ([17]). By following the procedure described in the latter work, the coupled-mode system of horizontal equations for the amplitudes of the incident wave field propagating over the variable bathymetry region is finally obtained:

$${\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$$

(13)

where the coefficients A_mn, B_mn, C_mn of the coupled-mode system (13) are defined in terms of the vertical modes $Z_n\left(z;x\right)$. The coefficients are dependent on $x$ through $h\left(x\right)$ and the corresponding expressions can be found in

Table 1. Optimum PML parameters.

ω < 2	${\tilde{\mathsf{\sigma }}}_{\mathsf{o}}$ = 1	$R/\lambda$ = 2	Ν/λ = 15	n = 5
2 ≤ ω < 7	${\tilde{\mathsf{\sigma }}}_{\mathsf{o}}$ = 1	$R/\lambda$ = 3	Ν/λ = 20	n = 3
7 ≤ ω < 8	${\tilde{\mathsf{\sigma }}}_{\mathsf{o}}$ = 1	$R/\lambda$ = 3	Ν/λ = 15	n = 3
8 ≤ ω ≤ 9	${\tilde{\mathsf{\sigma }}}_{\mathsf{o}}$ = 1	$R/\lambda$ = 3	Ν/λ = 10	n = 3

of [17]. The system is supplemented by appropriate boundary conditions specifying the incident waves and treating reflection, transmission and radiation of waves. It is worth mentioning here that if only the propagating mode (n = 0) is retained in the expansion (11) the above CMS reduces to an one-equation model which is exactly the modified mild-slope Equation derived in [13,28]. So, the present approach could be automatically reduced to mild-slope model in the subregions where such a simplification is permitted, saving a lot of computational cost. On the other hand in subregions where bottom variations are strong the extra (evanescent) modes are turned on and have substantial effects concerning the 3D wave field all over the water column, as illustrated in [17,18].

4. A Novel BEM for the Diffraction and Radiation Problems in 3D Environments

The corresponding problems on the diffraction and radiation potentials ${\phi }_D\left(x,z\right)\mathrm{and}{\phi }_k\left(x,z\right),$ associated with the operation of the floating WECs, are treated by means of low-order Boundary Element Method, based on simple singularity distributions and 4-node quadrilateral boundary elements ([29]), ensuring continuity of the geometry approximation of the various parts of the boundary. The potential and velocity fields are approximated by:

$$\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),$$

(14)

where the summation ranges over all panels and $F_p\left(r\right)$ and $U_p\left(r\right)$ denote induced potential and velocity from the p-th element with unit singularity distribution to the field point r; see, e.g., [30] and the references cited there. We mention here that a minimum number of 10–20 elements per wavelength is used in discretizing the free surface, in order to eliminate errors due to damping and dispersion associated with the above discrete scheme. In order to eliminate the infinite extent of the domain and treat the radiating behaviour of the diffraction and radiation fields at far distances from the bodies, an absorbing layer technique is used, based on a matched layer all around the fore and side borders of the computational domain on the free surface; see, e.g., [31]. The thickness of the absorbing layer is of the order of 1–2 characteristic wavelengths and its coefficient is taken increasing within the layer; see Figure 2. The efficiency of this technique to damp the outgoing waves with minimal reflection is dependent on the thickness of the layer.

Thus, the diffraction and radiation potentials are represented by integral formulations with support only on the wetted surface of the floating body (ies) $\partial D_C$, the bottom surface $\partial D_{\Pi }$ and free surface $\partial D_F$; see Figure 2. In accordance with the present absorbing layer model, the free surface boundary condition is modified as follows:

$$\frac{\partial \phi }{\partial n}-\mu \sigma \phi =0,r\in \partial D_F,$$

(15)

where $\mu ={\omega }^2/g$ and the coefficient $\sigma =1$ everywhere on $\partial D_F$, except in the absorbing layer (indicated in Figure 2), where this is given by:

$$\sigma =\left(1+i{\sigma }_0\frac{{\left(R-R_a\right)}^n}{{\lambda }^n}\right),R=\sqrt{x_1^2+x_2^2}>R_a$$

(16)

where λ is the local wavelength.

Also, it is assumed for the starting radius of the absorbing layer that $R_a>>\lambda$. The discrete solution is then obtained using collocation method, by satisfying the boundary conditions at the centroid of each panel on the various parts of the boundaries. Induced potentials and velocities from each panel to any collocation point are calculated by numerical quadratures, treating the self-induced quantities semi-analytically.

4.1. Investigation of the Optimal Parameters of the Absorbing Layer

The radiation condition expresses the weakening behavior of the outgoing waves at the far field, and it formulates the final solution. In complicated problems, where analytical or even semi-analytical solutions are unreachable, this condition cannot be formed a priori and further investigation is needed, in order to obtain its final form. One common way to overcome this obstacle, is the implementation of an absorbing layer from a specific length of activation and with defined characteristics, based on the Perfectly Matched Layer (PML) model [32,33]. In this specific approach, the wave absorbing is induced by an imaginary part of frequency, expressed by Equations (15) and (16), which operates as a damping filter for the waves without significant reflections. The formulation of the optimal PML is a multiparametric problem, mainly based on five parameters. The objective functions of this optimization procedure is the avoidance of any influence of the PML in the region before its appliance, due to reflections, and the progressive nullification of the wave field in the region after its activation. The effectiveness of the PML can be tested by comparing the numerical and the analytical solution in case of a cylindrical WEC body in steady depth regions [34]. The requirement for the PML not to disturb the solution in the computational domain before its activation is quantified with the usage of the Chebyshev Norm. Thus, the PML tunning parameters, discussed above, are:

• Dimensionless frequency $\left(\tilde{\mathsf{\omega }}=\mathsf{\omega }\sqrt{\frac{\mathrm{h}}{\mathrm{g}}}\right)$
• Coefficient $\left({\tilde{\mathsf{\sigma }}}_{\mathsf{o}}={\mathsf{\sigma }}_{\mathsf{o}}\frac{{\mathsf{\lambda }}^{\mathrm{n}}}{\mathsf{\omega }}\right)$$\widehat{\sigma }=\sigma {\lambda }^n/\omega$
• The activation length $R/\lambda$
• The exponent $n$
• The number of panels per wave length (N/λ)

Aiming to the minimization of the Chebyshev Norm, 64 different PMLs, corresponding to different sets of these parameters, are investigated. The final configuration of the optimum PML is described in Table 1. The efficient operation of the PML, especially in medium frequency bandwidths, where WEC devices operate most of the time and absorb the largest amount of energy, is illustrated in Figure 3, for different values of the non-dimensional frequency $\omega \sqrt{a/g}\tilde{\mathsf{\omega }}=\mathsf{\omega }\sqrt{\frac{\mathrm{a}}{\mathrm{g}}}$.

4.2. Power Output in the Case of Cylindrical WEC

A cylindrical heaving WEC is widely used in offshore installations of the devices for harnessing wave energy [35]. The numerical treatment of the wave-body interaction problem by means of BEM, described in this study, constitutes from three separate regions, namely the free surface, the body of the WEC and the bottom. Appropriate mesh generation in all these surfaces is crucial for obtaining reliable solution. For this purpose, the free surface is discretized in 4 × (N/λ) × 88 elements, expanding for 4 wavelengths, where the first number indicates discretization along the radial direction and the other along azimuthal direction, respectively, while the bottom mesh is 26 × 88 elements, spatially and azimuthally respectively. The WEC mesh is 10 × 88 elements, in depth and in azimuthal direction, as illustrated in Figure 4. It should be noticed the demand for consistency between the lengths of the elements, those of the WEC and these of the free surface, at the matching position of the body’s boundary. Very fine meshes only on the body and not on the free surface, which binds most of the computational capacity and therefore has its limitations, may cause worse approximation of the analytic solution on account of inconsistency.

Focusing on the power output coming from a single device, the first step for its calculation is the evaluation of the Froude-Krylov and total forces, which are the summation of Froude-Krylov and Diffraction forces, and the related hydrodynamic coefficients of added mass and damping. In Figure 5 and Figure 6 are illustrated the results for these aspects, as they calculated both from the analytical and the numerical treatment of the problem ([34,36]).

Furthermore, the WEC responses and the power output are evaluated and plotted in Figure 7, assuming typical PTO damping values, equal to 5, 10 and 20 times a mean value of hydrodynamic damping b_m. This value is estimated as $2\pi b_m/m{\omega }_R=0.12$, where the resonance frequency is ${\omega }_R\sqrt{a/g}=0.7$, also described in [22].

The output power of the WEC by this PTO is normalized with respect to the incident wave powerflux and is defined as $P/(0.5\rho C_g{H}^{2}a)$, for ${\eta }_{eff}=1$. It can be observed the fact that maximization of power output occurs at the resonance frequency. In addition, higher values of PTO damping are reasons for the observed decrease of peak values of heave Response Amplitude Operator (defined as $RAO={\xi }_{k3}/(H/2)$, where a is the amplitude of the incident wave). However, at the same time they are causing wider frequency spreads of energy productive function of the device.

5. Examination of Other Shapes of Axisymmetric Floaters

Regarding the examination of other WEC shapes, eight different axisymmetric geometries, including the cylinder, are tested. This is made with the conviction of efficiency improvement, in comparison with the reference cylindrical shape. Upon mesh generation, the elements used on the bodies, except cylinder, are 18 × 88, in order to achieve a better approximation of the shape, avoiding gaps and discontinuities of the geometry. For these shapes, there are no analytic solutions, and furthermore, not any prospect for validation by comparing this numerical model with analytical results. The reference cylindrical WEC has a ratio of radius to local depth, equal to 1/3.5 and a ratio of draft to radius equal to 1.5/1. In every other design test, the radius and the draft of each geometry are calculated with the assumption of constant mass. In other words, the area of the submerged vertical cross section of the tested geometry is equal to the area of the submerged vertical cross section of the cylindrical WEC, keeping with this approach the value of the mass unchangeable.

As referred previously, eight different shapes are put under investigation. Heave response and power output are evaluated by the BEM computational code. These geometries, illustrated in Figure 8, are strongly related with the current design trends of the industry and present similarities with many already installed WEC devices [35,37,38,39].

Using as an efficiency index, the area under the curve of the normalized power, which expresses the maximum values so as the functional frequency bandwidth, three of the above geometries are qualified and their heave response and power output are presented in Figure 9, Figure 10 and Figure 11. The qualified geometries are namely the nailhead-shaped, which also presents further interest due to its unconventional design for studies of multi-dof WECs, the conical and the floater-shaped.

On the assessment of these geometries, according to the figures above, despite the fact that the PTO with higher damping is responsible for lower values of heave RAO, the power output appears to be higher and with a higher frequency spreading. This dissimilar behavior of the device, in terms of heave RAO and power output, is very intense in case of the nailhead-shaped WEC, where RAOs are closely oriented, while the power output is far higher in the case of the “harder” PTO. A point of interest in the study of the conical WEC is that a switch of efficiency occurs in $\omega \sqrt{a/g}=1.25$ $\tilde{\mathsf{\omega }}=1.25$, when the medium ranked PTO is more efficient than higher ranked. Furthermore, the floater-shaped WEC is shown to be a little more efficient in lower frequencies, where the maximum of the power output curve is located, while the nailhead-shaped and the conical are more efficient in bandwidths of 11 < $\omega \sqrt{a/g}$ < 1.4. The conical WEC is far more efficient in frequencies of 0.9 < $\omega \sqrt{a/g}$ < 1.4, however, the final choice of the device and the PTO is depended on the sea climate and the dominant frequencies in the area of installation.

6. Extensions to Treat the WEC Arrays in Variable Bathymetry Regions

An important part of the present BEM implementation deals with the construction of the mesh on the various parts of the boundary. The details of the mesh generator are illustrated in Figure 12. More specifically, the mesh on the free surface around a single WEC is plotted. The latter consists of two subparts, the one close to the waterline of the floating unit and the far (outer) part. The near mesh is based on the cylindrical distribution of the panels around the waterline of the WEC that gradually deforms in order to end in a rectangular boundary. This permits the continuous junction of the near mesh around one floater with the adjacent one, as illustrated in Figure 12a. After the rectangular boundary, the mesh again deforms to become a cylindrical arrangement on the outer part. Taking into account that in 3D diffraction and radiation fields associated with floating bodies the far field behaves like essentially cylindrical outgoing waves [26], the cylindrical mesh in the outer part of the free surface boundary is considered to be optimum for the numerical solution of the studied problems. The discretization is accomplished by the incorporation corresponding meshes on each floating body and on the bottom variable bathymetry surface, as shown in Figure 12b–d. An important feature is the continuous junction of the various parts of the mesh, which, in conjunction with the quadrilateral elements, ensures global continuity of the geometry approximation of the boundary. It is remarked here that the present BEM is free of any kind of interior meshes or artificial intesection(s) of boundaries. Global continuity of geometry is important concerning the convergence of the numerical results in BEM.

As an example, we consider the array of 3 × 2 cylindrical heaving WECs of radius a = 10 m and draft d/a = 1.5, arranged as illustrated in Figure 12, in the middle of the variable bathymetry region (a smooth upslope with max bottom slope 7%), and operating in waves at the same frequency as before $\omega \sqrt{h_m/g}=1.5$, $\omega \sqrt{a/g}=0.8$. The horizontal spacing of the floaters along the $x_1$ and $x_2$ axes is $s_1/a=5,s_2/a=4$. In this case the ratio of the WEC spacing with respect to the characteristic wavelength in the area of the array is small (less than 50%) and thus, the interaction between the floaters is strong. In order to illustrate the applicability of the present BEM, a mesh is used, consisting of 6 × 40 elements on each WEC and 5 × 40 elements on the nearest part of the free surface around each WEC and 30 × 100 elements on the outer part. This includes the absorbing layer, and 14 × 20 elements on the bottom surface (see Figure 12). The total number of elements is 5920.

The propagating field over the shoaling region, for normally and for 45deg obliquely incident waves is shown in Figure 13. This field represents the available wave energy in the domain for possible extraction. The responses of the above array of cylindrical heaving WECs are then calculated, using, as before, the values of $B_S\stackrel{}{/b_m}$= 5, 10, 20 (where $2\pi b_m/m{\omega }_R=0.12$) to model the Power Take Off system for heving floaters. The results calculated by the present BEM approach, at the same frequency as before $\omega \sqrt{h_m/g}=1.5$, $\omega \sqrt{a/g}=0.8$, both for normal and 4deg incident waves over the variable bathymetry region are represented in Figure 14. In the specific arrangement and operating conditions the q-factor decreases with increasing PTO damping and ranges from 90–70%, for normal incident waves, and drops down to 75–60%, for 45 degrees obliquely incident waves.

7. Conclusions

In this work a numerical method is presented for the hydrodynamic analysis of the floating bodies over general seabed topography supporting the calculation of the wave power absorption by single WECs and the performance of arrays of devices in nearshore and coastal regions. As a first step, in order to subsequently formulate and solve 3D diffraction and radiation problems for floaters in the inhomogeneous domain, the present approach is based on the coupled-mode model for the calculation of the wave field propagating over the variable bathymetry region. The results are subsequently used for the hydrodynamic analysis of floating bodies over general bottom topography by means of a low-order BEM. Extensions of the present method supporting the estimation of single WEC and WEC array performance in variable bathymetry regions have been discussed. Future work will be focused on the validation of the present method by comparisons with other methods and experimental laboratory data. Moreover, phase-averaged models like SWAN can treat macroscopically WEC-array effects by including energy pumping in the energy balance equations by using sinks with specific intensity; see, e.g., [40]. Future work will examine the possibility of coupling phase-averaged with the present phase-resolving model, by using coupling schemes as the ones presented in [41,42] and present comparisons at the scale of an array in a nearshore region.