2.1. Governing equation and boundary conditions

We consider a thin elastic plate floating atop water. As illustrated in figure 1, the upper surface of the plate is periodically attached with resonators, forming the so-called floating metaplate. The plate is always in contact with the water surface and is thin enough such that Kirchhoff’s thin plate theory is valid. We consider a Cartesian coordinate system with the $x$ -axis coinciding with the water surface and the positive $z$ -axis pointing upward. We assume that the fluid is incompressible, homogeneous and inviscid, and the flow is irrotational. Under these assumptions, the linear potential theory is adopted. The velocity potential $\varPhi (x,t)$ takes the form of the following Laplace equation (Fox & Squire Reference Fox and Squire1994):

(2.1) \begin{align} {\nabla} ^2\varPhi =0\quad -H\leqslant z \leqslant 0, \end{align}

where $H$ is the depth of water. The seabed (no flow) boundary condition is given by

(2.2) \begin{equation} \frac {\partial \varPhi }{\partial z}=0\quad z=-H, -\infty \lt x\lt \infty . \end{equation}

The dynamic condition for the free surface reads

(2.3) \begin{equation} \frac {\partial \varPhi }{\partial t}+gw=0\quad z=0, \end{equation}

while that on the surface covered with the thin plate is

(2.4) \begin{equation} D\frac {\partial ^4w}{\partial x^4}+\rho _p h\frac {\partial ^2w}{\partial t^2}=p\quad z=0, \end{equation}

where $g$ is the acceleration due to gravity, $w(x,t)$ is the displacement of the fluid surface (which is also the vertical displacement of the plate), $D=Eh^3/12(1-\nu ^2)$ represents the flexural rigidity, with $E$ , $\nu$ , $h$ and $\rho _p$ being the Young’s modulus, Poisson’s ratio, thickness and density of the plate, respectively. The surface pressure $p$ in (2.4) satisfies

(2.5) \begin{equation} \rho \frac {\partial \varPhi }{\partial t}+\rho gw+p=0, \end{equation}

where $\rho$ is the density of water. It should be pointed out that the effect of the resonator is embodied in boundary conditions to be described below, but not in the governing equation (2.4). In addition, the kinematic boundary condition reads

(2.6) \begin{equation} \frac {\partial \varPhi }{\partial z}=\frac {\partial w}{\partial t}\quad z=0. \end{equation}

Figure 1. Schematic diagram of a floating plate, with periodic resonators attached on its surface, forming the so-called floating metaplate. The plate is assumed thin and elastic. One unit cell spans from $x=0$ to $x=a$ , in which the resonator is attached at $x=b$ . The velocity potential of the water wave and the deflection of the plate are labelled by $\phi$ and $\eta$ , respectively. The inset depicts the force diagram related to the interaction between the resonator and the plate.

Assuming that all motions are time-harmonic with angular frequency $\omega$ , the velocity potential of water waves  $\varPhi$ and the displacement of the plate $w$ can be expressed as

(2.7a) \begin{align}\varPhi (x,z,t)&={\rm{Re}}\big\{\phi (x,z)\textrm{e}^{-\textrm{i}\omega t}\big\},\end{align}

(2.7b) \begin{align}w(x,t)&={\rm{Re}}\big\{\eta (x)\textrm{e}^{-\textrm{i}\omega t}\big\},\end{align}

where the reduced velocity potential $\phi$ and plate deflection $\eta$ are both complex valued.

From (2.1)–(2.7), we can get the following boundary value problem in the frequency domain:

(2.8a) \begin{align}\frac {\partial ^2\phi }{\partial x^2}+\frac {\partial ^2\phi }{\partial z^2}&=0\quad -H\leqslant z\lt 0,\end{align}

(2.8b) \begin{align}\frac {\partial \phi }{\partial z}&=0\quad z=-H,\end{align}

together with the boundary conditions at $z=0$ , for the free surface, written as

(2.9) \begin{equation} \frac {\partial \phi }{\partial z}=\alpha \phi \quad z=0, \end{equation}

and for the surface covered with the plate as

(2.10) \begin{equation} \left(\beta \frac {\partial ^4}{\partial x^4}+1-\alpha \gamma _0\right)\frac {\partial \phi }{\partial z}=\alpha \phi \quad z=0, \end{equation}

where $\alpha$ , $\beta$ and $\gamma _0$ are

(2.11) \begin{equation} \alpha =\frac {\omega ^2}{g},\,\beta =\frac {D}{\rho g},\,\gamma _0=\frac {\rho _p h}{\rho }, \end{equation}

respectively. Note that $\beta$ and $\gamma _0$ , which represent the bending stiffness term and inertia term, are two key parameters to characterise the physical properties of the plate.

2.2. Eigenfunction matching method and Bloch’s theorem

We derive the solution by the eigenfunction matching method (Fox & Squire Reference Fox and Squire1994). Using separation of variables, we can write the general solution to the velocity potential $\phi (x,z)$ as

(2.12) \begin{equation} \phi (x,z)=\sum _{n=0}^{N}X_n(x)\varphi _n(z)\quad \text{under the free surface}, \end{equation}

or

(2.13) \begin{equation} \phi (x,z)=\sum _{m=-2}^{M}Z_m(x)\psi _m(z)\quad \text{under the plate}, \end{equation}

where $\varphi _n(z)$ and $\psi _m(z)$ correspond to the normalised vertical eigenfunction of the potential under the free surface and the plate region, respectively, and can be expressed as

(2.14) \begin{align}\varphi _n(z)=\frac {\cos \vartheta _n(z+H)}{\cos \vartheta _n H}\quad n\geqslant 0,\end{align}

(2.15) \begin{align}\psi _m(z)=\frac {\cos \kappa _m(z+H)}{\cos \kappa _m H}\quad m\geqslant -2,\end{align}

in which the separation constant $\vartheta _n$ and $\kappa _m$ are solutions of the following two dispersion equations, respectively:

(2.16) \begin{align}\vartheta \tan (\vartheta H)=-\alpha \quad \text{under the free surface},\end{align}

(2.17) \begin{align}\kappa \tan {\kappa H}=\frac {-\alpha }{\beta \kappa ^4+1-\gamma _0\alpha }\quad \text{under the plate}.\end{align}

Note that the solutions of (2.16) are $\vartheta _0$ , which is a purely negative imaginary number, and $\vartheta _n$ ( $n\geqslant 1$ ) which is a positive real number. While the solutions of (2.17) are $\kappa _m$ ( $m\geqslant -2$ ), in which $\kappa _{-2}$ and $\kappa _{-1}$ are complex with positive real parts, $\kappa _0$ is purely negative imaginary, and $\kappa _m$ ( $m\geqslant 1$ ) is purely positive real. Then, we can write $X_n(x)$ under the free surface in (2.12) and $Z_m(x)$ under the plate in (2.13) as

(2.18) \begin{equation} X_n(x)=l_n\textrm{e}^{\vartheta _n x}+r_n\textrm{e}^{-\vartheta _n x} , \end{equation}

and

(2.19) \begin{equation} Z_m(x)=a_m\textrm{e}^{\kappa _m x}+b_m\textrm{e}^{-\kappa _m x}, \end{equation}

respectively, where $l_n, r_n$ ( $n=0,1,2,\ldots ,N$ ) and $a_m, b_m$ ( $m=-2,-1,0,\ldots ,M$ ) are undetermined coefficients. It is worth noting that, while there are an infinite number of vertical modes (i.e. $\varphi (z)$ and $\psi (z)$ ), we truncate (2.12) and (2.13) at the $N$ th and $M$ th modes for numerical solutions.

From the kinematic boundary condition (2.6) and using (2.10), (2.13), (2.15) and (2.19), the deflection of the plate can be obtained as

(2.20) \begin{equation} \eta (x)=\frac {1}{-\textrm{i}\omega }\frac {\partial \phi }{\partial z}\Big |_{z=0}=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\left(a_m\textrm{e}^{\kappa _m x}+b_m\textrm{e}^{-\kappa _m x}\right)\!. \end{equation}

Since the resonators are periodically arranged on the floating plate with periodicity $a$ in the $x$ direction, it is possible to confine attention to one unit cell ranging from $x=0$ to $x=a$ that contains only a single resonator, as shown in figure 2(a). Bloch’s theorem states that, under a discrete lattice translation $na$ ( $n\in \mathbb{Z}$ ), any physical field, e.g. the velocity potential $\phi (x,z)$ or the plate deflection $\eta (x)$ , can be obtained from the one within this unit cell being modulated by a plane wave (see Chou Reference Chou1998), i.e.

