4.1 Leading-order problem $O(1)$

The $O(1)$ problem is homogeneous and unforced. Solution of the governing equations yields the well-known trapped natural modes of an array of $Q$ gates in an infinitely long channel (Mei et al. Reference Mei, Sammarco, Chan and Procaccini1994; Li & Mei Reference Li and Mei2003; Sammarco et al. Reference Sammarco, Michele and d’Errico2013). The boundary value problem at the leading order is governed by the following equations:

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{1}^{\pm }=0,\quad \text{in }\unicode[STIX]{x1D6FA}^{\pm }\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{1_{z}}^{\pm }=\unicode[STIX]{x1D719}_{1}^{\pm }\frac{\unicode[STIX]{x1D714}^{2}}{g},\quad z=0,\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{1}^{\pm }=\frac{\text{i}\unicode[STIX]{x1D714}}{g}\unicode[STIX]{x1D719}_{1}^{\pm },\quad z=0,\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{1_{x}}^{\pm }=\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D703}_{1}(z+h_{p})H(z+h-c),\quad x=x^{\pm }.\end{eqnarray}$$

By defining $\unicode[STIX]{x1D703}_{1}=\{r_{1q}\}\unicode[STIX]{x1D703}(t_{2})$ , with $\{r_{1q}\}=\{r_{11},\ldots ,r_{1Q}\}$ the modal shape, we obtain the solution of the velocity potential,

(4.8) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{1}^{\pm }=\mp \text{i}\unicode[STIX]{x1D703}\unicode[STIX]{x1D714}\mathop{\sum }_{q=1}^{Q}\mathop{\sum }_{m=1}^{\infty }\mathop{\sum }_{n=0}^{\infty }\frac{b_{mq}D_{n}}{C_{n}\unicode[STIX]{x1D6FC}_{mn}}\text{e}^{\mp \unicode[STIX]{x1D6FC}_{mn}x}\cos \frac{m\unicode[STIX]{x03C0}y}{b}\cosh k_{n}(h+z)\equiv \mp \text{i}\unicode[STIX]{x1D703}f_{1}^{\pm },\end{eqnarray}$$

and of the free-surface elevation,

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{1}^{\pm }=\pm \frac{\unicode[STIX]{x1D714}}{g}f_{1}^{\pm }\unicode[STIX]{x1D703},\end{eqnarray}$$

with $k_{n}$ being the roots of the dispersion relation

(4.10) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D714}^{2}=gk_{0}\tanh k_{0}h,\\ \unicode[STIX]{x1D714}^{2}=-g\overline{k}_{n}\tan \overline{k}_{n}h,\quad k_{n}=\text{i}\overline{k}_{n},n=1,\ldots ,\infty .\end{array}\right\}\end{eqnarray}$$

The real coefficients $b_{mq}$ , $\unicode[STIX]{x1D6FC}_{mn}$ , $C_{n}$ and $D_{n}$ are defined by

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle b_{mq}=r_{1q}\frac{2}{m\unicode[STIX]{x03C0}}\left[\sin \frac{qm\unicode[STIX]{x03C0}}{Q}-\sin \frac{(q-1)m\unicode[STIX]{x03C0}}{Q}\right], & \displaystyle\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}_{mn}=\sqrt{\left(\frac{m\unicode[STIX]{x03C0}}{b}\right)^{2}-k_{n}^{2}}, & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle C_{n}=\frac{1}{2}\left(h+\frac{g}{\unicode[STIX]{x1D714}^{2}}\sinh ^{2}k_{n}h\right), & \displaystyle\end{eqnarray}$$

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle D_{n}=\frac{1}{k_{n}^{2}}[\cosh k_{n}c-\cosh k_{n}h+k_{n}(h-c)\sinh k_{n}h]. & \displaystyle\end{eqnarray}$$

Note that the real function $f_{1}^{\pm }$ defined in (4.8) has the properties

(4.15a,b ) $$\begin{eqnarray}f_{1}^{+}(x)=f_{1}^{-}(-x),\quad f_{1_{x}}^{+}(x)=-f_{1_{x}}^{-}(-x),\end{eqnarray}$$

and that the wavenumber $k_{0}$ must satisfy the condition $k_{0}<\unicode[STIX]{x03C0}/b$ in order that the absence/existence of propagating waves/trapped modes is satisfied.

Conservation of angular momentum for each gate $G_{j}$ requires that

(4.16) $$\begin{eqnarray}r_{1j}(-\unicode[STIX]{x1D714}^{2}I+C)-\unicode[STIX]{x1D714}^{2}I_{j}=0,\end{eqnarray}$$

where $I_{j}$ is the added inertia defined as

(4.17) $$\begin{eqnarray}I_{j}=\mathop{\sum }_{q=1}^{Q}\mathop{\sum }_{m=1}^{\infty }\mathop{\sum }_{n=0}^{\infty }\frac{2b_{mq}b_{mj}D_{n}^{2}}{C_{n}\unicode[STIX]{x1D6FC}_{mn}},\end{eqnarray}$$

with

(4.18) $$\begin{eqnarray}b_{mj}=\frac{b}{m\unicode[STIX]{x03C0}}\left[\sin \frac{jm\unicode[STIX]{x03C0}}{Q}-\sin \frac{(j-1)m\unicode[STIX]{x03C0}}{Q}\right].\end{eqnarray}$$

As in Sammarco et al. (Reference Sammarco, Michele and d’Errico2013), solution of (4.16) yields $(Q-1)$ out-of-phase natural modes and related eigenfrequencies. The reader should refer to the work of Li & Mei (Reference Li and Mei2003) for a complete list of modal shapes in which $2\leqslant Q\leqslant 20$ .

4.2 The second-order problem $O(\unicode[STIX]{x1D716})$

Since the incident wave field is assumed to be at $O(\unicode[STIX]{x1D716})$ , the incident wave amplitude $A^{\prime }$ must be an order smaller than $A_{T}^{\prime }$ , thus $A^{\prime }/A_{T}^{\prime }=O(\unicode[STIX]{x1D716})$ . Hence, at $O(\unicode[STIX]{x1D716})$ we assume the coexistence of the second harmonic $2\unicode[STIX]{x1D714}$ with two components:

(4.19) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{2}^{\pm }=\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }}+\unicode[STIX]{x1D719}^{A^{\pm }},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }}$ is the second-order potential forced by the quadratic products (3.17) on the gate surface, while $\unicode[STIX]{x1D719}^{A^{\pm }}$ is the second-order potential forced by the incident wave field.

4.2.1 Radiated second harmonic – $\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }}$

The kinematic condition on the moving gates includes a forcing term,

(4.20) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{2_{x}}^{\unicode[STIX]{x1D703}^{\pm }}=\left[2\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D703}_{2}^{\unicode[STIX]{x1D703}}(z+h_{p})-\frac{\unicode[STIX]{x1D719}_{1_{z}}^{\pm }\unicode[STIX]{x1D703}_{1}}{\unicode[STIX]{x1D716}}\mp \text{i}\unicode[STIX]{x1D714}\,d\frac{\unicode[STIX]{x1D703}_{1}^{2}}{\unicode[STIX]{x1D716}}\right]H(z+h-c),\quad x=x^{\pm },\end{eqnarray}$$

while the governing equation and the remaining boundary conditions on the free surface remain homogeneous like the $O(1)$ problem:

(4.21) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }}=0,\end{eqnarray}$$

(4.22) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{2_{z}}^{\unicode[STIX]{x1D703}^{\pm }}=\frac{4\unicode[STIX]{x1D714}^{2}}{g}\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }},\quad z=0,\end{eqnarray}$$

(4.23) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{2}^{\unicode[STIX]{x1D703}^{\pm }}=\frac{2\text{i}\unicode[STIX]{x1D714}}{g}\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }},\quad z=0.\end{eqnarray}$$

As the forcing term in (4.20) contains only the second harmonic, we assume the solution of the potential $\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }}$ and of the angular motion $\unicode[STIX]{x1D703}_{2}$ in the form

(4.24a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }}=\text{i}\unicode[STIX]{x1D703}^{2}f_{2}^{\pm },\quad \unicode[STIX]{x1D703}_{2}^{\unicode[STIX]{x1D703}}=\text{i}\unicode[STIX]{x1D703}^{2}\unicode[STIX]{x1D703}_{2}.\end{eqnarray}$$

We get the following boundary value problem for the complex function $f_{2}^{\pm }$ :

(4.25) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}f_{2}^{\pm }=0,\end{eqnarray}$$

(4.26) $$\begin{eqnarray}f_{2_{z}}^{\pm }=\frac{4\unicode[STIX]{x1D714}^{2}}{g}f_{2}^{\pm },\quad z=0,\end{eqnarray}$$

(4.27) $$\begin{eqnarray}f_{2_{x}}^{\pm }=2\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D703}_{2}(z+h_{p})H(z+h-c)\pm \frac{1}{\unicode[STIX]{x1D716}}\mathop{\sum }_{p=0}^{\infty }\unicode[STIX]{x1D6E5}_{p}^{\pm }\cos \frac{p\unicode[STIX]{x03C0}y}{b},\quad x=x^{\pm },\end{eqnarray}$$

where

(4.28) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}_{p}^{\pm } & = & \displaystyle \frac{\unicode[STIX]{x1D6FF}_{p}}{b}\int _{0}^{b}\,\text{d}y(f_{1_{z}}^{\pm }r_{1q}-d\unicode[STIX]{x1D714}r_{1q}^{2})\cos \frac{p\unicode[STIX]{x03C0}y}{b}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2\unicode[STIX]{x1D6FF}_{p}}\mathop{\sum }_{q=1}^{Q}\left\{\mathop{\sum }_{m=1}^{\infty }b_{(m+p)q}\left[-d\unicode[STIX]{x1D714}b_{mq}+\mathop{\sum }_{n=0}^{\infty }\frac{\unicode[STIX]{x1D714}b_{mq}D_{n}}{C_{n}\unicode[STIX]{x1D6FC}_{mn}}k_{n}\sinh k_{n}(h+z)\right]\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\mathop{\sum }_{m=1}^{\infty }b_{mq}\left[-d\unicode[STIX]{x1D714}b_{(m+p)q}+\mathop{\sum }_{n=0}^{\infty }\frac{\unicode[STIX]{x1D714}b_{(m+p)q}D_{n}}{C_{n}\unicode[STIX]{x1D6FC}_{(m+p)n}}k_{n}\sinh k_{n}(h+z)\right]\right\}.\end{eqnarray}$$

Here $\unicode[STIX]{x1D6FF}_{p}$ denotes the Jacobi symbol, i.e.  $\unicode[STIX]{x1D6FF}_{0}=1$ and $\unicode[STIX]{x1D6FF}_{p}=2,p=1,\ldots .$ Note that the latter term is the same in both fluid regions: $\unicode[STIX]{x1D6E5}_{p}^{+}\equiv \unicode[STIX]{x1D6E5}_{p}^{-}$ .

