1. Introduction
Structures comprised of closely spaced parallel arrays of thin plates are useful devices in the bespoke manipulation of waves in several physical settings including acoustics (Zhu et al. Reference Zhu, Chen, Zhu, Garcia-Vidal, Yin, Zhang and Zhang2013; Jan & Porter Reference Jan and Porter2018; Porter Reference Porter2021; Bravo & Maury Reference Bravo and Maury2023), electromagnetics (Putley et al. Reference Putley, Guenneau, Porter and Craster2022, Reference Putley, Guenneau, Craster, Davies and Poulton2023), elasticity (Colombi et al. Reference Colombi, Colquitt, Roux, Guenneau and Craster2016; Colquitt et al. Reference Colquitt, Colombi, Craster, Roux, Guenneau and Craster2017; De Ponti, Iorio & Ardito Reference De Ponti, Iorio and Ardito2022) and water waves (Zheng, Porter & Greaves Reference Zheng, Porter and Greaves2020; Porter, Zheng & Liang Reference Porter, Zheng and Liang2022; Wilks, Montiel & Wakes Reference Wilks, Montiel and Wakes2022; Kucher et al. Reference Kucher, Koźluk, Petitjeans, Maurel and Pagneux2023; Zheng, Liang & Greaves Reference Zheng, Liang and Greaves2024). The key underpinning feature in all such applications is how flux is restricted by the narrow channels between adjacent plates in the device, compared with the isotropic nature of propagation in the surrounding medium. The wavelength is thus implicitly assumed to be much larger than the characteristic separation between adjacent plates. This contrast in length scales and the unusual wave phenomena, such as negative refraction (Porter Reference Porter2021), that can result from the anisotropy has led to such plate-array devices being classified as a type of metamaterial (Maier Reference Maier2017). Additionally, the finite length of the channels within compact devices means that they typically support local resonant modes thereby allowing small devices (less than a wavelength, say, in size) to have a disproportionately large effect on the external wave field (Zheng et al. Reference Zheng, Porter and Greaves2020).
Owing to the contrast in scales, several studies have investigated the effect of plate-array metastructures on waves by replacing the discrete structure of the plate array with an effective medium after implementing a low-frequency homogenisation approach. This allows wave interaction with plate-array devices having certain simple geometrical shapes to be analysed using established mathematical techniques for solving partial differential equations. For example, rectangular and cylindrical structures lend themselves to separation methods (e.g. as considered in Zheng et al. (Reference Zheng, Porter and Greaves2020) and Porter (Reference Porter2021)) and, in rare cases, mathematical methods can be applied to more complex geometries (e.g. Jan & Porter (Reference Jan and Porter2018) who considered a trapezoidal plate-array cavity in a waveguide wall). One of the restrictions of homogenisation, however, is that it does not apply close to internal channel resonance where local effects destroy the assumption of a contrast in scales. Thus, it has been shown in Putley et al. (Reference Putley, Guenneau, Porter and Craster2022) and Jan & Porter (Reference Jan and Porter2018) for example that the problems become ill-posed in frequency intervals where resonance is present on account of the assumptions of low-frequency homogenisation having been violated. Problems can be regularised by the introduction of a small amount of dissipation (as in Jan & Porter (Reference Jan and Porter2018) and Zheng et al. (Reference Zheng, Porter and Greaves2020)) into the effective field equations, but this ‘sticking-plaster approach’ overlooks the precise nature of the influence of the local channel scale.
In this paper, we present a methodology which allows us to investigate wave interaction with structures comprised of discrete plate-arrays; that is, without the homogenisation. Such an approach is not new: see Porter (Reference Porter2021) who used Fourier transform methods to compare wave scattering by an infinitely long rectangular strip filled with a periodic array of tilted plates with the equivalent homogenisation theory. Resonant amplification is not encountered in this problem and the discrete plate array description was shown to converge rapidly to the homogenised description with near-identical results for the far-field scattered amplitudes when the channel width to length ratio fell below 0.1. Experimental results of Kucher et al. (Reference Kucher, Koźluk, Petitjeans, Maurel and Pagneux2023) also supported this conclusion. The idea of using Fourier transforms also underpins the current work where the focus is on methods for determining wave scattering by more general, non-regular, metastructures. In particular, we focus on the effect on wave propagation of so-called graded plate-arrays in which the width of the channels in the device is non-constant (typically increasing linearly, and thus forming a wedge).
Graded metamaterials have been of interest to researchers in a range of different applications since they produce broadbanded effects. For example, in Colombi et al. (Reference Colombi, Colquitt, Roux, Guenneau and Craster2016) and Colquitt et al. (Reference Colquitt, Colombi, Craster, Roux, Guenneau and Craster2017) a graded array placed on the surface of an elastic half-space was shown to deflect surface Rayleigh waves into elastic body waves and it was later proposed (e.g. Brûlé, Enoch & Guenneau Reference Brûlé, Enoch and Guenneau2020) as a scheme for protecting infrastructure from earthquakes. In acoustics Zhu et al. (Reference Zhu, Chen, Zhu, Garcia-Vidal, Yin, Zhang and Zhang2013) have graded structures to provide broadbanded absorption of sound by a metasurface, and Jan & Porter (Reference Jan and Porter2018) and Bravo & Maury (Reference Bravo and Maury2023) showed that a metamaterial plate-array cavity could suppress acoustic transmission in waveguides over a wide range of frequencies. In water waves Wilks et al. (Reference Wilks, Montiel and Wakes2022); Wilks, Montiel & Wakes (Reference Wilks, Montiel and Wakes2023) have similarly shown the broadbanded reflective qualities of a graded array of plates submerged through the surface and also been proposed its extension as a wave energy harnessing device. So-called rainbow reflection and rainbow trapping and absorption by graded metamaterials have also featured in the work of Tsakmakidis, Boardman & Hess (Reference Tsakmakidis, Boardman and Hess2007), Jimenez et al. (Reference Jimenez, Romeo-Garcia, Pagneux and Groby2017), Bennetts, Peter & Craster (Reference Bennetts, Peter and Craster2018), Chaplain et al. (Reference Chaplain, Pajer, De Ponti and Craster2020) and De Ponti et al. (Reference De Ponti, Iorio and Ardito2022). Circular metacylinders comprised of a plate array are also graded, although not linearly, and have exhibited (e.g. Zheng et al. Reference Zheng, Porter and Greaves2020; Putley et al. Reference Putley, Guenneau, Craster, Davies and Poulton2023) similar features: a slowing wave speed and amplification of wave energy through the structure with a strong broadbanded reflective quality.
We consider three problems all set in the context of linearised water waves although the first two problems have analogues in other physical settings. In all three problems, oblique plane waves are scattered by metastructures consisting of a discrete plate array with elements which are arbitrary in separation and width allowing us to consider metastructures of general shape. In the first problem, described in § 2, we consider a single such device consisting of vertical plates extending fully through the water depth. In § 3 the second problem involves an infinite periodic array of these devices. In the final problem (§ 4) the plates extend only partially through the fluid depth, this problem being identical to that studied by Wilks et al. (Reference Wilks, Montiel and Wakes2022).
We propose a common method of solution based on transforms (infinite Fourier for the first problem, and finite transforms for the last two) in which the solution in the presence of 
$N+1$ plates of varying positions and lengths is shown to be expressed by the same simple characteristic formulation. This simplicity, an overlooked highlight of the related work of Noad & Porter (Reference Noad and Porter2015), is in contrast with, for example, Roy, De & Mandal (Reference Roy, De and Mandal2019) and Wilks et al. (Reference Wilks, Montiel and Wakes2022, Reference Wilks, Montiel and Wakes2023) who use separation solutions in each of the channel-based domains and then performed matching from one channel to the next using relatively convoluted methods.
Although there is a focus on the method of solution to these problems, the main emphasis is on the results which are presented in § 5. Here we compare discrete plate array results with existing results including those determined by homogenisation and present extensions to results inaccessible to homogenisation methods with a focus on resonance. This includes looking at the effects of graded arrays with a view to applications as sea defence systems. We conclude the work in § 6.
2. A plate array metastructure in an open domain
We consider waves on a fluid of constant depth 
$h$ with a free surface whose rest position is given by 
$z=0$, 
$z$ being the vertical coordinate, directed upwards out of the fluid. We suppose that a parallel array of 
$N+1$ thin vertical barriers occupy the surfaces 
$x=x_j$, 
$-h < z < 0$, 
$|y| < b_j$, for 
$j=0,\ldots,N$, as illustrated in figure 1. A surface wave of angular frequency 
$\omega$ is incident from infinity, heading at an anticlockwise angle 
$\theta _0$ with respect to the positive 
$x$-direction. On the assumptions of linearised water wave theory, its motion and the subsequent response of the fluid due to the interaction with the array of barriers may be described by a velocity potential (e.g. Linton & McIver Reference Linton and McIver2001)
\begin{equation} \varPhi(x,y,z,t) = {\mathrm{Re}} \{ \phi(x,y) \psi_0(z) {\mathrm{e}}^{-{\rm i} \omega t} \}, \end{equation}where the uniformity of the geometry through the depth allows us to factorise a depth dependence
\begin{equation} \psi_0(z) = N_0^{{-}1/2} \cosh k (z+h), \quad \text{and} \quad N_0 = \frac12 \left( 1 + \frac{\sinh 2 kh }{2 kh} \right) \end{equation}
is a normalising factor whilst 
$k$ is the positive real root of
\begin{equation} \omega^2/g \equiv K = k \tanh kh, \end{equation}
the usual dispersion relation for water waves with gravitational acceleration given by 
$g$. The wave elevation is proportional to 
$\phi (x,y)$. Consequently, the reduced two-dimensional complex velocity potential 
$\phi (x,y)$ satisfies the Helmholtz equation
\begin{equation} \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + k^2 \right) \phi = 0. \end{equation}Within this framework, the incident wave is described by the function
\begin{equation} \phi_{inc}(x,y) = {\mathrm{e}}^{{\rm i} \alpha_0 x} {\mathrm{e}}^{{\rm i} \beta_0 y}, \end{equation}
where 
$(\alpha _0 ,\beta _0) = k(\cos \theta _0, \sin \theta _0)$ and we require that 
$\phi (x,y) - \phi _{inc}(x,y)$ represents outgoing waves as 
$k r \to \infty$ where 
$r = (x^2 + y^2)^{1/2}$. Specifically, we write
\begin{equation} \phi(x,y) - \phi_{inc}(x,y) =\phi_{sca}(x,y)\sim A(\theta;\theta_0) \sqrt{ \frac{2}{{\rm \pi} k r}} {\mathrm{e}}^{{\rm i} k r - {\rm i} {\rm \pi}/4}, \end{equation}
where 
$(x,y)=r(\cos \theta,\sin \theta )$ and 
$A(\theta ;\theta _0)$ is defined as the diffraction coefficient, measuring the amplitude of circular waves scattered in the direction 
$\theta$ due to an incident wave heading 
$\theta _0$.
Figure 1. Sketch of wave interactions with a plate-array metastructure.
The scattering of waves is due to the presence of barriers on which the following conditions apply:
\begin{equation} \frac{\partial \phi}{\partial x} = 0, \quad x = x_j^\pm, \quad |y| < b_j, \ (j=0,\ldots,N). \end{equation} We remark that the boundary-value problem posed above can be interpreted in physical settings other than water waves including, for example, two-dimensional acoustics or transverse electrically polarised electromagnetics, in which the factorisation of the 
$z$-dependence and the dispersion relation will both differ.
The method of solution for this problem is described in the work of Noad & Porter (Reference Noad and Porter2015) but we include below a key simplification to the solution method which will be reused in later sections. Thus, we introduce the Fourier transform pair
\begin{equation} \bar{\phi}(x;\beta) = \int_{-\infty}^\infty [\phi(x,y)- \phi_{inc}(x,y)] {\mathrm{e}}^{- {\rm i} \beta y} \, {\rm d} y \end{equation}and
\begin{equation} \phi(x,y) = \phi_{inc}(x,y) + \frac{1}{2 {\rm \pi}} \int_{-\infty}^\infty \bar{\phi}(x;\beta) {\mathrm{e}}^{{\rm i} \beta y} \, {\rm d} \beta. \end{equation}Then the governing wave equation is transformed to
\begin{equation} \left(\frac{{\rm d}^2}{{\rm d}\kern0.7pt x^2} - \gamma^2 \right) \bar{\phi} = 0, \quad x \neq x_j \end{equation}
(
$\,j=0,\ldots,N$) where
\begin{equation} \gamma = \left\{ \begin{array}{@{}ll} \sqrt{\beta^2 - k^2}, & |\beta| \geq k,\\ -{\rm i} \alpha, & |\beta| < k, \end{array} \right. \end{equation}
where 
$\alpha = \sqrt {k^2 - \beta ^2}$ and the choice of the complex branch of the square root function is made to satisfy the radiation condition at infinity (this becomes clear only later on). We note the transformation of the barrier conditions leads to the jump conditions
\begin{equation} \bar{\phi}_x (x_j^+; \beta) - \bar{\phi}_x (x_j^-; \beta) = 0 \end{equation}and
\begin{equation} \bar{\phi} (x_j^+; \beta) - \bar{\phi} (x_j^-; \beta) = P_j(\beta), \end{equation}
for 
$j=0,\ldots, N$ where
\begin{equation} P_j(\beta) = \int_{{-}b_j}^{b_j} p_j(y) {\mathrm{e}}^{-{\rm i} \beta y}\, {\rm d} y \end{equation}using the definition
\begin{equation} \phi(x_j^+,y) - \phi(x_j^-,y) = \left\{ \begin{array}{@{}ll} p_j(y), & |y| < b_j,\\ 0, & |y| > b_j. \end{array} \right. \end{equation}
Rather than expand the solution in each of the 
$N+2$ domains 
$x < x_0$, 
$x_{j-1} < x < x_j$ (
$\,j=1,\ldots,N$) and 
$x > x_N$ and match using (2.12) and (2.13), as in Noad & Porter (Reference Noad and Porter2015), we adopt a much more elegant approach which results in the same final expression and is easy to adapt to other problems.
Let us define the canonical function 
$g(x,x_j;\beta )$ as the solution of
\begin{equation} \left(\frac{{\rm d}^2}{{\rm d}\kern0.7pt x^2} - \gamma^2 \right) g = 0, \quad x \gtrless x_j \end{equation}
satisfying jump conditions 
$g_x(x_j^+,x_j;\beta ) - g_x(x_j^-,x_j;\beta ) = 0$ and 
$g(x_j^+,x_j;\beta ) - g(x_j^-,x_j;\beta ) = 1$ such that 
$g$ is outgoing (when 
$|\beta |< k$) or exponentially decaying (when 
$|\beta | > k$) as 
$k|x-x_j| \to \infty$. It is straightforward to confirm that
\begin{equation} g(x,x_j;\beta) ={-} \frac{1}{2} \mbox{sgn}(x-x_j) {\mathrm{e}}^{-\gamma|x-x_j|}. \end{equation}The solution of (2.10), (2.12), (2.13), with outgoing waves at infinity is given by the weighted superposition
\begin{equation} \bar{\phi}(x;\beta) = \sum_{j=0}^N P_j(\beta) g(x,x_j;\beta) ={-}\frac12 \sum_{j=0}^N P_j(\beta) \mbox{sgn}(x-x_j) {\mathrm{e}}^{-\gamma|x - x_j|}. \end{equation}The general solution throughout the domain is given by inverting the transform, thus
\begin{equation} \phi(x,y) = \phi_{inc}(x,y) - \frac{1}{4 {\rm \pi}} \sum_{j=0}^N \mbox{sgn}(x-x_j) \int_{-\infty}^\infty {\mathrm{e}}^{-\gamma|x - x_j|} {\mathrm{e}}^{{\rm i} \beta y} \int_{{-}b_j}^{b_j} p_j(y') {\mathrm{e}}^{-{\rm i} \beta y'} \, {\rm d} y' \, {\rm d} \beta. \end{equation}We note that this representation of the general solution may also be obtained by distributing Green's functions over the barriers and applying the conditions on the barriers. The particular form expressed above requires that the integral representation of the Hankel function (representing the Green's function) given by (A2) is used and the ordering of integrals is interchanged. The advantage of using the representation (2.19) of the solution, rather than a Green's function representation, is that we encounter no technical issues relating to convergence. In contrast, the Green's function approach leads to integrals with hypersingular kernels having to be treated as Hadamard finite-part integrals (see Martin (Reference Martin1991) for example). Despite the complexity involved in handling the hypersingular kernel, methods based on boundary integral equations in conjunction with Green's function still remain widely used due to their flexibility and their ability to handle complex configurations (e.g. see Martin Reference Martin1991; Renzi & Dias Reference Renzi and Dias2012; Hariri Nokob & Yeung Reference Hariri Nokob and Yeung2015).
The particular solution is determined by applying the barrier conditions (2.7) which results in the coupled integral equations
\begin{equation} \frac{1}{4 {\rm \pi}} \sum_{j=0}^N \int_{-\infty}^\infty \gamma {\mathrm{e}}^{-\gamma|x_j - x_l|} {\mathrm{e}}^{{\rm i} \beta y} \int_{{-}b_j}^{b_j} p_j(y') {\mathrm{e}}^{-{\rm i} \beta y'} \, {\rm d} y' \, {\rm d} \beta ={-}{\rm i} \alpha_0 {\mathrm{e}}^{{\rm i} \alpha_0 x_l} {\mathrm{e}}^{{\rm i} \beta_0 y}, \quad |y| < b_l, \end{equation}
for 
$l=0,\ldots,N$ for the 
$N+1$ unknown functions 
$p_j(y)$. We approximate solutions to (2.20) by writing
\begin{equation} p_j(y) \approx \sum_{p=0}^{2Q+1} a_p^{(j)} w_p(y/b_j), \end{equation}
where 
$Q$ is a truncation parameter, 
$a_p^{(j)}$ are designated unknown expansion coefficients and
\begin{equation} w_p(u) = \frac{{\mathrm{e}}^{{\rm i} {\rm \pi}p/2} }{(p+1) {\rm \pi}} \sqrt{1-u^2} \mathrm{U}_p(u) \end{equation}
are expansion functions where 
$\mathrm {U}_p({\cdot })$ represents the Chebyshev polynomial of the second kind. We note the relation (see Gradshtyen & Ryhzik Reference Gradshtyen and Ryhzik1965, 10§ 3.715 (13), (18))
\begin{equation} D_p(\lambda) = \int_{{-}1}^1 w_p(u) {\mathrm{e}}^{-{\rm i} \lambda u} \, {\rm d} u = \left\{\begin{array}{@{}ll}
{\rm J}_{p+1}(\lambda)/\lambda, & \lambda \neq 0,\\
\dfrac12 \delta_{p0}, & \lambda = 0, \end{array} \right. \end{equation}
where 
${\rm J}_p({\cdot })$ is a Bessel function of order 
$p$ whilst 
$\delta$ represents the Kronecker delta. The representation (2.21) thus accounts explicitly for the anticipated square root behaviour in 
$p_j(y)$ as 
$|y| \to b_j^-$. We implement Galerkin's method which involves substituting (2.21) into (2.20) before multiplying by the conjugate function 
$w_q^*(y/b_k)$ and integrating over 
$|y| < b_k$, where the asterisk 
$*$ denotes the complex conjugate. This results in the following system of equations for the expansion coefficients:
\begin{equation} \sum_{p=0}^{2Q+1} \sum_{j=0}^N a_p^{(j)} K_{pq}^{(jl)} ={-}{\rm i} \alpha_0 b_l {\mathrm{e}}^{{\rm i} \alpha_0 x_l} D_q(\beta_0 b_l), \quad q=0,\ldots,2Q+1,\ l=0,\ldots,N, \end{equation}where
\begin{equation} K_{pq}^{(jl)} = \frac{b_j b_l}{4{\rm \pi}} \int_{-\infty}^\infty \gamma {\mathrm{e}}^{-\gamma|x_j - x_l|} D_p(\beta b_j) D_q(\beta b_l) \, {\rm d} \beta. \end{equation}
Computational savings are available by making further manipulations which, in part, reflect the symmetry about 
$y=0$ of the geometry and, in part, exploit the logarithmic singularity that is embedded in the formulation despite us having avoided the use of Green's functions. We note that 
$D_p(\lambda ) = (-1)^p D_p(-\lambda )$ whilst 
$\gamma$ is symmetric in 
$\beta$ with 
$\gamma \sim |\beta |$ as 
$\beta \to \pm \infty$. Furthermore, we note an orthogonality relation for Bessel functions (Gradshtyen & Ryhzik Reference Gradshtyen and Ryhzik1965, 10§ 6.5382(2))
\begin{equation} \int_0^\infty \frac{{\rm J}_{2p+1+\nu}(u) {\rm J}_{2q+1+\nu}(u)}{u} \, {\rm d} u = \frac{1}{4p+2\nu+2}\delta_{pq}, \end{equation}
for 
$\nu = 0,1$. Taken together, this allows the original system (2.24) to be decoupled into the pair of second-kind systems of equations
\begin{align}
&\frac{1}{2{\rm \pi}} \frac{a_{2q+\nu}^{(l)}}{4q+2 \nu + 2} + \sum_{p=0}^Q \sum_{j=0}^N a_{2p+\nu}^{(j)} \hat{K}_{2p+\nu, 2q+\nu}^{(jl)} \nonumber\\
&\quad={-}{\rm i} \alpha_0 b_l {\mathrm{e}}^{{\rm i} \alpha_0 x_l} D_{2q+\nu}(\beta_0 b_l), \quad \begin{cases} q=0,\ldots,Q, \\ l=0,\ldots,N, \end{cases} \end{align}
(
$\nu = 0,1$ encode symmetric and antisymmetric components) where, for 
$l \neq j$,
\begin{equation} \hat{K}_{2p+\nu, 2q+\nu}^{(jl)} = \frac{b_j b_l}{2 {\rm \pi}} \int_0^\infty \gamma {\mathrm{e}}^{-\gamma |x_j-x_l|} D_{2p+\nu}(\beta b_j) D_{2q+\nu}(\beta b_l) \, {\rm d} \beta \end{equation}
are dimensionless exponentially convergent integrals whilst, for 
$j=l$,
\begin{equation} \hat{K}_{2p+\nu, 2q+\nu}^{(jj)} = \frac{b_j^2}{2 {\rm \pi}} \int_0^\infty ( \gamma - \beta) D_{2p+\nu}(\beta b_j) D_{2q+\nu}(\beta b_j) \, {\rm d} \beta \end{equation}
contain oscillatory integrands whose amplitude decays as 
$O(1/\beta ^3)$ accelerated from a 
$O(1/\beta )$ decay in the original system (2.24) with (2.25). Furthermore, we have
\begin{equation} \hat{K}_{2p+\nu, 2q+\nu}^{(jl)} = \hat{K}_{2q+\nu, 2p+\nu}^{(lj)}. \end{equation} We note that in the special arrangement 
$x_j = j c$ and 
$b_j = b$, representative of a rectangular metastructure with regular spacing between array elements,
\begin{equation} \hat{K}_{2p+\nu, 2q+\nu}^{(jl)} = \frac{b^2}{2 {\rm \pi}} \int_0^\infty \gamma {\mathrm{e}}^{-\gamma |\,j-l|c} D_{2p+\nu}(\beta b) D_{2q+\nu}(\beta b) \, {\rm d} \beta \end{equation}
depends only on 
$|\,j-l| = 0,\ldots, N$ and requires only 
$N+1$ integrals for each 
$(p,q)$ pair, rather than 
$(N+1)(N+2)/2$ evaluations. Computation of the elements of the matrix system is thus an 
$O(N)$ task rather than 
$O(N^2)$ for this special case. For a matrix with 
$N\times N$, the inversion of a Toeplitz matrix, though reduced from 
$O(N^3)$ to 
$O(N^2)$, still remains a limiting factor as 
$N$ becomes very large.
The values of 
$a_{p}^{(j)}$ are numerically determined from the solution of (2.27) where, typically, a value of 
$Q = 5$ is sufficient for convergence to five or more decimal places unless the frequency is high when 
$Q$ must be increased. Subsequently, this allows 
$\phi$ to be determined everywhere by using
\begin{equation} \phi(x,y) = \phi_{inc}(x,y) + \sum_{l=0}^{N} \sum_{p=0}^{2Q+1} a_{p}^{(l)} \varLambda_{p}^{(l)}(x,y), \end{equation}
where 
$\varLambda _{p}^{(l)}(x,y)$ can be alternatively expressed as
\begin{equation} \varLambda_{p}^{(l)}(x,y) ={-}\frac{b_{l}}{4{\rm \pi}} \int_{-\infty}^{\infty} \mbox{sgn}(x-x_{l})\, D_{p}(\beta b_{l}) {\mathrm{e}}^{-\gamma|x-x_{l}|+{\rm i}\beta y} {\rm d} \beta \end{equation}or
\begin{equation} \varLambda_{p}^{(l)}(x,y) ={-}\frac{{\rm i}}{4} \int_{{-}b_{l}}^{b_{l}} \frac{k(x-x_l)}{\varrho} H_{1}(k\varrho) w_{p}(y'/b_{l}) {\rm d} y', \end{equation}
where the expression (2.34) has applied the integral representation of Hankel function, see the Appendix for details. In the computation of wave field, (2.33) is used when 
$|x-x_k|>\epsilon$ due to the exponential decay factor, and expression (2.34) is adopted otherwise.
We have a particular interest in the diffraction coefficient which may be calculated from (2.20) using 
$x = r \cos \theta$, 
$y = r \sin \theta$ and employing a stationary phase approximation following the parametrisation of 
$\beta \in (-\infty,\infty )$ as 
$\beta = k \sin \psi$ for 
$(-{\rm \pi} /2,{\rm \pi} /2)$ and 
$\beta = \pm k \cosh u$ for 
$u \in (0,\theta )$ via the relationship 
$\psi = \pm {\rm \pi}/2 \mp {\rm i} u$. In the limit 
$kr \to \infty$ the dominant contribution to the far field comes from the integral over 
$-{\rm \pi} /2 < \psi < {\rm \pi}/2$ at 
$\psi = \theta$ or 
$\psi = \theta +{\rm \pi}$ depending on the value of 
$\theta$. Within this branch, 
$\gamma = -{\rm i} \alpha = -{\rm i} \cos \psi$ and it is the negative sign of the branch, chosen earlier, that dictates that the scattered waves are outgoing. After some algebra we find
\begin{equation} A(\theta;\theta_0) \approx{-}\frac{k\cos\theta}{4}\sum_{l=0}^N {\mathrm{e}}^{-{\rm i} k x_l \cos \theta} \sum_{p=0}^{2Q+1} a_p^{(l)} b_{l} D_p (k b_l \sin \theta) \end{equation}
and the dependence on 
$\theta _0$ is embedded in the coefficients 
$a_p^{(l)}$ whose values are determined by the incident wave forcing in (2.27). We note that the diffraction coefficient satisfies the so-called optical theorem (Maruo Reference Maruo1960)
\begin{equation} \sigma = \frac{1}{2{\rm \pi}} \int_0^{2 {\rm \pi}} |A(\theta;\theta_0)|^2 \, {\rm d} \theta ={-}{\mathrm{Re}}\left[ A(\theta_0;\theta_0) \right] \end{equation}and represents the total scattering cross-section or scattering energy.
We are also interested in the total hydrodynamic force in the 
$x$-direction of the 
$j$th plate in the array which is proportional to
\begin{equation} F_x^{(j)} ={-}{\rm i} \omega \rho \int_{{-}h}^0 \psi_0(z) \int_{{-}b_j}^{b_j} p_j(y) \, {\rm d} y \, {\rm d} z \approx{-}{\rm i} \omega \rho \frac{N_0^{{-}1/2} \sinh k h}{2 k} a_0^{(j)} b_j. \end{equation}3. An infinite periodic array of plate array metastructures
We assume now that the metastructure considered in the previous section is repeated periodically in the 
$y$-direction with spacing between a reference point within adjacent identical structures given by 
$2d$. This is commonly referred to as the scattering of oblique waves by a periodic diffraction grating as described in the context of plate-array metastructures by Putley et al. (Reference Putley, Guenneau, Porter and Craster2022). When 
$\theta _0 = 0$ the periodicity allows the problem to be interpreted as geometrically equivalent to the reflection and transmission of incident waves by a single metastructure on the centreline of a uniform channel of width 
$2d$ with impermeable walls. However, we retain the generality of oblique incidence here and demonstrate that both the solution method and numerical procedure are very similar to that encountered in the open domain problem considered in the previous section. The usual arguments for plane wave scattering by a periodic grating follow. Thus, since 
$\phi _{inc}(x,y+2d) = {\mathrm {e}}^{2 {\rm i} \beta _0 d} \phi _{inc}(x,y)$ with 
$\beta _0 = k \sin \theta _0$ as before it also must follow that 
$\phi (x,y+2d) = {\mathrm {e}}^{2 {\rm i} \beta _0 d} \phi (x,y)$ and this allows one to consider the scattering problem in a fundamental cell, say 
$y \in [-d,d]$, 
$-\infty < x < \infty$ provided we also impose periodic boundary conditions on the lateral edges of the cell, these being (Porter & Evans Reference Porter and Evans1996)
\begin{equation} \phi(x,d) = {\mathrm{e}}^{2 {\rm i} \beta_0 d} \phi(x,-d), \quad \text{and} \quad \phi_y(x,d) = {\mathrm{e}}^{2 {\rm i} \beta_0 d} \phi_y(x,-d). \end{equation}
The extension to 
$y \not \in [-d,d]$ is provided by 
$\phi (x,y + 2m d) = {\mathrm {e}}^{2{\rm i} \beta _0 m d}\phi (x,y)$ for 
$m \in \mathbb {Z}$. As well as restricting the domain to a strip of width 
$2d$, the far-field conditions also change to
\begin{equation} \phi(x,y)- \phi_{inc}(x,y) \sim \sum_{n={-}n_-}^{n_+} R_n {\mathrm{e}}^{-{\rm i} \alpha_n x} {\mathrm{e}}^{{\rm i} \beta_n y}, \quad kx \to -\infty \end{equation}and
\begin{equation} \phi(x,y) \sim \sum_{n={-}n_-}^{n_+} T_n {\mathrm{e}}^{{\rm i} \alpha_n x} {\mathrm{e}}^{{\rm i} \beta_n y}, \quad kx \to \infty, \end{equation}
where 
$R_n$, 
$T_n$ are complex-valued reflection and transmission coefficients,
\begin{equation} \beta_n = \beta_0 + n {\rm \pi}/d, \quad n \in \mathbb{Z} \end{equation}and
\begin{equation} \alpha_n = \sqrt{k^2 - \beta_n^2}, \quad -n_-\leq n \leq n_+ \end{equation}
are real wavenumber components with 
$\alpha _0 = k \cos \theta _0$ as before and
\begin{equation} n_-= \lfloor k d (1+\sin\theta_0)/{\rm \pi} \rfloor, \quad n_{+}= \lfloor k d (1-\sin\theta_0)/{\rm \pi} \rfloor \end{equation}define the number of propagating diffracted modes (Porter & Evans Reference Porter and Evans1996). We choose to write
\begin{equation} \gamma_n = \sqrt{\beta_n^2 - k^2} \equiv{-} {\rm i} \alpha_n \end{equation}
such that 
$\gamma _n$ is real if 
$n \not \in [-n_-,n_+]$. The notation and definition mimic (2.11) and we are ready to follow the methods of the previous section. Thus, we define the Fourier transform over a finite interval
\begin{equation} \bar{\phi}_n(x) = \frac{1}{2d} \int_{{-}d}^d [\phi(x,y) - \phi_{inc}(x,y)] {\mathrm{e}}^{-{\rm i} \beta_n y} \, {\rm d} y, \end{equation}
for 
$n \in \mathbb {Z}$ and the inverse
\begin{equation} \phi(x,y) = \phi_{inc}(x,y) + \sum_{n={-}\infty}^\infty \bar{\phi}_n(x) {\mathrm{e}}^{{\rm i} \beta_n y} \end{equation}which follows from the orthogonality relation
\begin{equation} \frac{1}{2d} \int_{{-}d}^d {\mathrm{e}}^{{\rm i} \beta_m y} {\mathrm{e}}^{-{\rm i} \beta_n y} \, {\rm d} y = \delta_{mn}. \end{equation}The governing Helmholtz equation is reduced to
\begin{equation} \left( \frac{{\rm d}^2}{{\rm d}\kern0.7pt x^2} - \gamma_n^2 \right) \bar{\phi}_n = 0, \quad x \neq x_j, \ (j=0,\ldots,N) \end{equation}
and the transform of continuity of 
$\phi _x(x,y)$ at 
$x=x_j$ for all 
$y \in [-d,d]$ is expressed as
\begin{equation} \frac{\partial}{\partial x} \bar{\phi}_n(x_j^+) - \frac{\partial}{\partial x} \bar{\phi}_n(x_j^-) = 0, \quad j=0,\ldots,N. \end{equation}Likewise, we readily find that
\begin{equation} \bar{\phi}_n(x_j^+) - \bar{\phi}_n(x_j^-) = P_{n,j}, \quad j=0,\ldots,N \end{equation}where
\begin{equation} P_{n,j} = \frac{1}{2d} \int_{{-}b_j}^{b_j} p_j(y) {\mathrm{e}}^{-{\rm i} \beta_n y} \, {\rm d} y \end{equation}
and 
$\phi (x_j^+,y) - \phi (x_j^-,y) = p_j(y)$ for 
$|y| < b_j$ and is zero for 
$b_j < |y| < d$. With reference to the approach outlined in the previous section the transform solution can now clearly be written as
\begin{equation} \bar{\phi}_n(x) = \sum_{j=0}^N P_{n,j} g_n(x,x_j), \end{equation}
where 
$g_n(x,x_j)$ satisfies (3.11), has continuous 
$x$-derivative at 
$x=x_j$, has a jump of unity in its value from 
$x_j^+$ to 
$x_j^-$ and is outgoing at infinity for 
$n \in [-n_-,n_+]$ and exponentially decaying towards infinity otherwise. This gives
\begin{equation} g_n(x,x_j) ={-} \frac12 \mbox{sgn}(x-x_j) {\mathrm{e}}^{-\gamma_n |x-x_j|} \end{equation}and so the solution in physical space is
\begin{equation} \phi(x,y) = \phi_{inc}(x,y) - \frac{1}{4d} \sum_{j=0}^N \sum_{n={-}\infty}^\infty \mbox{sgn}(x-x_j) {\mathrm{e}}^{-\gamma_n |x-x_j|} {\mathrm{e}}^{{\rm i} \beta_n y} \int_{{-}b_j}^{b_j} p_j(y') {\mathrm{e}}^{- {\rm i} \beta_n y'} \, {\rm d} y'. \end{equation}
By comparing (3.17) with (3.2) and (3.3) in the limits 
$kx \to -\infty$ and 
$kx \to + \infty$, respectively, and we can deduce simply that
\begin{equation} R_n = \frac{1}{4d} \sum_{j=0}^N {\mathrm{e}}^{{\rm i} \alpha_n x_j} \int_{{-}b_j}^{b_j} p_j(y') {\mathrm{e}}^{- {\rm i} \beta_n y'} \, {\rm d} y' \end{equation}and
\begin{equation} T_n = \delta_{n,0} - \frac{1}{4d} \sum_{j=0}^N {\mathrm{e}}^{-{\rm i} \alpha_n x_j} \int_{{-}b_j}^{b_j} p_j(y') {\mathrm{e}}^{- {\rm i} \beta_n y'} \, {\rm d} y', \end{equation}
for 
$-n_- \leq n \leq n_+$.
Coupled integral equations for the unknowns 
$p_j(y)$ are constructed by applying the barrier conditions (2.7) at 
$x=x_l$, so that
\begin{equation} \frac{1}{4d} \sum_{j=0}^N \sum_{n={-}\infty}^\infty \gamma_n {\mathrm{e}}^{-\gamma_n|x_j - x_l|} {\mathrm{e}}^{{\rm i} \beta_n y} \int_{{-}b_j}^{b_j} p_j(y') {\mathrm{e}}^{-{\rm i} \beta_n y'} \, {\rm d} y' ={-}{\rm i} \alpha_0 {\mathrm{e}}^{{\rm i} \alpha_0 x_l} {\mathrm{e}}^{{\rm i} \beta_0 y}, \quad |y| < b_l \end{equation}
and 
$l = 0,\ldots,N$. This equation is the analogue of (2.20) in the open domain case: infinite integrals over continuous variables 
$\beta$ are replaced by infinite sums over discrete variables 
$\beta _n$. The approximation to the integral equations follows as in the previous section and the final system of equations that need to be solved in this problem remains (2.24) but with
\begin{equation} K_{pq}^{(jl)} = \frac{b_j b_l}{4d} \sum_{n={-}\infty}^\infty \gamma_n {\mathrm{e}}^{-\gamma_n |x_j-x_l|} D_{p}(\beta_n b_j) D_{q}(\beta_n b_l) \end{equation}
with 
$D_p(\lambda )$ still defined by (2.23).
It follows that
\begin{equation} R_n \approx \sum_{j=0}^N \frac{b_j}{4d} {\mathrm{e}}^{{\rm i} \alpha_n x_j} \sum_{p=0}^{2Q+1} a_p^{(j)} D_{p}(\beta_n b_j), \end{equation}and
\begin{equation} T_n \approx \delta_{n,0}- \sum_{j=0}^N \frac{b_j}{4d} {\mathrm{e}}^{-{\rm i} \alpha_n x_j} \sum_{p=0}^{2Q+1} a_p^{(j)} D_{p}(\beta_n b_j), \end{equation}
for 
$-n_- \leq n \leq n_+$. These reflection and transmission coefficients satisfy the conservation of energy condition (see, e.g. Porter & Evans Reference Porter and Evans1996)
\begin{equation} E_{R}+E_T=1 \quad\text{with}\ E_{R}= \sum_{n={-}n_-}^{n_+} \frac{\alpha_n}{\alpha_0}|R_{n}|^2 \quad\text{and}\quad E_{T}= \sum_{n={-}n_-}^{n_+}\frac{\alpha_n}{\alpha_0}|T_{n}|^2, \end{equation}
where 
$E_{R}$ and 
$E_T$ represent total reflected and transmitted energy, respectively.
4. Arrays of partially submerged surface-piercing barriers
In order to showcase the method further, we consider a different type of problem which is still geometrically two-dimensional. An array of 
$N+1$ vertical barriers is assumed to extend indefinitely and uniformly in the 
$y$-direction and, instead of extending fully through the depth of the fluid, are truncated. Thus, the barrier at 
$x = x_j$ occupies 
$-\infty < y < \infty$, and 
$-b_j < z < 0$, with 
$b_j < h$ (
$\,j=0,\ldots,N$), as in figure 2. We remark that 
$b_j$ now denotes the full length of the plate that has previously been represented by 
$2b_j$ a choice made to connect with earlier sections. We retain the generality of oblique incidence of incoming surface waves and, although we can no longer trivially factorise out the depth dependence, the uniformity of the barriers in 
$y$ allows us to write
\begin{equation} \varPhi(x,y,z,t) = {\mathrm{Re}} [ \phi(x,z) {\mathrm{e}}^{{\rm i} \beta_0 y} {\mathrm{e}}^{- {\rm i} \omega t} ], \end{equation}
where 
$\beta _0 = k \sin \theta _0$ is the component of the wavenumber aligned with the 
$y$-axis. Now the problem is given by
\begin{equation} \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial z^2} -\beta_0^2 \right) \phi = 0 \end{equation}with
\begin{equation} \phi_z = 0, \quad \text{on }z={-}h \end{equation}and
\begin{equation} \phi_z - K \phi = 0, \quad \text{on }z=0 \end{equation}along with
\begin{equation} \phi_x = 0, \quad \text{on }x=x_j^\pm, \quad -b_j < z < 0\ (j=0,\ldots,N). \end{equation}Within this revised framework an obliquely incident wave is described by the potential
\begin{equation} \phi_{inc}(x,z) = {\mathrm{e}}^{{\rm i} \alpha_0 x} \psi_0(z), \end{equation}
where 
$\alpha _0 = k \cos \theta _0$. The conditions in the far field are
\begin{equation} \phi(x,z) -\phi_{inc}(x,z) \sim \left\{ \begin{array}{@{}ll}
R {\mathrm{e}}^{-{\rm i} \alpha_0 x} \psi_0(z), & kx \to -\infty,\\
(T-1) {\mathrm{e}}^{{\rm i} \alpha_0 x} \psi_0(z), & kx \to \infty, \end{array} \right. \end{equation}
where 
$R$ and 
$T$ are reflection and transmission coefficients, respectively; 
$\phi - \phi _{inc}$ is outgoing of course. We solve the problem above by first defining orthonormal depth eigenfunctions for a domain without barriers as (e.g. Linton & McIver Reference Linton and McIver2001)
\begin{equation} \psi_n(z) = N_n^{{-}1/2} \cos k_n (z+h), \quad N_n = \frac12 \left( 1 + \frac{\sin 2 k_n h }{2 k_n h} \right), \end{equation}
for 
$n \geq 1$ and 
$k_n$ are an increasing sequence of real positive roots of
\begin{equation} K ={-} k_n \tan k_n h. \end{equation}
We can extend the definition to 
$n=0$ by letting 
$k_0 = -{\rm i} k$ and then
\begin{equation} \frac{1}{h} \int_{{-}h}^0 \psi_n(z) \psi_m(z) \, {\rm d} z = \delta_{mn}, \end{equation}
for all 
$m,n = 0,1,\ldots$.
Figure 2. Sketch of wave scattering by an array of surface-piercing barriers.
We write
\begin{equation} \bar{\phi}_n(x) = \frac{1}{h} \int_{{-}h}^0 [\phi(x,z) - \phi_{inc}(x,z)] \psi_n(z) \,{\rm d} z \end{equation}such that
\begin{equation} \phi(x,z) = \phi_{inc}(x,z) + \sum_{n=0}^\infty \bar{\phi}_n(x) \psi_n(z) \end{equation}follows from (4.11) and (4.10). It follows that
\begin{equation} \left( \frac{{\rm d}^2}{{\rm d}\kern0.7pt x^2} - \gamma_n^2 \right) \bar{\phi}_n(x) = 0, \quad x \neq x_j, \ (j=0,\ldots,N), \end{equation}where, now,
\begin{equation} \gamma_n = \sqrt{k_n^2 + \beta_0^2} \end{equation}
is real for 
$n \geq 1$ but, for 
$n=0$, 
$\gamma _0 = - {\rm i} \alpha _0$.
We note that 
$\phi _x$ is continuous everywhere including across 
$x=x_j$ for all 
$-h < z < 0$ and so it follows that
\begin{equation} \frac{\partial}{\partial x} \bar{\phi}_n(x_j^+) = \frac{\partial}{\partial x} \bar{\phi}_n(x_j^-). \end{equation}
Defining 
$p_j(z) = \phi (x_j^+,z) - \phi (x_j^-,z)$ which is zero for 
$-h < z < -b_j$ means that
\begin{equation} \bar{\phi}_n(x_j^+) - \bar{\phi}_n(x_j^-) = P_{n,j} \equiv \frac{1}{h} \int_{{-}b_j}^0 p_j(z) \psi_n(z) \, {\rm d} z \end{equation}
represents the ‘depth transform’ of the pressure jump across the 
$j$th barrier. With reference to the two preceding sections, we are immediately able now to write down the transform solution as
\begin{equation} \bar{\phi}_n(x) ={-} \frac12 \sum_{j=0}^N P_{n,j} \mbox{sgn}(x-x_j) {\mathrm{e}}^{-\gamma_n |x-x_j|}, \end{equation}and we can confirm this satisfies all the conditions above. Thus,
\begin{equation} \phi(x,z) = \phi_{inc}(x,z) - \frac{1}{2h} \sum_{j=0}^N \mbox{sgn}(x-x_j) \sum_{n=0}^\infty {\mathrm{e}}^{-\gamma_n |x-x_j|} \psi_n(z) \int_{{-}b_j}^0 p_j(z') \psi_n(z') \, {\rm d} z' \end{equation}
is the general solution, expressed in terms of the unknown functions 
$p_j(z')$. We take the limit 
$k x \to \pm \infty$ in the above, comparing with (4.7) to get
\begin{equation} R = \frac{1}{2h} \sum_{j=0}^N {\mathrm{e}}^{{\rm i} k x_j} \int_{{-}b_j}^0 p_j(z') \psi_0(z') \, {\rm d} z' \end{equation}and
\begin{equation} T = 1 - \frac{1}{2h} \sum_{j=0}^N {\mathrm{e}}^{-{\rm i} k x_j} \int_{{-}b_j}^0 p_j(z') \psi_0(z') \, {\rm d} z'. \end{equation}
The unknowns 
$p_j(z)$ are determined by imposing the remaining no-flow conditions (4.5) on 
$x = x_l$ to give
\begin{equation} \frac{1}{2h} \sum_{j=0}^N \sum_{n=0}^\infty \gamma_n {\mathrm{e}}^{-\gamma_n |x_j - x_l|} \psi_n(z) \int_{{-}b_j}^0 p_j(z') \psi_n(z') \, {\rm d} z' ={-}{\rm i} \alpha_0 {\mathrm{e}}^{{\rm i} \alpha_0 x_l} \psi_0(z), \quad -b_l < z < 0, \end{equation}
for 
$l=0,\ldots,N$. The coupled integral equations are solved using the method first described in Porter & Evans (Reference Porter and Evans1995) in which
\begin{equation} p_j(z) \approx \sum_{p=0}^Q a_p^{(j)} w_p(z/b_j) \end{equation}and
\begin{equation} \hat{w}_p(u) = w_p(u) - K b_j \int_{{-}1}^u w_p(s) \, {\rm d} s, \end{equation}where
\begin{equation} \hat{w}_p(u) = \frac{2 ({-}1)^p}{(2p+1) {\rm \pi}} \sqrt{1-u^2} \mathrm{U}_{2p}(u) \end{equation}is designed to ensure that the free surface condition (4.4) is satisfied as well as retaining the correct local square root behaviour of the pressure jump in the vicinity of the lower edge of the plates. It follows that (Porter & Evans Reference Porter and Evans1995)
\begin{equation} D_{np}^{(j)} = \int_{{-}b_j}^0 \psi_n(z) w_p(z/b_j) \, {\rm d} z = N_n^{{-}1/2} \cos (k_n h) \int_{{-}b_j}^0 \cos(k_n z) \hat{w}_p(z/b_j) \, {\rm d} z \end{equation}after integrating by parts, is given by
\begin{equation} D_{np}^{(j)} = N_n^{{-}1/2} \cos (k_n h) {\rm J}_{2p+1}(k_n b_j)/(k_n b_j) \end{equation}
which, for 
$n = 0$, is better expressed as
\begin{equation} D_{0p}^{(j)} = ({-}1)^p N_0^{{-}1/2} \cosh (k h) {\rm I}_{2p+1}(k b_j)/(k b_j), \end{equation}
where 
${\rm I}_p({\cdot })$ is a modified Bessel function of the first kind of order 
$p$. Substituting (4.22) into (4.21), and multiplying through by 
$w_{q}(z/b_l)$ before integrating over 
$-b_l < z < 0$ gives the system of equations
\begin{equation} \sum_{j=0}^N \sum_{p=0}^Q a_{p}^{(j)} K_{pq}^{(jl)} ={-} {\rm i} \alpha_0 {\mathrm{e}}^{{\rm i} \alpha_0 x_l} D_{0q}^{(k)}, \quad l=0,\ldots,N,\ q = 0,\ldots,Q, \end{equation}where
\begin{equation} K_{pq}^{(jl)} = \frac{b_{j}b_{l}}{2h} \sum_{n=0}^\infty \gamma_n {\mathrm{e}}^{-\gamma_n |x_j - x_l|} D_{np}^{(j)} D_{nq}^{(l)}. \end{equation}
For 
$j\neq l$ the series is exponentially convergent. When 
$j=l$, the series defining 
$K_{pq}^{(jj)}$ resembles that encountered in Porter & Evans (Reference Porter and Evans1995) for a plate in isolation in which terms decay like 
$O(1/n^2)$. It is possible to accelerate the convergence of the series defining 
$K_{pq}^{(jj)}$ by subtracting the leading-order asymptotic behaviour of each term in the series which can be deduced from 
$k_n h \sim n {\rm \pi}$, 
$N_n \sim \tfrac 12$, 
$\gamma _n h \sim n {\rm \pi}$ as 
$n \to \infty$. The infinite series which compensates for the subtraction can then be evaluated as a different infinite series (see Paris Reference Paris2018) which, for the present purposes, is not worth pursuing.
In the case that plates are positioned at regular intervals, 
$x_j = jc$, with spacing 
$c$ and submerged to the same depth, 
$b_j = b_0 = b$, which corresponds to the case considered by Huang & Porter (Reference Huang and Porter2023) then
\begin{equation} K_{pq}^{(jl)} = \frac{b^2}{2h} \sum_{n=0}^\infty \gamma_n {\mathrm{e}}^{-\gamma_n |\,j - l|c} D_{np}^{(0)} D_{nq}^{(0)} \end{equation}
depends only on 
$|\,j-l|$ and only needs 
$N+1$ evaluations for 
$|\,j-l| = 0,\ldots,N$.
Using (4.22) in (4.19) and (4.20) gives
\begin{equation} R \approx \sum_{j=0}^N \frac{b_j}{2h} {\mathrm{e}}^{{\rm i} k x_j} \sum_{p=0}^Q a_p^{(q)} D_{0p}^{(j)} \end{equation}and
\begin{equation} T \approx 1- \sum_{j=0}^N \frac{b_j}{2h} {\mathrm{e}}^{-{\rm i} k x_j} \sum_{p=0}^Q a_p^{(q)} D_{0p}^{(j)} \end{equation}
and these coefficients should satisfy 
$|R|^2 + |T|^2 = 1$.
5. Results in open domain
5.1. A circular cylinder
We first consider the scattering of waves by a circular metacylinder, as first studied by Zheng et al. (Reference Zheng, Porter and Greaves2020) and later by Putley et al. (Reference Putley, Guenneau, Porter and Craster2022). Both used homogenisation to replace the discrete plate array with an effective medium. The present work allows us to validate the numerical method described in this paper by demonstrating convergence to the homogenisation results as 
$N$, the number of plates in the discrete array, increases. Figure 3 depicts the scattering energy 
$\sigma$, defined in (2.36), as a function of the non-dimensional wavenumber 
$ka$ under the oblique wave excitation (
$\theta _{0}=45^{\circ }$), where 
$a$ denotes the radius of the metacylinder. We present curves associated with metacylinders having 
$N=10$, 
$15$ and 
$20$ channels of constant width which can be seen to converge to the results of Zheng et al. (Reference Zheng, Porter and Greaves2020) (the homogenisation results have been obtained by truncating their numerical system of equations at 
$20$ terms) as 
$N$ increases for 
$ka < {\rm \pi}/2$. The vertical line corresponds to 
$ka={\rm \pi} /2$ which signals the onset of fluid resonance in narrow channels and the homogenisation method fails for 
$ka$ beyond this value (Putley et al. Reference Putley, Guenneau, Craster, Davies and Poulton2023). Our method therefore allows us to consider results for 
$ka > {\rm \pi}/2$. A general observation is that larger 
$N$ are required for convergence as the frequency increases and that the scattering energy generally increases with the wavenumber and exhibits oscillations near integer multiples of 
${\rm \pi} /2$, representing the onset of new gap resonance modes in the central channel (Molin et al. Reference Molin, Remy, Kimmoun and Stassen2002). It is noteworthy that the wavenumbers 
$ka=n{\rm \pi} /2$ with 
$n\in \mathbb {Z}^+$ for gap resonance in the central channel are determined by the assumption of homogeneous Dirichlet conditions 
$\phi =0$ at the ends of the channel. However, this assumption holds true only if the gap width is very small (Liang et al. Reference Liang, Zheng, Shao, Cong and Greaves2023).
Figure 3. Scattering energy 
$\sigma$ by circular metacylinders with different number of channels 
$N$ under the quartering wave excitation 
$\theta _0=45^{\circ }$ as a function of non-dimensional wavenumber 
$ka$. Comparison is made with the homogenisation solution by Zheng et al. (Reference Zheng, Porter and Greaves2020) which is valid when 
$ka<{\rm \pi} /2$.
In figure 4 we compare the results of figure 3 for 
$N=20$ channels of uniform width with a distribution of the plates within the metacylinder which maintains a constant aspect ratio of channel width to (mean) length. This new scheme therefore concentrates plates towards the two extremes of the cylinder. Although there are only small differences, the uniform width case is found to marginally improve convergence to the 
$N = \infty$ limit.
Figure 4. Comparison of scattering energy by circular metacylinders composed of 
$N=20$ channels for different plate separations constrained by constant channel aspect ratio and equal spacing. Comparison is made with the homogenisation solution valid for 
$ka<{\rm \pi} /2$.
This observation is made clearer in figure 5 where a comparison of the effect of plate distribution and the value of 
$N$ on the free surface is presented. A wave incident from 
$\theta _0 = 45^\circ$ at frequencies determined by 
$ka = 1$ (figure 5a,d,g,j,m), 
$2$ (figure 5b,e,h,k,n) and 
$3$ (figure 5c,f,i,l,o). In figures 5(a–c) and 5(g–i), the channel spacing is uniform and there are 
$N=10$, 
$N=20$ channels, respectively. In figures 5(d–f) and 5(j–l) 
$N=10$, 
$N=20$ once again but the plate distribution maintains a constant channel aspect ratio. Figure 5(m–o) shows results from homogenisation. Note that the final two results for 
$ka=2$, 
$ka=3$ are invalid since there is resonance inside the cylinder which violates the homogenisation assumptions. The plot shows more significant differences in the results for different spacing schemes at higher frequencies. We also note the presence of large local resonance within the cylinder, and the wave amplitude displayed is saturated to 
$2.0$.
Figure 5. Modulus of wave patterns scattered by a circular metacylinder for different number of plates and separation strategies. The wave patterns associated with 10 channels uniform spacing (a–c), 10 channels constant aspect ratio (d–f), 20 channels uniform spacing (g–i), 20 channels constant aspect ratio (j–l) and homogenisation solution (m–o) are exhibited for 
$ka=1.0$ (a,d,g,j,m), 
$2.0$ (b,e,h,k,n) and 
$3.0$ (c,f,i,l,o).
5.2. Rectangular and graded metawedge
As a sequel to the study on circular metacylinders, we now investigate wave scattering by metarectangles and graded metawedges, which have been less explored in the literature. Figure 6 presents the instantaneous wave patterns at 
$t=0$ scattered by a metarectangle with a width of 
$2b$ for different values of aspect ratio (
$AR$), which is defined as the ratio of the length to the width of the metarectangle, including 
$AR=1.0$ and 
$AR=5.0$, shown in figures 6(a,b) and 6(c,d), respectively. The channel width for both metarectangles is 
$c/b=0.1$. Wave patterns for 
$kb={\rm \pi} /2$ and 
$kb={\rm \pi}$ are presented in figures 6(a,c) and 6(b,d).
Figure 6. Instantaneous wave patterns at 
$t=0$ scattered by a rectangular metacylinder for different aspect ratios at 
$kb={\rm \pi} /2$ (a,c) and 
$kb={\rm \pi}$ (b,d) under the quartering wave excitation 
$\theta _0=45^{\circ }$. Panels (a,b) and (c,d) show the results for 
${AR}=1.0$ and 
$5.0$, respectively.
For the metasquare (
$AR=1.0$), shown in figure 6(a,b), the symmetrical property with respect to 
$y=x$ is disrupted due to the presence of channels. Notably, wave resonance in the channel on the upwave side is observed at 
$kb={\rm \pi}$. In the case of an elongated metarectangle (
$AR=5.0$), depicted in figure 6(c,d), large free surface responses are observed in the first channel facing the wave incidence. Besides, there is a noticeable wave twisting within the metarectangle, similar to the phenomenon described by Porter (Reference Porter2021) for an infinite setting. Unlike the perfect transmission reported in Porter (Reference Porter2021), however, the presence of end effects leads to appreciable disturbances riding on the wave crest/trough.
In figure 7, we consider the diffraction energy 
$\sigma$ under the normal wave incidence 
$\theta _0=0^{\circ }$ for a metasquare and a metawedge, depicted in figure 7(a) and figure 7(b), respectively. The metasquare used here is identical to the one shown in figure 6. Both the metasquare and the metawedge share the same length and are composed of 20 channels. Here we define the base ratio of the metawedge as 
$\ell =b_{N}/b_{0}$, and the mean semiwidth 
$b_{m}=(b_{0}+b_{N})/2$. When the base ratio is unequal to unity, i.e. 
$\ell \ne 1$, the constant aspect ratio separation strategy is employed in the configuration of the metawedge. The results show a good agreement between the two alternative representations provided by (2.36), thereby confirming the accuracy of the computation. In both cases, the scattering energy exhibits a step-shaped increase. For the metasquare, depicted in figure 7(a), strong oscillations occur at the beginning of the step. Although the metawedge, shown in figure 7(b), also exhibits fluctuations in the scattering energy, the oscillation amplitude is much smaller.
Figure 7. Scattering energy 
$\sigma$ under the normal wave excitation (
$\theta _0=0^{\circ }$) as a function of non-dimensional wavenumber 
$kb_{m}$ for base ratios 
$\ell =1$ (a, metasquare) and 
$\ell =3$ (b, metawedge).
Figure 8 illustrates the free surface elevation along the centreline of the metasquare (
$\ell =1$) and metawedge (
$\ell =3$) considered in figure 7, shown in figure 7(a) and figure 7(b), respectively, as a function of the normalised wavenumber 
$kb_{m}$ ranging from 
$0$ to 
$10$. The white lines indicate the locations of the plates, and the layout is identical to the set-up in figure 7.
Figure 8. Free surface elevation along the centreline of the metasquare 
$\ell =1$ (a) and metawedge 
$\ell =3$ (b) varying with the normalised wavenumber 
$kb_{m}$.
Within the metastructure, significant wave resonance accompanied by large-amplitude wave responses is observed, see figure 9. For the metasquare, wave resonance occurs at discrete frequencies, whereas for the metawedge, waves are trapped over a broad range of frequencies, demonstrating a ‘rainbow reflection’ behaviour. In both cases, the downwave side of the metastructure experiences minimal disturbance, exhibiting shielding effects, see figure 8 for 
$x > b_m$. Notably, we see from figure 8 that the metawedge provides superior shielding effects compared with the metasquare because of rainbow reflection, resulting in a large quiet region over a wider range of frequencies.
Figure 9. Demonstration of rainbow trapping by a metawedge in the 1st, 6th, 16th and 20th channels at 
$kb_{m}=2.90$, 
$2.25$, 
$1.35$, 
$1.02$, respectively. The colourbar indicates the modulus of free surface elevation.
6. Results for periodic arrays
Following the physical findings of wave scattering by a single metastructure in the open domain considered in § 5, our focus now turns to the analysis of periodic array scenarios as studied in § 3. Specifically, we aim at delving into the underlying physics of wave patterns associated with nearly total reflection and nearly perfect transmission, as predicted by the energy relation given by (3.24).
6.1. Circular metacylinder
We first study the scattering of waves by a periodic array of circular metacylinders. Figure 10 illustrates the reflected energy 
$E_{R}$, defined in (3.24), by a periodic array of circular metacylinders, with each composed of 20 channels, as a function of the non-dimensional wavenumber 
$ka$, where 
$a$ represents the radius of metacylinder. Both normal incidence (
$\theta _{0}=0^{\circ }$) and oblique incidence (
$\theta _{0}=45^{\circ }$) are presented, displayed in figure 10(a) and figure 10(b), respectively. In this configuration, half the centre-to-centre distance between adjacent metacylinders is twice the radius (
$d=2a$). In this set-up, the lowest resonant wavenumber 
$ka={\rm \pi} /2$ in the metacylinder coincides with the crossing mode wavenumber 
$kd={\rm \pi}$.
Figure 10. Reflected energy for a periodic array of circular metacylinders with 
$a/d=0.5$ for 
$\theta _{0}=0^{\circ }$ (a) and 
$\theta _{0}=45^{\circ }$ (b). The vertical line corresponds to 
$ka={\rm \pi} /2$, where 
$a$ denotes the radius of the circular metacylinder.
In figure 10(a) depicting normal incidence, we observe a sharp transition in the reflected energy. As the wavenumber approaches 
$ka={\rm \pi} /2$, the reflection changes from nearly perfect transmission (
$E_{R}\rightarrow 0$) to nearly total reflection (
$E_{R}\rightarrow 1$) occurred at 
$ka\approx 1.5036$ and 
$ka\approx 1.5707$, respectively. On the other hand, under oblique wave excitation, as in figure 10(b), specific wavenumbers exist where reflection is negligible, whereas complete reflection does not occur in this set-up.
To further elucidate the underlying physics governing the phenomena of nearly total transmission and nearly perfect reflection described in figure 10, we examine the free surface responses at these wavenumbers.
Figure 11 presents the wave patterns scattered by a circular metacylinder under the action of normal incidence (
$\theta _{0}=0^{\circ }$) at 
$ka=1.5036$ corresponding to nearly total transmission. Figures 11(a) and 11(b) show modulus and instantaneous wave patterns, respectively. It is notably observed that waves are trapped within the gaps of the plate arrays constituting the circular metacylinder, resulting in large free surface responses. Furthermore, at significant distances from the metacylinder, the wave field maintains the profile of the incident waves, indicating the occurrence of perfect transmission.
Figure 11. Wave patterns scattered by a periodic array of circular metacylinders under normal wave incidence (
$\theta _{0}=0^{\circ }$) at 
$ka=1.5036$ with a normalised radius of 
$a/d=0.5$, illustrating nearly perfect wave transmission. Panels (a,b) exhibit the modulus and real part of the wave pattern, respectively.
Figure 12 illustrates the diffraction wave field at 
$ka=1.5707$ under the head wave excitation 
$\theta _{0}=0^{\circ }$, at which waves are nearly totally reflected. On the downwave side, the flow field still remains disturbed, and the crossing mode 
$\cos ({\rm \pi} y/d)$ is predominantly exhibiting standing wave behaviours. Considering the wavenumber 
$ka=1.5707$, slightly less than 
${\rm \pi} /2$, it can be expressed as 
$kd=2ka = {\rm \pi}-\epsilon$, where 
$\epsilon \ll 1$. The characteristic wavenumber 
$\gamma _{1}$ is approximated as
\begin{equation} \gamma_1=\sqrt{{\rm \pi}^2/d^2 - ({\rm \pi}-\epsilon)^2/d^2} \approx \sqrt{2\epsilon {\rm \pi}/d^2}. \end{equation}
The smallness of the characteristic wavenumber 
$\gamma _{1}$ leads to a slow decay of the associated evanescent mode. Although this mode will eventually diminish at a significant distance from the metacylinder, it persists within a fairly large region surrounding the metacylinder.
Figure 12. Wave patterns scattered by a periodic array of circular metacylinders with a normalised radius of 
$a/d=0.5$ under normal wave incidence (
$\theta _{0}=0^{\circ }$) at 
$ka=1.5707$ close to crossing mode wavenumber 
$ka={\rm \pi} /2$, exhibiting nearly total reflection. Panels (a,b) exhibit the modulus and real part of the wave pattern, respectively.
In the case of oblique wave excitation, we focus on the wavenumber 
$ka=1.5025$, characterised by minimal energy reflection. Figure 13 showcases the wave patterns scattered by a periodic array of circular metacylinders at 
$ka=1.5025$, where the energy reflection is minimal, leading to nearly total transmission. Notably, the transmitted waves propagate at a different angle compared with the incident waves. Specifically, at 
$ka=1.5025$, the far-field transmitted waves are dominated by the two components 
$T_{-1}$ and 
$T_{0}$ based on (3.6a,b). Specifically, the computation indicates 
$|T_{-1}|>|T_{0}|$. As a consequence, the propagation of transmitted waves is primarily governed by the angle 
$\theta _{-1}=\arctan (\beta _{-1}/\alpha _{-1})\approx -19.78^{\circ }$. Therefore, if the component 
$T_{0}$ is smaller than other components, the transmitted waves will propagate at an angle different from the incident waves, resulting in wave-bending effects. This feature of metagratings was also discussed by Putley et al. (Reference Putley, Guenneau, Porter and Craster2022).
Figure 13. Wave pattern scattered by a periodic array of circular metacylinders with a normalised radius 
$a/d=0.5$ under the oblique wave excitation (
$\theta _{0}=45^{\circ }$) at a wavenumber 
$ka=1.5025$, showing nearly perfect wave transmission and wave bending effects on the downwave side.
6.2. Metasquare
We turn our attention to wave scattering by a periodic array of metasquares with each composed of 20 channels, where the plate width is 
$b/d=0.5$. Figure 14 depicts the variation of reflected energy 
$E_R$ with respect to the non-dimensional wavenumber 
$kb$ considering both head wave incidence (
$\theta _{0}=0^{\circ }$) and oblique wave incidence (
$\theta _{0}=45^{\circ }$) displayed in figures 14(a) and 14(b), respectively. Under the normal wave incidence as in figure 14(a), the reflected energy experiences strong oscillations near 
$kb={\rm \pi} /2$, rapidly alternating between total transmission and perfect reflection. The same oscillatory behaviours were also observed in the scattering of acoustic waves by a rectangular metamaterial cavity (Jan & Porter Reference Jan and Porter2018) due to complex interference. In the oblique wave excitation as in figure 14(b), the strong oscillations near 
$kb={\rm \pi} /2$ are also observed, and there exist dense discrete wavenumbers at which the nearly perfect wave transmission occurs. However, the value of reflected energy 
$E_R$ does not exceed 
$0.5$ within the considered wavenumber range, and thus perfect reflection is not achieved.
Figure 14. Reflected energy for a periodic array of metasquares with 
$b/d=0.5$ under head wave incidence 
$\theta _{0}=0^{\circ }$ (a) and oblique incidence 
$\theta _{0}=45^{\circ }$ (b). The vertical line corresponds to 
$kb={\rm \pi} /2$, where 
$b$ denotes the semiwidth of the plate constituting the metasquare.
To illustrate the total reflection 
$E_R\rightarrow 1$ under the normal wave incidence by a metasquare, we examine the wave patterns at 
$kb=1.5350$, where the wave transmission is minimised, as shown in figure 15. Unlike the scenario of perfect reflection by a periodic array of circular metacylinders in figure 12, where the wavenumber 
$ka=1.5707$ closely aligns with the crossing mode wavenumber 
$ka={\rm \pi} /2$, the current wavenumber deviates from the crossing mode wavenumber. As a consequence, the evanescent mode, associated with the characteristic wavenumber 
$\gamma _{1}$, decays rapidly with distance from the metasquare, resulting in a quiescent flow field on the downwave side of the structure.
Figure 15. Wave pattern scattered by a periodic array of metasquares with a semiwidth ratio of 
$b/d=0.5$, under head wave excitation (
$\theta _{0}=0^{\circ }$) at 
$kb=1.5350$, illustrating nearly total wave reflection. Panels (a,b) exhibit the modulus and real part of the wave pattern, respectively.
To showcase the perfect wave transmission predicted by the reflected energy plot, shown in figure 14(b), for the scattering of an array of metasquares by oblique waves (
$\theta _{0}=45^{\circ }$), wave patterns at a wavenumber 
$kb=1.3975$ are presented in figure 16. It is observed that the upwave flow field is minimally disturbed, indicating nearly perfect transmission of wave energy. Additionally, the wave field downstream aligns closely with the incident wave pattern, different from the scenario of oblique wave interactions with an array of circular metacylinders shown in figure 13, where wave propagation bends. In the current set-up, however, the transmitted wave associated with 
$T_0$ predominates over the component with 
$T_{-1}$, i.e. 
$T_0\gg T_{-1}$. Therefore, wave propagation remains unchanged, with only a phase shift occurring.
Figure 16. Wave pattern scattered by a periodic array of metasquares with a semiwidth ratio of 
$b/d=0.5$, under the action of oblique waves (
$\theta _{0}=45^{\circ }$) at 
$kb=1.3975$, illustrating nearly perfect wave transmission. Panels (a,b) exhibit the modulus and real part of the wave pattern, respectively.
6.3. Metawedge
For a periodic array of metawedges, we consider the set-up with an averaged semiwidth of 
$b_{m}/d=0.5$ and a ratio of longer base to shorter base 
$\ell =3.0$. Again the metawedge is composed of 20 channels. Figure 17 presents the reflected energy under the head wave incidence (
$\theta _{0}=0^{\circ }$) and oblique wave incidence (
$\theta _{0}=45^{\circ }$), displayed in figures 17(a) and 17(b), respectively. One notable feature in figure 17(a) is the nearly total reflection of waves across a wide spectrum of wavenumbers, exhibiting ‘rainbow reflections.’ Therefore, this device can act as a ‘broadband wave reflector.’ Under the quartering wave excitation as in figure 17(b), neither total wave reflection nor perfect wave transmission occurs within the considered range of wavenumbers.
Figure 17. Reflected energy for a periodic array of metawedges with the averaged semiwidth 
$b_{m}/d=0.5$ and base ratio 
$\ell =3.0$ under the actions of head waves 
$\theta _{0}=0^{\circ }$ (a) and oblique waves 
$\theta _{0}=45^{\circ }$ (b). The vertical line corresponds to 
$kb_{m}={\rm \pi} /2$.
To illustrate the near-perfect reflection achieved by the metawedge array, figure 18 presents the modulus, real part and imaginary part of the wave pattern corresponding to 
$kb_{m}=1.1980$ under head sea excitation. The set-up of the metawedge is identical to the one considered in figure 17. In this case, the wave energy experiences complete reflection resulting in a quiet flow field on the downwave side. On the upwave side, the real part is predominant whereas the imaginary part is negligible. As a consequence, the wave pattern on the upwave side manifests standing wave characteristics. Moreover, the wave crest lines are straight except for the flow region in the vicinity of the metawedge, then exhibiting two-dimensional behaviours.
Figure 18. Wave pattern scattered by a periodic array of metawedges, with an averaged semiwidth of 
$b_{m}/d=0.5$ and longer-to-shorter base ratio 
$\ell =3$, under the excitation of head waves (
$\theta _{0}=0^{\circ }$) at 
$kb_{m}=1.1980$ illustrating nearly perfect reflection. Panels (a,b) exhibit the modulus and real part of the wave pattern, respectively.
7. Results for surface-piercing plate-arrays
Finally, we investigate the scattering of waves by an array of two-dimensional partially submerged surface-piercing barriers.
7.1. Verification
For verification purposes, we show in figure 19 the modulus of the reflection coefficient, 
$|R|$, for an array of vertical barriers with uniform truncated depth 
$b$. The results presented in figures 19(a,c) and 19(b,d) correspond to a gap width of 
$c/b=0.5$ and 
$c/b=0.05$, and figures 19(a,b) and 19(c,d) exhibit the results for 
$N=1$ and 
$N=10$ cavities, respectively. Good agreement is made with the solutions obtained from the discrete model developed in Huang & Porter (Reference Huang and Porter2023).
Figure 19. Modulus of the reflection coefficient 
$|R|$ by an array of vertical identical barriers for gaps 
$c/b=0.5$ (a,c) and 
$c/b=0.05$ (b,d) at 
$b/h=0.2$, where 
$c$ denotes the distance between adjacent barriers and 
$b$ is the truncated depth. Panels (a,b) and (c,d) are for 
$N=1$ and 
$N=10$ cavities, respectively. Comparison is made with the discrete model by Huang & Porter (Reference Huang and Porter2023).
In the case of a single cavity, depicted in figure 19(a,b) , the reflection coefficient experiences a transition from total transmission 
$|R|=0$ to perfect reflection 
$|R|=1$. This transition becomes sharp as the cavity gap 
$c/b$ decreases, and it occurs in the vicinity of the resonance frequency 
$\omega \approx \sqrt {g/b}$ corresponding to 
$Kb \approx 1$ (Newman Reference Newman1974). For multiple cavities as shown in figure 19(c,d), the solution exhibits increasingly rapid oscillations as the frequency approaches the resonant frequency for a single cavity and practically no transmission for frequencies beyond. As discussed in Huang & Porter (Reference Huang and Porter2023), oscillations arise from constructive/destructive interference effects from the ends of the array compounded with a retardation of the effective wave speed through the array (exemplified in the subsequent subsection) as resonance is approached.
7.2. Uniform and graded plate-arrays
We continue by making a comparison between uniform arrays of Huang & Porter (Reference Huang and Porter2023) and the graded arrays considered in Wilks et al. (Reference Wilks, Montiel and Wakes2022) and Wilks et al. (Reference Wilks, Montiel and Wakes2023). Figure 20 presents reflection coefficient 
$|R|$ for both uniform and graded surface-piercing plate-arrays under normal wave incidence (
$\theta _{0}=0^{\circ }$). The metastructure is composed of 
$N=20$ cavities, spanning the interval 
$x/h \in [-0.5, +0.5]$, with an average plate immersion of 
$b_m/h = 0.5$. For the graded plate-array, we adopted a constant aspect ratio strategy, with a base length ratio of 
$b_N/b_0 = 3.0$.
Figure 20. The modulus of the reflection coefficient 
$|R|$ by an array of uniform and graded vertical barriers for 
$\theta _{0}=0^{\circ }$ with panel (b) highlighting the area where the reflection curve for the uniform array touch the zero. The vertical grey line at 
$Kb_m=2/3$ corresponds to the lowest resonant wavenumber for the graded plate-array over which perfect reflection occurs.
As already described, 
$|R|$ for the uniform plate-array exhibits rapid oscillations between 
$|R| = 0$ and peaks approaching 
$|R|=1$ at resonance. The region of strong oscillations is magnified in figure 20(b). In contrast, the reflection curve for the graded plate-array is smooth, free of oscillatory behaviours, transitioning to 
$|R|=1$ at 
$Kb_N = 1$, corresponding to 
$Kb_{m}=2/3$ plotted by the grey vertical line in the figure.
Figure 21 exhibits the imaginary part of spatial potential distribution 
${\mathrm {Im}}[\phi (x,z)]$ within the flow field for wave scattering by a surface-piercing plate-array. Figures 21(a) and 21(b) illustrate the potential distribution for a uniform plate-array at 
$Kb_m=0.9777$ and 
$Kb_m=0.9784$, respectively. Despite slight variation in wavenumber, the reflection coefficient undergoes a sharp transition from 
$|R|=0$ to 
$|R|=1$ corresponding to complete transmission and perfect reflection, respectively, indicating a dramatic shift in the flow field dynamics. Figure 21(a) shows a multiple interference effect from the ends of the array with a large fluid response within the cavities and figure 21(b) shows an exponential decay through the array. In contrast, figure 21(c) exhibits the scenario where the plate-array is graded, where perfect reflection is observed for all 
$Kb_m \gtrsim 0.66$. In this configuration, a wave is trapped within the middle cavity where the group velocity has slowed to zero and hardly any fluid motion is observed downwave of this. A careful analysis of the rainbow reflection characteristics of graded arrays including a discussion relating the evolution of the wave field within the array to local Bloch wavenumbers for the corresponding infinite periodic array is given in the work of Wilks et al. (Reference Wilks, Montiel and Wakes2023).
Figure 21. Distribution of the imaginary part of the velocity potential in the flow field for wave scattering by a surface-piercing plate-array under normal incidence 
$\theta _{0}=0^{\circ }$: (a) uniform plate-array at 
$Kb_m=0.977698$; (b) uniform plate-array at 
$Kb_m=0.978375$; (c) graded plate-array at 
$Kb_m=0.977698$.
The theory developed in the paper allows for oblique wave incidence, but we found that the results did not change too much in character after replacing 
$k$ by 
$k \cos \theta _0$, being the 
$x$-component of the wavenumber.
7.3. Semicircular plate-array
Finally, we consider the wave scattering by a semicircular profiled plate-array. Figure 22 depicts the reflection curve as a function of non-dimensional wavenumber 
$ka$, where 
$a$ denotes the radius of the semicircle. This is also a graded array with the onset of resonance associated with the longest channel, therefore at 
$Ka = 1$. We observe a similar type of behaviour in 
$|R|$ and the plot for the potential field as for graded arrays. That is, we transition to 
$|R|=1$ for 
$Ka > 1$ preceded by a small number of oscillations in the reflection before 
$Ka =1$; and the fluid motion dies downwave of the cavity at which resonance occurs.
Figure 22. The modulus of the reflection coefficient 
$|R|$ by an array of vertical barriers subject to semicircular profile for 
$\theta _{0}=0^{\circ }$.
Similar to figure 21, figure 23 presents the imaginary components of the potential distribution, 
${\mathrm {Im}}[\phi (x,z)]$, within the flow field for wave scattering by a semicircular profiled plate array. Figures 23(a) and 23(b) illustrate the cases of total transmission and perfect reflection at 
$Ka=0.958022$ and 
$Ka=1.092743$, respectively, corresponding to 
$|R|=0$ and 
$|R|=1$ as in figure 22. Due to the graded nature of the semicircular metastructure, the physical properties are analogous to those of the wedge-shaped plate-array.
Figure 23. Distribution of the imaginary part of the velocity potential in the flow field for wave scattering by a semicircular profiled surface-piercing plate-array under normal incidence 
$\theta _{0}=0^{\circ }$ at 
$Ka=0.958022$ (a) and 
$Ka=1.092743$ (b).
8. Conclusions
In this paper, we have considered a variety of settings in which water waves interact with metastructures consisting of dense plate arrays. These settings include the scattering of plane waves by isolated vertical metacylinders extending uniformly through the depth in an open ocean, scattering of plane waves by periodic arrays of vertical metacylinders and oblique wave scattering by horizontal surface-piercing metacylinders. The metacylinders are formed by closely spaced parallel arrays of thin barriers whose variable length defines the shape of the structure. We have concentrated on square, rectangular, wedge and circular structures in this paper. In each setting, local fluid resonance in the cavities between the plates produces a global effect on the wave field which produces an unorthodox behaviour.
The key novelty of the work is that we have used an exact description of the plate array rather than replacing it with an effective medium. This has allowed us to consider wave frequencies above resonance where the effective medium theory breaks down and where the most interesting results are found. The method of solution that has been used is also novel and has been crucial in simplifying the otherwise complicated interaction between the multiple plate elements of the metastructures. We have shown how to apply a transform-based approach in each of the three settings to reduce the problem to a canonical type meaning that all three problems, though superficially quite different, are resolved as solutions to almost identical systems of equations.
A range of results have been produced across the three settings which have been shown to compare favourably with existing results (where that is possible) but showing new results, especially highlighting the role that resonance plays. Arguably, the most interesting results involve graded arrays in which the length of the plates in the array increases with distance into the structure (forming a wedge-shaped metacylinder). This produces a dense spectrum of resonance frequencies associated with the variable length of the cavities in the array and allows for broadbanded ‘rainbow reflection’ effects. We imagine these results will be of interest to coastal engineers developing defence schemes or devices with the potential to manufacture bespoke wave control or harness wave energy. The problems in this paper are set in the context of water waves but the methodology developed herein can be applied to problems in the areas of acoustics, elasticity and electromagnetics.
Acknowledgements
Authors warmly thank Dr J. Huang for providing the raw data to verify the algorithm.
Funding
H.L. was supported by A*STAR Science and Engineering Research Council, Singapore, grant no. 172 19 00089 under the Marine & Offshore Strategic Research Programme (M&O SRP). S.Z. gratefully acknowledges the financial support from the Open Research Fund Program of the State Key Laboratory of Ocean Engineering (Shanghai Jiao Tong University) (grant no. 1916).
Declaration of interests
The authors report no conflict of interest.
Appendix. Far-field scattering waves
The potential in (2.19) indicates that the scattering potential 
$\phi _{sca} = \phi -\phi _{inc}$ is written as
\begin{equation} \phi_{sca}(x,y) \approx{-} \frac{1}{4 {\rm \pi}} \sum_{l=0}^N \sum_{p=0}^{2Q+1} a_p^{(l)} \mbox{sgn}(x-x_l) \int_{-\infty}^\infty \int_{{-}b_l}^{b_l} {\mathrm{e}}^{-\gamma|x - x_l|+{\rm i} \beta (y-y')} w_p(y/b_l) \, {\rm d} y' \, {\rm d} \beta. \end{equation}By using the integral form of the zeroth-order Hankel function (Twersky Reference Twersky1962)
\begin{equation} H_{0}(k \varrho)= \frac{1}{{\rm \pi} {\rm i}} \int_{-\infty}^{\infty} \frac{{\mathrm{e}}^{-\gamma|x-x'|+{\rm i}\beta (y-y')}}{\gamma} {\rm d} \beta, \quad \varrho = \sqrt{(x-x')^2+(y-y')^2}, \end{equation}
where 
$\gamma$ has been defined in (2.11), the scattering potential can be rewritten as
\begin{align} \phi_{sca}(x,y) & \approx{-}\frac{{\rm i}}{4} \sum_{l=0}^{N} \sum_{p=0}^{2Q+1} a_{p}^{(l)} \int_{{-}b_{l}}^{b_{l}} \left[\frac{\partial}{\partial x'} H_{0}(k\varrho)\right]_{x'=x_l} w_{p}(y'/b_{l}) \, {\rm d} y'\nonumber\\ & ={-}\frac{{\rm i}}{4} \sum_{l=0}^{N} \sum_{p=0}^{2Q+1} a_{p}^{(l)} \int_{{-}b_{l}}^{b_{l}} \left[ \frac{k(x-x')}{\varrho} H_{1}(k \varrho) \right]_{x' = x_l} w_{p}(y'/b_{l}) \, {\rm d} y'. \end{align}
In the limit that 
$k r = k \sqrt {x^2 + y^2} \to \infty$, 
$\varrho \to r$ and 
$x-x_l \to \varrho \cos \theta$, 
$\theta = \tan ^{-1}(y/x)$ and using the asymptotic representation of first-order Hankel function for large argument (Abramowitz & Stegun Reference Abramowitz and Stegun1964)
\begin{equation} H_{1}(kr) \sim \sqrt{\frac{2}{{\rm \pi} kr}} {\mathrm{e}}^{{\rm i}(kr-3{\rm \pi}/4)}, \end{equation}
the scattering potential in the far field 
$kr \to \infty$ is approximated as
\begin{equation} \phi_{sca}(x,y) \sim \sqrt{\frac{2}{{\rm \pi} kr}} A(\theta;\theta_{0}) {\mathrm{e}}^{{\rm i}(kr-{\rm \pi}/4)} \end{equation}
such that the scattering amplitude 
$A(\theta ;\theta _{0})$ is approximated numerically by
\begin{equation} A(\theta;\theta_{0}) \approx{-}\frac{k\cos\theta}{4} \sum_{l=0}^{N} b_k {\mathrm{e}}^{-{\rm i} kx_{l}\cos\theta} \sum_{p=0}^{2Q+1} a_{p}^{(l)} D_{p}(kb_{l}\sin\theta). \end{equation}