(2.21) \begin{eqnarray} \phi (x+na,z)=\phi (x,z)\textrm{e}^{\textrm{i}kna},\, \eta (x+na)=\eta (x)\textrm{e}^{\textrm{i}kna}, \end{eqnarray}

where $k$ is the Bloch wavenumber.

We segment the unit cell at the point $x=b$ where the resonator is located, and express the velocity potential underneath the two segments (see figure 1) as

(2.22a) \begin{align}\phi _1(x,z)&=\sum _{m=-2}^{M}\left(a_m\textrm{e}^{\kappa _m x}+b_m\textrm{e}^{-\kappa _m x}\right)\psi _m(z)\quad 0\leqslant x\leqslant b,\end{align}

(2.22b) \begin{align}\phi _2(x,z)&=\sum _{m=-2}^{M}\left(c_m\textrm{e}^{\kappa _m x}+d_m\textrm{e}^{-\kappa _m x}\right)\psi _m(z)\quad b\leqslant x\leqslant a,\end{align}

where $\kappa _m$ is the $m$ th root of (2.17) and $a_m$ , $b_m$ , $c_m$ and $d_m$ are undetermined coefficients. According to (2.21), we can write out the velocity potential that belongs to the left (i.e. $x\lt 0$ ) and right (i.e. $x\gt a$ ) segments adjacent to this unit cell as

(2.23a) \begin{align}\phi _l(x-a,z)&=\phi _2(x,z)\textrm{e}^{-\textrm{i}ka}\quad b\leqslant x\leqslant a,\end{align}

(2.23b) \begin{align}\phi _r(x+a,z)&=\phi _1(x,z)\textrm{e}^{\textrm{i}ka}\quad 0\leqslant x\leqslant b. \\[9pt] \nonumber\end{align}

From (2.20) and applying Bloch’s theorem again, the deflection of the plate that belongs to distinct segments (see figure 1) can be expressed as

(2.24a) \begin{align}\eta _1(x)&=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\left(a_m\textrm{e}^{\kappa _m x}+b_m\textrm{e}^{-\kappa _m x}\right)\quad 0\leqslant x\leqslant b,\end{align}

(2.24b) \begin{align}\eta _2(x)&=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\left(c_m\textrm{e}^{\kappa _m x}+d_m\textrm{e}^{-\kappa _m x}\right)\quad b\leqslant x\leqslant a,\end{align}

(2.24c) \begin{align}&\qquad\qquad\eta _l(x-a)=\eta _2(x)\textrm{e}^{-\textrm{i}ka}\quad b\leqslant x\leqslant a,\end{align}

(2.24d) \begin{align}&\qquad\qquad\,\eta _r(x+a)=\eta _1(x)\textrm{e}^{\textrm{i}ka}\quad 0\leqslant x\leqslant b.\end{align}

Identical spatial relations hold for other physical quantities such as the horizontal velocity of the fluid $\partial \phi /\partial x$ , rotation angle $\partial \eta /\partial x$ , bending moment $D\partial ^2\eta /\partial x^2$ and shear force $-D\partial ^3\eta /\partial x^3$ of the plate (Iida & Umazume Reference Iida and Umazume2020). It is worth noting that the wavenumbers $k$ in (2.23) and (2.24) need only be considered between $-\pi /a \leqslant k \leqslant \pi /a$ , which is the first Brillouin zone (BZ), as solutions with all other real values $k$ can be obtained by adding or subtracting a multiple of $2\pi /a$ .

2.3. Generalised eigenvalue problem

We seek non-trivial solutions of the unknown coefficients mentioned in the previous section by considering the matching conditions associated with the potential and the boundary conditions of the plate. Specifically, we solve for values of Bloch wavenumbers $k$ for a given angular frequency $\omega$ .

For the unit cell we selected, we consider the matching conditions of the velocity potential and its normal derivative across the boundary between distinct regions separated by $x=0$ , $x=b$ and $x=a$ . This gives

(2.25) \begin{align}&\phi _1(b,z)=\phi _2(b,z), \text{i.e.,}\nonumber \\ &\sum _{m=-2}^{M}\big(a_m\textrm{e}^{\kappa _m b}+b_m\textrm{e}^{-\kappa _m b}\big)\psi _m(z)=\sum _{m=-2}^{M}\big(c_m\textrm{e}^{\kappa _m b}+d_m\textrm{e}^{-\kappa _m b}\big)\psi _m(z),\end{align}

(2.26) \begin{align}&\quad\phi _l(0,z)=\phi _1(0,z), \text{i.e.,}\nonumber \\ &\quad\sum _{m=-2}^{M}\left(c_m\textrm{e}^{\kappa _m a}+d_m\textrm{e}^{-\kappa _m a}\right)\psi _m(z)\textrm{e}^{-\textrm{i}ka}=\sum _{m=-2}^M(a_m+b_m)\psi _m(z),\end{align}

and

(2.27) \begin{align}&\frac {\partial \phi _1(x,z)}{\partial x}\Big |_{x=b}=\frac {\partial \phi _2(x,z)}{\partial x}\Big |_{x=b}, \text{i.e.,}\nonumber \\ &\sum _{m=-2}^{M}\kappa _m\big(a_m\textrm{e}^{\kappa _m b}-b_m\textrm{e}^{-\kappa _m b}\big)\psi _m(z)=\sum _{m=-2}^{M}\kappa _m\big(c_m\textrm{e}^{\kappa _m b}-d_m\textrm{e}^{-\kappa _m b}\big)\psi _m(z),\end{align}

(2.28) \begin{align}&\quad\frac {\partial \phi _l(x,z)}{\partial x}\Big |_{x=0}=\frac {\partial \phi _1(x,z)}{\partial x}\Big |_{x=0}, \text{i.e.,}\nonumber \\ &\quad\sum _{m=-2}^{M}\kappa _m\left(c_m\textrm{e}^{\kappa _m a}-d_m\textrm{e}^{-\kappa _m a}\right)\psi _m(z)\textrm{e}^{-\textrm{i}ka}=\sum _{m=-2}^M\kappa _m(a_m-b_m)\psi _m(z).\end{align}

It is worth noting that the matching conditions across the right boundary of the unit cell (i.e. $x=a$ ) are exactly equivalent to (2.26) and (2.28). Since these relations must hold at every depth $z$ , we solve these equations by multiplying both sides of (2.25), (2.26), (2.27) and (2.28) by the vertical eigenfunction $\varphi _n(z)$ (see (2.14)) and integrating from $z=-H$ to $z=0$ to obtain

(2.29) \begin{align}\sum _{m=-2}^{M}\big(a_m\textrm{e}^{\kappa _m b}+b_m\textrm{e}^{-\kappa _m b}\big)P_{\textit{nm}}=\sum _{m=-2}^{M}\big(c_m\textrm{e}^{\kappa _m b}+d_m\textrm{e}^{-\kappa _m b}\big)P_{\textit{nm}},\end{align}

(2.30) \begin{align}&\quad \sum _{m=-2}^{M}\left(c_m\textrm{e}^{\kappa _m a}+d_m\textrm{e}^{-\kappa _m a}\right)P_{\textit{nm}}\textrm{e}^{-\textrm{i}ka}=\sum _{m=-2}^M(a_m+b_m)P_{\textit{nm}},\end{align}

and

(2.31) \begin{align}& \sum _{m=-2}^{M}\kappa _m\big(a_m\textrm{e}^{\kappa _m b}-b_m\textrm{e}^{-\kappa _m b}\big)P_{\textit{nm}}=\sum _{m=-2}^{M}\kappa _m\big(c_m\textrm{e}^{\kappa _m b}-d_m\textrm{e}^{-\kappa _m b}\big)P_{\textit{nm}},\end{align}

(2.32) \begin{align}&\quad \sum _{m=-2}^{M}\kappa _m\left(c_m\textrm{e}^{\kappa _m a}-d_m\textrm{e}^{-\kappa _m a}\right)P_{\textit{nm}}\textrm{e}^{-\textrm{i}ka}=\sum _{m=-2}^M\kappa _m(a_m-b_m)P_{\textit{nm}},\end{align}