Solution of the problem can be found by separation of variables:

(4.29) $$\begin{eqnarray}f_{2}^{\pm }=-\text{i}\mathop{\sum }_{p=0}^{\infty }\mathop{\sum }_{l=0}^{\infty }\frac{1}{\unicode[STIX]{x1D6FC}_{pl}}\left(\frac{\unicode[STIX]{x1D6E5}_{pl}}{\unicode[STIX]{x1D716}}\pm \mathop{\sum }_{q=1}^{Q}\frac{2\text{i}b_{pq}\unicode[STIX]{x1D714}D_{l}}{C_{l}}\right)\text{e}^{\pm \text{i}\unicode[STIX]{x1D6FC}_{pl}x}\cos \frac{p\unicode[STIX]{x03C0}y}{b}\cosh \unicode[STIX]{x1D705}_{l}(h+z),\end{eqnarray}$$

where

(4.30) $$\begin{eqnarray}b_{mp}=r_{2q}\frac{2}{p\unicode[STIX]{x03C0}}\left[\sin \frac{qp\unicode[STIX]{x03C0}}{Q}-\sin \frac{(q-1)p\unicode[STIX]{x03C0}}{Q}\right],\end{eqnarray}$$

$\unicode[STIX]{x1D705}_{l}$ are the real roots of the dispersion relation

(4.31) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}4\unicode[STIX]{x1D714}^{2}=g\unicode[STIX]{x1D705}_{0}\tanh \unicode[STIX]{x1D705}_{0}h,\\ \unicode[STIX]{x1D705}_{l}=\text{i}\overline{\unicode[STIX]{x1D705}}_{l},\quad l=1,\ldots ,\infty ,\\ 4\unicode[STIX]{x1D714}^{2}=-g\overline{\unicode[STIX]{x1D705}}_{l}\tan \overline{\unicode[STIX]{x1D705}}_{l}h,\end{array}\right\}\end{eqnarray}$$

while the coefficient $\unicode[STIX]{x1D6E5}_{pl}$ is given by

(4.32) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}_{pl} & = & \displaystyle \frac{1}{C_{l}}\int _{-h}^{0}\,\text{d}z\unicode[STIX]{x1D6E5}_{p}\cosh \unicode[STIX]{x1D705}_{l}(h+z)\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2\unicode[STIX]{x1D6FF}_{p}C_{l}}\mathop{\sum }_{q=1}^{Q}\mathop{\sum }_{m=1}^{\infty }\mathop{\sum }_{n=0}^{\infty }\left\{b_{(m+p)q}\left[-d\unicode[STIX]{x1D714}b_{mq}E_{l}+\frac{\unicode[STIX]{x1D714}b_{mq}D_{n}C_{ln}}{C_{n}\unicode[STIX]{x1D6FC}_{mn}}\right]\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,b_{mq}\left[-d\unicode[STIX]{x1D714}b_{(m+p)q}E_{l}+\frac{\unicode[STIX]{x1D714}b_{(m+p)q}D_{n}C_{ln}}{C_{n}\unicode[STIX]{x1D6FC}_{(m+p)n}}\right]\right\},\end{eqnarray}$$

in which

(4.33) $$\begin{eqnarray}\displaystyle C_{ln} & = & \displaystyle \frac{1}{\unicode[STIX]{x1D705}_{l}^{2}-k_{n}^{2}}\left[\cosh k_{n}h\cosh \unicode[STIX]{x1D705}_{l}h\left(\frac{4\unicode[STIX]{x1D714}^{2}}{g^{2}}-k_{n}^{2}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,k_{n}^{2}\cosh k_{n}c\cosh \unicode[STIX]{x1D705}_{l}c-\unicode[STIX]{x1D705}_{l}k_{n}\sinh k_{n}h\sinh \unicode[STIX]{x1D705}_{l}h\right],\end{eqnarray}$$

(4.34) $$\begin{eqnarray}E_{l}=\frac{1}{\unicode[STIX]{x1D705}_{l}}(\sinh \unicode[STIX]{x1D705}_{l}h-\sinh \unicode[STIX]{x1D705}_{l}c).\end{eqnarray}$$

Because of the symmetry properties (4.15), the equation of motion reads

(4.35) $$\begin{eqnarray}r_{2j}(-4\unicode[STIX]{x1D714}^{2}I+C)=2\text{i}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D714}\int _{(j-1)a}^{ja}\,\text{d}y\int _{-h_{p}}^{0}\text{i}\unicode[STIX]{x1D703}\unicode[STIX]{x0394}f_{2}(z+h_{p})\,\text{d}z,\end{eqnarray}$$

however,

(4.36) $$\begin{eqnarray}\unicode[STIX]{x0394}f_{2}=\mathop{\sum }_{p=0}^{\infty }\mathop{\sum }_{l=0}^{\infty }\frac{4}{\unicode[STIX]{x1D6FC}_{pl}}\mathop{\sum }_{q=1}^{Q}\frac{b_{pq}\unicode[STIX]{x1D714}D_{l}}{C_{l}}\cos \frac{p\unicode[STIX]{x03C0}y}{b}\cosh \unicode[STIX]{x1D705}_{l}(h+z),\quad x=x^{\pm },\end{eqnarray}$$

so that substitution of the above expression in (4.35) yields a system of homogeneous equations in $r_{2j}$ , $j=1,2$ , with non-zero determinant of the coefficient matrix. As a result, the only solution is $r_{2q}=0$ , $q=1,\ldots ,Q$ , i.e. the angular motion is zero. Finally, the complete solution of the second-order velocity potential forced by the trapped mode is

(4.37) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{2}^{\unicode[STIX]{x1D703}^{\pm }}=\unicode[STIX]{x1D703}^{2}\mathop{\sum }_{p=0}^{\infty }\mathop{\sum }_{l=0}^{\infty }\frac{\unicode[STIX]{x1D6E5}_{pl}}{\unicode[STIX]{x1D6FC}_{pl}\unicode[STIX]{x1D716}}\text{e}^{\pm \text{i}\unicode[STIX]{x1D6FC}_{pl}x}\cos \frac{p\unicode[STIX]{x03C0}y}{b}\cosh \unicode[STIX]{x1D705}_{l}(h+z),\end{eqnarray}$$

while the corresponding free-surface displacement is

(4.38) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{2}^{\unicode[STIX]{x1D703}^{\pm }}=-\frac{2\unicode[STIX]{x1D714}}{g}f_{2}^{\pm }\unicode[STIX]{x1D703}^{2}.\end{eqnarray}$$

4.2.2 Response to incident waves – $\unicode[STIX]{x1D719}^{A^{\pm }}$

This is a simple diffraction problem where the incident wave field forces the gates to move at unison in phase. Decompose $\unicode[STIX]{x1D719}^{A^{\pm }}$ as follows:

(4.39) $$\begin{eqnarray}\unicode[STIX]{x1D719}^{A^{\pm }}=\unicode[STIX]{x1D719}^{I}+\unicode[STIX]{x1D719}^{S}+\unicode[STIX]{x1D719}^{R^{\pm }},\end{eqnarray}$$

where

(4.40) $$\begin{eqnarray}\unicode[STIX]{x1D719}^{I}=-\frac{\text{i}Ag}{4\unicode[STIX]{x1D714}}\frac{\cosh \unicode[STIX]{x1D705}_{0}(h+z)}{\cosh \unicode[STIX]{x1D705}_{0}h}\text{e}^{-\text{i}\unicode[STIX]{x1D705}_{0}x}\end{eqnarray}$$

is the velocity potential of the incident waves of amplitude $A$ ,

(4.41) $$\begin{eqnarray}\unicode[STIX]{x1D719}^{S}=-\frac{\text{i}Ag}{4\unicode[STIX]{x1D714}}\frac{\cosh \unicode[STIX]{x1D705}_{0}(h+z)}{\cosh \unicode[STIX]{x1D705}_{0}h}\text{e}^{\text{i}\unicode[STIX]{x1D705}_{0}x}\end{eqnarray}$$

is the scattered wave potential, while

(4.42) $$\begin{eqnarray}\unicode[STIX]{x1D719}^{R^{\pm }}=\pm \mathop{\sum }_{l=0}^{\infty }\frac{2\unicode[STIX]{x1D714}\unicode[STIX]{x1D703}^{A}D_{l}}{\unicode[STIX]{x1D705}_{l}C_{l}}\cosh \unicode[STIX]{x1D705}_{0}(h+z)\text{e}^{\pm \text{i}\unicode[STIX]{x1D705}_{0}x}\end{eqnarray}$$

is the radiation potential. The angular response $\unicode[STIX]{x1D703}^{A}$ is given by

(4.43) $$\begin{eqnarray}\unicode[STIX]{x1D703}^{A}=\frac{\unicode[STIX]{x1D70C}aAgD_{0}/\cosh \unicode[STIX]{x1D705}_{0}h}{-4\unicode[STIX]{x1D714}^{2}I+C-8\text{i}\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D70C}a\displaystyle \mathop{\sum }_{l=0}^{\infty }\frac{D_{l}^{2}}{\unicode[STIX]{x1D705}_{l}C_{l}}}.\end{eqnarray}$$

4.3 The third-order problem $O(\unicode[STIX]{x1D716}^{2})$

At the order $O(\unicode[STIX]{x1D716}^{2})$ , the boundary conditions on the free surface and on the gate array surface are inhomogeneous. Since $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D719}_{1}$ solve the homogeneous problem at the leading order we invoke the solvability condition applying Green’s theorem to $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{3}$ over the entire fluid domain $\unicode[STIX]{x1D6FA}^{\pm }$ :

(4.44) $$\begin{eqnarray}\iiint _{\unicode[STIX]{x1D6FA}^{\pm }}(\unicode[STIX]{x1D719}_{1}^{\pm }\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{3}^{\pm }-\unicode[STIX]{x1D719}_{3}^{\pm }\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{1}^{\pm })\,\text{d}\unicode[STIX]{x1D6FA}^{\pm }=\iint _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{\pm }}\left(\unicode[STIX]{x1D719}_{1}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}n}-\unicode[STIX]{x1D719}_{3}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}n}\right)\,\text{d}S^{\pm },\end{eqnarray}$$

where the normal $n$ points outward the volume boundaries $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ and the volume integral on the left-hand side is null. Performing straightforward algebra we obtain

