1. Introduction
A large number of concepts of wave energy converters (WECs) have been proposed since the 1790s (Clément et al. Reference Clément2002; Drew, Plummer & Sahinkaya Reference Drew, Plummer and Sahinkaya2009), and most of them are composed of rigid bodies. There are theoretical limits for wave power absorption of the rigid-body-based WECs. For example, the maximum time-averaged power generated by a heaving WEC, which is composed of an axisymmetric rigid body, in monochromatic water waves of wavelength 
$\lambda$, is equivalent to that contained in a 
$\lambda /2{\rm \pi}$ length of an incident wave crest, i.e. the maximum wave power capture width is 
$\lambda /2{\rm \pi}$ (Budal & Falnes Reference Budal and Falnes1975; Evans Reference Evans1976; Newman Reference Newman1976). A notable feature is that the axisymmetric heaving WEC can be reduced in size without affecting the capture width, provided that the amplitude of heaving motion increases simultaneously. The capture width limit (
$\lambda /2{\rm \pi}$) for the heaving WEC also applies for an offshore axisymmetric oscillating water column device, which harnesses energy from the oscillation of the water inside a chamber caused by the action of waves (Martins-Rivas & Mei Reference Martins-Rivas and Mei2009). For attenuator devices, which operate parallel to the wave direction and capture energy from the relative motion of interconnecting floats (e.g. Pelamis; Yemm et al. Reference Yemm, Pizer, Retzler and Henderson2012), the theoretical limits are less clear, although a similar principle applies: the wave power from a much greater length of incident wave crest than the slender width of the device is possibly harnessed, and the maximum wave power capture width is proportional to wavelength 
$\lambda$ (Newman Reference Newman1979). The above-mentioned theoretical maximum wave capture width may be unobtainable for WECs in a small size and/or oscillating in long waves since the amplitude of the device motions can be substantial and must be constrained. Some practical considerations – e.g. a global constraint (Evans Reference Evans1981; Pizer Reference Pizer1993), a constraint on even and odd modes (Newman Reference Newman1979), and a relative rotation constraint (Zheng, Zhang & Sheng Reference Zheng, Zhang and Sheng2016) – can be applied to the motion to ensure that predictions remain within the limitations of the underlying theory.
In the recent decade, there has been a growing trend towards flexible-body-based WECs, e.g. elastic plate WECs (Renzi Reference Renzi2016), SQ devices (Kurniawan et al. Reference Kurniawan, Chaplin, Greaves and Hann2017), and bulge wave devices (Farley, Rainey & Chaplin Reference Farley, Rainey and Chaplin2012). The flexible WECs are believed to offer improved reliability/survivability and reduced cost superior to steel/concrete alternatives (Collins et al. Reference Collins, Hossain, Dettmer and Masters2021; Renzi et al. Reference Renzi, Michele, Zheng, Jin and Greaves2021). Moreover, they are advantageous for the larger potential of wave power absorption. Porter, Zheng & Greaves (Reference Porter, Zheng and Greaves2021) demonstrated that the capture width for axisymmetric WECs could be extended theoretically without bound through the use of generalised (non-rigid-body) modes of motion. Their concept was applied to two flexible cylinders whose surface is surrounded by an array of narrow vertical absorbing paddles, and capture widths in excess of 
$4\lambda /{\rm \pi}$ – which is eight times as large as the maximum capture width of an axisymmetric heaving WEC – were reported in their computations. The advantages of flexible attenuator WECs, which operate in ‘generalised modes’ of motion, in wave power absorption were also reported by Newman (Reference Newman1979, Reference Newman1994), Mei (Reference Mei2014) and Ancellin et al. (Reference Ancellin, Dong, Jean and Dias2020).
To make full use of the ‘generalised modes’ of the flexible WECs in harnessing wave power, a continuum distributed power take-off (PTO) system, or alternatively, a discrete PTO system consisting of a series of PTO units distributed on the flexible components of the WECs, should be adopted. Piezoelectric materials, which are able to establish a voltage upon application of an external stress, can work as a continuum PTO system in wave power conversion. A flexible substrate with both faces perfectly bonded with piezoelectric layers can work as a simple flexible plate WEC. As water waves propagate through the device, the tension variations at the plate–water interface are converted into a voltage, owing to the piezoelectric effect. Renzi (Reference Renzi2016) proposed a hydroelectromechanical-coupled two-dimensional (2-D) model to evaluate the wave power extraction from the dynamics of a submerged piezoelectric plate WEC, which is made of a 2-D flexible bimorph plate with both ends clamped. The device was found to extract a significant amount of energy and could be used for low-power applications such as supplying LEDs and sensors. His model was later extended to study the performance of a 2-D piezoelectric plate WEC mooring on a seabed-mounted/pile-supported breakwater (Buriani & Renzi Reference Buriani and Renzi2017; Zheng et al. Reference Zheng, Meylan, Zhang, Iglesias and Greaves2021a,Reference Zheng, Meylan, Zhang, Iglesias and Greavesb) and an offshore three-dimensional (3-D) piezoelectric circular plate WEC (Zheng et al. Reference Zheng, Greaves, Meylan and Iglesias2020a). Soft dielectric elastomer generators may be embedded into the tube material for continuous energy extraction of the attenuator WECs, e.g. SBM Offshore's S3 WEC (Jean et al. Reference Jean, Wattez, Ardoise, Melis, Kessel, Fourmon, Barrabino, Heemskerk and Queau2012) and the Electric Eel proposed by AWS (AWS Ocean Energy Ltd 2016), which are at the early stages of development. Compared with the continuum PTO system, the discrete PTO system has a lower technical threshold and is closer to engineering application. The ‘Wave Carpet’ device is a flexible WEC with the employment of a discrete PTO system (Alam Reference Alam2012) – it consists of a long, submerged elastic plate equipped with equally spaced PTO units, which are composed of springs and dampers to harvest energy. To evaluate the performance of the Wave Carpet, Desmars et al. (Reference Desmars, Tchoufag, Younesian and Alam2018) developed a wave-flexible structure interaction model through modal decomposition of the structure deformations, and determined the hydrodynamic coefficients of each deformation mode with a boundary element method. The wave power absorbed by the device was found to be more sensitive to the location of the PTO units than to their damping coefficients. Michele et al. (Reference Michele, Buriani, Renzi, van Rooij, Jayawardhana and Vakis2020) considered a 2-D flexible plate floating on the sea surface connected to a series of linear PTO units under the water. To investigate the wave power extraction by the device, a 2-D theoretical model was proposed based on the dry modes and eigenfunction matching method, and each PTO unit was modelled as a linear damper. It was revealed that the system's overall efficiency improved with an increasing number of PTO units.
As one of the simplest flexible WECs, the flexible plate WEC is expected to have a range of potential applications. Indeed, many scholars have investigated and reported the water wave interaction with floating or submerged elastic plates with applications in different fields. Meylan & Squire (Reference Meylan and Squire1996) developed two independent methods – an expansion in the eigenfunctions of a disk, and a more general method of eigenfunctions adopted to construct a Green's function for the disk, to predict the behaviour of a circular flexible ice floe in long-crested sea waves. Zilman & Miloh (Reference Zilman and Miloh2000) obtained a 3-D closed-form solution based on the angular eigenfunction matching method for wave scattering an elastic disk floating in shallow water. Their model was later extended by Peter, Meylan & Chung (Reference Peter, Meylan and Chung2004) to study water wave interaction with an elastic disk floating in finite-depth water. Montiel et al. (Reference Montiel, Bennetts, Squire, Bonnefoy and Ferrant2013a,Reference Montiel, Bonnefoy, Ferrant, Bennetts, Squire and Marsaultb) carried out experimental and theoretical studies to examine the flexural response of one or two floating elastic disks in regular waves. Meylan (Reference Meylan2019) computed the time-dependent vibration of floating elastic plates subject to a transient force. An eigenfunction matching method was employed first to solve the forced harmonic motion of the elastic plate, based on which the time-domain solution was then obtained using the Fourier transform. In addition to the hydroelastics of the floating elastic plates, the scattering from submerged elastic plates has also been investigated by some researchers. Mahmood-Ul-Hassan, Meylan & Peter (Reference Mahmood-Ul-Hassan, Meylan and Peter2009) proposed a theoretical model to solve hydroelastics of a submerged elastic plate by using the eigenfunction matching method. The properties of the wavenumbers of the eigenfunctions above and below the elastic plate were discussed. Later, Williams & Meylan (Reference Williams and Meylan2012) adopted the Wiener–Hopf and residue calculus techniques to solve the water scattering problem of a 2-D submerged semi-infinite elastic plate. Additionally, the wave scattering of porous elastic plates was studied, but not extensively, e.g. see Behera & Sahoo (Reference Behera and Sahoo2015), Koley, Mondal & Sahoo (Reference Koley, Mondal and Sahoo2018) and Selvan et al. (Reference Selvan, Gayathri, Behera and Meylan2021).
In this paper, an offshore floating elastic disk-shaped WEC is considered. The device is composed of a disk-shaped floating elastic plate moored to the seabed through a series of PTO units. The present work is also motivated by the potential integration of wave power with offshore floating flexible solar photovoltaic (PV) farms (Bjørneklett Reference Bjørneklett2018) (see figure 1). The integration can be achieved easily by deploying PTO systems below the floating flexible solar PV farms. The advantages of integrating wave power into floating solar power include but are not limited to: (1) it enables the device to not only capture power from solar radiation but also absorb wave power from water waves, which would efficiently improve the overall power output; (2) the wave attenuation due to wave power absorption weakens the dynamic response of the solar panels, providing benefits for solar power absorption and structural survivability; (3) the wave power device and floating solar share the floating foundation, bringing cost-sharing benefits, including construction, installation and maintenance; (4) the floating solar cannot work during the night, whereas the wave power device is able to work continuously provided that there are incoming waves, which would effectively ensure the continuous and stable power output of the whole power generation system.
Figure 1. Floating solar: (a) 0.5 MWp floating PV system made by Ocean Sun, owned by Statkraft on their 72 MWp hydro power dam at Banja, Albania (Photo: Ocean Sun); (b) 0.22 MWp floating PV system made by Ocean Sun, owned by SNAP (SN–Aboitiz Power) located at the Magat hydro power reservoir at Isabela, Philippines (Photo: Ocean Sun).
To evaluate the hydrodynamics of a floating flexible circular WEC, a theoretical model was developed by Michele, Zheng & Greaves (Reference Michele, Zheng and Greaves2022) based on free-edge dry mode expansion of plate motion. The PTO units in their model were distributed symmetrically about the incident wave direction, and the response of the plate was expanded into heave mode, pitch mode, and a series of flexible modes, with each treated as a generalised mode in solving the wave radiation problems. The flexible WEC was found to have larger efficiency than rigid devices. Yet there are still many fundamental questions to be addressed, e.g. how the incident wave direction affects wave power absorption of the elastic disk-shaped WEC, what if there are a large number of PTO units, and how to optimise the PTO system. In this paper, a theoretical model based on the dispersion relation and eigenfunction matching method is developed to study the hydroelastic characteristics and wave power absorption of the device, and meanwhile to address unanswered questions including but not limited to the above ones.
The remainder of this paper is structured as follows. The mathematical model is outlined in § 2. The convergence analysis and model validation are supplied in § 3. A multiparameter study is carried out with the validated model, the results of which can be found in § 4. Finally, § 5 outlines the conclusions.
2. Mathematical model
The elastic disk with radius 
$R$ is floating on water of depth 
$h$ (see figure 2). A local cylindrical coordinate system 
$Or \theta z$ is chosen, with the mean free surface coinciding with the plane 
$z=0$, and 
$Oz$ coinciding with the vertical axis of the disk and pointing upwards. The device is subjected to a train of regular waves propagating in the direction 
$\beta$ relative to the positive 
$Ox$ axis with amplitude 
$A$ and angular frequency 
$\omega$. The PTO system is composed of 
$N$ discrete PTO units, which are distributed all over the circle 
$r=r_0$ below the disk. Each PTO unit consists of a linear damper and a linear stiffness/mass. The damping and stiffness/mass of the 
$n$th PTO unit is denoted as 
$c_n\in \mathbb {R}^++\mathrm {i}\mathbb {R}$ (units of N s m
$^{-1}$), where the real part is positive and associated with the damping, whereas the imaginary part is related to the stiffness/mass. The stiffness and mass play opposite roles in affecting the motions response of the device. For simplicity, it is assumed that either stiffness or mass is vanishing in the PTO system. If 
${\rm Im}(c_n)>0$, then there is no mass in the PTO system, and the stiffness is 
$\omega \,{\rm Im}(c_n)$; if 
${\rm Im}(c_n)<0$, then there is no stiffness in the PTO system, and the mass is 
$-{\rm Im}(c_n)/\omega$. The 
$n$th PTO unit is located at 
$(r,\theta )=(r_0,\theta _n)$, with 
$\theta _n=2{\rm \pi} (n-1)/N$, unless specified otherwise. When 
$N$ is large enough, it is assumed that the discrete PTO system may be represented by a continuum PTO system, the damping and stiffness/mass of which can be denoted as 
$c(\theta )$ (units of N s m
$^{-2}$). This assumption will be tested later, in § 3.
Figure 2. Sketch of a floating elastic disk-shaped WEC with a discrete PTO system (
$N=4$): (a) bird's-eye view, (b) top view, (c) side view.
The fluid domain is divided into three regions (figure 2a): Region 1 is the interior region below the disk with 
$r\in [0,r_0]$ and 
$z\in [-h,0]$; Region 2 is the annular region below the disk with 
$r\in [r_0,R]$ and 
$z\in [-h,0]$; and Region 3 is the exterior region extending from Region 2 to infinity horizontally, with 
$r\in [R,\infty )$ and 
$z\in [-h,0]$.
2.1. Governing equation and boundary conditions
The velocity potential is defined as 
$\varPhi (x,y,z,t)=\mathrm {Re}[\phi (x,y,z)\,\mathrm {e}^{-\mathrm {i}\omega t}]$, where 
$\phi$ is the spatial potential, 
$\omega$ is the angular frequency of oscillation, and 
$t$ is the time. The spatial velocity potential 
$\phi$ is a solution of the governing equations
\begin{equation} (\partial_x^2+\partial_y^2+\partial_z^2)\phi=0 \end{equation}in the fluid domain, with
\begin{equation} \partial_z\phi=0,\quad z={-}h, \end{equation}on the seabed,
\begin{equation} {-K}\phi+\frac{\partial\phi}{\partial z}=0,\quad z=0,\ r>R, \end{equation}
at the water surface of the exterior region, where 
$K={\omega ^2}/{g}$, in which 
$g$ represents the acceleration of gravity, and
\begin{equation} \frac{\partial\phi}{\partial z}+\mathrm{i} \omega \eta=0,\quad z=0,\ r< R, \end{equation}
at the the lower surface of the disk, in which 
$\eta$ denotes the complex deflection of the disk.
Additionally, the dynamic equation of the elastic disk should also be satisfied.
For the discrete PTO system, the dynamic equation may be written as
\begin{equation} g\left[\chi\varDelta^2+1-K\gamma\right]\eta-\mathrm{i}\omega\phi|_{z=0}=\frac{1}{\rho r} \sum_{n=1}^N F_n\,\delta(r-r_0)\,\delta(\theta-\theta_n), \quad r< R, \end{equation}
where 
$\gamma$ and 
$\chi$ denote the mass per unit area and the flexural rigidity of the disk, respectively, scaled with respect to the water density; 
$\varDelta$ is the Laplacian operator in the horizontal plane; 
$\rho$ denotes the water density; and 
$\delta$ denotes the Dirac delta function. Note that in the cylindrical coordinates, we divide the Dirac delta function by 
$r$ because we multiply by 
$r$ when we integrate in cylindrical coordinates. Here, 
$F_n$ denotes the unknown PTO force due to the 
$n$th PTO unit acting on the disk at 
$(r,\theta )=(r_0,\theta _n)$, and satisfies
\begin{equation} F_n=\mathrm{i} \omega c_n \eta, \quad (r,\theta)=(r_0,\theta_n), \end{equation}
for 
$n=1,2, \ldots, N$.
When 
$N$ is large and the width of the interval, 
$2{\rm \pi} r_0/N$, is small with respect to the wavelength and the radius of the disk, we assume that the discrete PTO system may be represented by a continuum PTO system, hence 
$F_n$ may be replaced by a discrete evaluation 
$2{\rm \pi} r_0\,F(\theta _n)/N$ of the continuous function 
$F(\theta )$ (units of N m
$^{-1}$), and 
$c_n=2{\rm \pi} r_0\,c(\theta _n)/N$, resulting in
\begin{equation} \sum_{n=1}^{N}\frac{2{\rm \pi}}{N}\,\delta(\theta-\theta_n)\,F(\theta_n)\rightarrow \int_{0}^{2{\rm \pi}}\delta(\theta-\theta_n)\,F(\theta_n)\,{\rm d}\theta_n=F(\theta), \end{equation}and allowing (2.5) and (2.6) to be approximated by
\begin{equation} g\left[\chi\varDelta^2+1-K\gamma\right]\eta-\mathrm{i}\omega\phi|_{z=0}=\frac{1}{\rho}\, F(\theta)\,\delta(r-r_0), \quad r< R, \end{equation}and
\begin{equation} F(\theta)=\mathrm{i} \omega\,c(\theta)\,\eta, \quad r=r_0, \end{equation}
respectively. Here, 
$c$ may be non-dimensionalised as 
$\bar {c}=c/(\rho R \sqrt {gh})$. Functions 
$F(\theta )$ and 
$c(\theta )$ can be further expanded into
\begin{equation} F(\theta)=\sum_{m={-}\infty}^{\infty}f_m\,\mathrm{e}^{\mathrm{i} m\theta},\quad c(\theta)=\sum_{m={-}\infty}^{\infty}c^{(m)}\,\mathrm{e}^{\mathrm{i} m\theta}, \end{equation}
where 
$f_m$ are unknown coefficients to be determined, and 
$c^{(m)}$ are the Fourier coefficients of 
$c(\theta )$ that can be obtained from 
$c^{(m)}=({1}/{2 {\rm \pi}})\int _{-{\rm \pi} }^{{\rm \pi} }c(\theta )\, \mathrm {e}^{-\mathrm {i} m \theta } \,{\rm d}\theta$.
With the employment of the Laplace equation as given in (2.1), the kinematic and dynamic conditions as given in (2.4) and (2.5) or (2.8) for 
$r\in [0,r_0)\cup (r_0,R]$ can be combined into
\begin{equation} [\chi\partial_z^4+1-K\gamma]\,\partial_z\phi- K\phi=0. \end{equation} Additionally, in the far field horizontally, the scattered wave potential 
$\phi _{S}=\phi -\phi _{I}$, where 
$\phi _{I}$ is the velocity potential of the undisturbed incident waves whose expression will be given in § 2.2, is subject to the Sommerfeld radiation condition.
The boundary conditions at the free edge of the disk should be satisfied as well. For a free edge, both moment and shearing stresses vanish at the edge, providing
\begin{equation} \left[\varDelta-\frac{1-\upsilon}{r}\left(\frac{\partial}{\partial r}+\frac{1}{r}\,\frac{\partial^2}{\partial \theta^2}\right)\right]\eta=0 \end{equation}and
\begin{equation} \left[\frac{\partial}{\partial r}\varDelta-\frac{1-\upsilon}{r^2}\left(-\frac{\partial}{\partial r}+\frac{1}{r}\right)\frac{\partial^2}{\partial \theta^2}\right]\eta=0, \end{equation}
in which 
$\upsilon$ denotes the Poisson ratio, and a typo of the item 
${(1-\upsilon )\, \partial ^3\eta }/{ r^2\,\partial r\,\partial ^2\theta }$ as derived by Meylan (Reference Meylan2019) has been corrected.
2.2. Expression of the velocity potential
In Region 
$1$, the velocity potential may be expressed as
\begin{equation} \phi_1(r,\theta,z)=\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty A_{m,l}\,{\rm J}_m(\kappa_lr)\,Y_l(z)\, \mathrm{e}^{\mathrm{i} m\theta}, \end{equation}
where 
$A_{m,l}$ are the unknown coefficients to be determined, 
${\rm J}_m$ denotes the Bessel function of the 
$m$th order, 
$Y_l={\cosh [\kappa _l(z+h)]}/{\cosh (\kappa _lh)}$, and 
$\kappa _l$ for 
$l=-2, -1,0,1,2,\ldots$ are the roots of the dispersion relation for the interior region,
\begin{equation} \left[\chi\kappa_l^4+1-K\gamma\right]\kappa_l\tanh(\kappa_lh)=K. \end{equation} Here, 
$\kappa _0\in \mathbb {R}^+$ and 
$\kappa _l\in \mathrm {i} \mathbb {R}^+$ for 
$l =1, 2,3, \ldots$ can be obtained, which support the propagating waves and evanescent waves, respectively. The remaining two roots, 
$\kappa _{-2}$ and 
$\kappa _{-1}$, support damped propagating waves, and satisfy 
$\kappa _{-1}\in \mathbb {R}^++\mathrm {i}\mathbb {R}^+$ and 
$\kappa _{-2}=-\kappa _{-1}^*$, in which * denotes the complex conjugate.
In Region 
$2$, the velocity potential may read
\begin{equation} \phi_2(r,\theta,z)=\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty [B_{m,l}\,{\rm J}_m(\kappa_lr)+C_{m,l}\, H_m(\kappa_lr)]\,Y_l(z)\,\mathrm{e}^{\mathrm{i} m\theta}, \end{equation}
where 
$B_{m,l}$ and 
$C_{m,l}$ are the unknown coefficients to be determined, and 
$H_m={\rm J}_m +\mathrm {i} Y_m$ is the 
$m$th-order Hankel function of the first kind.
In Region 
$3$, the velocity potential can be written as
\begin{equation} \phi_3=\phi_{I}+\sum_{m={-}\infty}^\infty\sum_{l=0}^\infty D_{m,l}\,H_m(k_lr)\,Z_l(z)\,\mathrm{e}^{\mathrm{i} m\theta}, \end{equation}
where the accumulative term denotes the scattered wave potential 
$\phi _{S}$, 
$D_{m,l}$ are the unknown coefficients to be determined, 
$Z_l(z)={\cosh [k_l(z+h)]}/{\cosh (k_lh)}$. Here, 
$k_0\in \mathbb {R}^+$ and 
$k_l\in \mathrm {i}\mathbb {R}^+$ for 
$l=1, 2, 3, \ldots$ support the propagating waves and evanescent waves, respectively, for the exterior region, and they are the roots of the dispersion relation in the open water:
\begin{equation} K=k_l\tanh(k_lh). \end{equation}
For simplicity, we define 
$k=k_0$. Here, 
$\phi _{I}$ denotes the undisturbed incident wave velocity potential, which can be expressed as
\begin{equation} \phi_{I} (x,y,z)={-}\frac{\mathrm{i} gA}{\omega}\,Z_0(z)\exp({\mathrm{i} k(x\cos \beta+y \sin \beta)}) \end{equation}and
\begin{equation} \phi_{I} (r,\theta,z)={-}\frac{\mathrm{i} gA}{\omega}\,Z_0(z) \sum_{m={-}\infty}^\infty \mathrm{i}^m\,\mathrm{e}^{-\mathrm{i} m \beta}\, {\rm J}_m(kr)\,\mathrm{e}^{\mathrm{i} m\theta}, \end{equation}
in the general Cartesian coordinate system 
$Oxyz$ and the local cylindrical coordinate system 
$Or\theta z$, respectively.
2.3. Expressions of the deflection of the disk and the free surface response
After inserting (2.14) and (2.16) into (2.4), we have the deflection of the disk expressed as
\begin{equation} \eta= \begin{cases} \displaystyle \frac{\mathrm{i} \omega}{g}\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{A_{m,l}\,{\rm J}_m(\kappa_lr)}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}, & r< r_0,\\ \displaystyle \frac{\mathrm{i} \omega}{g}\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{B_{m,l}\,{\rm J}_m(\kappa_lr)+C_{m,l}\,H_m(\kappa_lr)}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}, & r_0< r< R. \end{cases} \end{equation} The response of the free surface approximately 
$z = 0$, 
$\eta _w$, can be expressed as
\begin{equation} \eta_w=\frac{\mathrm{i} \omega}{g}\,\phi_3|_{z=0}=A \exp({\mathrm{i} k(x\cos \beta+y \sin \beta)})+\frac{\mathrm{i} \omega}{g} \sum_{m={-}\infty}^\infty\sum_{l=0}^\infty D_{m,l}\,H_m(k_lr)\,\mathrm{e}^{\mathrm{i} m\theta}. \end{equation}2.4. Solution of the unknown coefficients
The velocity potentials and the deflection of the disk should satisfy the following continuity conditions.
(i) Continuity of pressure at the boundary
$r=r_0$:
(2.23)\begin{equation} \phi_1=\phi_2, \quad r=r_0. \end{equation}(ii) Continuity of radial velocity at the boundary
$r=r_0$:
(2.24)\begin{equation} \frac{\partial\phi_1}{\partial r}=\frac{\partial\phi_2}{\partial r}, \quad r=r_0. \end{equation}(iii) Continuity of pressure at the boundary
$r=R$:
(2.25)\begin{equation} \phi_2=\phi_3, \quad r=R. \end{equation}(iv) Continuity of radial velocity at the boundary
$r=R$:
(2.26)\begin{equation} \frac{\partial\phi_2}{\partial r}=\frac{\partial\phi_3}{\partial r}, \quad r=R. \end{equation}(v) Continuity of the deflection of the disk at the boundary
$r=r_0$:
(2.27)\begin{equation} \eta|_{r=r_0^-}=\eta|_{r=r_0^+}. \end{equation}(vi) Continuity of the first derivative of the disk deflection at the boundary
$r=r_0$:
(2.28)\begin{equation} \left.\frac{\partial \eta}{\partial r}\right|_{r=r_0^-}=\left.\frac{\partial \eta}{\partial r}\right|_{r=r_0^+}. \end{equation}(vii) Continuity of the second derivative of the disk deflection at the boundary
$r=r_0$:
(2.29)\begin{equation} \left.\frac{\partial^2 \eta}{\partial r^2}\right|_{r=r_0^-}=\left.\frac{\partial^2 \eta}{\partial r^2}\right|_{r=r_0^+}. \end{equation}(viii) A jump in
$\varDelta ^2\eta$ due to the PTO forces at the boundary $r=r_0$:
(2.30)and\begin{equation} \varDelta^2 \eta \left|_{r=r_0^+}-\varDelta^2 \eta \right|_{r=r_0^-}= \frac{1}{\rho g r \chi} \sum_{n=1}^N F_n\,\delta(r-r_0)\,\delta(\theta-\theta_n) \end{equation}(2.31)for discrete and continuum PTO systems, respectively, which are derived by using the dynamic equations, i.e. (2.5) and (2.8), together with the continuity conditions of the hydrodynamic pressure and disk deflection at the PTO systems, i.e. (2.23) and (2.27).\begin{equation} \varDelta^2 \eta \left|_{r=r_0^+}- \varDelta^2 \eta \right|_{r=r_0^-}= \frac{1}{\rho g\chi }\,F(\theta)\,\delta(r-r_0), \end{equation}
The continuity conditions, i.e. (2.23)–(2.31), together with the free-edge condition, i.e. (2.12) and (2.13), and the PTO system related equation, i.e. (2.6) or (2.9) depending on discrete or continuum PTO system, can be used to derive a complex linear matrix equation by using the orthogonality characteristics of 
$Z_l(z)$ and 
$\mathrm {e}^{\mathrm {i}m\theta }$, and the eigenfunction matching method. The unknown coefficients 
$A_{m,l}$, 
$B_{m,l}$, 
$C_{m,l}$, 
$D_{m,l}$, and 
$F_n$ or 
$f_m$, can then be calculated by solving the complex linear matrix equation. Detailed derivation and calculations for the unknown coefficients are given in Appendix A.
2.5. Wave power absorption
2.5.1. Direct method
The time-averaged wave power absorption can be evaluated in a straightforward manner, which is named the ‘direct method’ in this paper, by accumulating/integrating the power absorbed by each PTO unit.
For a discrete PTO system, the absorbed wave power can be written as the sum of the power absorbed by the 
$N$ PTO units (Michele et al. Reference Michele, Zheng and Greaves2022):
\begin{equation} P=\frac{\omega^2}{2}\sum_{n=1}^{N}{\rm Re}(c_n)\left|\eta|_{r=r_0,\theta=\theta_n}\right|^2=\frac{1}{2}\sum_{n=1}^{N}\frac{{\rm Re}(c_n)}{|c_n|^2}\,|F_n|^2. \end{equation}The wave power captured by a continuum PTO system may be expressed as
\begin{equation} P=\frac{r_0 \omega^2}{2}\int_0^{2{\rm \pi}}{\rm Re} (c(\theta))\left|\eta|_{r=r_0}\right|^2 {\rm d}\theta=\frac{r_0}{2}\int_0^{2{\rm \pi}}\frac{{\rm Re} (c(\theta))}{|c(\theta)|^2}\,|F(\theta)|^2 \,{\rm d}\theta. \end{equation} For an angular-independent PTO system, i.e. 
$c(\theta )=c^{(0)}$, (2.33) gives
\begin{equation} P=\frac{{\rm Re} (c^{(0)})\,r_0}{2 |c^{(0)}|^2}\int_0^{2{\rm \pi}}\left|\sum_{m={-}\infty}^{\infty}f_m\, \mathrm{e}^{\mathrm{i} m\theta}\right|^2 \,{\rm d}\theta=\frac{{\rm Re} (c^{(0)})\,r_0 {\rm \pi}}{|c^{(0)}|^2}\sum_{m={-}\infty}^{\infty}|f_m|^2. \end{equation}The power absorption of the device can be measured using the dimensionless wave power capture factor, defined as
\begin{equation} \eta_{e}=\frac{k P}{P_{in}}, \end{equation}
where 
$P_{in}=({\rho g A^2}/{2})c_g$ denotes the incoming wave power per unit width of the wave front, in which 
$c_g$ represents the group velocity,
\begin{equation} c_g=\frac{\omega}{2 k}\left[1+\frac{2kh}{\sinh(2kh)}\right]. \end{equation}2.5.2. Indirect method
At the far field with 
$kr$ large enough, only propagating modes of the velocity potential exist, and the velocity potential can be expressed in terms of incoming waves, 
$\phi _{in}$, and outgoing waves, 
$\phi _{out}$, as (Porter et al. Reference Porter, Zheng and Greaves2021)
\begin{align} \phi_3&={-}\frac{\mathrm{i} gA}{\omega}\,Z_0(z) \sum_{m={-}\infty}^\infty \mathrm{i}^m\,\mathrm{e}^{-\mathrm{i} m \beta}\,{\rm J}_m(kr)\, \mathrm{e}^{\mathrm{i} m\theta}+\sum_{m={-}\infty}^\infty D_{m,0}\,H_m(k r)\,Z_0(z)\,\mathrm{e}^{\mathrm{i} m\theta}\nonumber\\ & =\underbrace{-\frac{\mathrm{i} gA}{2 \omega}\,Z_0(z) \sum_{m={-}\infty}^\infty \mathrm{i}^m\,\mathrm{e}^{-\mathrm{i} m \beta}\, H_m^{(2)}(kr)\,\mathrm{e}^{\mathrm{i} m\theta}}_{\phi_{in}\ \mathrm{incoming waves}}\nonumber\\ &\quad \underbrace{{}-\frac{\mathrm{i} gA}{2 \omega}\,Z_0(z) \sum_{m={-}\infty}^\infty \mathrm{i}^m\left( \mathrm{e}^{-\mathrm{i} m \beta}+\frac{2\omega \mathrm{i}^{1-m}}{gA}\,D_{m,0}\right) H_m(kr)\,\mathrm{e}^{\mathrm{i} m\theta}}_{\phi_{out}\ \mathrm{outgoing\ waves}}, \end{align}
in which 
$H_m^{(2)}={\rm J}_m - \mathrm {i} Y_m$ is the 
$m$th-order Hankel function of the second kind, with the asymptotic forms of 
$H_m^{(2)}(kr)\rightarrow \sqrt {2/{\rm \pi} } \exp ({\mathrm {i} (m{\rm \pi} /2+{\rm \pi} /4)})(kr)^{-1/2}\, \mathrm {e}^{-\mathrm {i} kr}$ when 
$kr\rightarrow \infty$.
The power lost to the device can be expressed as the difference between the incoming wave power coming towards the device and the outgoing wave power propagating away from the device (Zheng et al. Reference Zheng, Greaves, Meylan and Iglesias2020a):
\begin{align} P&=\frac{1}{2}\,{\rm Re}\int_{{-}h}^0\int_0^{2{\rm \pi}}\mathrm{i} \omega \rho \left. \left[\left(\phi_{in}\,\frac{\partial \phi_{in}^*}{\partial r}-\phi_{out}\,\frac{\partial \phi_{out}^*}{\partial r}\right) r\right]\right|_{r\rightarrow \infty} \,{\rm d}z \,{\rm d} \theta\nonumber\\ &=\frac{\rho g A^2 c_g}{2k}\sum_{m={-}\infty}^{\infty}\left(1-\left|\mathrm{e}^{-\mathrm{i} m \beta}+\frac{2\omega \mathrm{i}^{1-m}}{gA}\,D_{m,0}\right|^2\right), \end{align}which is named the ‘indirect method’ in this paper. This expression can also be derived by using Green's theorem (e.g. see Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005).
The corresponding wave power capture factor can be written as
\begin{equation} \eta_{e}=\sum_{m={-}\infty}^{\infty}\left(1-\left|\mathrm{e}^{-\mathrm{i} m \beta}+\frac{2\omega \mathrm{i}^{1-m}}{gA}\,D_{m,0}\right|^2\right), \end{equation}regardless of the type of PTO system.
Compared with the ‘direct method’ for the angular-dependent continuum PTO system, i.e. (2.33), which includes an angular integral that needs to be solved numerically, the ‘indirect method’ as given in (2.38) gives a more straightforward and accurate evaluation of the wave power absorption.
2.6. Optimisation of the angular-independent continuum PTO system
The performance of the device applying an angular-independent continuum PTO system (i.e. 
$c(\theta )=c^{(0)}$) is independent of incident wave direction. For the sake of convenience, here we assume 
$\beta =0$, hence the response of the device and the wave field are symmetric about the plane of 
$Oxz$, giving 
$D_{-m,0}=(-1)^m\,D_{m,0}$. Correspondingly, (2.39) can be rewritten as
\begin{equation} \eta_{e}=\sum_{m=0}^{\infty}\eta_m=\sum_{m=0}^{\infty}\epsilon_m\left(1-\left|1+\frac{2\omega \mathrm{i}^{1-m}}{gA}\,D_{m,0}\right|^2\right), \end{equation}
where 
$\epsilon _m=1$ and 2 for 
$m=0$ and 
$m>0$, respectively. Here, 
$\eta _m$ denotes the contribution of the wave power extracted from the 
$m$th circular component of the wave field:
\begin{equation} \eta_m=\epsilon_m\left(1-\left|1+\frac{2\omega \mathrm{i}^{1-m}}{gA}\,D_{m,0}\right|^2\right)=\frac{\epsilon_m\,{\rm Re} (c^{(0)})\,k r_0 {\rm \pi}}{|c^{(0)}|^2 P_{in}}\,|f_m|^2, \end{equation}in which the right-hand side term after the second equals sign follows from (2.34).
The maximum dimensionless wave power that could be captured from the 
$m$th circular component of the wave field is 
$\eta _{MAX}^{(m)}=\epsilon _m$ if its dynamics can be orchestrated to meet the condition 
$D_{m,0}=({-gA})/{2\omega \mathrm {i} ^{1-m}}$. For a rigid heaving disk-shaped WEC, the device radiates waves in the zeroth circular mode (
$\eta _{e}=\eta _0$), and the maximum theoretical wave power capture factor is 
$\eta _{e}=\eta _{MAX}^{(0)}=1$. For a rigid pitching disk-shaped WEC, the device radiates waves in the first circular mode (
$\eta _{e}=\eta _1$), and the maximum theoretical wave power capture factor is 
$\eta _{e}=\eta _{MAX}^{(1)}=2$. Combined heave and pitch provides a maximum of 
$\eta _{MAX}^{(0)}+\eta _{MAX}^{(1)}=3$.
The elastic disk-shaped WEC is expected to have the capacity to absorb energy in excess of these limits due to its ability to radiate in multiple circular modes.
The value of 
$c$ that maximises the power absorbed in any individual circular mode can be predicted by using the theoretical framework developed in § 2.4. This can be done by imposing 
$D_{m,0}=({-gA})/{2\omega \mathrm {i} ^{1-m}}$ as a known condition for maximum power absorption from the 
$m$th circular wave component, and treating 
$c^{(0)}$ as an unknown coefficient to be determined. An alternative way to determine 
$\eta _{MAX}^{(m)}$ is to express 
$\eta _m$ in terms of 
$c^{(0)}$ and then set 
${\partial \eta _m}/{\partial {\rm Re}(c^{(0)})}=0$ and 
${\partial \eta _m}/{\partial {\rm Im}(c^{(0)})}=0$. The latter method can be applied to evaluate the optimised PTO damping of a device with a fixed PTO stiffness/mass for maximising wave power absorption from the 
$m$th circular wave component by setting 
${\partial \eta _m}/{\partial {\rm Re}(c^{(0)})}=0$ and 
${\rm Im}(c^{(0)})={\rm constant}$, resulting in the corresponding dimensionless wave power denoted as 
$\eta _{{\it max}}^{(m)}$. To express the optimised PTO coefficients and the maximum wave power absorption, the key is to separate 
$c^{(0)}$ from the complicated terms of the expression for 
$f_m$.
After setting 
$\tau =m$ and rearranging (A12)–(A15), (A18)–(A20), (A22)–(A24) and (A28), we have
\begin{equation} \left[
\begin{array}{@{}cc@{}} \boldsymbol{\mathsf{A}}^{(m)}_0 &
\boldsymbol{\mathsf{B}}^{(m)}_0 \\
\boldsymbol{\mathsf{C}}^{(m)}_0 & {1}/{c^{(0)}} \end{array}
\right]_{(4L+11)\times(4L+11)} \left\{
\begin{array}{@{}c@{}} \boldsymbol{\mathsf{X}}_m \\ f_m \end{array}
\right\}_{(4L+11)\times1}= \left\{ \begin{array}{@{}c@{}}
\boldsymbol{\mathsf{F}}_0^{(m)} \\ 0 \end{array}
\right\}_{(4L+11)\times1},
\end{equation}
where we truncate all infinite series of vertical eigenfunctions at 
$L$, i.e. 
$(L+1)$ terms (
$l=0,1,\ldots,L$) for 
$D_{m,l}$, and 
$(L+3)$ terms (
$l = -2, -1, 0, 1,\ldots, L$) for 
$A_{m,l}$, 
$B_{m,l}$ and 
$C_{m,l}$. Here, 
$\boldsymbol{\mathsf{X}}_m$ denotes the vector of the unknown coefficients 
$A_{m,l}$, 
$B_{m,l}$, 
$C_{m,l}$ and 
$D_{m,l}$; the subscript ‘0’ means the matrix/vector is independent of 
$c^{(0)}$.
Using the formulae for the block matrix inversion, we have
\begin{equation}
f_m=\frac{\boldsymbol{\mathsf{C}}^{(m)}_0
\left(\boldsymbol{\mathsf{A}}^{(m)}_0\right)^{{-}1}
\boldsymbol{\mathsf{F}}_0^{(m)} c^{(0)}} {\boldsymbol{\mathsf{C}}^{(m)}_0
(\boldsymbol{\mathsf{A}}^{(m)}_0)^{{-}1} \boldsymbol{\mathsf{B}}^{(m)}_0
c^{(0)}-1}, \end{equation}
with which (2.41) can be rewritten as
\begin{equation} \eta_m=\frac{\epsilon_m\,{\rm Re} (c^{(0)})\,k r_0 {\rm \pi}}{P_{in}}\,\frac{\left|b_0^{(m)}\right|^2} {\left|a_0^{(m)} c^{(0)} -1\right|^2}, \end{equation}
where 
$a_0^{(m)}=\boldsymbol{\mathsf{C}}^{(m)}_0 (\boldsymbol{\mathsf{A}}^{(m)}_0)^{-1} \boldsymbol{\mathsf{B}}^{(m)}_0$ and 
$b_0^{(m)}=\boldsymbol{\mathsf{C}}^{(m)}_0 (\boldsymbol{\mathsf{A}}^{(m)}_0)^{-1} \boldsymbol{\mathsf{F}}_0^{(m)}$.
For any certain value of PTO stiffness/mass (
${\rm Im}(c^{(0)})$) there is a corresponding maximum of absorbed power from the 
$m$th component when 
$\partial {\eta _m}/{\partial {\rm Re}(c^{(0)})}=0$, which occurs if
\begin{equation} {\rm Re}(c^{(0)})= \frac{\left|a_0^{(m)}\,{{\rm Im}(c^{(0)})}+\mathrm{i}\right|}{\left|a_0^{(m)}\right|}\equiv c_{{\it opt}}^{(m)}, \end{equation}resulting in
\begin{equation} \eta_{{\it max}}^{(m)}=\frac{\epsilon_m k r_0 {\rm \pi}}{2 P_{in}}\,\frac{\left|b_0^{(m)}\right|^2} {\left|a_0^{(m)}\right|\left|a_0^{(m)}\,{\rm Im}(c^{(0)})+\mathrm{i}\right|-{\rm Re}(a_0^{(m)})}. \end{equation} This 
$\eta _{{\it max}}^{(m)}$ can be further enlarged with an optimal PTO stiffness/mass, satisfying 
${\partial \eta _m}/{\partial {\rm Re}(c^{(0)})}=0$ and 
${\partial \eta _m}/{\partial {\rm Im}(c^{(0)})}=0$, which gives
\begin{equation} {\rm
Re}(c^{(0)})=\frac{-{\rm Re}(a _0)}{|a_0|^2} \equiv {\rm
Re}(c_{OPT}^{(m)}),\quad {\rm
Im}(c^{(0)})=\frac{-{\rm Im}(a _0)}{|a_0|^2} \equiv {\rm
Im}(c_{OPT}^{(m)}),
\end{equation}
and
\begin{equation} \eta_{MAX}^{(m)}={-}\frac{\epsilon_m k r_0 {\rm \pi}}{4 P_{in}}\,\frac{|b_0^{(m)}|^2}{{\rm Re}(a_0^{(m)})}= \epsilon_m, \end{equation}as expected (e.g. see Porter et al. Reference Porter, Zheng and Greaves2021).
3. Convergence analysis and model validation
Prior to case studies with the present analytical model, a convergence analysis is carried out, and the model is validated by the results documented in the literature. Hereinafter, 
$\chi /h^4=\gamma /h=0.01$ is adopted following Meylan, Bennetts & Peter (Reference Meylan, Bennetts and Peter2017) and Zheng et al. (Reference Zheng, Meylan, Zhu, Greaves and Iglesias2020b), and for the sake of simplification, 
$N=1,2,3, \ldots$ denotes the device with a discrete PTO system consisting of 
$N$ PTO units, whereas 
$N=\infty$ represents the device applying a continuum PTO system. We also assume that the PTO units are identical, and the continuum PTO system is angular-independent, unless otherwise specified, i.e. 
$c_n=2{\rm \pi} r_0 c^{(0)}/N$, and 
$c^{(0)}$ may be non-dimensionalised as 
$\bar {c}=c^{(0)}/(\rho R\sqrt {gh})$.
3.1. Convergence analysis
Figure 3 illustrates the impact of the angular and vertical truncated cut-offs (i.e. terms of 
$M$ and 
$L$) on the frequency response of the wave power capture factor of a flexible floating WEC with 
$N=4$ or 
$N=\infty$. In order to obtain the converged results, 
$M\ge 20$ and 
$L\ge 10$ are suggested. Hereinafter, 
$M= 20$ and 
$L=10$ are adopted.
Figure 3. Frequency response of the wave power capture factor 
$\eta _e$ for 
$R/h=2.0$, 
$r_0/R=0.5$, 
$\beta ={\rm \pi} /6$ and 
$\bar {c}=0.2$: (a) impact of the angular cut-offs in terms of 
$M$ for 
$N=4$ and 
$L=10$; (b) impact of the vertical cut-offs in terms of 
$L$ for 
$N=4$ and 
$M=20$; (c) impact of the angular cut-offs in terms of 
$M$ for 
$N=\infty$ and 
$L=10$; (d) impact of the vertical cut-offs in terms of 
$L$ for 
$N=\infty$ and 
$M=20$.
3.2. Model validation
If the PTO damping coefficient is rather small, i.e. 
$\bar {c} \rightarrow 0$, then the constraint that is induced by the PTO system, regardless of whether it is discrete or continuum, can be neglected. Therefore, the response of the device will be close to that of a free-floating elastic disk. Figure 4 illustrates the theoretical and experimental (Montiel et al. Reference Montiel, Bonnefoy, Ferrant, Bennetts, Squire and Marsault2013b) deflection of four markers (M1–M4) placed on a free-floating elastic disk. The present model is found to be able to well predict the frequency response of the disk motion, although a slight overestimation of the deflection is observed due to the ignorance of the fluid viscosity in the potential flow theory. The difference may also be caused by the presence of constraints and possible viscous-elastic behaviour of real materials in the experimental model. We have also predicted the deflection of a floating flexible WEC for 
$R/h=2.0$, 
$\beta ={\rm \pi}$ and 
$h\omega ^2/g=2.0$ with either a discrete PTO system (
$N=4$) or a continuum PTO system (
$N=\infty$) for 
$\bar {c}=10^{-5}$, the results of which show excellent agreement with those of Meylan et al. (Reference Meylan, Bennetts and Peter2017). In addition, we have evaluated the wave power capture factor of a specific system by using both the direct and indirect methods, and their results are found to agree perfectly with one another, as expected.
Figure 4. Deflection of four markers (M1–M4) placed on a free-floating elastic disk with 
$R= 0.72$ m, 
$h= 1.9$ m, 
$\bar {\chi }=3.55 \times 10^{-4}$ and 
$\bar {\gamma }= 2.79 \times 10^{-3}$ as a function of frequency. Lines represent the present theoretical results with 
$N=4$, 
$r_0/R=0.5$, 
$\bar {c}=10^{-5}$; symbols represent the experimental data (Montiel et al. Reference Montiel, Bonnefoy, Ferrant, Bennetts, Squire and Marsault2013b).
Figure 5 compares the variations of the wave power capture factor with 
$kh$ and 
$\bar {c}$ for devices with a discrete PTO system consisting of a different number of PTO units and with a continuum PTO system, respectively. As expected, the larger the number (
$N$) of the PTO units composed in the discrete PTO system, the closer the wave power capture factors are to those of the continuum PTO. In the examined cases, the results of the discrete PTO system with 
$N=16$ nearly overlap with those of the continuum ones.
Figure 5. Comparison of the wave power capture factors of the devices with discrete and continuum PTO systems, with 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta ={\rm \pi} /6$: (a) frequency response of 
$\eta _e$ for 
$\bar {c}=0.2$; (b) variation of 
$\eta _e$ with 
$\bar {c}$ for 
$kh=4.0$.
The excellent agreement of the results, together with those plotted in figures 4 and 5, gives confidence in the present theoretical model for solving water wave interaction with a floating flexible WEC.
4. Results and discussions
The validated model is applied to a series of case studies to investigate the influence of multiparameters, e.g. incident wave direction, wave frequency and the position of the PTO system, on wave power absorption of the device.
4.1. Effect of incident wave direction $\beta$
Figure 6 presents the contour of the wave power capture factor (
$\eta _e$) as a function of incident wave direction (
$\beta$) and PTO damping coefficient (
$\bar {c}$) for 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$kh=4.0$. The panels represent the device with a discrete PTO system consisting of a different number of PTO units and the device with a continuum PTO system. Since the PTO system is distributed symmetrically about the 
$x$ axis, the examined 
$\beta$ are limited to the range 
$[0,{\rm \pi} ]$. For the discrete PTO system consisting of 
$N$ PTO units, the property of the device is symmetric about 
$\theta =\beta =0,{\rm \pi} /N, 2{\rm \pi} /N, \ldots\,$, hence the results for 
$\beta \in [0,{\rm \pi} /N]$ can be used to generate the results for the remaining range of 
$\beta$ by placing a mirror on each of the lines of symmetry, which is also observed in figure 6.
Figure 6. Contour plots for the variation of 
$\eta _e$ as a function of incident wave direction 
$\beta$ and PTO damping coefficient 
$\bar {c}$, for 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$kh=4.0$: (a) 
$N=1$, (b) 
$N=2$, (c) 
$N=3$, (d) 
$N=4$, (e) 
$N=5$, and (f) 
$N=\infty$.
For any specified value of 
$\beta$, there is an optimal 
$\bar {c}$ to maximise wave power absorption of the device, and the optimised value of 
$\bar {c}$ is insensitive to 
$\beta$; for any specified value of 
$\bar {c}$, there are one or two optimised 
$\beta$ in the interval range 
$\beta \in [0,{\rm \pi} /N]$ such that the wave power can be captured most efficiently. For 
$N=\infty$, the performance of the device is independent of 
$\beta$.
The peak values of 
$\eta _e$ and the peak positions in the interval range 
$\beta \in [0,{\rm \pi} /N]$ for different 
$N$ are listed in table 1. In each case, except 
$N=1$, for which there are two peaks of 
$\eta _e$ observed, there is merely one peak for 
$\beta \in [0,{\rm \pi} /N]$. As 
$N$ increases from 1 to 
$\infty$, the area of the field of 
$\eta _e>1.0$ plotted in figure 6, including the peak value of 
$\eta _e$, presents an overall increase trend. The 
$\eta _e$ value for 
$N=4$ is the most sensitive one to the change of 
$\beta$ among the six examined cases, making the peak value of 
$\eta _e$ even larger than that for 
$N=5$.
Table 1. Peak values of 
$\eta _e$ and the peak positions in the contours plotted in figure 6.
4.2. Effect of wavenumber $k$
The performance of the device with different values of 
$\bar {c}$ for different wave frequencies in terms of 
$kh$ ranging from 0.05 to 10.0 and 
$\beta ={\rm \pi} /6$ is also examined, the result of which is plotted in figure 7. In the computed range of wave conditions and PTO damping coefficients, the peak value of wave power capture factor and the corresponding wavenumber and PTO damping coefficient in terms of 
$(\eta _e, kh, \bar {c})$ for 
$N=1,2, \ldots, 5, \infty$ are 
$(1.239, 4.28, 0.04)$, 
$(1.826, 3.91, 0.06)$, 
$(3.114, 7.64, 0.08)$, 
$(3.677, 5.51, 0.20)$, 
$(3.695, 4.65, 0.12)$ and 
$(5.397, 5.03, 0.24)$, respectively. The device deflection and the near-field wave response for those optimised circumstances are plotted in figure 8.
Figure 7. Contour plots for the variation of 
$\eta _e$ as a function of wavenumber 
$kh$ and PTO damping coefficient 
$\bar {c}$, for 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta ={\rm \pi} /6$: (a) 
$N=1$, (b) 
$N=2$, (c) 
$N=3$, (d) 
$N=4$, (e) 
$N=5$, and (f) 
$N=\infty$.
Figure 8. Deflection of the floating elastic disk, 
${\rm Re}(\eta \,\mathrm {e}^{-\mathrm {i} \omega t})/A$, and the near-field wave motion, 
${\rm Re}(\eta _w\,\mathrm {e}^{-\mathrm {i} \omega t})/A$, at 
$t=0$ for 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta ={\rm \pi} /6$: (a) 
$N=1$, 
$(kh, \bar {c})=(4.28, 0.04)$; (b) 
$N=2$, 
$(kh, \bar {c})=(3.91, 0.06)$; (c) 
$N=3$, 
$(kh, \bar {c})=(7.64, 0.08)$; (d) 
$N=4$, 
$(kh, \bar {c})=(5.51, 0.20)$; (e) 
$N=5$, 
$(kh, \bar {c})=(4.65, 0.12)$; and (f) 
$N=\infty$, 
$(kh, \bar {c})=(5.03, 0.24)$.
As 
$N$ varies from 1 to 5, apart from the peak of the 
$\eta _e$ contour, the overall contour shape changes as well, which, to some extent, may be explained with the flexural wave interaction induced by the relative motion between the PTO units. Due to the water wave interaction with the floating elastic disk, a train of propagating flexural waves, which correspond to the root 
$\kappa _0$, are excited on the elastic disk by the incident waves. The wavelength of the propagating flexural waves (
$\lambda _0$) may be evaluated by 
$\lambda _0=2{\rm \pi} /\kappa _0$. When the distance between two PTO units is an integer multiple of 
$\lambda _0/2$, a local resonance of the flexural waves between the two PTO units may occur, stimulating the response of the PTO units and providing benefits to wave power absorption. The wavenumber of the propagating flexural waves can be evaluated with 
$\kappa _0=2{\rm \pi} /\lambda _0$, and the corresponding non-dimensionalised wavenumber 
$kh$ can be determined by using the dispersion relation for the interior region, i.e. (2.15). Table 2 reports the wavenumbers in the computed range of wave conditions associated with the local resonant flexural waves between the PTO units for the devices with 
$N=2$, 3, 4 and 5. The values of 
$\lambda _0$ smaller than 
$1.7h$ are not included in the table because the corresponding 
$kh$ are larger than 10.0, which is out of the computed range of wave frequencies. The predicted wavenumbers are found to agree well with the horizontal positions, i.e. 
$kh$, where the ‘crests’ of 
$\eta _e$ generally happen (see figure 7), indicating that the propagating component of the flexural waves travelling on the disk plays a significant role on the behaviour of the device. When a continuum PTO system is adopted, the device is found to capture wave power efficiently in a large range of wave frequencies.
Table 2. Evaluation of wavenumbers associated with the local resonant flexural waves between two PTO units.
4.3. Effect of PTO position $r_0$
The effect of the PTO position 
$r_0$ is examined for 
$R/h=2.0$, 
$\beta ={\rm \pi} /6$, and 
$\bar {c}$ varying from 0.01 to 0.8. The results in terms of the contours of 
$\eta _e$ as a function of 
$r_0/R$ and 
$\bar {c}$ for the device consisting of different number of PTO units are plotted in figure 9. For the device with any specified 
$N$, the optimised PTO damping coefficient generally presents an overall declining trend as 
$r_0/R$ increases from 0.05 to 0.95. For 
$N=1$, there is only one peak of 
$\eta _e$ observed at 
$r_0/R=0.44$. As 
$N$ increases, two or more peaks of the 
$\eta _e$ contour are excited in the computed range of 
$r_0/R$, and the contour shape changes accordingly. This may also be explained by the local resonant flexural waves between the PTO units. The non-dimensionalised wavenumber and wavelength of the propagating flexural waves travelling on the disk for 
$kh=4.0$ are 
$\kappa _0h=2.70$ and 
$\lambda _0=2.33h$, respectively. As discussed in § 4.2, a local resonance of the flexural waves between two PTO units may happen when the distance between them is 
$1.16h$, 
$2.33h$ or 
$3.49h$, satisfying 
$j\lambda _0/2$ (
$j=1, 2, 3, \ldots$). Table 3 lists the values of 
$r_0/R$ for which the local resonance of the propagating flexural waves is likely to happen. For some examined cases, e.g. 
$N=3$, the values of 
$r_0/R$ reported in table 3 and the corresponding results plotted in figure 9 agree well with one another. However, for some of the other circumstances, a deviation between them is observed, e.g. 
$N=2$, for which the predicted 
$r_0/R$ are slightly larger than those observed in figure 9. Despite this, the data reported in table 3 may help to bring some insight into the performance of the floating flexible device. The peak value of 
$\eta _e$ is found to increase with the increase of 
$N$. For 
$N=\infty$, the peak of 
$\eta _e$ occurs at 
$(r_0/R, \bar {c})=(0.80,0.14)$, and the peak value can be as large as 8.90.
Table 3. Evaluation of PTO positions satisfying the local resonant condition of the propagating flexural waves.
Figure 9. Contour plots for the variation of 
$\eta _e$ as a function of PTO position 
$r_0/R$ and PTO damping coefficient 
$\bar {c}$, for 
$R/h=2.0$, 
$\beta ={\rm \pi} /6$ and 
$kh=4.0$: (a) 
$N=1$; (b) 
$N=2$, (c) 
$N=3$, (d) 
$N=4$, (e) 
$N=5$, and (f) 
$N=\infty$.
For a traditional single WEC consisting of an axisymmetric rigid heaving body, the theoretical maximum wave power capture factor is 
$\eta _e=1.0$. If the rigid body is also allowed to oscillate in pitch/surge mode, i.e. the axisymmetric rigid body captures wave power in both heave and surge/pitch, then the theoretical limit of 
$\eta _e$ increases further to 3.0. The present floating elastic WEC is demonstrated to be able to extend the limit for a large range of circumstances (e.g. see figures 6, 7 and 9).
4.4. Optimisation of the angular-independent continuum PTO system
As reported in the previous subsections, wave power absorption of the device is affected significantly by the PTO damping coefficient. The maximum wave power capture factor and the corresponding optimised PTO damping coefficient (
$\bar {c}_{{\it opt}}$) can be determined by using a trial-and-error method.
Instead of applying the trial-and-error method, in this subsection, we implement the optimisation outlined in § 2.6, which provides a recipe for equal damper settings of the angular-independent continuum PTO system to extract the available power from any given circular mode component, 
$m$ in the incident wave as much as possible. Results are illustrated in figure 10, where figure 10(b) shows the variation of the value of 
$\bar {c}_{{\it opt}}^{(m)}$ with frequency, resulting in the capture factor in figure 10(a). Although the PTO damping coefficient is designed to capture the wave power contained in a specific circular mode component, the capture factor can be much larger than the available power in this circular mode, since power is absorbed from all the circular wave components in the incident wave other than the one being targeted. In the computed range of 
$kh$, the envelope of the 
$\eta _e\unicode{x2013}kh$ curves with 
$\bar {c}=\bar {c}_{{\it opt}}^{(m)}$ for 
$m=0, 1, 2, 3$ is found to agree well with that for 
$\bar {c}=\bar {c}_{{\it opt}}$.
Figure 10. (a) Capture factor against dimensionless wavenumber for 
$R/h=2.0$, 
$r_0/R=0.5$, 
$\beta =0$ and 
${\rm Im}(\bar {c})=0$, with (b) corresponding damper values optimised in order to capture the available power in the 
$m$th circular mode component as much as possible.
Figure 11 presents the contribution of the capture factor from different circular wave components (
$n$, different curves) when PTO dampers have been tuned to extract the available power in a particular mode 
$m$ as much as possible. There is significant wave power absorption across multiple circular modes, indicating and explaining how the elastic disk-shaped WEC can capture more wave power compared to that of the rigid axisymmetric WECs.
Figure 11. The partition of capture factor into contributions from different circular mode components for operation tuned to be optimal for the 
$m$th mode, with 
$R/h=2.0$, 
$r_0/R=0.5$, 
$\beta =0$ and 
${\rm Im}(\bar {c})=0$: (a) 
$m=0$, (b) 
$m=1$, (c) 
$m=2$, (d) 
$m=3$, (e) 
$m=4$, and (f) 
$m=5$.
It is noted from figure 10(a) that there is a trough of the 
$\eta _e\unicode{x2013}kh$ curve occurring at 
$kh\approx 3$. This can be explained from the point of view of natural vibration mode shapes. At 
$kh$ approaches 3.1, for the dominant axisymmetric natural vibration mode of the disk excited by the incident waves, the angular nodal line could be rather close to the circle 
$r=r_0$, where the PTO units are located, resulting in a trough of the wave power absorption related to that mode (i.e. see 
$\eta _0$ in figure 11) around 
$kh=3.1$ regardless of the value of 
$\bar {c}$. This in turn affects the frequency response of 
$\eta _e$ (see figure 10a), leading to a trough around 
$kh=3.1$, except in the case with 
$\bar {c}=\bar {c}^{(0)}_{{\it opt}}$, for which the trough occurs at 
$kh=2.9$. This is because the damping coefficient to maximise the power captured wave power from the zeroth mode at 
$kh=2.9$ is too large (see figure 10b), and it makes the response of the high-order modes repressed, which is reflected from the wave power absorption related to those modes as shown in figure 11(a), and ultimately results in a sudden drop of the overall wave power capture factor 
$\eta _e$.
Apart from the pure damping continuum PTO system, the one consisting of both damper and stiffness/mass is also considered. Figure 12 illustrates the frequency response of 
$\eta _e$ and the corresponding optimised setting of the continuum PTO system in terms of PTO damping (
${\rm Re}(\bar {c}_{OPT}^{(m)})$) and PTO stiffness/mass (
${\rm Im}(\bar {c}_{OPT}^{(m)})$) to extract all the available power from any given circular mode component 
$m$ in the incident wave. The contribution of the capture factor from different circular wave components (
$n$, different curves) is plotted in figure 13. As expected, 
$\eta _n=\epsilon _n$ is achieved when 
$n=m$ (see figure 13). Although the device is able to capture all the available power from the 
$m$th circular mode component, its performance in capturing wave power from some of the other circular mode components is low, resulting in a lower envelope of the 
$\eta _e\unicode{x2013}kh$ curves (figure 12a) compared to that for the pure damper PTO system (figure 10a) at 
$kh>6.0$.
Figure 12. (a) Capture factor against dimensionless wavenumber for 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta =0$, with (b,c) corresponding damper and stiffness/mass values optimised in order to capture all the available power in the 
$m$th circular mode.
Figure 13. The partition of capture factor into contributions from different circular modes for operation tuned to be optimal for the 
$m$th mode, with 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta =0$: (a) 
$m=0$, (b) 
$m=1$, (c) 
$m=2$, (d) 
$m=3$, (e) 
$m=4$, and (f) 
$m=5$.
The theoretical optimisation method of the PTO damping coefficient is applied to examine the scalability of the device in terms of 
$R/h$ with 
$N=\infty$, 
${\rm Im}(\bar {c})=0$ and 
$r_0/R=0.5$ on the wave power extraction. Figure 14 illustrates the frequency response of the maximum wave power capture factor 
$\eta _{{\it max}}$ for 
$R/h=2.0$, 2.5, 3.0, 3.5 and 4.0. For each case, the 
$\eta _{{\it max}}\unicode{x2013}kh$ curve is represented by the corresponding envelope of the 
$\eta _e\unicode{x2013}kh$ curves with 
$\bar {c}=\bar {c}_{{\it opt}}^{(m)}$ for 
$m=0,1, \ldots, 5$. It is observed that in the computed range of wave conditions, more power can generally be captured from water waves for a larger diameter of the device.
Figure 14. Frequency response of 
$\eta _{{\it max}}$ for different values of 
$R/h$, with 
$N=\infty$, 
$r_0/R=0.5$ and 
${\rm Im}(\bar {c})=0$.
4.5. Non-uniformly distributed discrete PTO system
The results presented in the previous subsections all concern cases with the PTO units uniformly distributed over a circle below the disk. The model proposed in this paper can also be applied to cases with the PTO units uniformly distributed within a sector of the circle. In this subsection, the device consisting of 
$N=5$ units is selected as a case to demonstrate the effect of the non-uniform distribution on wave power absorption. Six configurations, as illustrated in figure 15, are considered. The five PTO units in the configurations in figures 15(a,d) are both distributed on a quarter of a circle, in which the PTO units for figure 15(a) are located at the sea-side of the disk whereas those for figure 15 (d) are positioned at the lee-side. Similar arrangements are used in the configurations in figures 15(b,e), except that the PTO units are deployed on a half-circle. The configurations in figures 15(c,f) are the control cases, with the five PTO units uniformly distributed around the full circle, in which 
$\theta _1={\rm \pi}$ and 0, respectively. Each device is subjected to incident waves propagating from the left-hand side to the right-hand side (
$\beta =0$). Figure 16 plots the variation of the wave power capture width ratio with 
$\bar {c}$ and 
$kh$ for these six configurations correspondingly.
Figure 15. Different distributions of the case with 
$N=5$ for 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta =0$: (a) PTO units distributed over 
$\theta \in [0.75{\rm \pi},1.25{\rm \pi} ]$; (b) PTO units distributed over 
$\theta \in [0.5{\rm \pi},1.5{\rm \pi} ]$; (c) PTO units uniformly distributed all over the circle with 
$\theta _1={\rm \pi}$; (d) PTO units distributed over 
$\theta \in [-0.25{\rm \pi},0.25{\rm \pi} ]$; (e) PTO units distributed over 
$\theta \in [-0.5{\rm \pi},0.5{\rm \pi} ]$; (f) PTO units distributed all over the circle with 
$\theta _1=0$.
Figure 16. Contour plots for the variation of 
$\eta _e$ as a function of wavenumber 
$kh$ and PTO damping coefficient 
$\bar {c}$ for the six configurations of the PTO units as shown in figure 15, i.e. 
$N=5$, 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta =0$: (a) PTO units distributed over 
$\theta \in [0.75{\rm \pi},1.25{\rm \pi} ]$; (b) PTO units distributed over 
$\theta \in [0.5{\rm \pi},1.5{\rm \pi} ]$; (c) PTO units uniformly distributed all over the circle with 
$\theta _1={\rm \pi}$; (d) PTO units distributed over 
$\theta \in [-0.25{\rm \pi},0.25{\rm \pi} ]$; (e) PTO units distributed over 
$\theta \in [-0.5{\rm \pi},0.5{\rm \pi} ]$; (f) PTO units uniformly distributed all over the circle with 
$\theta _1=0$.
For the PTO units clustered at the sea-side of the device, the device presents a better performance at 
$kh\in [7.0,10.0]$ in terms of a more significant wave power capture factor and less sensitivity to the change of 
$kh$ compared to those with the PTO units located at the lee-side. For the PTO units clustered at the lee-side of the device, a dramatic peak of 
$\eta _e$ is observed at 
$kh\in [4.0,6.0]$. The optimised PTO damping coefficient for the device with the PTO units clustered at the lee-side is larger than that for the case with the PTO units distributed at the sea-side.
The control cases (figures 16c,f) perform rather differently from one another at 
$kh\in [6.0,10.0]$, where the configuration in figure 15(c) absorbs the least wave power, whereas that in figure 15(f) captures the most, among the six examined cases, demonstrating the significant effect of PTO units distribution (or incident wave direction) on wave power absorption of the device. It should be noted that for wave conditions around 
$kh=9.0$, the maximum value of 
$\eta _e$ for the configuration in figure 15(f) is even larger than that for 
$N=\infty$ (see figure 7f), meaning that for some specific circumstances, the discrete PTO system could be more advantageous than the continuum one in harnessing wave power.
4.6. Angular-dependent continuum PTO system
In this section, we take six cases of the angular-dependent continuum PTO system, the angular PTO damping distribution of which is plotted in figure 17, as an example to examine the performance of the device in capturing wave power with different angular-dependent continuum PTO systems. The incident waves propagate along the 
$Ox$ axis (i.e. 
$\beta =0$). The PTO system for each case is located at 
$r_0/R=0.5$, and the distribution of the PTO damping is symmetric about the 
$Ox$ axis. For cases 1, 2 and 3, the PTO damping at the lee-side is smaller than that at the wind-side, and as 
$\theta$ increases from 0 to 
${\rm \pi}$, 
$c(\theta )$ increases monotonically from 0 to the maximum, 
$c_{{\it max}}$, following linear, trigonometric and quadratic functions, respectively (see figure 17). When the PTO system is reflected over the 
$Oy$ axis, cases 1, 2 and 3 turn into cases 4, 5 and 6, respectively.
Figure 17. Angular distribution of the PTO damping coefficient for six angular-dependent continuum PTO systems.
Figure 18 illustrates the variation of 
$\eta _e$ as a function of wavenumber and PTO damping coefficient in terms of 
$\bar {c}^*$, which denotes the non-dimensional value of 
$c_{{\it max}}$ (
$\bar {c}^*=c_{{\it max}}/(\rho R \sqrt {gh})$). When the lee-side PTO damping is smaller than that at the wind-side (figures 18a–c), there are two peaks of 
$\eta _e$ in the computed range of 
$kh$ and 
$\bar {c}^*$. The main peak occurs around 
$kh=7.5$, and the other one around 
$kh=4.0$, which is nearly a half as large as that of the main peak. The corresponding optimised PTO damping in terms of 
$\bar {c}^*$ for cases 2 and 3 are smaller and larger, respectively, than that of case 1. This may be explained with the angular distribution shape of 
$c(\theta )$: for any specified 
$\bar {c}^*$, although the 
$c(\theta )$ values are the same for these three cases at 
$\theta =0$ and 
$\pm {\rm \pi}$, 
$c(\theta )$ values for cases 2 and 3 are the largest and smallest, respectively, among them for any other angle. To have a same equivalent value of the full-angle PTO damping, smaller and larger values of 
$c_{max}$ for cases 2 and 3 would be expected, respectively. A similar phenomenon is observed for cases 4, 5 and 6 (see figures 18d–f) when the lee-side PTO damping is larger than that at the wind-side. The wavenumbers where the peaks of 
$\eta _e$ occur shift towards smaller wave frequencies, whereas the optimised 
$\bar {c}^*$ become larger. Compared with the angular-independent continuum PTO case (see figure 7f), wave power absorption of the device in terms of both the peak value and bandwidth can be improved by applying an angular-dependent continuum PTO with larger PTO damping placed at the lee-side (figures 18d–f). On the contrary, wave power absorption can be inhibited should the larger PTO damping be placed at the wind-side.
Figure 18. Contour plots for the variation of 
$\eta _e$ as a function of wavenumber 
$kh$ and PTO damping coefficient 
$\bar {c}^*$ for six angular-dependent continuum PTO systems as shown in figure 17, with 
$N=\infty$, 
$R/h=2.0$, 
$r_0/R=0.5$ and 
$\beta =0$: (a) case 1, i.e. 
$c(\theta )=c_{max}|\theta |/{\rm \pi}$; (b) case 2, i.e. 
$c(\theta )=c_{max}\sin (|\theta |/2)$; (c) case 3, i.e. 
$c(\theta )=c_{max}(\theta /{\rm \pi} )^2$; (d) case 4, i.e. 
$c(\theta )=c_{max}(1-|\theta |/{\rm \pi} )$; (e) case 5, i.e. 
$c(\theta )=c_{max}\cos (\theta /2)$; (f) case 6, i.e. 
$c(\theta )=c_{max}(1-|\theta |/{\rm \pi} )^2$.
5. Conclusions
A concept of a floating elastic disk-based WEC is proposed in this paper. The device captures wave power by driving a series of PTO units deployed on a circle below the disk. To evaluate the wave power extraction from the device, a theoretical model is developed based on the eigenfunction matching method within the framework of the linear potential flow theory. As well as the discrete PTO system, a continuum PTO system is considered in the model to represent the case when the discrete PTO system is composed of a large number of PTO units. The wave power absorption is evaluated using two approaches: one is the straightforward method by accumulating/integrating the power absorbed by each PTO unit; the other is the indirect method expressing the wave power by the difference between the far-field incoming wave power and outgoing wave power. The indirect method was shown to produce the same wave power absorption as the direct method. A multiparameter impact analysis was carried out by applying the theoretical model.
The main findings are as follows.
(i) It is confirmed that the continuum PTO system converges as
$N$ of the discrete PTO system increases. In the examined cases, the results of the discrete PTO system with $N=16$ are found to be commensurate with those of the continuum PTO system.(ii) For any specified value of the incident wave direction
$\beta$, there is an optimal PTO damping coefficient $\bar {c}$ to maximise wave power absorption, and the optimised value of $\bar {c}$ is insensitive to $\beta$.(iii) A train of propagating flexural waves is excited on the floating elastic disk. When the distance between two PTO units is an integer multiple of
$\lambda _0/2$ (where $\lambda _0$ denotes the wavelength of the propagating flexural waves), a local resonance of the flexural waves between the two PTO units may occur, stimulating the response of the PTO units and providing benefits to wave power absorption.(iv) The device adopting a continuum PTO system is found to capture wave power efficiently in a wide range of wave frequencies. More power can generally be captured from water waves for a larger diameter of the device.
(v) There is an optimised PTO position in terms of
$r_0/R$ to maximise wave power absorption of the device. The optimised $r_0/R$ changes with the variation of the number ($N$) of the PTO units included in the PTO system. The maximised wave power absorption increases with the increase of $N$.(vi) The present floating elastic WEC is demonstrated to achieve well in excess of the standard limit of a capture factor of
$\eta _e=3$ for an axisymmetric rigid body. The wave power capture factors above $\eta _e=8.5$ are reported in computations in this paper.(vii) It is theoretically possible to adopt optimised PTO damper and stiffness/mass for the continuum PTO system to guarantee the absorption of 100 % of the energy flux available in one circular component of the plane incident wave.
(viii) Wave power absorption of the device in terms of both the peak value and bandwidth can be improved by applying an angular-dependent continuum PTO with larger PTO damping placed at the lee-side.
The present work applied the direct approach to solving the combined diffraction and radiation potential. An alternative model could be based on the standard dry-mode approach, with which the response of the plate can be expanded into heave mode, pitch mode, and a series of flexible modes, and each mode is treated as a generalised mode in solving the wave diffraction and wave radiation problems, respectively (e.g. see Michele et al. Reference Michele, Zheng and Greaves2022). A comparison between the performance of the present model and the standard dry-mode approach in terms of the accuracy and the cost of computing is an aspect of interest. This evaluation is beyond the scope of this paper and is therefore left for future work.
Funding
S.Z. and S.M. gratefully acknowledge support from the Supergen ORE Hub ECR Research Fund project ECRRF2021-12 – ‘Analytical and experimental modelling of a floating/submerged elastic disk’. S.Z. also acknowledges the State Key Laboratory of Hydroscience and Engineering (Tsinghua University) for supporting part of this work through the Open Research Fund Program (grant no. sklhse-2021-E-02). D.G. gratefully acknowledges the EPSRC for supporting part of this work through the Supergen ORE Hub, EP/S000747/1. H.L. was supported by A*STAR under its RIE 2020 Industry Alignment Fund, Grant no. A19F1a0104.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation process of the formulae and calculation for the unknown coefficients
After inserting the expressions of the velocity potentials at Regions 1, 2 and 3 into (2.23)–(2.29), we have
\begin{gather} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty [(A_{m,l}-B_{m,l})\,{\rm J}_m(\kappa_lr_0)-C_{m,l}\, H_m(\kappa_lr_0)]\,Y_l(z)\,\mathrm{e}^{\mathrm{i} m\theta}=0, \end{gather}
\begin{gather} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \kappa_l[(A_{m,l}-B_{m,l})\,{\rm J}'_m(\kappa_lr_0)-C_{m,l}\,H'_m(\kappa_lr_0)]\,Y_l(z)\,\mathrm{e}^{\mathrm{i} m\theta}=0, \end{gather}
\begin{gather} \hspace{-13pc}\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty [B_{m,l}\,{\rm J}_m(\kappa_lR)+C_{m,l}\, H_m(\kappa_lR)]\,Y_l(z)\,\mathrm{e}^{\mathrm{i} m\theta}\nonumber\\ \quad=\sum_{m={-}\infty}^\infty\sum_{l=0}^\infty D_{m,l}\,H_m(k_lR)\,Z_l(z)\,\mathrm{e}^{\mathrm{i} m\theta}-\frac{\mathrm{i} gA}{\omega}\,Z_0(z) \sum_{m={-}\infty}^\infty \mathrm{i}^m\,\mathrm{e}^{-\mathrm{i} m \beta}\,{\rm J}_m(kR)\,\mathrm{e}^{\mathrm{i} m\theta}, \end{gather}
\begin{gather} \hspace{-11.5pc}\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \kappa_l[B_{m,l}\,{\rm J}'_m(\kappa_lR)+C_{m,l}\, H'_m(\kappa_lR)]\,Y_l(z)\,\mathrm{e}^{\mathrm{i} m\theta}\nonumber\\ \hspace{2pc} =\sum_{m={-}\infty}^\infty\sum_{l=0}^\infty D_{m,l}\,k_l\,H'_m(k_lR)\,Z_l(z)\, \mathrm{e}^{\mathrm{i} m\theta}-\frac{\mathrm{i} gkA}{\omega}\,Z_0(z) \sum_{m={-}\infty}^\infty \mathrm{i}^m\, \mathrm{e}^{-\mathrm{i} m \beta}\,{\rm J}'_m(kR)\,\mathrm{e}^{\mathrm{i} m\theta}, \end{gather}
\begin{gather} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{(A_{m,l}-B_{m,l})\,{\rm J}_m(\kappa_lr_0)-C_{m,l}\,H_m(\kappa_lr_0)}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}=0, \end{gather}
\begin{gather} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{\kappa_l[(A_{m,l}-B_{m,l})\,{\rm J}'_m(\kappa_lr_0)-C_{m,l}\,H'_m(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}=0, \end{gather}
\begin{gather} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{\kappa_l^2[(A_{m,l}-B_{m,l})\,{\rm J}''_m(\kappa_lr_0)-C_{m,l}\,H''_m(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}=0. \end{gather} After inserting (2.21) into (2.30) and (2.31), producing the radial delta function, and taking the continuous 
$\eta$, 
$\eta '$ and 
$\eta ''$ (i.e. (A5)–(A7)) into consideration, a jump in 
$\eta '''$ can be obtained, which gives
\begin{equation} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{\kappa_l^3[(A_{m,l}-B_{m,l})\,{\rm J}'''_m(\kappa_lr_0)-C_{m,l}\,H'''_m(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}=\frac{-1}{\mathrm{i} \omega \rho r_0 \chi} \sum_{n=1}^N F_n\,\delta(\theta-\theta_n) \end{equation}and
\begin{equation} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{\kappa_l^3[(A_{m,l}-B_{m,l})\,{\rm J}'''_m(\kappa_lr_0)-C_{m,l}\,H'''_m(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}=\frac{-1}{\mathrm{i} \omega \rho \chi}\sum_{m={-}\infty}^{\infty}f_m\,\mathrm{e}^{\mathrm{i} m\theta} \end{equation}for the discrete PTO and continuum PTO, respectively.
Using the fact that 
$\varDelta ({\rm J}_m(\kappa _l r)\,\mathrm {e}^{\mathrm {i} m \theta })=-\kappa _l^2\,{\rm J}_m(\kappa _l r)\,\mathrm {e}^{\mathrm {i} m \theta }$ and 
$\varDelta (H_m(\kappa _l r)\,\mathrm {e}^{\mathrm {i} m \theta })= $ 
$-\kappa _l^2\, H_m(\kappa _l r)\,\mathrm {e}^{\mathrm {i} m \theta }$, the free-floating edge conditions (2.12) and (2.13) can be rewritten as
$$\begin{gather} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{1}{\chi \kappa_l^4+1-K\gamma} \left\{B_{m,l}\left[-\kappa_l^2\,{\rm J}_m(\kappa_l R)-\frac{1-\nu}{R}\left(\kappa_l\,{\rm J}'_m(\kappa_l R)-\frac{m^2}{R}\,{\rm J}_m(\kappa_l R)\right)\right]\right.\nonumber\\ \left. {}+C_{m,l}\left[-\kappa_l^2\,H_m(\kappa_l R)-\frac{1-\nu}{R}\left(\kappa_l\,H'_m(\kappa_l R)-\frac{m^2}{R}\,H_m(\kappa_l R)\right)\right]\right\} \mathrm{e}^{\mathrm{i} m\theta}=0, \end{gather}$$
$$\begin{gather} \sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{1}{\chi \kappa_l^4+1-K\gamma} \left\{B_{m,l}\left[-\kappa_l^3\,{\rm J}'_m(\kappa_l R)+\frac{m^2(1-\nu)}{R^2}\left(-\kappa_l\,{\rm J}'_m(\kappa_l R)+\frac{1}{R}\,{\rm J}_m(\kappa_l R)\right)\right]\right.\nonumber\\ \left.{}+C_{m,l}\left[-\kappa_l^3\,H'_m(\kappa_l R)+\frac{m^2(1-\nu)}{R^2}\left(-\kappa_l\,H'_m(\kappa_l R)+\frac{1}{R}\,H_m(\kappa_l R)\right)\right]\right\} \mathrm{e}^{\mathrm{i} m\theta}=0, \end{gather}$$respectively.
After multiplying both sides of (A1)–(A4) by 
$Z_\zeta (z)\, \mathrm {e}^{-\mathrm {i}\tau \theta }$, integrating in 
$z\in [-h, 0]$ and 
$\theta \in [0, 2{\rm \pi} ]$, and using their orthogonality characteristics, (A1)–(A4) can be rewritten as
\begin{gather} \sum_{l={-}2}^\infty [(A_{\tau,l}-B_{\tau,l})\,{\rm J}_\tau(\kappa_lr_0)-C_{\tau,l}\, H_\tau(\kappa_lr_0)]Y_{l,\zeta}=0, \end{gather}
\begin{gather} \sum_{l={-}2}^\infty \kappa_l[(A_{\tau,l}-B_{\tau,l})\,{\rm J}'_\tau(\kappa_lr_0)-C_{\tau,l}\, H'_\tau(\kappa_lr_0)]\,Y_{l,\zeta}=0, \end{gather}
\begin{gather} \sum_{l={-}2}^\infty [B_{\tau,l}\,{\rm J}_\tau(\kappa_lR)+C_{\tau,l}\, H_\tau(\kappa_lR)]\,Y_{l,\zeta}- D_{\tau,\zeta}\,H_\tau(k_{\zeta}R)\,A_{\zeta}={-}\delta_{\zeta,0}\,\frac{\mathrm{i} gA}{\omega}\,A_0\,\mathrm{i}^{\tau}\,\mathrm{e}^{-\mathrm{i} \tau \beta}\,{\rm J}_\tau(kR), \end{gather}
\begin{gather} \sum_{l={-}2}^\infty \kappa_l[B_{\tau,l}\,{\rm J}'_{\tau}(\kappa_lR)+C_{\tau,l}\, H'_{\tau}(\kappa_lR)]\,Y_{l,\zeta}- D_{\tau,\zeta}\,k_{\zeta}\,H'_{\tau}(k_{\zeta}R)\,A_{\zeta}\notag\\ \quad={-}\delta_{\zeta,0}\,\frac{\mathrm{i} gkA}{\omega}\,A_0\,\mathrm{i}^{\tau}\,\mathrm{e}^{-\mathrm{i} \tau \beta}\,{\rm J}'_{\tau}(kR), \end{gather}where
\begin{gather} A_{\zeta}=\int_{{-}h}^0Z_{\zeta}^2(z)\,\mathrm dz=\frac{\sinh(k_{\zeta}h)\cosh(k_{\zeta}h)+k_{\zeta}h}{2 k_{\zeta}\cosh^2(k_{\zeta}h)}, \end{gather}
\begin{gather} Y_{l,\zeta}=\int_{{-}h}^0Y_l(z)\,Z_\zeta(z)\,\mathrm dz=\frac{\kappa_l\sinh(\kappa_lh)\cosh(k_\zeta h)-k_\zeta\cosh(\kappa_l h)\sinh(k_\zeta h)}{(\kappa_l^2-k_\zeta^2)\cosh(\kappa_lh)\cosh(k_\zeta h)}, \end{gather}
and 
$\delta _{\zeta,0}$ denotes the Kronecker delta function, which is equal to 1 if 
$\zeta =0$ and 0 otherwise.
After multiplying both sides of (A5)–(A11) by 
$\mathrm {e}^{-\mathrm {i} \tau \theta }$, integrating in 
$\theta \in [0, 2{\rm \pi} ]$ and using their orthogonality characteristics, (A5)–(A11) can be rewritten as
\begin{gather} \sum_{l={-}2}^\infty \frac{(A_{\tau,l}-B_{\tau,l})\,{\rm J}_{\tau}(\kappa_lr_0)-C_{\tau,l}\,H_{\tau}(\kappa_lr_0)}{\chi \kappa_l^4+1-K\gamma} =0, \end{gather}
\begin{gather} \sum_{l={-}2}^\infty \frac{\kappa_l[(A_{\tau,l}-B_{\tau,l})\,{\rm J}'_{\tau}(\kappa_lr_0)-C_{\tau,l}\,H'_{\tau}(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma} =0, \end{gather}
\begin{gather} \sum_{l={-}2}^\infty \frac{\kappa_l^2[(A_{\tau,l}-B_{\tau,l})\,{\rm J}''_{\tau}(\kappa_lr_0)-C_{\tau,l}\,H''_{\tau}(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma} =0, \end{gather}
\begin{gather} \sum_{l={-}2}^\infty \frac{\kappa_l^3[(A_{\tau,l}-B_{\tau,l})\,{\rm J}'''_{\tau}(\kappa_lr_0)-C_{\tau,l}\,H'''_{\tau}(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma} +\frac{1}{2{\rm \pi} \mathrm{i} \omega \rho r_0 \chi}\sum_{n=1}^N F_n\,\mathrm{e}^{-\mathrm{i} \tau \theta_n}=0 \end{gather}and
\begin{equation} \sum_{l={-}2}^\infty \frac{\kappa_l^3[(A_{\tau,l}-B_{\tau,l})\,{\rm J}'''_{\tau}(\kappa_lr_0)-C_{\tau,l}\,H'''_{\tau}(\kappa_lr_0)]}{\chi \kappa_l^4+1-K\gamma} +\frac{f_{\tau}}{\mathrm{i} \omega \rho \chi}=0, \end{equation}for the discrete and continuum PTO units, respectively. Also,
$$\begin{gather} \sum_{l={-}2}^\infty \frac{1}{\chi \kappa_l^4+1-K\gamma} \left\{B_{\tau,l}\left[-\kappa_l^2\,{\rm J}_{\tau}(\kappa_l R)-\frac{1-\nu}{R}\left(\kappa_l\,{\rm J}'_{\tau}(\kappa_l R)-\frac{{\tau}^2}{R}\,{\rm J}_{\tau}(\kappa_l R)\right)\right]\right.\nonumber\\ \left.{}+C_{\tau,l}\left[-\kappa_l^2\,H_{\tau}(\kappa_l R)-\frac{1-\nu}{R}\left(\kappa_l\,H'_{\tau}(\kappa_l R)-\frac{{\tau}^2}{R}\,H_{\tau}(\kappa_l R)\right)\right]\right\}=0, \end{gather}$$
$$\begin{gather} \sum_{l={-}2}^\infty \frac{1}{\chi \kappa_l^4+1-K\gamma} \left\{B_{\tau,l}\left[-\kappa_l^3\,{\rm J}'_{\tau}(\kappa_l R)+\frac{{\tau}^2(1-\nu)}{R^2}\left(-\kappa_l\,{\rm J}'_{\tau}(\kappa_l R)+\frac{1}{R}{\rm J}_{\tau}(\kappa_l R)\right)\right]\right.\nonumber\\ \left.{}+C_{\tau,l}\left[-\kappa_l^3\,H'_{\tau}(\kappa_l R)+\frac{{\tau}^2(1-\nu)}{R^2}\left(-\kappa_l\,H'_{\tau}(\kappa_l R)+\frac{1}{R}\,H_{\tau}(\kappa_l R)\right)\right]\right\}=0. \end{gather}$$For the discrete PTO system, the relationship between the disk deflection and the PTO forcing, i.e. (2.6), gives
\begin{equation} {-K} c_0\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{A_{m,l}\,{\rm J}_m(\kappa_lr_0)}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta_n}-F_n=0. \end{equation}For the continuum PTO system, (2.9) can be rewritten as
\begin{equation} {-K}\sum_{m'={-}\infty}^{\infty}c^{(m')}\,\mathrm{e}^{\mathrm{i} m'\theta}\sum_{m={-}\infty}^\infty\sum_{l={-}2}^\infty \frac{A_{m,l}\,{\rm J}_m(\kappa_lr_0)}{\chi \kappa_l^4+1-K\gamma}\,\mathrm{e}^{\mathrm{i} m\theta}-\sum_{m={-}\infty}^{\infty}f_m\,\mathrm{e}^{\mathrm{i} m\theta}=0, \end{equation}
which, after multiplying both sides by 
$\mathrm {e}^{-\mathrm {i} \tau \theta }$, integrating in 
$\theta \in [0, 2{\rm \pi} ]$ and using their orthogonality characteristics, gives
\begin{equation} {-K}\sum_{m={-}\infty}^{\infty}\sum_{l={-}2}^\infty \frac{A_{m,l}\,{\rm J}_{m}(\kappa_lr_0)\,c^{(\tau-m)}}{\chi \kappa_l^4+1-K\gamma} -f_{\tau}=0. \end{equation}For an angular-independent PTO system, (A27) gives
\begin{equation} {-K} c^{(0)}\sum_{l={-}2}^\infty \frac{A_{\tau,l}\,{\rm J}_{\tau}(\kappa_lr_0)}{\chi \kappa_l^4+1-K\gamma} -f_{\tau}=0. \end{equation}In order to evaluate the unknown coefficients 
$A_{m,l}$, 
$B_{m,l}$, 
$C_{m,l}$, 
$D_{m,l}$, and 
$F_n$ or 
$f_m$, which depend on the PTO types, we truncate all infinite series of vertical eigenfunctions at 
$L$, i.e. 
$(L+1)$ terms (
$l=0,1,\ldots,L$) for 
$D_{m,l}$, and 
$(L+3)$ terms (
$l = -2, -1, 0, 1,\ldots, L$) for 
$A_{m,l}$, 
$B_{m,l}$ and 
$C_{m,l}$, and we take 
$(2M+1)$ angular terms (
$m = -M, \ldots, 0, \ldots, M$), resulting in 
$2(2M + 1)(2L +5)+N$ and 
$(2M + 1)(4L +11)$ unknown coefficients to be determined for the discrete PTO and continuum PTO, respectively. After taking 
$\tau = -M,\ldots, 0,\ldots, M$ and 
$\zeta = 0, 1, \ldots, L$ in (A12)–(A27), 
$2(2M + 1)(2L +5)+N$ and 
$(2M + 1)(4L +11)$ equations are obtained for the discrete PTO and continuum PTO, respectively, which can be used to determine the exact same number of unknown coefficients.
Open access