where $0\leqslant n\leqslant N$ , and $P_{\textit{nm}}$ is the inner product of $\varphi _n(z)$ and $\psi _m(z)$ (Meylan Reference Meylan2019), i.e.

(2.33) \begin{eqnarray} P_{\textit{nm}}&=&\int _{-H}^0\varphi _n(z)\psi _m(z)\textrm{d}z\nonumber \\ &=&\frac {\vartheta _n\sin \vartheta _nH\cos \kappa _mH-\kappa _m\cos \vartheta _nH\sin \kappa _mH}{\cos \vartheta _nH\cos \kappa _mH(\vartheta _n^2-\kappa _m^2)}. \end{eqnarray}

Note that the total number of equations contained in (2.29)–(2.32) is $4(N+1)$ . It is also worth noting that, by taking advantage of the completeness, the vertical eigenfunctions $\varphi _n(z)$ of the free surface waves, instead of those $\psi _m(z)$ underneath the elastic plate, are selected for the projection operations (Kohout et al. Reference Kohout, Meylan, Sakai, Hanai, Leman and Brossard2007; Meng & Lu Reference Meng and Lu2017).

Then we consider the matching and boundary conditions of the plate. At $x=0$ , the deflection, the rotation angle, the bending moment and the shear force are continuous, i.e.

(2.34) \begin{align}\eta _l(0)&=\eta _1(0), \text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\left(c_m\textrm{e}^{\kappa _m a}+d_m\textrm{e}^{-\kappa _m a}\right)\textrm{e}^{-\textrm{i}ka}\nonumber \\ &=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }(a_m+b_m),\end{align}

(2.35) \begin{align}\frac {\partial \eta _l}{\partial x}\Big |_{x=0}&=\frac {\partial \eta _1}{\partial x}\Big |_{x=0}, \text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m\left(c_m\textrm{e}^{\kappa _m a}-d_m\textrm{e}^{-\kappa _m a}\right)\textrm{e}^{-\textrm{i}ka}\nonumber \\ &=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m(a_m-b_m),\end{align}

(2.36) \begin{align}D\frac {\partial ^2\eta _l}{\partial x^2}\Big |_{x=0}&=D\frac {\partial ^2\eta _1}{\partial x^2}\Big |_{x=0}, \text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m^2\left(c_m\textrm{e}^{\kappa _m a}+d_m\textrm{e}^{-\kappa _m a}\right)\textrm{e}^{-\textrm{i}ka}\nonumber \\ &=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m^2(a_m+b_m),\end{align}

(2.37) \begin{align}-D\frac {\partial ^3\eta _l}{\partial x^3}\Big |_{x=0}&=-D\frac {\partial ^3\eta _1}{\partial x^3}\Big |_{x=0}, \text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m^3\left(c_m\textrm{e}^{\kappa _m a}-d_m\textrm{e}^{-\kappa _m a}\right)\textrm{e}^{-\textrm{i}ka}\nonumber \\ &=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m^3(a_m-b_m).\end{align}

Note that the continuity conditions at $x=a$ are equivalent to those at $x=0$ . Finally, at the point where the resonator is located ( $x=b$ ), the deflection, the rotation angle and the bending moment are continuous, but there is a jump of $F$ in the shear force, where $F$ comes from the dynamic force in the spring of the resonator (see the inset in figure 1) and can be expressed as

(2.38) \begin{equation} F=s_r(\eta |_{x=b}-\varXi )=-m_r\omega ^2\varXi , \end{equation}

where the assumption $\xi =\varXi \textrm{e}^{-\textrm{i}\omega t}$ was used, and $s_r$ and $m_r$ are the stiffness of the spring and the mass of the resonator, respectively. Note that the rightmost term in (2.38) arises from applying Newton’s second law to the resonator. Equation (2.38) implies $\varXi =s_r\eta |_{x=b}/(s_r-m_r\omega ^2)$ , so we have

(2.39) \begin{equation} F=s_r\eta |_{x=b}\frac {-\tilde {\omega }^2}{1-\tilde {\omega }^2}, \end{equation}

where $\tilde {\omega }$ is the normalised angular frequency by the natural frequency of the resonator $\omega _0$ , defined as

(2.40) \begin{equation} \tilde {\omega }=\frac {\omega }{\omega _0},\,\omega _0=\sqrt {\frac {s_r}{m_r}}. \end{equation}

Consequently, the matching conditions at $x=b$ read

(2.41) \begin{align}\eta _1(b)&=\eta _2(b), \text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\big(a_m\textrm{e}^{\kappa _m b}+b_m\textrm{e}^{-\kappa _m b}\big)\nonumber \\ &=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\big(c_m\textrm{e}^{\kappa _m b}+d_m\textrm{e}^{-\kappa _m b}\big),\end{align}

(2.42) \begin{align}\frac {\partial \eta _1}{\partial x}\Big |_{x=b}&=\frac {\partial \eta _2}{\partial x}\Big |_{x=b}, \text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m\big(a_m\textrm{e}^{\kappa _m b}-b_m\textrm{e}^{-\kappa _m b}\big)\nonumber \\ &=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m\big(c_m\textrm{e}^{\kappa _m b}-d_m\textrm{e}^{-\kappa _m b}\big),\end{align}

(2.43) \begin{align}D\frac {\partial ^2\eta _1}{\partial x^2}\Big |_{x=b}&=D\frac {\partial ^2\eta _2}{\partial x^2}\Big |_{x=b}, \text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m^2\big(a_m\textrm{e}^{\kappa _m b}+b_m\textrm{e}^{-\kappa _m b}\big)\nonumber \\ &=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m^2\big(c_m\textrm{e}^{\kappa _m b}+d_m\textrm{e}^{-\kappa _m b}\big),\end{align}

(2.44) \begin{align}&\kern-50pt D\frac {\partial ^3\eta _1}{\partial x^3}\Big |_{x=b}-D\frac {\partial ^3\eta _2}{\partial x^3}\Big |_{x=b}=F,\text{i.e.,} \sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\nonumber \\ &\qquad \times \left [ \left(\kappa _m^3+\frac {s_r}{D}\frac {\tilde {\omega }^2}{1-\tilde {\omega }^2}\right)a_m\textrm{e}^{\kappa _m b}- \left(\kappa _m^3-\frac {s_r}{D}\frac {\tilde {\omega }^2}{1-\tilde {\omega }^2}\right) b_m\textrm{e}^{-\kappa _m b}\right ] \nonumber \\ &\quad=\sum _{m=-2}^{M}\frac {\textrm{i}}{\omega }\frac {\alpha }{\beta \kappa _m^4+1-\gamma _0\alpha }\kappa _m^3\big(c_m\textrm{e}^{\kappa _m b}-d_m\textrm{e}^{-\kappa _m b}\big).\end{align}

Combining (2.29)–(2.37) and (2.41)–(2.44), we have $4(N+3)$ equations and $4(M+3)$ unknowns in total. Letting $M=N$ and after some recombination, we get the following homogeneous equation:

(2.45) \begin{equation} \unicode{x1D646}\boldsymbol{x}=\boldsymbol{0}, \end{equation}

where $\boldsymbol{x}$ is the vector containing all the unknowns, i.e.

(2.46) \begin{equation} \boldsymbol{x}=[ a_{-2}, a_{-1}, \ldots , a_{M}, b_{-2}, b_{-1}, \ldots , b_{M}, c_{-2}, c_{-1}, \ldots , c_{M}, d_{-2}, d_{-1}, \ldots , d_{M}]^{\textrm{T}}, \end{equation}

and $\unicode{x1D646}$ is the coefficient matrix, written as

(2.47) \begin{equation} \unicode{x1D646}(\omega )=\left [ \begin{array}{cc} \boldsymbol{{\mathcal{A}}}(\omega ) & \textrm{e}^{-\textrm{i}ka}\boldsymbol{\mathcal{B}}(\omega )\\ \displaystyle \boldsymbol{\mathcal{C}}(\omega ) & \boldsymbol{\mathcal{D}}(\omega ) \\ \end{array} \right ]_{4(M+3)\times 4(M+3)}, \end{equation}

where submatrices $\boldsymbol{\mathcal{A}}, \boldsymbol{\mathcal{B}}, \boldsymbol{\mathcal{C}}$ and $\boldsymbol{\mathcal{D}}$ can be found in Appendix A.

The condition of non-zero solutions requires $\mbox{Det}[\unicode{x1D646}\,]=0$ , which reduces to the following generalised eigenvalue problem