(4.45) $$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \frac{\unicode[STIX]{x1D714}^{2}}{g}\int _{0}^{b}\,\text{d}y\left(\int _{-\infty }^{0}f_{1}^{-}{\mathcal{F}}_{3}^{-}\,\text{d}x-\int _{0}^{\infty }f_{1}^{+}{\mathcal{F}}_{3}^{+}\,\text{d}x\right)\nonumber\\ \displaystyle & & \displaystyle +\,2\int _{0}^{b}\,\text{d}y\int _{-h}^{0}\overline{f_{1}{\mathcal{G}}_{3}}\,\text{d}z+\mathop{\sum }_{q=1}^{Q}\frac{1}{\unicode[STIX]{x1D70C}}\left(\text{i}{\mathcal{D}}_{3}r_{1q}-\frac{\unicode[STIX]{x1D708}_{pto}r_{1q}^{2}\unicode[STIX]{x1D714}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D716}^{2}}\right),\end{eqnarray}$$

where the forcing terms ${\mathcal{F}}_{3}$ , ${\mathcal{G}}_{3}$ , ${\mathcal{D}}_{3}$ are given by

(4.46) $$\begin{eqnarray}{\mathcal{F}}_{3}^{\pm }=\pm \frac{2f_{1}^{\pm }\unicode[STIX]{x1D703}_{t}}{\unicode[STIX]{x1D714}},\end{eqnarray}$$

(4.47) $$\begin{eqnarray}\displaystyle {\mathcal{G}}_{3}^{\pm } & = & \displaystyle [\!-r_{1q}\unicode[STIX]{x1D703}^{\ast }(\text{i}\unicode[STIX]{x1D703}^{2}f_{2_{z}}^{\pm }+A\unicode[STIX]{x1D719}_{z}^{A^{\pm }})\mp \text{i}\unicode[STIX]{x1D703}^{A}f_{1_{z}}^{\pm }\unicode[STIX]{x1D703}^{\ast }-r_{1q}\unicode[STIX]{x1D703}_{t}(z+h_{p})\nonumber\\ \displaystyle & & \displaystyle \pm \,d(\unicode[STIX]{x1D703}_{1}^{\ast }\unicode[STIX]{x1D703}_{2})_{t}\! ]H(z+h-c),\end{eqnarray}$$

(4.48) $$\begin{eqnarray}\displaystyle {\mathcal{D}}_{3}^{\pm } & = & \displaystyle -\unicode[STIX]{x1D70C}\int _{(q-1)a}^{qa}\,\text{d}y\left\{\int _{-h_{p}}^{0}f_{1}(-6r_{1q}^{2}\unicode[STIX]{x1D703}^{\ast }\unicode[STIX]{x1D703}^{2}\unicode[STIX]{x1D714}-2\text{i}\unicode[STIX]{x1D703}_{t})(z+h_{p})\,\text{d}z\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\,\frac{\unicode[STIX]{x1D714}^{2}f_{1}}{g}\unicode[STIX]{x1D703}^{\ast }\left[-\frac{\unicode[STIX]{x1D714}f_{1}^{2}}{g}\unicode[STIX]{x1D703}^{2}+4h_{p}f_{2}\unicode[STIX]{x1D703}^{2}-4\text{i}h_{p}\unicode[STIX]{x1D719}^{D}\right]\right\}+2\unicode[STIX]{x1D714}Ir_{1q}\text{i}\unicode[STIX]{x1D703}_{t}-\frac{1}{2}gSr_{1q}^{3}\unicode[STIX]{x1D703}^{2}\unicode[STIX]{x1D703}^{\ast }\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D70C}\unicode[STIX]{x1D714}\,d\unicode[STIX]{x1D703}^{\ast }\int _{(q-1)a}^{qa}r_{1q}\,\text{d}y\left\{\int _{-h_{p}}^{0}(-4f_{2}\unicode[STIX]{x1D703}^{2}+2\text{i}\unicode[STIX]{x1D719}^{D})\,\text{d}z-3\unicode[STIX]{x1D714}\frac{f_{1}^{2}}{g}\unicode[STIX]{x1D703}^{2}\right\},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}^{D}=\unicode[STIX]{x1D719}^{I}+\unicode[STIX]{x1D719}^{S}$ is the diffraction potential and $\unicode[STIX]{x1D703}^{\ast }$ denotes the complex conjugate of $\unicode[STIX]{x1D703}$ . Solution of the integrals in (4.45) gives the evolution equation of the Stuart–Landau type:

(4.49) $$\begin{eqnarray}-\text{i}\unicode[STIX]{x1D703}_{t}=\unicode[STIX]{x1D703}^{2}\unicode[STIX]{x1D703}^{\ast }(c_{N}+\text{i}c_{R})+A\unicode[STIX]{x1D703}^{\ast }c_{F}+\text{i}\unicode[STIX]{x1D703}\unicode[STIX]{x1D708}_{pto}c_{L}.\end{eqnarray}$$

If there is a small detuning $2\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ such that $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}\sim O(\unicode[STIX]{x1D716}^{2})$ , the above equation modifies as follows:

(4.50) $$\begin{eqnarray}-\text{i}\unicode[STIX]{x1D703}_{t}=\unicode[STIX]{x1D703}^{2}\unicode[STIX]{x1D703}^{\ast }(c_{N}+\text{i}c_{R})+A\text{e}^{-2\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}t}\unicode[STIX]{x1D703}^{\ast }c_{F}+\text{i}\unicode[STIX]{x1D703}\unicode[STIX]{x1D708}_{pto}c_{L}.\end{eqnarray}$$

Through the transformation

(4.51) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\overline{\unicode[STIX]{x1D703}}\text{e}^{-\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}t},\end{eqnarray}$$

equation (4.50) becomes

(4.52) $$\begin{eqnarray}-\text{i}\overline{\unicode[STIX]{x1D703}}_{t}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}\overline{\unicode[STIX]{x1D703}}+\overline{\unicode[STIX]{x1D703}}^{2}\overline{\unicode[STIX]{x1D703}}^{\ast }(c_{N}+\text{i}c_{R})+A\overline{\unicode[STIX]{x1D703}}^{\ast }c_{F}+\text{i}\overline{\unicode[STIX]{x1D703}}\unicode[STIX]{x1D708}_{pto}c_{L}.\end{eqnarray}$$

The coefficient $c_{N}$ represents the shift of the trapped mode eigenfrequency from the incident wave frequency, $c_{R}$ is the radiation damping due to wave radiation at $O(\unicode[STIX]{x1D716})$ , $c_{F}$ represents the energy influx by the incident waves while $c_{L}$ represents damping due to the PTO mechanism.

The coefficients of the evolution equation (4.52) $c_{L}$ , $c_{F}$ , $c_{N}+\text{i}c_{R}$ , are given respectively by

(4.53) $$\begin{eqnarray}\displaystyle & \displaystyle c_{L}=\frac{\unicode[STIX]{x1D714}}{c_{T}\unicode[STIX]{x1D70C}}\mathop{\sum }_{q=1}^{Q}r_{1q}^{2}, & \displaystyle\end{eqnarray}$$

(4.54) $$\begin{eqnarray}\displaystyle & \displaystyle c_{F}=\frac{1}{\text{i}c_{T}}\left(\int _{0}^{b}\,\text{d}y\int _{-h_{p}}^{0}(-f_{1}r_{1q}\unicode[STIX]{x1D719}_{z}^{D}+2\unicode[STIX]{x1D714}\,dr_{1q}^{2}\unicode[STIX]{x1D719}^{D})\,\text{d}z-\int _{0}^{b}2r_{1q}\unicode[STIX]{x1D714}^{2}h_{p}\unicode[STIX]{x1D719}^{D}\frac{f_{1}}{g}\,\text{d}y\right),\qquad & \displaystyle\end{eqnarray}$$

(4.55) $$\begin{eqnarray}\displaystyle c_{N}+\text{i}c_{R} & = & \displaystyle \frac{1}{c_{T}}\left\{\int _{0}^{b}\,\text{d}y\int _{-h_{p}}^{0}[-2\overline{f_{1}^{\pm }r_{1q}f_{2_{z}}^{\pm }}+6\unicode[STIX]{x1D714}f_{1}r_{1q}^{3}(z+h_{p})-4\text{i}dr_{1q}^{2}f_{2}\unicode[STIX]{x1D714}]\,\text{d}z\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\int _{0}^{b}r_{1q}\unicode[STIX]{x1D714}^{2}\left(\unicode[STIX]{x1D714}\frac{f_{1}^{3}}{g^{2}}-4h_{p}\frac{f_{1}}{g}f_{2}+3dr_{1q}\frac{f_{1}^{2}}{g}\right)\,\text{d}y-\frac{gSr_{1q}^{4}}{2\unicode[STIX]{x1D70C}}\right\},\end{eqnarray}$$

with

(4.56) $$\begin{eqnarray}c_{T}=\frac{4\unicode[STIX]{x1D714}}{g}\int _{0}^{b}\,\text{d}y\int _{0}^{\infty }f_{1}^{2}\,\text{d}x+2\int _{0}^{b}\,\text{d}y\int _{-h_{p}}^{0}f_{1}r_{1q}(z+h_{p})\,\text{d}z+\mathop{\sum }_{q=1}^{Q}2I\unicode[STIX]{x1D714}r_{1q}^{2}.\end{eqnarray}$$

Computation of $c_{F}$ , $c_{N}$ , $c_{R}$ and $c_{T}$ can be evaluated by numerical calculation of the integrals with the known expressions of the integrands.

4.4 Dynamical system analysis for uniform incident waves

By making use of action-angle variables $R$ and $\unicode[STIX]{x1D713}$ expressed by $\overline{\unicode[STIX]{x1D703}}=\text{i}\sqrt{R}\text{e}^{\text{i}\unicode[STIX]{x1D713}}$ we obtain the following autonomous dynamical system in $R$ and $\unicode[STIX]{x1D713}$ :

(4.57) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}R_{t}=-2R(-2c_{R}R+c_{F}A\sin 2\unicode[STIX]{x1D713}+c_{L}\unicode[STIX]{x1D708}_{pto})\\ \unicode[STIX]{x1D713}_{t}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}+c_{N}R-c_{F}A\cos 2\unicode[STIX]{x1D713}.\end{array}\right\}\end{eqnarray}$$

With the requirement $R_{t}=0$ , $\unicode[STIX]{x1D713}_{t}=0$ we find at least three fixed points. The trivial fixed point is located at

(4.58a,b ) $$\begin{eqnarray}R=0,\quad \unicode[STIX]{x1D713}=(1/2)\cos ^{-1}\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}}{c_{F}A}=\unicode[STIX]{x1D713}^{\ast }.\end{eqnarray}$$

It is easily shown that this point is an unstable saddle for $|\unicode[STIX]{x0394}\unicode[STIX]{x1D714}|<\sqrt{A^{2}c_{F}^{2}-\unicode[STIX]{x1D708}_{pto}^{2}c_{L}^{2}}$ , and a stable fixed point otherwise. The latter condition fixes the threshold of resonance,

(4.59) $$\begin{eqnarray}A>\frac{c_{L}}{c_{F}}\unicode[STIX]{x1D708}_{pto},\end{eqnarray}$$