(2.48) \begin{equation} \textrm {Det}[\boldsymbol{\mathcal{D}}(\omega )-\lambda \boldsymbol{\mathcal{M}}(\omega )]=0, \end{equation}

where $\lambda =\textrm{e}^{-\textrm{i}ka}$ and $\mathcal{M}(\omega )=\boldsymbol{\mathcal{C}}\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\mathcal{B}}$ . For any given angular frequency $\omega$ , we can solve for eigenvalues $\lambda$ , and the non-trivial solutions which represent the propagating waves correspond to the real values of $k$ such that $|\lambda |=1$ . This results in the dispersion relation $\omega (k)$ of the periodic floating metaplate.

3.1. Analytical model: uniform floating plate with effective medium

We refer to the model solved by the aforementioned numerical method as the real model, shown in figure 2(a), in which each unit cell contains a single spring–mass resonator. Its equivalent model is illustrated in figure 2(b), where a uniformly distributed spring–mass layer across the unit cell is used to replace the single resonator in the unit cell of the real model. It is ensured that the total mass and the total stiffness of the distributed spring–mass layer in the equivalent model are identical to those of the resonator in the real model. The underlying mathematical principle of this approximation is the use of a rectangular pulse of width $a$ to approximate the $\delta$ -function within one unit cell (Haberman Reference Haberman2013). After this discrete-to-continuous approximation, (2.4) should be modified as

(3.1) \begin{equation} D\frac {\partial ^4w}{\partial x^4}+\rho _p h\frac {\partial ^2w}{\partial t^2}=p+\frac {s_r}{a}(\tilde {w}-w)\quad z=0, \end{equation}

where $\tilde {w}(x,t)$ represents the vertical displacement of the distributed mass layer. Note that the last term in (3.1) characterises the traction (or pressure) acting on the upper surface of the plate, which reflects Hooke’s law. On the other hand, the governing equation with respect to the displacement $\tilde {w}(x,t)$ of the distributed mass layer reads

(3.2) \begin{equation} \frac {m_r}{a}\frac {\partial ^2\tilde {w}}{\partial t^2}=\frac {s_r}{a}(w-\tilde {w}), \end{equation}

which reflects Newton’s second law. Assuming the harmonic solution, eliminating the term $\tilde {w}$ by substituting (3.2) into (3.1), and combining with (2.5) and (2.6), we get

(3.3) \begin{equation} \left \{ D\frac {\partial ^4 }{\partial x^4}+\rho g-\left [\frac {m_r}{a}\frac {1}{1-\left(\frac {\omega }{\omega _0}\right)^2}+\rho _p h\right ]\omega ^2 \right \}\frac {\partial \phi }{\partial z}=\rho \omega ^2\phi , \end{equation}

where $\omega _0$ can be found in (2.40). Plugging the basic solution $\phi (x,z)=\textrm{e}^{\textrm{i}kx}\cosh k(z+H)$ , which satisfies (2.1) and boundary condition (2.2), into (3.3), we obtain the dispersion equation expressed as

(3.4) \begin{equation} \left \{\beta k^4+1-\left [\gamma _1\frac {1}{1-\left(\frac {\omega }{\omega _0}\right)^2}+\gamma _0\right ]\alpha \right \}k\tanh {kH}=\alpha , \end{equation}

where $\alpha$ , $\beta$ and $\gamma _0$ are given in (2.11), and

(3.5) \begin{equation} \gamma _1=\frac {m_r}{a\rho }. \end{equation}

Through (3.4) and (3.5), the dispersion relation $\omega (k)$ can be explicitly solved. For the sake of simplicity, it is listed in Appendix B.

Comparing (3.4) with the well-known dispersion equation for a bare plate floating atop water, which can be recast from (2.17) and written as (Fox & Squire Reference Fox and Squire1994)

(3.6) \begin{equation} \big(\beta k^4+1-\gamma _0\alpha \big)k\tanh {kH}=\alpha , \end{equation}

one can immediately recognise that (3.4) is nothing but the dispersion equation for a floating plate with the same value of $\beta$ but a modified $\gamma _0$ term, expressed as

(3.7) \begin{equation} \big(\beta k^4+1-\gamma _{\textit{eff}}\,\alpha \big)k\tanh {kH}=\alpha , \end{equation}

where

(3.8) \begin{equation} \gamma _{\textit{eff}}=\gamma _1\frac {1}{1-\left(\frac {\omega }{\omega _0}\right)^2}+\gamma _0. \end{equation}

From this perspective, the floating metaplate shown in figure 2(a) can be regarded as a uniform floating plate, represented by the effective model displayed in figure 2(c), with the parameter $\beta$ unchanged but $\gamma _0$ replaced by $\gamma _{\textit{eff}}$ . We can notice that $\gamma _{\textit{eff}}$ characterises the effective mass density of this effective-medium plate, embodying the inertia term of the original bare plate, i.e. $\gamma _0$ , and an additional mass density term, i.e. $\gamma _1$ , contributed from the resonator. This expression indicates that $\gamma _{\textit{eff}}$ of the effective plate is highly dependent on the ratio between the forcing frequency $\omega$ and the resonant frequency $\omega _0$ , closely resembling the effective dynamic mass defined in the literature on acoustic and elastic metamaterials (Huang, Sun & Huang Reference Huang, Sun and Huang2009; Fedele, Suryanarayana & Yavari Reference Fedele, Suryanarayana and Yavari2023). In the following, we show that the analytical dispersion relation (3.7), derived from this effective model, provides valuable insights into the band structure characteristics of such a complex system.

3.2. Band structure analysis

Unless otherwise specified, we will consider the values for the relevant parameters listed in table 1. Note that the values of $\beta$ and $\gamma _0$ are chosen according to those reported earlier (e.g. Meylan Reference Meylan2019; Zheng et al. Reference Zheng, Meylan, Zhu, Greaves and Iglesias2020; Iida et al. Reference Iida, Zareei and Alam2023).

Table 1. Values of the relevant parameters for solving the problem.

We first examine the dispersion curve for a bare plate floating on the water surface. To numerically derive the dispersion relation $\omega (k)$ , using the numerical scheme described in § 2, we set the spring stiffness of the resonator, $s_r$ , to approach zero, causing the resonant frequency $\omega _0$ to vanish. This limiting case of the metaplate exactly corresponds to a bare plate. In figure 3(a), we plot the dispersion curve as $s_r\to 0$ , alongside the analytical dispersion relation for a floating bare plate given by (3.6). The results obtained from these two approaches exhibit excellent agreement, validating our numerical scheme.

Figure 3. (a) Comparison of dispersion curve for a floating bare plate, obtained from numerical method and analytical formula. (b) Comparison of band structure for a floating metaplate, obtained from numerical method and analytical formulas enabled by the equivalent model shown in figure 2. The grey shading indicates the locally resonant bandgap obtained by numerical method. The inset shows a magnified view of the band structure at $k=\pi /a$ .

Then we consider the floating metaplate with reasonably selected parameter values shown in table 1. The band structures obtained from the numerical method are plotted in figure 3(b), together with those obtained from the analytical model described by (3.7) and (3.8) (see Appendix B for the explicit expression of $\omega (k)$ ). It is evident that a bandgap forms between frequencies around 9.5–10 rad s⁻¹, highlighted in grey, close to the resonant frequency of the resonator, i.e. 10 rad s⁻¹. This bandgap is induced by the local resonance rather than the Bragg scattering. We draw this conclusion based on the following two considerations. First, for a periodic floating metaplate with periodicity $a$ , the frequency at which a bandgap arises due to the Bragg scattering can be estimated using the condition $2\pi n/k=2a$ ( $n\in \mathbb{N}$ ) (Lorenzo et al. Reference Lorenzo, Pezzutto, De Lillo, Ventrella, De Vita, Bosia and Onorato2023), i.e.

(3.9) \begin{equation} \omega \approx \sqrt {g\frac {\beta k^5+k}{\gamma _0k +1}}\Bigg |_{k=\frac {n\pi }{a}}, \end{equation}

where the dispersion equation (3.6) has been used. For $n=1$ , the lowest Bragg bandgap occurs at around 13.2 rad s⁻¹, which is higher than the bandgap presented. Second, this bandgap exhibits the characteristic profile of a locally resonant bandgap (i.e. the lower bound of the bandgap is dictated by the lower branch of the dispersion curve which is flat at the first BZ edge $k=\pi /a$ , while the upper bound is governed by the higher branch which is flat at the first BZ centre $k=0$ ). In the following, we provide explanations of the resonance-induced bandgap obtained here by directly examining the behaviour as $\omega \to \omega _0$ in the original dispersion (3.4) (i.e. (3.7)), which can be fully expressed as

(3.10) \begin{equation} \left \{\beta k^4+1-\left [\gamma _1\frac {1}{1-\left(\frac {\omega }{\omega _0}\right)^2}+\gamma _0\right ]\frac {\omega ^2}{g} \right \}k\tanh {kH}=\frac {\omega ^2}{g}. \end{equation}

On the one hand, if $\omega \to \omega _{0-}$ , i.e. $\omega =\omega _0-\Delta \omega$ where $\Delta \omega \to 0$ , the term

(3.11) \begin{align} \gamma _1\frac {1}{1-\left (\frac {\omega _0-\Delta \omega }{\omega _0}\right )^2}\sim \frac {\gamma _1\omega _0}{2\Delta \omega }\to \infty , \end{align}

which dominates the terms in $[\boldsymbol{\cdot }]$ , then, the dispersion (3.10) can be approximated as

(3.12) \begin{equation} \left (\beta k^4+1-\frac {\gamma _1\omega _0^3}{2g\Delta \omega } \right )k\tanh {kH}=\frac {\omega _0^2}{g}. \end{equation}

Because the Bloch wavenumber $k$ is confined within the first BZ (i.e. $k\in [-\pi /a, \pi /a]$ ), the term $\beta k^4+1$ is positive and bounded by $\beta (\pi /a)^4+1$ . Since $k\tanh {kH}$ is non-negative, if $\beta (\pi /a)^4+1\lt \gamma _1\omega _0^3/2g\Delta \omega$ , the left-hand side of (3.12) is non-positive, while its right-hand side is positive, indicating no real solutions can be found for (3.12). This means that when $\omega$ is smaller than but very close to $\omega _0$ , there is a bandgap.

On the other hand, if $\omega \to \omega _{0+}$ , i.e. $\omega =\omega _0+\Delta \omega$ where $\Delta \omega \to 0$ , the term

(3.13) \begin{align} \gamma _1\frac {1}{1-\left (\frac {\omega _0+\Delta \omega }{\omega _0}\right )^2}\sim -\frac {\gamma _1\omega _0}{2\Delta \omega }\to -\infty , \end{align}

then, the dispersion (3.10) becomes

(3.14) \begin{equation} \left (\beta k^4+1+\frac {\gamma _1\omega _0^3}{2g\Delta \omega } \right )k\tanh {kH}=\frac {\omega _0^2}{g}. \end{equation}

Both sides of (3.14) are always positive if $k\neq 0$ , indicating it is possible to find real solutions and therefore propagating bands exist when $\omega$ is slightly larger than $\omega _0$ .

When $\omega =\omega _0$ , to make (3.10) valid, $k$ must be zero. To check if $(\omega =\omega _0, k=0)$ is a solution, we re-examine the expression of $\omega (k)$ derived from (3.4) (shown in Appendix B) and found that it is indeed a solution.

These findings are clearly manifested in the band structure shown in figure 3(b), where the bandgap occurs right below the resonator’s natural frequency $\omega _0$ , and the dispersion branch above the bandgap starts from $(\omega =\omega _0,k=0)$ . These features will be the basis in the design a broadband water wave reflector with a customisable frequency range that will be discussed in § 5.

3.3. Effect of key parameters on bandgap evolution

To demonstrate the generality of the features revealed by the band structure analysis in § 3.2, we investigate how the bandgap evolves as several key physical parameters are varied. Specifically, we examine the influence of the normalised stiffness parameter $\beta$ , the normalised mass term $\gamma _1$ , the depth of the water $H$ and the periodicity $a$ of the metaplate on the location and width of the bandgap.

Figure 4. Evolution of the bandgap on varying parameters (a) $\beta$ , (b) $\gamma _1$ , (c) $H$ and (d) $a$ . The four figures were obtained by changing their respective parameters while keeping the other parameters unchanged as original ones taken in figure 3, where $\beta =0.05, \gamma _1=0.01, H=10$ and $a=1$ .

First, by varying $\beta$ , the variation of the bandgap is shown in figure 4(a). It can be seen that the gap width decreases as $\beta$ increases. Note that the bandgap’s upper edge always occurs at $\omega _0= 10\,\textrm{rad s}^{-1}$ .

Second, we consider the variation of the parameter $\gamma _1=m_r/a\rho$ , which characterises the mass of the resonator. Increasing $\gamma _1$ means increasing the mass of the resonator $m_r$ , and consequently decreasing the resonant frequency as $\omega _0=\sqrt {s_r/m_r}$ . This is shown in figure 4(b), where the upper edge of the bandgap, determined by $\omega _0$ , indeed decreases as $\gamma _1$ increases.

Third, we change $H$ from its original value of 10 to 60 m, which means a deeper water region is considered. The variation of the bandgap is shown in figure 4(c), indicating that the bandgap is almost unchanged.

Last but not least, we change $a$ from its original value of 1 to 0.5 m. The variation of the bandgap is shown in figure 4(d), indicating that the bandgap size decreases as the periodicity $a$ decreases.

In a nutshell, although the frequency and the width of the bandgap may change when the parameters are altered, the overall features of the band structure of the decorated floating plate remain the same, i.e. there exists a resonance-induced bandgap located right below the resonant frequency $\omega _0$ .

3.4. Error analysis and applicable criteria for the analytical model

As shown in figure 3(b), the band structure derived from the simplified analytical model closely matches the numerical results across almost the entire wavenumber $k$ range in the first BZ, with noticeable discrepancies only near the first BZ edge. We will provide an explanation for the origin of these discrepancies and give an brief error-analysis–based discussion of when the analytical (equivalent) model remains reliably applicable, with a focus on the unit cell periodicity $a$ .

It is worth emphasising that the analytical model is derived by replacing discrete resonators with uniformly distributed masses and springs (Chen et al. Reference Chen, Li, Nassar, Norris, Daraio and Huang2019), as shown in figure 2. Mathematically, this averaging process can be regarded as using a rectangular pulse function to approximate the Dirac delta function (Haberman Reference Haberman2013), which results in the discrepancies observed in the band structures obtained from analytical and numerical models. In this approximation, the rectangular pulse function better approximates the $\delta$ -function as $a$ decreases, meaning the smaller the periodicity, the better the equivalent model. The discrepancies at large wave vectors stem from this approximation.

Figure 5. Variation of the relative error ${\textit{RE}}$ (3.15) between the frequency at $k=\pi /a$ from the two models with respect to the periodicity $a$ .

We define the relative error ( ${\textit{RE}}$ ) between the frequency at $k=\pi /a$ from the real model (numerical) and the equivalent model (analytical) as

(3.15) \begin{equation} RE=\frac {\left|\omega _{{re}}\left(\frac {\pi }{a}\right)-\omega _{{eq}}\left(\frac {\pi }{a}\right)\right|}{\omega _{{re}}\left(\frac {\pi }{a}\right)}\times 100\, \%, \end{equation}

where $\omega _{{re}}( {\pi }/{a})$ denotes the low-branch frequency at $\pi /a$ from the real model, while $\omega _{{eq}}( {\pi }/{a})$ is that from the equivalent model (referring to the inset in figure 3(b) to clearly visualise the discrepancy between them).

We evaluate the relative error ${\textit{RE}}$ for various values of $a$ while keeping all other parameters fixed as in table 1. Figure 5 presents the result, illustrating that, as the periodicity $a$ decreases, the relative error between the two models diminishes. This indicates that the equivalent model can more accurately represent the real system when the unit cell size of the periodic floating metaplate is relatively small, in agreement with the conclusion derived above from mathematical viewpoint.

4.1. Frequency-domain analysis

We can separate the region into three parts: free surface on the left-hand side for $x\lt 0$ (labelled by subscript ‘ $l$ ’, with transmitted wave), plate-covered region for $0\lt x\lt N_ca$ (labelled by subscript ‘ $c$ ’, with gravity flexural wave) and free surface on the right-hand side for $x\gt N_c a$ (labelled by subscript ‘ $r$ ’, with incident and reflected waves). By considering all equations and conditions mentioned in § 2.1 and employing the eigenfunction matching method described in § 2.2, we obtain the velocity potential underneath the three parts as follows (Kohout et al. Reference Kohout, Meylan, Sakai, Hanai, Leman and Brossard2007; Meylan Reference Meylan2019; Iida & Umazume Reference Iida and Umazume2020):