i.e. for large values of the PTO coefficient $\unicode[STIX]{x1D708}_{pto}$ there is no subharmonic resonance from the rest position and only in-phase motion occurs. Now we seek the non-trivial fixed points which correspond to the roots of the quadratic equation:

(4.60) $$\begin{eqnarray}R^{2}(c_{N}^{2}+c_{R}^{2})+2R(c_{R}\unicode[STIX]{x1D708}_{pto}c_{L}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714}c_{N})+\unicode[STIX]{x1D708}_{pto}^{2}c_{L}^{2}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714}^{2}-A^{2}c_{F}^{2}=0,\end{eqnarray}$$

i.e.

(4.61) $$\begin{eqnarray}R^{\pm }=\frac{-\unicode[STIX]{x1D708}_{pto}c_{L}c_{R}-c_{N}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}\pm \sqrt{Ac_{F}^{2}(c_{N}^{2}+c_{R}^{2})-\unicode[STIX]{x1D708}_{pto}c_{L}c_{N}+c_{R}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}}}{c_{N}^{2}+c_{R}^{2}}.\end{eqnarray}$$

The latter corresponds to two branches of an ellipse in the plane $\unicode[STIX]{x0394}\unicode[STIX]{x1D714},R$ . Since the square root of (4.61) must be real and $R^{\pm }>0$ , we obtain a single non-trivial fixed point $R^{+}$ for

(4.62) $$\begin{eqnarray}-\sqrt{A^{2}c_{F}^{2}-\unicode[STIX]{x1D708}_{pto}^{2}c_{L}^{2}}<\unicode[STIX]{x0394}\unicode[STIX]{x1D714}<\sqrt{A^{2}c_{F}^{2}-\unicode[STIX]{x1D708}_{pto}^{2}c_{L}^{2}},\end{eqnarray}$$

and the coexistence of two non-trivial fixed points $R^{\pm }$ for

(4.63) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D708}_{pto}c_{L}c_{N}-Ac_{F}\sqrt{c_{N}^{2}+c_{R}^{2}}}{c_{R}}<\unicode[STIX]{x0394}\unicode[STIX]{x1D714}<-\sqrt{A^{2}c_{F}^{2}-\unicode[STIX]{x1D708}_{pto}^{2}c_{L}^{2}}.\end{eqnarray}$$

The latter holds if

(4.64) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{pto}<\frac{c_{N}Ac_{F}}{c_{L}\sqrt{c_{N}^{2}+c_{R}^{2}}},\end{eqnarray}$$

this is because $\unicode[STIX]{x1D708}_{pto}$ alters the threshold of instability and determines the vanishing of $R^{-}$ . Analysis of the Jacobian matrix at the fixed points reveals that $R^{+}$ is stable while $R^{-}$ correspond to an unstable saddle. It is worth evaluating the maximum value of the equilibrium amplitude which corresponds to the maximum of the stable state $R^{+}$ . Taking the derivative of $R^{+}$ with respect to $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ we obtain

(4.65a,b ) $$\begin{eqnarray}R_{max}=\frac{Ac_{F}-\unicode[STIX]{x1D708}_{pto}c_{L}}{c_{R}},\quad \text{at }\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{max}=\frac{c_{N}(\unicode[STIX]{x1D708}_{pto}c_{L}-Ac_{F})}{c_{R}},\end{eqnarray}$$

or in terms of the modal amplitude:

(4.66) $$\begin{eqnarray}|\overline{\unicode[STIX]{x1D703}}|_{max}=\sqrt{\frac{Ac_{F}-\unicode[STIX]{x1D708}_{pto}c_{L}}{c_{R}}}.\end{eqnarray}$$

Once $|\overline{\unicode[STIX]{x1D703}}|$ is evaluated, the average generated power by the array due to subharmonic resonance of the natural mode is given by

(4.67) $$\begin{eqnarray}P_{s}(\unicode[STIX]{x1D708}_{pto},R^{+})=\frac{\unicode[STIX]{x1D714}}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}}\unicode[STIX]{x1D708}_{pto}\mathop{\sum }_{q=1}^{Q}\left(\frac{\text{d}\unicode[STIX]{x1D6E9}_{q}}{\text{d}t}\right)^{2}\,\text{d}t=2\unicode[STIX]{x1D708}_{pto}(\unicode[STIX]{x1D714}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714})^{2}\mathop{\sum }_{q=1}^{Q}r_{1q}^{2}R^{+}.\end{eqnarray}$$

The maximum of $P_{s}$ depends on $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D708}_{pto}$ and cannot be found analytically. Numerical calculation of (4.67) reveals that a good approximation consists in finding the optimal value of $\unicode[STIX]{x1D708}_{pto}$ which maximizes the generated power at $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{max}$ where the gate oscillation is maximum. So we impose the derivative of (4.67) with respect to $\unicode[STIX]{x1D708}_{pto}$ evaluated at $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{max}$ to be zero:

(4.68) $$\begin{eqnarray}\left.\frac{\text{d}P_{s}}{\text{d}\unicode[STIX]{x1D708}_{pto}}\right|_{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{max}}=0.\end{eqnarray}$$

Solution of the latter equation yields a criterion to find the optimal value of the PTO coefficient which maximizes the power output:

(4.69) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{pto}=\frac{c_{T}\unicode[STIX]{x1D70C}}{8c_{N}\unicode[STIX]{x1D714}\displaystyle \mathop{\sum }_{q=1}^{Q}r_{1q}^{2}}\left(5Ac_{F}c_{N}-2c_{R}\unicode[STIX]{x1D714}+\sqrt{9A^{2}c_{F}^{2}c_{N}^{2}-4Ac_{F}c_{N}c_{R}\unicode[STIX]{x1D714}+4c_{R}^{2}\unicode[STIX]{x1D714}^{2}}\right).\end{eqnarray}$$

Finally we assess the maximum efficiency of the system by considering the capture width ratio ${\mathcal{C}}^{F}$ (Mei et al. Reference Mei, Stiassnie and Yue2005) defined as the ratio between $P_{s}$ and the incident wave energy flux per array width $b$ (Michele et al. Reference Michele, Sammarco and d’Errico2016b ):

(4.70) $$\begin{eqnarray}{\mathcal{C}}^{F}=\frac{P_{s}}{EC_{g}b},\end{eqnarray}$$

where

(4.71) $$\begin{eqnarray}EC_{G}=\frac{\unicode[STIX]{x1D70C}gA^{2}(\unicode[STIX]{x1D714}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714})}{2k}\left(1+\frac{2kh}{\sinh 2kh}\right).\end{eqnarray}$$

In the latter, $C_{G}$ is the group velocity and $k$ the wavenumber related to the frequency $2(\unicode[STIX]{x1D714}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714})$ .

4.5 The non-autonomous dynamical system for modulated incident waves

In this section we investigate the effects of a modulation of the incident wave envelope on the dynamics and power generated by the array. Let us assume

(4.72) $$\begin{eqnarray}A=\overline{A}+\tilde{A}\cos \unicode[STIX]{x1D6FA}t,\end{eqnarray}$$

where $\overline{A}$ is the constant amplitude of the short waves, while $\tilde{A}$ and $\unicode[STIX]{x1D6FA}$ are respectively the amplitude and frequency of the modulation with order $O(\unicode[STIX]{x1D716}^{2})$ . The evolution equation now includes a forcing term in time and becomes

(4.73) $$\begin{eqnarray}-\text{i}\overline{\unicode[STIX]{x1D703}}_{t}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}\overline{\unicode[STIX]{x1D703}}+\overline{\unicode[STIX]{x1D703}}^{2}\overline{\unicode[STIX]{x1D703}}^{\ast }(c_{N}+\text{i}c_{R})+(\overline{A}+\tilde{A}\cos \unicode[STIX]{x1D6FA}t)\overline{\unicode[STIX]{x1D703}}^{\ast }c_{F}+\text{i}\overline{\unicode[STIX]{x1D703}}c_{L}.\end{eqnarray}$$

In action-angle coordinates the latter equation reads

(4.74) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}R_{t}=-2R[-2Rc_{R}+(\overline{A}+\tilde{A}\cos \unicode[STIX]{x1D6FA}t)c_{F}\sin 2\unicode[STIX]{x1D713}+c_{L}]\\ \unicode[STIX]{x1D713}_{t}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}+Rc_{N}-(\overline{A}+\tilde{A}\cos \unicode[STIX]{x1D6FA}t)c_{F}\cos 2\unicode[STIX]{x1D713}.\end{array}\right\}\end{eqnarray}$$

The latter non-autonomous system is similar to that for Venice gates (Sammarco et al. Reference Sammarco, Tran, Gottlieb and Mei1997b ) and can exhibit chaos for a certain range of the long wave amplitude $\tilde{A}$ . There are different theoretical criteria to determine under what conditions the response of a dynamical system becomes chaotic. In this paper we make use of the Melnikov method. This mathematical technique allows us to predict global chaos and is based on the search of horseshoe maps and homoclinic/heteroclinic orbits of the associated undamped Hamiltonian system (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983). Let us assume the following variables:

(4.75) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D6FC}=\frac{c_{N}}{c_{R}},\quad \unicode[STIX]{x1D6FD}=\frac{c_{L}}{\overline{A}c_{F}},\quad W=\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}}{\overline{A}c_{F}},\quad \unicode[STIX]{x1D717}=\sqrt{\frac{c_{N}}{\overline{A}c_{F}}}\unicode[STIX]{x1D703},\\ \displaystyle T=\overline{A}c_{F}t,\quad a=\frac{\tilde{A}}{\overline{A}},\quad \unicode[STIX]{x1D70E}=\frac{\unicode[STIX]{x1D6FA}}{c_{F}\overline{A}},\end{array}\right\}\end{eqnarray}$$

and the damping and forcing terms to be $O(\unicode[STIX]{x1D6FF})$ with $\unicode[STIX]{x1D6FF}\ll 1$ . By denoting with $\unicode[STIX]{x1D717}=\text{i}\sqrt{R^{\prime }}\text{e}^{\text{i}\unicode[STIX]{x1D713}^{\prime }}$ , the system (4.74) can be rewritten in the form

(4.76) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}R_{t}^{\prime }=-2R\sin 2\unicode[STIX]{x1D713}^{\prime }-2R^{\prime }\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D6FC}R^{\prime }+a\sin 2\unicode[STIX]{x1D713}^{\prime }\cos \unicode[STIX]{x1D70E}T+\unicode[STIX]{x1D6FD})\\ \unicode[STIX]{x1D713}_{t}^{\prime }=W+R^{\prime }-\cos 2\unicode[STIX]{x1D713}^{\prime }-\unicode[STIX]{x1D6FF}a\cos 2\unicode[STIX]{x1D713}^{\prime }\cos \unicode[STIX]{x1D70E}T.\end{array}\right\}\end{eqnarray}$$