(4.1) \begin{equation} \varPhi (x,z,t)={\textit{Re}}\left [\frac {g}{\textrm{i}\omega }\phi _{l,c,r}\textrm{e}^{-\textrm{i}\omega t}\right ]\!, \end{equation}

where

(4.2) \begin{align}&\qquad\qquad\,\, \phi _l(x,z)=\sum _{n=0}^{N}A_n\textrm{e}^{\vartheta _nx}\varphi _n(z),\end{align}

(4.3) \begin{align}& \phi _c^{(i)}(x,z)=\sum _{m=-2}^{M}\big[B_m^{(i)}\textrm{e}^{\kappa _mx}+C_m^{(i)}\textrm{e}^{-\kappa _mx}\big]\psi _m(z),\end{align}

(4.4) \begin{align}&\quad\,\,\, \phi _r(x,z)=\textrm{e}^{\vartheta _0x}\varphi _0(z)+\sum _{n=0}^{N}D_n\textrm{e}^{-\vartheta _nx}\varphi _n(z).\end{align}

Figure 6. Schematic of a finite-size floating metaplate, consisting of $N_c$ unit cells. Waves are incident from the positive $x$ -direction. Resonators are located at $x=b+(i-1)a$ , where $i$ ranges from $1$ to $N_c$ .

Here, $\varphi _n(z)(n\geqslant 0)$ and $\psi _m(z)(m\geqslant -2)$ are given in (2.14) and (2.15), respectively; $\vartheta _n$ and $\kappa _m$ are determined from (2.16) and (2.17); $A_n, B_m^{(i)}, C_m^{(i)}$ and $D_n$ are unknown coefficients. In (4.3), the superscript $(i)$ ( $i=1,2,{\cdots} ,N+1$ ) indicates the $i$ th segment of the thin plate separated by the $(i-1)$ th and $i$ th resonators on it. Specifically, for $i=1$ , $0\lt x\lt b$ ; for $i=N_c+1$ , $b+(N_c-1)a\lt x\lt N_ca$ ; and for any other $i$ except $1$ and $N_c+1$ , $b+(i-2)a\lt x\lt b+(i-1)a$ . Note that (4.2) and (4.4) are obtained by taking into account the radiation conditions (Iida & Umazume Reference Iida and Umazume2020), which stipulate that the scattered waves (i.e. reflected and transmitted waves here) are outgoing to the infinite far field. In addition, the first term of (4.4) represents incident waves, i.e. $\varPhi _I={\textit{Re}} [(g/\textrm{i}\omega )\textrm{e}^{\vartheta _0x}\varphi _0(z)\textrm{e}^{-\textrm{i}\omega t} ]$ . The factor $g/\textrm{i}\omega$ is a normalising value to unify the wave elevation of the free surface, which is calculated by the dynamic condition as

(4.5) \begin{equation} w=-\frac {1}{g}\frac {\partial \varPhi }{\partial t}\Big |_{z=0}={\textit{Re}}\big[\phi _{l, r}(x,0)\textrm{e}^{-\textrm{i}\omega t}\big]. \end{equation}

And the harmonic deflection $\eta$ of the plate (referring to (2.20)) can be written as

(4.6) \begin{equation} \eta ^{(i)}=\left .\frac {1}{-\textrm{i}\omega }\frac {g}{\textrm{i}\omega }\frac {\partial \phi _c^{(i)}}{\partial z}\right |_{z=0}\quad i=1,2,\ldots ,N_c+1. \end{equation}

Other physical quantities, such as the horizontal velocity of the fluid $\partial \phi /\partial x$ , rotation angle $\partial \eta /\partial x$ , bending moment $D\partial ^2\eta /\partial x^2$ and shear force $-D\partial ^3\eta /\partial x^3$ of the plate, can also be expressed by the velocity potential, as we did in § 2.

We use truncated summations of finite terms to approximate the solution of the velocity potential. According to (4.2)–(4.4), there are $(N+1)+2(M+3)\times (N_c+1)+(N+1)=2(N+1)+2(N_c+1)(M+3)$ unknowns in total. They are determined by matching and boundary conditions, which will be explained in the following.

First, we consider the matching conditions for the velocity potential and its normal derivative across the boundary between distinct regions underneath the free surface or plate segments. This implies that

(4.7) \begin{align}&\qquad\qquad\qquad\qquad\qquad \phi _l(0,z)=\phi _c^{(1)}(0,z),\end{align}

(4.8) \begin{align}& \phi _c^{(i)}(b+(i-1)a,z)=\phi _c^{(i+1)}(b+(i-1)a,z)\quad (i=1,2,3,\ldots ,N_c),\end{align}

(4.9) \begin{align}&\qquad\qquad\qquad\qquad \phi _c^{(N_c+1)}(N_ca,z)=\phi _r(N_ca,z),\end{align}

and

(4.10) \begin{align}&\qquad\qquad\qquad\qquad\qquad \frac {\partial \phi _l}{\partial x}\Big |_{x=0}=\frac {\partial \phi _c^{(1)}}{\partial x}\Big |_{x=0},\end{align}

(4.11) \begin{align}& \frac {\partial \phi _c^{(i)}}{\partial x}\Big |_{x=b+(i-1)a}=\frac {\partial \phi _c^{(i+1)}}{\partial x}\Big |_{x=b+(i-1)a}\quad (i=1,2,3,\ldots ,N_c),\end{align}

(4.12) \begin{align}&\qquad\qquad\qquad\qquad \frac {\partial \phi _c^{(N_c+1)}}{\partial x}\Big |_{x=N_ca}=\frac {\partial \phi _r}{\partial x}\Big |_{x=N_ca}.\end{align}

The detailed forms of the matching conditions listed above are given in Appendix C. Following the approach used in § 2.3, we solve these equations by multiplying both sides by the vertical eigenfunction $\varphi _j(z)$ and integrating over the vertical domain from $z=-H$ to $z=0$ . This process yields a total of $(N+1)\times [2(N_c+2)]$ equations, which are elaborated in Appendix C.

Then we consider the boundary and matching conditions for the plate. Free boundary conditions are set for the two ends of the plate, i.e. the bending moments and the shear forces at $x=0$ and $x=N_ca$ are zero. This implies that

(4.13) \begin{align}& \frac {\partial ^2\eta ^{(1)}}{\partial x^2}\Big |_{x=0}=0,\end{align}

(4.14) \begin{align}& \frac {\partial ^3\eta ^{(1)}}{\partial x^3}\Big |_{x=0}=0,\end{align}

and

(4.15) \begin{align}& \frac {\partial ^2\eta ^{(N_C+1)}}{\partial x^2}\Big |_{x=N_ca}=0,\end{align}

(4.16) \begin{align}& \frac {\partial ^3\eta ^{(N_c+1)}}{\partial x^3}\Big |_{x=N_ca}=0.\end{align}

At the point where the resonator is located (i.e. $x=b+(i-1)a,i=1,2,\ldots ,N_c$ ), the deflection, the rotation angle and the bending moment of the plate are continuous, but there is a jump of force $F$ in the shear force, where the dynamic force $F$ was detailed in (2.38)–(2.40). This gives that

(4.17) \begin{align}\eta ^{(i)}\Big |_{x=b+(i-1)a}&=\eta ^{(i+1)}\Big |_{x=b+(i-1)a},\end{align}

(4.18) \begin{align}\frac {\partial \eta ^{(i)}}{\partial x}\Big |_{x=b+(i-1)a}&=\frac {\partial \eta ^{(i+1)}}{\partial x}\Big |_{x=b+(i-1)a},\end{align}

(4.19) \begin{align}D\frac {\partial ^2\eta ^{(i)}}{\partial x^2}\Big |_{x=b+(i-1)a}&=D\frac {\partial ^2\eta ^{(i+1)}}{\partial x^2}\Big |_{x=b+(i-1)a},\end{align}