The Hamiltonian function of the unforced system is then given by

(4.77) $$\begin{eqnarray}H=\frac{R^{\prime 2}}{2}-R^{\prime }\cos 2\unicode[STIX]{x1D713}^{\prime }+WR^{\prime }.\end{eqnarray}$$

At $O(1)$ the system (4.76) admits three fixed points:

(4.78a-c ) $$\begin{eqnarray}O=\{0,{\textstyle \frac{1}{2}}\cos ^{-1}W\},\quad s=\{1-W,0\},\quad u=\{-1-W,0\},\end{eqnarray}$$

where $O$ is an unstable saddle for $|W|<1$ and a centre if $|W|>1$ , the fixed point $s$ is a centre and exist for $W<1$ while $u$ is an unstable saddle and exist for $W<-1$ . Because heteroclinic tangle is difficult to excite experimentally (Sammarco et al. Reference Sammarco, Tran, Gottlieb and Mei1997b ), we focus our attention to the bifurcation of homoclinic paths assuming $W=0$ . In this case the homoclinic orbit has the equation

(4.79a,b ) $$\begin{eqnarray}R^{h}=\frac{2}{\cosh 2(T-T_{0})},\quad \unicode[STIX]{x1D713}^{h}=\tan ^{-1}\{\tanh (T-T_{0})\},\end{eqnarray}$$

and the corresponding Melnikov function is given by

(4.80) $$\begin{eqnarray}M(T_{0})=\int _{-\infty }^{\infty }2R^{h}[a_{1}R^{h}\sin 2\unicode[STIX]{x1D713}^{h}\cos \unicode[STIX]{x1D70E}T+(R^{h}-\cos 2\unicode[STIX]{x1D713}^{h})(\unicode[STIX]{x1D6FC}_{1}R^{h}+\unicode[STIX]{x1D6FD}_{1})]\,\text{d}T,\end{eqnarray}$$

where

(4.81) $$\begin{eqnarray}(\unicode[STIX]{x1D6FC}_{1},\unicode[STIX]{x1D6FD}_{1},a1)=\frac{(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},a)}{\unicode[STIX]{x1D707}}.\end{eqnarray}$$

Homoclinic tangle occurs when $M(T_{0})$ has a simple zero. The integral (4.80) can be evaluated analytically, and gives the lower threshold in terms of the modulated amplitude $\tilde{A}$ to give rise to horseshoe structures in the neighbourhood of the homoclinic orbit $R^{h},\unicode[STIX]{x1D713}^{h}$ :

(4.82) $$\begin{eqnarray}\tilde{A}=\frac{4\overline{A}^{3}c_{F}^{2}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FA}^{2}}\left(\frac{\unicode[STIX]{x03C0}c_{R}}{c_{N}}+\frac{2c_{L}}{\overline{A}c_{F}}\right)\sinh \frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FA}}{4\overline{A}c_{F}}.\end{eqnarray}$$

Note that the smaller the damping due to $O(\unicode[STIX]{x1D716})$ radiated waves and PTO coefficient, the lower is the threshold. Moreover, taking the limit of (4.82) for $\unicode[STIX]{x1D6FA}\rightarrow (0,\infty )$ we obtain that $\tilde{A}$ tends to infinity and chaos is not generated.

4.6 Results for uniform incident waves and $Q=2$

Here we discuss the theoretical results of the previous section for the case of $Q=2$ gates with channel and array dimensions listed in table 3. Let us assume the amplitude of the incident wave $A=0.1$ m. Solution of the $O(1)$ problem yields a trapped mode with eigenvector $r_{1q}=\{1,-1\}$ . The values of the coefficients $c_{L}$ , $c_{N}$ , $c_{R}$ and $c_{F}$ with respect to the eigenfrequency $\unicode[STIX]{x1D714}$ are shown in figure 2. All the coefficients show similar behaviour, i.e. they increase with the eigenfrequency of the mode.

Table 3. Flap and channel dimensions.

Figure 2. Behaviour of the coefficients of the evolution equation versus the eigenfrequency $\unicode[STIX]{x1D714}$ .

In figure 3 we show the dependence of the maximum value of the capture factor ${\mathcal{C}}_{max}^{F}$ (4.70) and of the optimal value of the PTO coefficient (4.69) on the eigenfrequency $\unicode[STIX]{x1D714}$ . The larger the eigenfrequency the larger the efficiency of the system, hence gates with large buoyancy (large eigenfrequency $\unicode[STIX]{x1D714}$ ) are more efficient than heavier gates with large inertia. Moreover, figure 3(a) reveals that the capture factor can be larger than $0.5$ . This latter value corresponds to the maximum of a two-dimensional wave absorber working in synchronous resonance conditions with PTO coefficient equal to the radiation damping (Mei et al. Reference Mei, Stiassnie and Yue2005). Hence, nonlinear resonance has a significant beneficial effect on power extraction.

Figure 3. (a) Maximum of the capture factor ${\mathcal{C}}_{max}^{F}$ versus the eigenfrequency of the array $\unicode[STIX]{x1D714}$ . The capture factor reaches values larger than 0.5, i.e. the maximum that can be reached with synchronous motion only. (b) The optimal value of $\unicode[STIX]{x1D708}_{pto}$ versus eigenfrequency $\unicode[STIX]{x1D714}$ . The power take-off increases with $\unicode[STIX]{x1D714}$ .

We now examine the behaviour of the bandwidth of instability expressed by (4.62) on the eigenfrequency $\unicode[STIX]{x1D714}$ for the optimal values of $\unicode[STIX]{x1D708}_{pto}$ (see figure 4). For large values of the inertia, i.e. small eigenfrequencies $\unicode[STIX]{x1D714}$ , the frequency band of instability approaches zero. Hence heavier gates are more difficult to resonate. Conversely, gates characterized by large values of the buoyancy can be easily resonated because the bandwidth of instability tends to increase. In summary, from a nonlinear point of view maximization of the performance of an array of OWCSs can be reached by choosing gates with large hydrodynamic stiffness.

Figure 4. Bandwidth of instability versus the eigenfrequency of the array $\unicode[STIX]{x1D714}$ . Gates with small inertia $I$ and large buoyancy $C$ increase the eigenfrequency of the mode and render the subharmonic resonance more observable and efficient to extract power from water waves.

Let us fix the eigenfrequency $\unicode[STIX]{x1D714}=1.5~\text{rad}~\text{s}^{-1}$ . The corresponding values of the coefficients of the evolution equation are (see also figure 2):

(4.83a-d ) $$\begin{eqnarray}c_{L}=1.6\times 10^{-4},\quad c_{N}=3.81,\quad c_{R}=0.24,\quad c_{F}=0.91,\end{eqnarray}$$

while the value of $\unicode[STIX]{x1D708}_{pto}$ is chosen following (4.69) (see figure 3 b), i.e.  $\unicode[STIX]{x1D708}_{pto}=423~\text{kg}~\text{m}^{2}~\text{s}^{-1}$ . Figure 5 shows the bifurcation diagram describing the dependence of the unstable and stable equilibria $R^{\pm }$ on $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ . The ellipse is inclined to the left hence the system behaves like a Duffing oscillator with soft spring. Similar behaviour has been already obtained for the subharmonic resonance of Venice gates analysed by Sammarco et al. (Reference Sammarco, Tran and Mei1997a ,Reference Sammarco, Tran, Gottlieb and Mei b ).

Figure 5. Bifurcation diagram of stable/unstable branches for the eigenfrequency of the natural mode $\unicode[STIX]{x1D714}=1.5~\text{rad}~\text{s}^{-1}$ : – – –, unstable branch; ——, stable branch.

Figure 6 shows the behaviour of the capture factor ${\mathcal{C}}^{F}$ given by (4.70) versus the detuning $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ . The maximum capture factor that can be reached is ${\mathcal{C}}^{F}\sim 1$ (see also figure 3 a for $\unicode[STIX]{x1D714}\sim 1.5~\text{rad}~\text{s}^{-1}$ ), i.e. twice the size of the maximum possible for linear synchronous excitation in a channel (Mei et al. Reference Mei, Stiassnie and Yue2005).

Figure 6. Behaviour of the capture factor ${\mathcal{C}}^{F}$ versus detuning $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ .

4.7 Results for modulated incident waves and $Q=2$

In this section we show the behaviour of the array in modulated incident waves with $\overline{A}=0.1$ m. To investigate the effects of the gate stiffness on the minimum of the Melnikov function, consider the geometrical characteristics listed in table 3 and two different values of the eigenfrequency, that is, $\unicode[STIX]{x1D714}_{1}=1.5$ and $\unicode[STIX]{x1D714}_{2}=1.0~\text{rad}~\text{s}^{-1}$ . The respective values of the coefficients are

(4.84) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}c_{L1}=1.6\times 10^{-4},\quad c_{N1}=3.81,\quad c_{R1}=0.24,\quad c_{F1}=0.91,\\ c_{L2}=1.2\times 10^{-4},\quad c_{N2}=1.38,\quad c_{R2}=0.02,\quad c_{F2}=0.17.\end{array}\right\}\end{eqnarray}$$

For each system we choose the PTO coefficient that satisfies (4.69) and maximizes power output. We obtain

(4.85a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{1}=423,\quad \unicode[STIX]{x1D708}_{2}=111~\text{kg}~\text{m}^{2}~\text{s}^{-1}.\end{eqnarray}$$

Once the coefficients are evaluated, we assume $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=0$ and derive the threshold for global chaos that can lead to homoclinic tangle. Figure 7 shows the dependence of the amplitude $\tilde{A}$ on the long-wave frequency $\unicode[STIX]{x1D6FA}$ for both systems. The minimum value shifts to the right with increasing eigenfrequency of the array. Note that the bandwidth for gates with small buoyancy is narrower than that for gates characterized by small inertia. These results demonstrate that with respect to the case of heavier gates, systems with larger eigenfrequencies (more buoyant) can exhibit chaos for a broader range of long-wave forcing frequencies.

Figure 7. Melnikov criterion for determining lower bound of chaos in the plane $\tilde{A}{-}\unicode[STIX]{x1D6FA}$ . Location of the minima depends on the eigenfrequency $\unicode[STIX]{x1D714}$ and tends to zero with the gate inertia. Moreover note that $\tilde{A}\rightarrow \infty$ for $\unicode[STIX]{x1D6FA}\rightarrow (0,\infty )$ and chaotic response becomes impracticable.

Figure 8. Gate array response at $\tilde{A}=0.05~\text{m}$ . Synchronous response. (a) Phase plane $X{-}Y$ and Poincaré map for the evolution of $\overline{\unicode[STIX]{x1D703}}$ ; (b) power spectrum of the gate oscillation.

Figure 9. Subharmonic response at $\tilde{A}=0.14~\text{m}$ .

Figure 10. Chaos at $\tilde{A}=0.2~\text{m}$ .

Now integrate system (4.74) numerically and investigate the occurrence of bifurcation scenarios with increasing amplitude $\tilde{A}$ . Let us assume the eigenfrequency of the array to be $\unicode[STIX]{x1D714}=1.5~\text{rad}~\text{s}^{-1}$ and fix the frequency of the long wave which corresponds to the minimum of the threshold shown in figure 7, i.e. $\unicode[STIX]{x1D6FA}=0.225~\text{rad}~\text{s}^{-1}$ . In order to follow bifurcation patterns we use the Poincaré map sampling time series of the gate envelope response $\overline{\unicode[STIX]{x1D703}}$ at times $t=2\unicode[STIX]{x03C0}n/\unicode[STIX]{x1D6FA}$ , $n=1,2,3\ldots ,$ once transient motion disappears. To better analyse $\overline{\unicode[STIX]{x1D703}}$ we use the definition $\overline{\unicode[STIX]{x1D703}}=X+\text{i}Y$ separating real and imaginary parts. Moreover, the output of (4.74) includes several harmonics, hence we further investigate the time series of the gate response $\unicode[STIX]{x1D6E9}(t)$ looking at its power spectrum (Jordan & Smith Reference Jordan and Smith2011). For $\tilde{A}<0.118$ m the response is synchronous and a stable $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FA}$ periodic solution exists. Figure 8 shows the case for $\tilde{A}=0.05~\text{m}$ . The phase paths of the modal amplitude envelope together with its Poincaré map are shown in figure 8(a), while the power spectrum of the gate response is shown in figure 8(b). At $\tilde{A}=0.118$ m a period doubling occurs and the response of the system becomes subharmonic. This case is shown in figure 9. A further increase of $\tilde{A}$ to values larger than $\tilde{A}=0.151$ m reveals a stable period-four response, while chaos is fully developed for $0.158<\tilde{A}<0.173$  m. In the neighbourhood of $\tilde{A}=0.173$ m a stable period-three response appears. For $\tilde{A}=0.175$ m a chaotic regime reappears with a strange attractor of S shape, as shown in figure 10. This shape follows from the homoclinic path described by (4.79). The forcing term breaks the homoclinic connection and intersections between stable and unstable manifolds, causing chaos to occur. Moreover, the corresponding broad banded power spectrum (figure 10 b) shows that a large number of frequencies is present. Note also the presence of peaks corresponding to the forcing frequency $\unicode[STIX]{x1D714}=1.5~\text{rad}~\text{s}^{-1}$ and its subharmonics. Chaotic motion disappears again when the amplitude of the long wave is further increased up to values $\tilde{A}>0.236~\text{m}$ . Figure 11 shows subharmonic response with frequency downshift of the peak power spectrum for $\tilde{A}=0.25~\text{m}$ (Trulsen & Dysthe Reference Trulsen and Dysthe1997). Indeed there is no peak located at $\unicode[STIX]{x1D714}=1.5~\text{rad}~\text{s}^{-1}$ . Figure 12 shows the bifurcation diagram for the values of $X$ at each time step. The same figure highlights period doubling routes to chaos, periodic windows (in this case around $\tilde{A}=0.173~\text{m}$ ), chaotic bands and periodic responses with frequency downshift. Note also that values of $\tilde{A}$ corresponding to the bifurcation points almost satisfy the universal law found by Feigenbaum (Jordan & Smith Reference Jordan and Smith2011):

(4.86) $$\begin{eqnarray}\lim _{n\rightarrow \infty }\frac{\tilde{A}_{n-1}-\tilde{A}_{n-2}}{\tilde{A}_{n}-\tilde{A}_{n-1}}=\unicode[STIX]{x1D6FF}=4.669\ldots .\end{eqnarray}$$

Recalling that the first bifurcation points are located at $\tilde{A_{1}}=0.118,\tilde{A_{2}}=0.151$ and $\tilde{A_{3}}=0.158$ (see also figure 12), we obtain

(4.87) $$\begin{eqnarray}\frac{\tilde{A}_{2}-\tilde{A}_{1}}{\tilde{A}_{3}-\tilde{A}_{2}}=4.7142\ldots\end{eqnarray}$$

i.e. a value near the universal number $\unicode[STIX]{x1D6FF}$ found by Feigenbaum.

Figure 11. Period-two response and downshift at $\tilde{A}=0.25~\text{m}$ .

Figure 12. Period doubling scenarios as $\tilde{A}$ increases from 0 to 0.35 m.

Finally we evaluate the effects of the chaotic response on the performance of the system. The incident wave amplitude is $\overline{A}+\tilde{A}\cos \unicode[STIX]{x1D6FA}t$ , hence the modulated incident wave potential $\unicode[STIX]{x1D6F7}_{m}^{I}$ is

(4.88) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{m}^{I}=-\frac{\text{i}g(\overline{A}+\tilde{A}\cos \unicode[STIX]{x1D6FA}t)}{4\unicode[STIX]{x1D714}\cosh k_{0}h}\cosh k_{0}(h+z)\text{e}^{-\text{i}(k_{0}x+2\unicode[STIX]{x1D714}t)}+\ast ,\end{eqnarray}$$

and the corresponding averaged energy influx

(4.89) $$\begin{eqnarray}\lim _{\unicode[STIX]{x1D70F}\rightarrow \infty }-\frac{\unicode[STIX]{x1D70C}w}{\unicode[STIX]{x1D70F}}\int _{0}^{\unicode[STIX]{x1D70F}}\int _{-h}^{0}\unicode[STIX]{x1D6F7}_{m_{t}}^{I}\unicode[STIX]{x1D6F7}_{m_{x}}^{I}\,\text{d}z\,\text{d}t=\frac{\unicode[STIX]{x1D6FE}w\unicode[STIX]{x1D714}(2\overline{A}^{2}+\tilde{A})}{4k_{0}}\left(1+\frac{2k_{0}h}{2\sinh 2k_{0}h}\right).\end{eqnarray}$$

The averaged generated power by the array can be evaluated numerically by solving the following integral:

(4.90) $$\begin{eqnarray}P_{m}=\lim _{\unicode[STIX]{x1D70F}\rightarrow \infty }\frac{1}{\unicode[STIX]{x1D70F}}\int _{0}^{\unicode[STIX]{x1D70F}}(\unicode[STIX]{x1D6E9}_{t})^{2}\unicode[STIX]{x1D708}_{pto}\,\text{d}t.\end{eqnarray}$$

Figure 13 shows the capture factor of the array in modulated waves, defined as the ratio between $P_{m}$ (4.90) and the energy influx (4.89) for the same range of $\tilde{A}$ used before. Note that the minimum of the efficiency of the system is located in the range of $\tilde{A}$ where chaos occurs, while the maximum corresponds to $\tilde{A}=0~\text{m}$ . When chaos disappears and regular motion returns with downshift, ${\mathcal{C}}^{F}$ spikes to ${\sim}0.2$ and starts again to decrease. In other words chaos and period doubling is detrimental in terms of wave energy extraction efficiency. We thus obtained a deterministic confirmation of the previous findings of Michele et al. (Reference Michele, Sammarco and d’Errico2016a ,Reference Michele, Sammarco and d’Errico b ) for gate energy production under stochastic incident wave spectra. Finally, note that the value of the capture factor for $\tilde{A}=0$ m is ${\mathcal{C}}^{F}\sim 0.35$ and corresponds to the same value shown in figure 6 for $\unicode[STIX]{x1D714}=3~\text{rad}~\text{s}^{-1}$ and uniform incident waves.

Figure 13. Behaviour of the capture factor of the array in modulated incident waves. Minima correspond to the range of $\tilde{A}$ where chaos occurs while the maximum is located at $\tilde{A}=0~\text{m}$ .

5.1 Leading-order problem $O(1)$

This is the same problem already explained in § 4.1. The solution consists of trapped modes whose related potentials $\unicode[STIX]{x1D719}_{11}$ , $\unicode[STIX]{x1D719}_{12}$ have the form

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{1j}^{\pm }=\mp \text{i}\unicode[STIX]{x1D703}_{j}\mathop{\sum }_{q=1}^{Q}\mathop{\sum }_{m=1}^{\infty }\mathop{\sum }_{n=0}^{\infty }\frac{\unicode[STIX]{x1D714}_{j}b_{mqj}D_{nj}}{C_{n}\unicode[STIX]{x1D6FC}_{mnj}}\text{e}^{\mp \unicode[STIX]{x1D6FC}_{mnj}x}\cos \frac{m\unicode[STIX]{x03C0}y}{b}\cosh k_{nj}(h+z)\equiv \mp \text{i}\unicode[STIX]{x1D703}_{j}f_{1j}^{\pm },\end{eqnarray}$$

where the terms including the subscript $j=1,2$ , are relative to the $j$ th mode with eigenfrequency $\unicode[STIX]{x1D714}_{j}$ and modal form $\unicode[STIX]{x1D703}_{1j}=\{r_{1jq}\}\unicode[STIX]{x1D703}_{j}(t_{2})$ , $q=1,\ldots ,Q$ . To investigate the interactions between two modes it is necessary to consider the number of the gates $Q>2$ because the number of trapped modes is $Q-1$ .

5.2 The second-order problem $O(\unicode[STIX]{x1D716})$

As mentioned before, let us assume at $O(1)$ the existence of a couple of trapped modes $N_{1}$ and $N_{2}$ and the incident wave with frequency $(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})$ . At the order $O(\unicode[STIX]{x1D716})$ we have eight wave components generated by quadratic interactions of the $O(1)$ solution, i.e. the four harmonics (5.1) plus the respective complex conjugates.

We get the governing equation of each problem,

(5.6) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{2j}^{\pm }=0,\quad j=1,\ldots ,4,\end{eqnarray}$$

the relative boundary conditions on the moving gates,