(4.20) \begin{align}D\frac {\partial ^3\eta ^{(i)}}{\partial x^3}\Big |_{x=b+(i-1)a}-D\frac {\partial ^3\eta ^{(i+1)}}{\partial x^3}&\Big |_{x=b+(i-1)a}=s_r\eta ^{(i)}\Big |_{x=b+(i-1)a}\frac {-\tilde {\omega }^2}{1-\tilde {\omega }^2},\end{align}

where $i=1,2,\ldots ,N_c$ ; $\tilde {\omega }$ can be found in (2.40). The detailed forms of (4.13)–(4.20) are also presented in Appendix C.

Finally, we get a total of $2(N+1)(N_c+2) + 4 + 4N_c = (N_c+1)(2N+6) + 2(N+1)$ equations. For the system to be well determined, this must exactly match the number of unknowns, $2(N+1) + 2(N_c+1)(M+3)$ . This condition is ensured by letting $M = N$ .

4.2. Time-domain response

Using the procedure outlined above, the surface displacement in the frequency domain can be determined. The time-dependent solution is then obtained through superposition, leveraging the linearity of the system and applying the Fourier transform to transition from the frequency domain to the time domain. The surface displacement is a function of $\omega$ , therefore, we denote the complex frequency-domain surface displacement by $\eta (x,\omega )$ . We assume that the incident wave is a Gaussian at $t=0$ . The time-dependent displacement is given by the following Fourier integral (Meylan Reference Meylan2019):

(4.21) \begin{eqnarray} w(x,t)={\textit{Re}}\left \{\frac {1}{\pi }\int _0^{+\infty }\hat {f}(\omega )\eta (x,\omega )\textrm{e}^{-\textrm{i}\omega t}\textrm{d}\omega \right \}\!, \end{eqnarray}

where $\hat {f}(\omega )$ is

(4.22) \begin{eqnarray} \hat {f}(\omega )=\sqrt {\frac {s}{\pi }}\textrm{e}^{-s(\omega -\omega _c)^2}. \end{eqnarray}

Here, $s$ is a scaling factor that controls the overall width of the Gaussian wave packet, while $\omega _c$ represents the central frequency of the Gaussian.

4.3. Numerical results

We consider a finite-size floating metaplate consisting of $N_c=30$ unit cells. All other parameters remain identical to those listed in table 1. In the following, we will first present the reflection and transmission coefficients obtained from the frequency-domain solutions for each individual incident frequency. Then, we will provide animations demonstrating the propagation and reflection of incident waves at various frequencies based on the time-domain results.

Figure 7. (a) Transmittance of the finite-size floating metaplate with resonant frequency $\omega _0=10$ rad s⁻¹. The grey shading indicates the locally resonant bandgap as shown in figure 3(b). (b) Value of $|T|^2+|R|^2$ for verifying the law of conservation of energy.

In the frequency-domain response analysis, the characteristics of reflection and transmission are evaluated based on the amplitude of the velocity potential in the far field from the plate. The reflection and transmission coefficients, $R$ and $T$ , are given as

(4.23) \begin{eqnarray} R=D_0,\quad T=A_0, \end{eqnarray}

which need to satisfy

(4.24) \begin{eqnarray} |R|^2+|T|^2=1, \end{eqnarray}

as required by the energy conservation principle.

Figure 7(a) illustrates the transmittance, $20\log |T|$ , of the finite-size floating metaplate as a function of the incident wave frequency $\omega$ . A sharp dip is observed within the grey-shaded frequency range, corresponding to the locally resonant bandgap identified in figure 3(b). This highlights the distinctive characteristic of the locally resonant bandgap, namely, its pronounced attenuation effect on incident waves within this frequency range. Figure 7(b) presents $|R|^2+|T|^2$ as a function of $\omega$ , providing a clear verification of the energy conservation principle while also validating the accuracy and reliability of our numerical calculations.

On the other hand, the time-domain response results can provide a clearer demonstration of the blocking effect of the locally resonant bandgap on incident waves. We consider incident waves centred at different frequencies, corresponding to the passband below the bandgap, the stop band within the bandgap and the passband above the bandgap. The time-dependent responses for these three cases are shown as animations in Movies 1–3 which are given as supplementary material. From these results, we observe that for incident waves with a central frequency within the passbands (both below and above the bandgap), a significant portion of the wave energy propagates through the metaplate region. In contrast, for incident waves with a central frequency within the bandgap, almost all the energy is reflected. These results not only further validate our bandgap analysis but also highlight the potential of the floating metaplate as an effective water wave reflector.

5.1. Rainbow reflection mechanism for the graded metaplate

We consider a finite-size metaplate consisting of $N_c=30$ unit cells, where the resonant frequency $\omega _0$ of the resonator mounted in each unit cell gradually increases linearly from 8 rad s⁻¹ at the right end to 11 rad s⁻¹ at the left end, as illustrated in figure 8(a). To enhance the practical feasibility of this set-up, we assume that the mass of the resonators remains constant, while the gradient variation in the resonant frequency of the resonators is achieved by altering the stiffness, i.e. $s_r$ , of the springs connecting the resonators to the plate.

Unlike the conventional technique, which typically relies on the lower branch of the dispersion curve below the bandgap to achieve rainbow reflection (e.g. Bennetts et al. Reference Bennetts, Peter and Craster2018; Wilks et al. Reference Wilks, Montiel and Wakes2022), here, we utilise the second branch of the dispersion curve located above the bandgap in the metaplate structure. This choice is motivated by two key factors. First, the second branch features a broad flat-band region (referring to figure 3 b), which facilitates the phenomenon where different frequency components of the incident wave are effectively stopped at different positions within the structure. Second, as discussed in § 3.2, the lower edge of the second branch is precisely determined by the resonant frequencies of the resonators. This property enables highly customisable operating frequencies for the wave reflector, as will be elaborated later.

In order to reveal the rainbow reflection mechanism of this graded floating metaplate, the second branch of dispersion curves corresponding to the 1st, 5th, 10th, 15th, 20th, 25th and 30th unit cells are plotted together, as shown in figure 8(b). Considering an incident wave with a frequency of, for example, around 10 rad s⁻¹, indicated by the black dashed arrow, as it propagates from the right end of the host plate into the graded region, it can pass the first few unit cells since it lies within their passing band until it reaches the unit cell whose band edge corresponds to 10 rad s⁻¹. The wave modes are gradually transferred along the arrow from the dark blue curve to the yellow one, while the group velocity of the wave gradually decreases until it reaches zero. Then the wave will be stopped here, and, because the frequency of the incident wave is close to the resonant frequency of the resonator at that location, a significant portion of the incident wave energy will be transferred into the vibration of the resonator. This energy is subsequently converted into the bending vibration of the plate and reflected back. In general, incident waves of different frequencies reach zero group velocity at different unit cells. Since the lower edge of the flat band is precisely determined by the resonant frequency of the resonators, waves are stopped and reflected at the position where the resonant frequency of the local resonator matches or is close to the incident frequency (beyond this position, the waves enter the bandgap region, as shown in figure 8 c), resulting in a spatial separation of incident water waves with different frequencies, akin to the effect of a rainbow. Besides, owing to the locally resonant mechanism, energy will be accumulated at the corresponding resonator, leading to the energy localisation in this resonator. This mechanism explains the reflection process of incident waves and provides an intuitive theoretical foundation for the implementation of a broadband water wave reflector discussed in the following sections.

5.2. Frequency-domain response

To demonstrate the rainbow reflection phenomenon, we perform a frequency response analysis on the graded metaplate. The numerical analysis procedure is almost the same as that described in § 4.1, with the exception that the boundary conditions related to the resonators must be modified to account for the fact that the resonant frequencies are no longer identical but gradually vary. This requires rewriting (4.20) as follows:

(5.1) \begin{eqnarray} D\frac {\partial ^3\eta ^{(i)}}{\partial x^3}\Big |_{x=b+(i-1)a}-D\frac {\partial ^3\eta ^{(i+1)}}{\partial x^3}\Big |_{x=b+(i-1)a}=s_r^{(i)}\eta ^{(i)}\Big |_{x=b+(i-1)a}\frac {-[\tilde {\omega }^{(i)}]^2}{1-{[\tilde {\omega }^{(i)}}]^2}, \end{eqnarray}

where $i=1,2,\ldots ,N_c$ , $\tilde {\omega }^{(i)}=\omega /\omega _0^{(i)}$ , in which $\omega _0^{(i)}$ is the value of the resonant frequency of the $i$ th resonator (counting from the left end), and $s_r^{(i)}=m_r{ [\omega _0^{(i)} ]}^2$ is the spring stiffness of the $i$ th resonator, where $m_r$ is the mass of the resonator. In the numerical calculations, all other related parameters are consistent with those adopted in previous sections.