(5.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{21_{x}}^{\pm }=\left[2\text{i}\unicode[STIX]{x1D714}_{1}\unicode[STIX]{x1D703}_{21}(z+h_{p})-\frac{\unicode[STIX]{x1D719}_{1_{z}}^{\pm }\unicode[STIX]{x1D703}_{1}}{\unicode[STIX]{x1D716}}\mp \text{i}\unicode[STIX]{x1D714}_{1}\,d\frac{\unicode[STIX]{x1D703}_{1}^{2}}{\unicode[STIX]{x1D716}}\right]H(z+h-c), & \displaystyle\end{eqnarray}$$

(5.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{22_{x}}^{\pm }=\left[2\text{i}\unicode[STIX]{x1D714}_{2}\unicode[STIX]{x1D703}_{22}(z+h_{p})-\frac{\unicode[STIX]{x1D719}_{2_{z}}^{\pm }\unicode[STIX]{x1D703}_{2}}{\unicode[STIX]{x1D716}}\mp \text{i}\unicode[STIX]{x1D714}_{2}\,d\frac{\unicode[STIX]{x1D703}_{2}^{2}}{\unicode[STIX]{x1D716}}\right]H(z+h-c), & \displaystyle\end{eqnarray}$$

(5.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{23_{x}}^{\pm } & = & \displaystyle \left[\text{i}(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})\unicode[STIX]{x1D703}_{23}(z+h_{p})-\frac{\unicode[STIX]{x1D719}_{2_{z}}^{\pm }\unicode[STIX]{x1D703}_{1}+\unicode[STIX]{x1D719}_{1_{z}}^{\pm }\unicode[STIX]{x1D703}_{2}}{\unicode[STIX]{x1D716}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.\mp \,\text{i}d(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})\frac{\unicode[STIX]{x1D703}_{1}\unicode[STIX]{x1D703}_{2}}{\unicode[STIX]{x1D716}}\right]H(z+h-c),\end{eqnarray}$$

(5.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{24_{x}}^{\pm } & = & \displaystyle \left[\text{i}(\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{1})\unicode[STIX]{x1D703}_{24}(z+h_{p})-\frac{\unicode[STIX]{x1D719}_{2_{z}}^{\pm }\unicode[STIX]{x1D703}_{1}^{\ast }+\unicode[STIX]{x1D719}_{1_{z}}^{\pm ^{\ast }}\unicode[STIX]{x1D703}_{2}}{\unicode[STIX]{x1D716}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.\mp \,\text{i}d(\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{1})\frac{\unicode[STIX]{x1D703}_{1}^{\ast }\unicode[STIX]{x1D703}_{2}}{\unicode[STIX]{x1D716}}\right]H(z+h-c),\quad x=x^{\pm },\end{eqnarray}$$

the mixed boundary conditions on the free surface,

(5.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{21_{z}}^{\pm }=\unicode[STIX]{x1D719}_{21}^{\pm }\frac{4\unicode[STIX]{x1D714}_{1}^{2}}{g}, & \displaystyle\end{eqnarray}$$

(5.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{22_{z}}^{\pm }=\unicode[STIX]{x1D719}_{22}^{\pm }\frac{4\unicode[STIX]{x1D714}_{2}^{2}}{g}, & \displaystyle\end{eqnarray}$$

(5.13) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{23_{z}}^{\pm }=\unicode[STIX]{x1D719}_{23}^{\pm }\frac{(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})^{2}}{g}, & \displaystyle\end{eqnarray}$$

(5.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{24_{z}}^{\pm }=\unicode[STIX]{x1D719}_{24}^{\pm }\frac{(\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{1})^{2}}{g},\quad z=0, & \displaystyle\end{eqnarray}$$

and the corresponding free-surface elevations,

(5.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{21}^{\pm }=\frac{2\text{i}\unicode[STIX]{x1D714}_{1}}{g}\unicode[STIX]{x1D719}_{21}^{\pm }, & \displaystyle\end{eqnarray}$$

(5.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{22}^{\pm }=\frac{2\text{i}\unicode[STIX]{x1D714}_{2}}{g}\unicode[STIX]{x1D719}_{22}^{\pm }, & \displaystyle\end{eqnarray}$$

(5.17) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{23}^{\pm }=\frac{\text{i}(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})}{g}\unicode[STIX]{x1D719}_{23}^{\pm }, & \displaystyle\end{eqnarray}$$

(5.18) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{24}^{\pm }=\frac{\text{i}(\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{1})}{g}\unicode[STIX]{x1D719}_{24}^{\pm },\quad z=0. & \displaystyle\end{eqnarray}$$

As in § 4.2, we obtain the second-order solution of each $j$ th wave field component,

(5.19) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{2j}^{\pm }=\tilde{\unicode[STIX]{x1D703}}_{j}\mathop{\sum }_{p=0}^{\infty }\mathop{\sum }_{l=0}^{\infty }\frac{\unicode[STIX]{x1D6E5}_{plj}}{\unicode[STIX]{x1D6FC}_{plj}\unicode[STIX]{x1D716}}\text{e}^{\pm \text{i}\unicode[STIX]{x1D6FC}_{plj}x}\cos \frac{p\unicode[STIX]{x03C0}y}{b}\cosh k_{lj}(h+z),\quad j=1,\ldots ,4,\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D703}}(t)=\{\unicode[STIX]{x1D703}_{1}^{2},\unicode[STIX]{x1D703}_{2}^{2},\unicode[STIX]{x1D703}_{1}\unicode[STIX]{x1D703}_{2},\unicode[STIX]{x1D703}_{1}\unicode[STIX]{x1D703}_{2}^{\ast }\}$ is the vector that contains the combinations of the resonated mode amplitudes, $k_{lj}$ are the roots of the dispersion relations,

(5.20) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}4\unicode[STIX]{x1D714}_{1}^{2}=gk_{01}\tanh k_{01}h,\quad 4\unicode[STIX]{x1D714}_{1}^{2}=-g\overline{k}_{l1}\tan \overline{k}_{l1}h,\\ 4\unicode[STIX]{x1D714}_{2}^{2}=gk_{02}\tanh k_{02}h,\quad 4\unicode[STIX]{x1D714}_{2}^{2}=-g\overline{k}_{l2}\tan \overline{k}_{l2}h,\\ (\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})^{2}=gk_{03}\tanh k_{03}h,\quad (\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2})^{2}=-g\overline{k}_{l3}\tan \overline{k}_{l3}h,\\ (\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{1})^{2}=gk_{04}\tanh k_{04}h,\quad (\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D714}_{1})^{2}=-g\overline{k}_{l4}\tan \overline{k}_{l4}h,\\ k_{lj}=\text{i}\overline{\unicode[STIX]{x1D705}}_{l},\quad j=1,\ldots ,4,\quad l=1,\ldots ,\infty ,\end{array}\right\}\end{eqnarray}$$

while $\unicode[STIX]{x1D6E5}_{plj}$ is a complex coefficient evaluated in appendix A. Similarly to the previous case, the gate motion at this order is null ( $\unicode[STIX]{x1D703}_{2j}=0$ , $j=1,\ldots ,4$ ), however there are propagating waves that radiate energy to infinity and damp mode motion through cubic interactions at the third order.

5.3 The coupled evolution equations

At this order we apply the solvability condition which is found by Green’s formula:

(5.21) $$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \frac{\unicode[STIX]{x1D714}_{j}^{2}}{g}\int _{0}^{b}\,\text{d}y\left(\int _{-\infty }^{0}f_{1j}^{-}{\mathcal{F}}_{3j}^{-}\,\text{d}x-\int _{0}^{\infty }f_{1j}^{+}{\mathcal{F}}_{3j}^{+}\,\text{d}x\right)\nonumber\\ \displaystyle & & \displaystyle +\,2\int _{0}^{b}\,\text{d}y\int _{-h}^{0}\overline{f_{j}{\mathcal{G}}_{3j}}\,\text{d}z+\mathop{\sum }_{q=1}^{Q}\left(\frac{\text{i}{\mathcal{D}}_{3j}r_{1jq}}{\unicode[STIX]{x1D70C}}-\frac{\unicode[STIX]{x1D708}_{pto}r_{1jq}^{2}\unicode[STIX]{x1D714}_{j}\unicode[STIX]{x1D703}_{j}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D716}^{2}}\right),\quad j=1,2.\end{eqnarray}$$

Details of the forcing terms in the latter equation can be found in appendix B. Carrying out the integrals of (5.21) and grouping the terms according to $\unicode[STIX]{x1D703}_{11}$ , $\unicode[STIX]{x1D703}_{12}$ , $\unicode[STIX]{x1D703}_{11}\left|\unicode[STIX]{x1D703}_{11}\right|^{2}$ , $\unicode[STIX]{x1D703}_{11}\left|\unicode[STIX]{x1D703}_{12}\right|^{2}$ , $\unicode[STIX]{x1D703}_{12}\left|\unicode[STIX]{x1D703}_{12}\right|^{2}$ , $\unicode[STIX]{x1D703}_{12}\left|\unicode[STIX]{x1D703}_{11}\right|^{2}$ , $\unicode[STIX]{x1D703}_{11}^{\ast }$ , $\unicode[STIX]{x1D703}_{12}^{\ast }$ , we obtain two coupled evolution equations for both modal oscillations $\unicode[STIX]{x1D703}_{11}$ and $\unicode[STIX]{x1D703}_{12}$ :

(5.22) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}\unicode[STIX]{x1D703}_{11_{t}}=\unicode[STIX]{x1D703}_{11}|\unicode[STIX]{x1D703}_{11}|^{2}(c_{N}+\text{i}c_{R})+\unicode[STIX]{x1D703}_{11}|\unicode[STIX]{x1D703}_{12}|^{2}(c_{S}+\text{i}c_{U})+A\unicode[STIX]{x1D703}_{12}^{\ast }c_{F}+\text{i}\unicode[STIX]{x1D703}_{11}\unicode[STIX]{x1D708}_{pto}c_{L}, & \displaystyle\end{eqnarray}$$

(5.23) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}\unicode[STIX]{x1D703}_{12_{t}}=\unicode[STIX]{x1D703}_{12}|\unicode[STIX]{x1D703}_{12}|^{2}(d_{N}+\text{i}d_{R})+\unicode[STIX]{x1D703}_{12}|\unicode[STIX]{x1D703}_{11}|^{2}(d_{S}+\text{i}d_{U})+A\unicode[STIX]{x1D703}_{11}^{\ast }d_{F}+\text{i}\unicode[STIX]{x1D703}_{12}\unicode[STIX]{x1D708}_{pto}d_{L}, & \displaystyle\end{eqnarray}$$

where the coupling coefficients are evaluated in a similar manner as explained in § 4.3. Differently to the case of single-mode subharmonic excitation we have two new pairs of coefficients as in the edge wave theories, $c_{S},d_{S},c_{U},d_{U}$ , that represent the coupling between the modes. Note also that in (5.22) the incident wave amplitude $A$ is multiplied by the complex conjugate of $\unicode[STIX]{x1D703}_{12}$ , while in (5.23) the same term is multiplied by $\unicode[STIX]{x1D703}_{11}^{\ast }$ . This depends on the fact that the incident wave frequency is assumed to be equal to $\left(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2}\right)$ and not equal to twice the eigenfrequency of a single mode. Now consider the detuning of the incident wave frequency $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ with $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}\sim O(\unicode[STIX]{x1D716}^{2})$ and let us assume the following change of variables:

(5.24) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{j}=\overline{\unicode[STIX]{x1D703}}_{j}\text{e}^{-\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{j}t},\quad j=1,2,\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{2}$ . Then the evolution equations modify as follows:

(5.25) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}\overline{\unicode[STIX]{x1D703}}_{11_{t}}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{1}\overline{\unicode[STIX]{x1D703}}_{11}+\overline{\unicode[STIX]{x1D703}}_{11}|\overline{\unicode[STIX]{x1D703}}_{11}|^{2}(c_{N}+\text{i}c_{R})+\overline{\unicode[STIX]{x1D703}}_{11}|\overline{\unicode[STIX]{x1D703}}_{12}|^{2}(c_{S}+\text{i}c_{U})+A\overline{\unicode[STIX]{x1D703}}_{12}^{\ast }c_{F}+\text{i}\overline{\unicode[STIX]{x1D703}}_{11}\unicode[STIX]{x1D708}_{pto}c_{L}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(5.26) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}\overline{\unicode[STIX]{x1D703}}_{12_{t}}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{2}\overline{\unicode[STIX]{x1D703}}_{12}+\overline{\unicode[STIX]{x1D703}}_{12}|\overline{\unicode[STIX]{x1D703}}_{12}|^{2}(d_{N}+\text{i}d_{R})+\overline{\unicode[STIX]{x1D703}}_{12}|\overline{\unicode[STIX]{x1D703}}_{11}|^{2}(d_{S}+\text{i}d_{U})+A\overline{\unicode[STIX]{x1D703}}_{11}^{\ast }d_{F}+\text{i}\overline{\unicode[STIX]{x1D703}}_{12}\unicode[STIX]{x1D708}_{pto}d_{L}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

5.4 Analysis of the coupled equations and theoretical results

Consider an array of $Q=4$ gates with characteristics listed in table 3. We obtain two odd natural modes and a single even mode (see also Sammarco et al. Reference Sammarco, Michele and d’Errico2013). Interaction between odd–even modes yields the coefficients $c_{F},d_{F}$ equal to zero and the system of (5.25)–(5.26) becomes unforced. As a consequence, we can consider the interaction between the odd modes only, with shape $r_{11q}=\{1,0.41,-0.41,-1\}$ and $r_{12q}=\{1,-2.41,2.41,-1\}$ . Sample values of the coefficients of (5.25)–(5.26) are shown in figure 14. Note that $c_{U}$ and $d_{U}$ (figure 14 e) are very small if compared to the others and tend to increase for large eigenfrequencies while $c_{F}$ is larger than $d_{F}$ although $\unicode[STIX]{x1D714}_{1}<\unicode[STIX]{x1D714}_{2}$ for fixed gate characteristics. This has an important consequence because $N_{1}$ can be forced more easily by the incident waves.

Figure 14. Behaviour of the coefficients of the evolution equations versus the eigenfrequency of the mode $N_{2}$ , $\unicode[STIX]{x1D714}_{2}$ .

Now define $\overline{\unicode[STIX]{x1D703}}_{1j}$ in action-angle variables form, i.e.  $\overline{\unicode[STIX]{x1D703}}_{1j}=\text{i}\sqrt{R_{j}}\text{e}^{\text{i}\unicode[STIX]{x1D713}_{j}}$ , and neglect $c_{U},d_{U}$ with respect to the other coefficients to reduce the algebra without simplifying the physics. Then, from (5.25)–(5.26), we obtain a system of four real differential equations:

(5.27) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle R_{1_{t}}=-2R_{1}\left[R_{1}c_{R}+Ac_{F}\sqrt{\frac{R_{2}}{R_{1}}}\sin \unicode[STIX]{x1D713}+\unicode[STIX]{x1D708}_{pto}c_{L}\right]\\ \displaystyle R_{2_{t}}=-2R_{2}\left[R_{2}d_{R}+Ad_{F}\sqrt{\frac{R_{1}}{R_{2}}}\sin \unicode[STIX]{x1D713}+\unicode[STIX]{x1D708}_{pto}d_{L}\right]\\ \displaystyle \unicode[STIX]{x1D713}_{1_{t}}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{1}-Ac_{F}\sqrt{\frac{R_{2}}{R_{1}}}\cos \unicode[STIX]{x1D713}+R_{1}c_{N}+R_{2}c_{S}\\ \displaystyle \unicode[STIX]{x1D713}_{2_{t}}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{2}-Ad_{F}\sqrt{\frac{R_{1}}{R_{2}}}\cos \unicode[STIX]{x1D713}+R_{1}d_{S}+R_{2}d_{N},\end{array}\right\}\end{eqnarray}$$

in which $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}_{1}+\unicode[STIX]{x1D713}_{2}$ denotes the sum of the modal phases. After a long time the system reaches equilibrium, hence we now focus our attention on determining the fixed points of (5.27). The trivial fixed point is at

(5.28a,b ) $$\begin{eqnarray}R_{1}=R_{2}=0,\quad \unicode[STIX]{x1D713}^{\ast }=\cos ^{-1}\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}}{A(c_{F}+d_{F})},\end{eqnarray}$$

while non-trivial fixed points related to unstable and stable equilibria correspond to the roots of the system

(5.29) $$\begin{eqnarray}R_{2}=\frac{1}{2d_{R}}\left[-\unicode[STIX]{x1D708}_{pto}d_{L}+\sqrt{\unicode[STIX]{x1D708}_{pto}^{2}d_{L}^{2}+\frac{4d_{F}d_{R}R_{1}(\unicode[STIX]{x1D708}_{pto}c_{L}+c_{R}R_{1})}{c_{F}}}\right],\end{eqnarray}$$

(5.30) $$\begin{eqnarray}\displaystyle & & \displaystyle [\unicode[STIX]{x0394}\unicode[STIX]{x1D714}+R_{1}(c_{N}+d_{S})+R_{2}(d_{N}+c_{S})]^{2}+(R_{1}c_{R}+\unicode[STIX]{x1D708}_{pto}c_{L}+R_{2}d_{R}+\unicode[STIX]{x1D708}_{pto}d_{L})^{2}\nonumber\\ \displaystyle & & \displaystyle \quad -\,A^{2}\left(\sqrt{\frac{R_{2}}{R_{1}}}c_{F}+\sqrt{\frac{R_{1}}{R_{2}}}d_{F}\right)^{2}=0,\end{eqnarray}$$

which can be found numerically. Note that if $R_{1}=0$ also $R_{2}=0$ and triad resonance does not occur. Stability analysis of the linearized system (5.25)–(5.26) at the origin allows us to find the threshold of resonance triad,

(5.31) $$\begin{eqnarray}A>\unicode[STIX]{x1D708}_{pto}\sqrt{\frac{c_{L}d_{L}}{c_{F}d_{F}}},\end{eqnarray}$$

hence the amplitude of the incident wave must be larger than the latter quantity to trigger mode–mode interactions. Once the equilibrium states $R_{1}$ and $R_{2}$ are evaluated we determine the sum of the modal phases by manipulating $\unicode[STIX]{x1D713}_{1_{t}}$ and $\unicode[STIX]{x1D713}_{2_{t}}$ :

(5.32) $$\begin{eqnarray}\unicode[STIX]{x1D713}^{0}=\cos ^{-1}\left[\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D714}+R_{1}(c_{N}+d_{S})+R_{2}(c_{S}+d_{N})}{A\displaystyle \left(c_{F}\sqrt{\frac{R_{2}}{R_{1}}}+d_{F}\sqrt{\frac{R_{1}}{R_{2}}}\right)}\right].\end{eqnarray}$$

Substitution of $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}^{0}$ inside (5.27) gives $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{1}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{2}$ . Physically, the latter terms are the speed along the limit cycles in the complex plane $\text{Re}\{\unicode[STIX]{x1D703}\},\text{Im}\{\unicode[STIX]{x1D703}\}$ . The generated power due to triad resonance is finally given by

(5.33) $$\begin{eqnarray}\displaystyle P_{t} & = & \displaystyle \lim _{\unicode[STIX]{x1D70F}\rightarrow \infty }\frac{1}{\unicode[STIX]{x1D70F}}\int _{0}^{\unicode[STIX]{x1D70F}}[(\unicode[STIX]{x1D703}_{11}\text{e}^{-\text{i}(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{1})t}+\unicode[STIX]{x1D703}_{12}\text{e}^{-\text{i}(\unicode[STIX]{x1D714}_{2}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{2})t}+\ast )_{t}]^{2}\unicode[STIX]{x1D708}_{pto}\nonumber\\ \displaystyle & = & \displaystyle 2\unicode[STIX]{x1D708}_{pto}\left[(\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{1})^{2}\mathop{\sum }_{q=1}^{Q}r_{11q}^{2}R_{1}^{0}+(\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{2}+\unicode[STIX]{x1D714}_{2})^{2}\mathop{\sum }_{q=1}^{Q}r_{12q}^{2}R_{2}^{0}\right].\end{eqnarray}$$

Now focus the attention on a fixed OWSC configuration and analyse the effects of triad resonance on the performances of the array. Let us consider the case where the system has the eigenfrequency $\unicode[STIX]{x1D714}_{2}=1.5~\text{rad}~\text{s}^{-1}$ . Solution of the eigenvalue condition at the first order gives the eigenfrequency of the mode $N_{1}$ , i.e.  $\unicode[STIX]{x1D714}_{1}=0.9~\text{rad}~\text{s}^{-1}$ . The corresponding evolution equation coefficients are (see also figure 14):

(5.34) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}c_{N}=1.84,\quad c_{R}=0.1,\quad c_{F}=1.34,\quad c_{S}=5.68,\\ d_{N}=12.96,\quad d_{R}=0.06,\quad d_{F}=0.11,\quad d_{S}=7.93.\end{array}\right\}\end{eqnarray}$$

Choose the $\unicode[STIX]{x1D708}_{pto}$ which maximizes power output for perfect resonance, i.e.  $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{1}=\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{2}=0$ and $A=0.05$ m. Maximization of (5.33) yields $\unicode[STIX]{x1D708}_{pto}=120~\text{kg}~\text{m}^{2}~\text{s}^{-1}$ . The corresponding equilibrium branches of $R_{1}$ and $R_{2}$ are plotted in figure 15. Simple numerical evaluation of the eigenvalues reveals results similar to the case of single-mode resonance. The continuous lines correspond to stable fixed points, while the dot lines are related to unstable saddles.

Figure 15. Equilibrium branches of $R_{1}$ ad $R_{2}$ versus detuning of the incident wave $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ . The solid lines represent the stable equilibrium branch while the dotted lines represent the unstable branch.

At this point, we evaluate the efficiency of the system excited through mode-mode interactions. Figure 16 shows the behaviour of ${\mathcal{C}}^{F}$ versus the frequency of the incident waves. The maximum of the capture factor is ${\sim}0.25$ , hence the effect of triad interaction is not as significant as the subharmonic case (see figure 6) but positive anyway.

Figure 16. Behaviour of the capture factor due to triad resonance. Mode–mode interaction yields smaller values of ${\mathcal{C}}^{F}$ than the subharmonic resonance of a single mode.