Figure 9. Amplitude of the reflection coefficient $|R|$ for the finite-size floating metaplate with the graded resonant frequencies ranging from 8 to 11 rad s⁻¹.

Figure 9 illustrates the amplitude of the reflection coefficient $|R|$ for the graded floating metaplate depicted in figure 8(a). The light cyan shaded region represents the frequency range of the resonant frequencies, which spans from $8$ to 11 rad s⁻¹. It is evident that total wave reflection occurs for incident waves with frequencies within this range. This result is consistent with the previously explained rainbow effect: as the wave propagates along the metaplate, it reaches the resonator whose resonant frequency matches that of the incident wave, leading to reflection. Therefore, this graded floating metaplate can function as a ‘broadband wave reflector’, with its working frequency range entirely determined by the graded frequency range of the resonators.

Figure 10. Displacement field of (a) the plate and (b) the resonators at different frequencies to demonstrate the rainbow reflection effect. In (a), A and B represent two specific cases where incident waves coming from the positive $x$ -axis at different frequencies are completely stopped at distinct positions.

To illustrate the rainbow reflection and the broadband reflection phenomenon, figure 10 depicts the frequency response field for the graded metaplate. Figure 10(a) displays the normalised displacement field of the plate. The dark black triangular region indicates the behaviour within the graded resonant frequency range ( $8$ to 11 rad s⁻¹), where the incident wave propagates a certain distance through the array before being stopped at a specific location. As the incident frequency $\omega$ increases linearly within this range, the stopping point shifts correspondingly in a linear manner, progressing from the right end toward the left end. Figure 10(b) illustrates the normalised absolute displacement of all resonators. One can observe that within the graded resonant frequency range of the resonators, the amplitudes of the resonators near the stopping point of the incident wave are significantly enhanced, and for different incident frequencies, the positions of the enhanced resonators vary accordingly. This unambiguously demonstrates the rainbow trapping effect. Moreover, it validates the previously described rainbow reflection mechanism, wherein the incident wave stops at a location where the resonant frequency of the local resonator matches the frequency of the incident wave and, subsequently, part of the wave’s energy is converted into the vibration of the resonator.

5.3. Time-dependent response

Following the procedure outlined in § 4.2 and based on the frequency-domain response field given above, we now present the time-dependent response of the graded floating metaplate under the excitation of an incident wave centred at a specific frequency.

Figure 11 shows the results in the form of waterfall plots for the transient displacement field of the plate and resonators. For an incident Gaussian wave with central frequency $\omega _c=9\,\textrm{rad s}^{-1}$ , it impinges into the right end of the plate ( $x=30\,\textrm{m}$ ) and propagates to a position approximately one third of the plate’s total length from the right end ( $x\approx 20\,\textrm{m}$ ). At this point, the group velocity of the wave decreases to zero, causing the wave to stop (see figure 11 b). A portion of the incident wave’s energy is then converted into the vibrational energy of the resonators near this location (see figure 11 c). Over time, the wave is reflected back, and this back-scattering process repeats multiple times before gradually reaching a state of equilibrium. This continuous radiation of the wave signal is a characteristic phenomenon observed in rainbow structures, akin to the time spreading of reflected pulses in acoustics (Cebrecos et al. Reference Cebrecos, Picó, Sánchez-Morcillo, Staliunas, Romero-García and Garcia-Raffi2014). To visualise this result clearly, some snapshots of the time-dependent responses of the plate and resonators are presented in Appendix D, and the full animation can be found in Movie 4 as supplementary material. For comparison, figure 11(d,e, f) displays the results for an incident wave centred at $\omega _c= 10\,\textrm{rad s}^{-1}$ . In this case, the wave is trapped and reflected at a position approximately one third of the plate’s total length from the left end ( $x\approx 10\,\textrm{m}$ ) (an animation illustrating this phenomenon is provided in Movie 5). These two cases align precisely with markers A and B in the frequency-domain wavefield shown in figure 10(a), which clearly indicates that rainbow reflection occurs around $x=20\,\textrm{m}$ and $x=10\,\textrm{m}$ , respectively. In a nutshell, the characteristic behaviour of the rainbow trapping phenomenon manifests for incident waves centred at frequencies within the range of the resonant frequencies of the graded resonators, facilitating broadband total wave reflection.

Figure 11. (a) An incident Gaussian wave centred at $\omega _c= 9\,\textrm{rad s}^{-1}$ , corresponding to marker A in figure 10(a). Waterfall plot representing the time-domain responses of (b) the plate and (c) the resonators to the incident Gaussian wave coming from the positive $x$ -axis. Snapshots at some different times are shown in figure 14 in Appendix D. Panels (d), (e) and ( f) are the same as (a), (b) and (c), respectively, except that now the incident wave is centred at $\omega _c=10\,\textrm{rad s}^{-1}$ , corresponding to marker B in figure 10(a). Hereafter, ‘FFT’ refers to the magnitude of the fast Fourier transform of the Gaussian wave packet.

Figure 12. (a) A broadband Gaussian wave centred at $\omega _c=9.5\,\textrm{rad s}^{-1}$ , covering the total-reflection frequency range. (b) A broadband Gaussian wave centred at $\omega _c=8\,\textrm{rad s}^{-1}$ , spanning frequencies both inside and outside the range of total reflection. (c) Waterfall plot showing the time-domain responses of the water–metaplate–water region to the incident Gaussian wave in (a), coming from the positive $x$ -axis. (d) Same as (c) but corresponding to the incident Gaussian wave in (b).

To further substantiate our findings, we performed two additional time-domain simulations using broadband Gaussian wave packets. One with central frequency $\omega _c=9.5\,\textrm{rad s}^{-1}$ , shown as figure 12(a), covers the entire frequency range of total reflection (i.e. 8–11 rad s⁻¹). The corresponding time-domain displacement field of the water–metaplate–water region, as shown in figure 12(c), indicates that the packet is almost completely reflected, which confirms the broadband reflection of the graded floating metaplate. An animation illustrating this result is provided in Movie 6. The other one with spectrum centred at 8 rad s⁻¹, as shown in figure 12(b), extends frequencies from 6 to 10 rad s⁻¹, covering the range both inside and outside total reflection. From the time-domain displacement field shown in figure 12(d), it can be seen that a portion of the incident wave is transmitted through the metaplate, confirming that the waves with frequencies outside the total-reflection band cannot be blocked. The animation illustrating this phenomenon can be found in Movie 7.

5.4. Customisable operating frequency range

As previously mentioned, thanks to the fact that the upper bound (i.e. cutoff frequency of the graded structure) of the locally resonant bandgap is precisely determined by the resonant frequency of the resonators (see § 3.2), rainbow reflection occurs for incident waves with central frequencies falling within the range that exactly matches the resonant frequency range of the graded resonators. In the following, we will briefly demonstrate that this frequency range, in which the graded-metaplate-enabled wave reflector operates, can be precisely customised by presetting the resonant frequency range of the graded resonators.

Figure 13. (a) Amplitude of the reflection coefficient $|R|$ for the finite-size floating metaplate with the graded resonant frequencies ranging from $5$ to 8 rad s⁻¹. (b) Displacement field of the plate and (c) the resonators at different frequencies, illustrating the rainbow reflection effect. Panels (b) and (c) are the same as figure 10, execpt that the frequency range associated with the rainbow reflection shifted to the new frequency range of the graded resonators.

Now, we set the resonant frequency of the resonator to change, without the loss of generality, from 5 rad s⁻¹ at the right-most end unit cell to 8 rad s⁻¹ at the left-most end unit cell. Following the numerical procedure outlined previously, we obtain the reflection coefficient and the frequency-domain wavefield for this case, as shown in figure 13. One can clearly notice that the frequency range over which total reflection occurs (i.e. $|R| \to 1$ ) shifts to the new frequency range of the graded resonators, as indicated by the light cyan region in figure 13(a). Once again, figures 13(b) and 13(c) illustrate that rainbow reflection is observed for incident waves with frequencies precisely within the resonant range of the resonators. In summary, the operating frequency range of this wave reflector, facilitated by the graded floating metaplate, can be accurately tuned by adjusting the predefined resonant frequency range of the resonators.