Rainbow reflection and broadband energy absorption of water waves by graded arrays of vertical barriers

Abstract

The rainbow reflection effect describes the broadband spatial separation of wave spectral components caused by a spatially graded array of resonators. Although mainly studied in optics and acoustics, this phenomenon has recently been demonstrated both theoretically and experimentally for water waves travelling through an array of vertical cylinders. The scattering of two-dimensional linear water waves by an array of vertical, surface-piercing barriers is considered here, in which both the submergence and spacing between the barriers are spatially graded. The rainbow reflection effect arises naturally as wave energy temporarily becomes amplified at different locations depending on frequency. Band diagram calculations are used to demonstrate that this is a consequence of the wave gradually slowing down throughout the array. The wave/barriers scattering problem is then augmented by positioning heave-restricted, rectangular floating bodies equipped with a linear damping mechanism between each adjacent pair of barriers. A solution to the resulting boundary value problem is obtained using an integral equation/Galerkin method. Using constrained optimisation, passive rainbow absorbers are designed that achieve near-perfect absorption (i) over a discrete set of frequencies and (ii) over an octave. This suggests potential applications of rainbow absorbers in the design of smart coastal technologies.

1. Introduction

There is growing interest throughout the wave sciences in metamaterials that have a subwavelength structure which is spatially varying in some way. This spatial variation can give rise to so-called rainbow phenomena, through which broadband wave signals are slowed and spatially separated into their spectral components. These phenomena can be classified into rainbow trapping, where the separated energy becomes permanently localised, and rainbow reflection, where the energy is temporarily amplified before being reflected (He, He & Jin 2012; Chaplain et al. 2020). First demonstrated by Tsakmakidis, Boardman & Hess (2007) for optical waves, these concepts have since inspired the development of analogues in several fields, including elastic waves (Skelton et al. 2018; Arreola-Lucas et al. 2019), seismic waves (Colombi et al. 2016), water waves (Bennetts, Peter & Craster 2018) and acoustics. In this latter discipline, rainbow reflection has been prominent in the design of acoustic filters and sensors, where one-dimensional waveguides with graded resonant side cavities have been shown to spatially separate frequency components (Zhu et al. 2013; Ni et al. 2014; Zhao & Zhou 2019). The local energy amplification of acoustic waves has also been studied in two-dimensional arrays of (i) rigid cylinders with chirped spacing (Romero-García et al. 2013; Cebrecos et al. 2014) and (ii) C-shaped resonators with graded radii (Bennetts, Peter & Craster 2019). Jiménez et al. (2017) further extended the study of graded acoustic metamaterials by accounting for viscothermal losses. Theoretical and numerical methods were used to study energy absorption by arrays of planar Helmholtz resonators, which were optimised numerically to create devices which achieve (i) perfect absorption over a discrete set of frequencies and (ii) near-perfect absorption over a frequency interval. The performance of the device was later validated experimentally using a three-dimensionally printed prototype. Our primary goal is to investigate the feasibility of achieving this rainbow absorption for water waves in a model wave energy converter with an appropriate power take-off mechanism.

Recently, several other remarkable capabilities of metamaterials, originally studied for optical or acoustic waves, have been predicted in the context of water waves, where research is motivated by applications that include the development of coastal protection technology or wave energy parks. For instance, two-dimensional periodic arrays of uniform rigid cylinders immersed in water exhibit negative refraction, ultra-refraction and bandgaps – the latter characterising frequency intervals in which waves cannot propagate without a change in amplitude (McIver 2000; Hu et al. 2004; Farhat et al. 2010). The lowest-frequency bandgap of such arrays is shifted to lower frequencies when the cylinders are replaced with C-shaped cylindrical resonators, which suggests that they could be used as effective breakwaters (Hu et al. 2011; Dupont et al. 2017). Bennetts et al. (2018) demonstrated rainbow reflection in arrays of C-shaped cylindrical resonators with graded radii. In these arrays, the group velocity of the waves gradually approaches zero through the array at a location that depends on frequency, which results in local energy amplification. Arrays of rigid cylinders with chirped spacing also give rise to the rainbow reflection effect and associated local energy amplifications, which was verified experimentally by Archer et al. (2020). Here, the absence of a resonant substructure means that this effect is limited to higher-frequency intervals when compared with arrays of C-shaped cylinders. Although local energy amplification in graded arrays is now fairly well understood, energy absorption within such systems remains mostly unexplored.

In the last few decades, wave energy converters have been the subject of an intensive research effort (Falcão 2010). Several devices have been proposed which are capable of absorbing 100 % of wave energy at a particular frequency, but their efficiency is usually sensitive to a small change in frequency (Salter 1974; Evans et al. 1979; Evans & Porter 1995). Because real ocean waves occur over a broad range of frequencies, the limited band of effectiveness of single absorbers is typically overcome using optimal control theory (see Coe et al. (2017) and references therein). Recent research has also focused on the design of arrays of multiple absorbers, or wave energy parks, which optimally extract energy from realistic sea states (Babarit 2013; Göteman et al. 2020). Using shallow water theory, Porter (2021) recently demonstrated perfect absorption at all frequencies by a semi-infinite array of damped buoys. The buoys were modelled using a homogenised boundary condition at the upper surface of the two-dimensional fluid domain. To achieve perfect absorption, the first segment of the absorbing surface must have a negative damping coefficient and thus add energy to the fluid. Using full linear water wave theory, a similar device of finite length was shown to achieve near-perfect absorption for sufficiently large frequencies.

In § 2, we explore rainbow reflection in a one-dimensional waveguide consisting of a graded array of parallel, surface-piercing rigid vertical barriers. The limited complexity of this problem compared with two-dimensional arrays allows us to reveal a clear relationship between the structural geometry and the local energy amplifications. We build upon the results of Ursell (1947) and Porter & Evans (1995), who studied the scattering of water waves by a single vertical barrier, and those of Evans & Morris (1972), Newman (1974) and McIver (1985), who studied resonators consisting of two vertical barriers. In the present study, we demonstrate how these resonators can be combined into an array, and how local energy amplification within this array relates back to the local geometry, as well as to the band structure of the corresponding periodic infinite array.

The local energy amplification in graded arrays of vertical barriers suggests that highly efficient broadband absorption could be achieved by positioning suitable energy-absorbing mechanisms within the array. In § 3 we modify the problem so that the surface between each successive pair of barriers is occupied by a damped rectangular piston that is restricted to move in heave. Our solution to the corresponding linear boundary value problem is obtained by reformulating it as a coupled system of integral equations, which we solve using a Galerkin method. In § 4, we briefly discuss absorption by a single piston oscillating between two barriers. Then, we design rainbow absorbers (RAs) which achieve near-perfect absorption (i) over a discrete set of frequencies and (ii) over an octave – the latter highlighting the importance of non-absorbing resonant structures in transferring and localising energy. Our results extend those of Porter (2021) by using devices with a resonant substructure to investigate broadband absorption in the context of a more realistic model. They also extend those concerning the rainbow absorption of acoustic waves obtained by Jiménez et al. (2017) to water waves, which demonstrates the effect of employing an energy loss mechanism relevant to wave energy conversion.

2. Graded arrays of vertical barriers

2.1. Preliminaries

We consider the scattering of linear water waves by $N+1$ parallel surface-piercing vertical barriers in a fluid domain of infinite horizontal extent and constant, finite depth $H$. A two-dimensional Cartesian coordinate system $(x,z)$ is used, where the $z$ axis points vertically upwards and the horizontal coordinate $x$ corresponds to the direction of propagation of the waves. The fluid domain is bounded vertically by a flat seabed at $z=-H$ and a free surface located at $z=0$ when at equilibrium.

Each vertical barrier is of infinitesimal thickness and is situated at $x=x^{(n)}$, where $x^{(n)}\in \mathbb {R}$ for all $n\in \{0,\ldots,N\}$, and $x^{(n)}>x^{(n-1)}$ for $n\geq 1$. Each barrier is defined by the set of points $\varGamma _b^{(n)}=\{(x^{(n)},z)|z\geq -d^{(n)}\}$, where $d^{(n)}\in (-H,0)$ is the submergence depth of the $n$th barrier. A schematic of this geometry is given in figure 1.

Figure 1. Schematic of the multiple vertical barrier problem.

We denote by $\varOmega$ the region occupied by the fluid at equilibrium, that is

(2.1)\begin{equation} \varOmega=\left\{(x,z)\in\mathbb{R}\times({-}H,0)\left|(x,z)\notin\bigcup_{n=0}^{N}\varGamma_b^{(n)} \right.\right\}. \end{equation}

We also define the following subregions:

(2.2)\begin{equation} \left.\begin{array}{c} \varOmega^{(0)} = \{(x,z)|x< x^{(0)}\text{ and }z\in({-}H,0)\},\\[3pt] \varOmega^{(n)} = \{(x,z)|x\in({-}x^{(n-1)},x^{(n)})\text{ and } z\in({-}H,0)\}\quad \text{for all }n\in\{1,\ldots,N\},\\[3pt] \varOmega^{(N+1)} = \{(x,z)|x>x^{(N)}\text{ and }z\in({-}H,0)\}. \end{array}\right\} \end{equation}

In order to apply linear water wave theory to this problem, we assume that the steepness of the free-surface elevation of the incident waves is small. Thus, we consider the fluid to be inviscid, incompressible and undergoing irrotational, time-harmonic motion with angular frequency $\omega$. We assume that there is no normal flow across the seabed and the barriers. The velocity field of the fluid can therefore be expressed as the spatial gradient of some potential $\varPhi :\varOmega \times \mathbb {R}\to \mathbb {R}$, which is of the form

(2.3)\begin{equation} \varPhi(x,z,t) = {{\rm Re}}(\phi(x,z)\exp(-\mathrm{i}\omega t)). \end{equation}

The complex-valued spatial function $\phi$ then satisfies the following boundary value problem:

(2.4a)$$\begin{gather} \Delta\phi=0 \quad \forall\, (x,z)\in \varOmega, \end{gather}$$

(2.4b)$$\begin{gather}\partial_z\phi=0 \quad \text{for }z={-}H, \end{gather}$$

(2.4c)$$\begin{gather}\partial_z\phi=\frac{\omega^{2}}{g}\phi \quad \text{for } z=0, \end{gather}$$

(2.4d)$$\begin{gather}\partial_x\phi= 0 \quad \text{for all }(x,z)\in\bigcup_{n=0}^{N} \varGamma_b^{(n)}, \end{gather}$$

where (2.4c) is the linear free-surface condition with $g$ denoting acceleration due to gravity.

Since $\phi$ is twice differentiable, it must be continuous and have a continuous first partial derivative with respect to $x$ across the gaps beneath each barrier. This can be expressed in terms of one-sided limits as

(2.5a,b)\begin{equation} \lim_{x\to x^{(n)-}}\phi(x,z)=\lim_{x\to x^{(n)+}}\phi(x,z) \quad\mbox{and}\quad \lim_{x\to x^{(n)-}}\partial_x\phi(x,z)=\lim_{x \to x^{(n)+}}\partial_x\phi(x,z), \end{equation}

for all $n=0,\ldots,N$, where $z\in (-H,-d^{(n)})$. We require that $|\phi |$ remains finite as ${x\to \pm \infty }$. The solution must also obey an appropriate Sommerfeld radiation condition, which enforces that the scattered plane waves must propagate away from the array of barriers (see (2.12)).

The submerged edges of the barriers located at $(x^{(n)},-d)$ create singularities of order $-1/2$ in the fluid velocity (Evans & Morris 1972; Mosig 2018). Thus we require that for all $n\in \{0,\ldots,N\}$

(2.6)\begin{equation} ((x-x^{(n)})^{2}+(z+d^{(n)})^{2})^{1/2}\nabla \phi\to \boldsymbol{0}\mbox{ as }((x-x^{(n)})^{2}+(z+d^{(n)})^{2})^{1/2}\to 0. \end{equation}

The general solution to (2.4a)–(2.4c) is obtained using separation of variables:

(2.7)\begin{equation} \phi(x,z) = \left\{\begin{array}{ll} \begin{aligned} & \displaystyle \sum_{m=0}^{\infty} (A_m^{(n)}\exp(\mathrm{i}k_m(x-x^{(n)}))\\ & \quad +B_m^{(n)}\exp(-\mathrm{i}k_m(x-x^{(n)})))\psi_m^{(0)}(z) \end{aligned} & \begin{aligned} & \mbox{for } (x,z)\in\varOmega^{(n)}\\ & \quad \mbox{where }n\in\{0,\ldots,N\},\end{aligned}\\ \begin{aligned} & \displaystyle \sum_{m=0}^{\infty} (A_m^{(N+1)}\exp(\mathrm{i}k_m(x-x^{(N)}))\\ & \quad +B_m^{(N+1)}\exp(-\mathrm{i}k_m(x-x^{(N)})))\psi_m^{(0)}(z) \end{aligned} & \mbox{for } (x,z)\in\varOmega^{(N+1)}. \end{array}\right. \end{equation}

The quantities $k_m$ are the solutions to the free-surface dispersion relation

(2.8)\begin{equation} k\tanh kH=\frac{\omega^{2}}{g}.\end{equation}

We use the convention that ${k_0=2{\rm \pi} /\lambda \in \mathbb {R}_+}$, where $\lambda$ is the wavelength of the propagating waves, and that ${-\mathrm {i}k_n\in ((n-1){\rm \pi} /H,n{\rm \pi} /H)}$ for all ${n\in \mathbb {N}}$. The vertical eigenfunctions are defined as

(2.9)\begin{equation} \psi_m^{(0)}(z)=\beta_m^{{-}1/2}\cosh(k_m(z+H)), \end{equation}

where the superscript $(0)$ has been included for consistency with § 3. The normalisation constants $\beta _m=\sinh (2k_m H)/(4k_m H)+\frac {1}{2}$ are chosen so that the vertical eigenfunctions satisfy

(2.10)\begin{equation} \int_{{-}H}^{0}\psi_m^{(0)}(\xi)\psi_p^{(0)}(\xi)\, {{\rm d}} \xi = H\delta_{m p}, \end{equation}

for all $m,p\in \{0\}\cup \mathbb {N}$, where $\delta _{mp}$ is the Kronecker delta.

The fluid motion is forced by incident plane waves with angular frequency $\omega$. This forcing determines the incident potential on the exterior regions of the fluid, which is of the form

(2.11)\begin{equation} \phi_{{in}}(x,z)= \begin{cases} A_0^{(0)}\exp(\mathrm{i}k_0(x-x^{(0)}))\psi_0^{(0)}(z) & \mbox{for }(x,z)\in\varOmega^{(0)},\\[3pt] B_0^{(N+1)} \exp(-\mathrm{i}k_0(x-x^{(N)}))\psi_0^{(0)}(z) & \mbox{for }(x,z)\in\varOmega^{(N+1)}. \end{cases} \end{equation}

We will mostly consider problems where the forcing originates from the left. In this case, $B_0^{(N+1)}=0$ and the transmission and reflection coefficients are $T=A_0^{(N+1)}/A_0^{(0)}$ and $R=B_0^{(0)}/A_0^{(0)}$, respectively. Since $|\phi |$ remains finite as $x\to \pm \infty$, we must have $A_m^{(0)}=B_m^{(N+1)}=0$ for all $m\geq 1$. The Sommerfeld radiation condition is thus

(2.12)\begin{equation} \phi(x,z)-\phi_{{in}}(x,z)\sim \begin{cases} B_0^{(0)}\exp(-\mathrm{i}k_0(x-x^{(0)}))\psi_0^{(0)}(z) & \mbox{as }x\to-\infty,\\[3pt] A_0^{(N+1)}\exp(\mathrm{i}k_0(x-x^{(N)}))\psi_0^{(0)}(z) & \mbox{as }x\to\infty, \end{cases} \end{equation}

for all $z\in (-H,0)$.

2.2. Overview of the solution technique

We invoke the continuity conditions (2.5a,b) and the no-flow condition (2.4d) for each barrier in order to map the amplitudes of waves incident on the $n$th barrier to the amplitudes of the waves scattered by that barrier. This mapping is approximated by the scattering matrix relation

(2.13)\begin{equation} {\boldsymbol{\mathsf{S}}}^{(n,n+1)}\begin{bmatrix} \boldsymbol{A}^{(n)}\\\boldsymbol{B}^{(n+1)}\end{bmatrix}= \begin{bmatrix} \boldsymbol{A}^{(n+1)}\\\boldsymbol{B}^{(n)}\end{bmatrix} \end{equation}

for $n\in \{0,\ldots,N\}$. Here, $\boldsymbol {A}^{(n)}$ and $\boldsymbol {B}^{(n)}$ are vectors with entries $A_m^{(n)}$ and $B_m^{(n)}$, for $0\leq m \leq M$, where $M$ is a truncation parameter.

The entries of each scattering matrix $\boldsymbol {S}^{(n,n+1)}$ are determined as a function of the barrier submergence $d^{(n)}$ as follows: (i) express the scattered wave amplitudes in terms of the horizontal fluid velocity $u(z)=\partial _x\phi$ on the interval $(-H,-d^{(n)})$, (ii) invoke the continuity conditions to construct an integral equation satisfied by $u$ and (iii) solve the integral equation via a Galerkin method which involves expanding $u$ over a basis of weighted Chebyshev polynomials, which encapsulate the singularities at $(x^{(n)},-d^{(n)})$ (see Porter & Evans (1995) for more details).

The scattering matrices $\{\boldsymbol {S}^{(0,1)},\ldots,\boldsymbol {S}^{(N,N+1)}\}$ and the horizontal positions of the barriers $\{x^{(0)},\ldots,x^{(N)}\}$ are used to assemble the global scattering matrix for the entire array of barriers, which satisfies

(2.14)\begin{equation} {\boldsymbol{\mathsf{S}}}^{(0,N+1)}\begin{bmatrix} \boldsymbol{A}^{(0)}\\ \boldsymbol{B}^{(N+1)}\end{bmatrix}=\begin{bmatrix} \boldsymbol{A}^{(N+1)}\\\boldsymbol{B}^{(0)}\end{bmatrix}. \end{equation}

The global scattering matrix is used to compute the amplitudes of the scattered waves in $\varOmega ^{(0)}$ and $\varOmega ^{(N+1)}$. Finally, the wave amplitudes in the interior regions $\varOmega ^{(n)}$ are computed recursively for $n\in \{1,\ldots,N\}$. The algorithm used to construct the global scattering matrix and compute the unknown coefficients is based on the method developed by Ko & Sambles (1988) to study electromagnetic wave propagation in layered media. This method has been adapted for problems involving the multiple scattering of linear water waves, yielding numerically stable results (Bennetts & Squire 2009; Montiel, Squire & Bennetts 2015). To the best of our knowledge, ours is the first application of this method to the multiple vertical barriers problem.

2.3. The fundamental resonance of a pair of barriers

We first consider the response of an isolated resonator consisting of a pair of barriers, i.e. when $N=1$. This problem has been extensively studied when $d^{(0)}=d^{(1)}$. Such a pair of identical barriers is associated with an infinite collection of resonant frequencies (Evans & Morris 1972). Here, we are exclusively concerned with the fundamental resonance (i.e. the lowest-frequency resonance), which is analogous to the Helmholtz resonance in acoustic resonators. The corresponding fundamental frequency, which we denote $\omega _{{fund}}$, can be approximated by $\sqrt {g/d^{(0)}}$. This approximation is valid when the barriers are close to each other and submerged in deep water, i.e. $k_0 l^{(1)}\ll 1$ and $d^{(0)}/H\ll 1$. Under these assumptions, the mass of fluid between the barriers can be treated as a heaving rectangular body oscillating about its hydrostatic equilibrium (Newman 1974). As $l^{(1)}$ increases, this approximation loses validity and the fundamental frequency decreases (McIver 1985). When $d^{(0)}\neq d^{(1)}$, the fundamental frequency of the corresponding resonator is similar to that of a pair of identical barriers with submergence $(d^{(0)}+d^{(1)})/2$, provided $|d^{(1)}/d^{(0)}-1|\ll 1$. This result is shown in figure 2(a). Throughout this paper, we fix the depth of the fluid to be $H=20$ m and the amplitude of the incident wave to be $0.1$ m. Linear water wave theory has considerable experimental validation at the resulting ratios of amplitude to wavelength considered in this paper (Mei, Stiassnie & Yue 2005; Montiel et al. 2013a,b).

Figure 2. (a) A comparison of the fundamental frequency of two unequal barriers (solid curve) with the fundamental frequency of two equal barriers with submergence $(d^{(0)}+d^{(1)})/2$ (dotted curve). Here, $l^{(1)}=1$ is fixed. (b) A plot of $1/Q$ as a function of $l^{(1)}/d^{(0)}$, which suggests that $Q$ and $l^{(1)}/d^{(0)}$ are inversely proportional. In both panels we use $d^{(0)}=10$ m, and the depth of the fluid is fixed to be $H=20$ m here and throughout this paper.

To further characterise the fundamental resonance, we introduce the following analogue for the average power of the oscillations between two barriers:

(2.15)\begin{equation} \tilde{P}(\omega)=\omega^{2}\int_{x^{(0)}}^{x^{(1)}}|\phi(x,0)|\, {\rm d} x, \end{equation}

where $\phi$ is the velocity potential of the fluid when excited by a plane wave incident from the left. The bandwidth of the fundamental resonance, which we denote ${\rm \Delta} \omega _{{fund}}$, is described by the full width at half the maximum of the corresponding peak of $\tilde {P}$. Specifically, ${\rm \Delta} \omega _{{fund}}=\delta _-+\delta _+$, where $\delta _\pm >0$ are the minimum values satisfying $\tilde {P}(\omega _{{fund}}\pm \delta _\pm )=\tilde {P}(\omega _{{fund}})/2$. Finally, the ‘sharpness’ of the peak is quantified by the dimensionless $Q$-factor, which is defined as $Q =\omega _{{fund}} /{\rm \Delta} \omega _{{fund}}$. As shown in figure 2(b), $Q$ and $l^{(1)}/d^{(0)}$ are approximately inversely proportional; namely, there is a very strong agreement between the data used to generate the plot and the model $1/Q=a(l^{(1)}/d^{(0)})$ for some coefficient $a\in \mathbb {R}$. This extends the assertion of McIver (1985) that the fundamental resonance ‘detunes’ as the spacing between the barriers increases, which was quantified by the distance between the zeros of the reflection and transmission coefficients associated with the resonance. To the best of our knowledge, the apparent inverse relationship between $Q$ and $l^{(1)}/d^{(0)}$ has not been reported before.

2.4. Local energy amplifications in graded arrays

The observations from the previous section motivate the design of a rainbow-reflecting device consisting of multiple vertical barriers, where the submergence and spacing of the barriers increase throughout the array. We aim to adjust the parameters $d^{(n)}$ and $l^{(n)}$ so that the fundamental resonances of each cavity are evenly spaced over the octave $\omega \in [0.8,1.6]$ s$^{-1}$. These frequencies are typical of wind-forced coastal waves. The resonant frequencies of the cavities must decrease in the direction of wave propagation, in order to prevent the deep barriers of the low-frequency cavities from sheltering the high-frequency cavities. This is because the transmission coefficient of the $n$th barrier in isolation is a decreasing function of $\omega ^{2}d^{(n)}$ (Ursell 1947). We design the array so as to excite local energy amplifications when the waves are incident from the left. Thus for $N=9$ the fundamental resonance of the $n$th cavity is given by $\omega ^{(n)}=(17-n)0.1$ s$^{-1}$.

First, we find the submergence $\tilde {d}^{(n)}$ for pairs of identical barriers so that the fundamental resonance is $\omega ^{(n)}$, while fixing the separation of the barriers to be $\tilde {l}^{(n)}=\tilde {d}^{(n)}/5$. This is done by maximising $\tilde {P}(\omega ^{(n)})$ as a function of $\tilde {d}^{(n)}\in (0,H)$, which yields the optimised submergences of the barriers presented in table 1. Second, we use these optima to choose parameters for the array of barriers. We want the geometry of the $n$th cavity to be similar to the pair of identical barriers which produces a resonance at $\omega ^{(n)}$. Thus, we select $l^{(n)}=\tilde {l}^{(n)}$. In addition, we require that the average submergence of the two barriers forming each cavity is equal to the depth of the optimised pair of barriers, that is

(2.16)\begin{equation} \frac{d^{(n-1)}+d^{(n)}}{2}=\tilde{d}^{(n)} \end{equation}

for all $n\in \{1,\ldots,N\}$. This is an underdetermined system for the $N+1$ unknown parameters $d^{(n)}$. We choose the solution which satisfies the constraints ${0< d^{(1)}<\cdots < d^{(N)}< H}$ and minimises the grading of the array (see Appendix A for further details).

Table 1. The submergences $\tilde {d}^{(n)}$ of pairs of identical barriers whose fundamental resonance occurs at $\omega ^{(n)}$, which were obtained by maximising the time-average power of free-surface oscillations. The optimised submergence values occur at approximately 0.9 times the estimate $g/(\omega ^{(n)})^{2}$.

Colour plots of $|\phi (x,z)|$ for the array parametrised by the submergences reported in table 1 are given in figure 3 at all the frequencies $\omega ^{(n)}$. As expected, the strongest local energy amplification in the graded array for a given frequency occurs in the cavity associated with that fundamental frequency. In results not shown here, this does not occur when waves are incident from the right, which is a consequence of the low transmission of high-frequency waves discussed earlier. Furthermore, the differing scales of the colourbars in figure 3 vary in a non-uniform way, which indicates that the resonant frequencies of the graded array do not perfectly coincide with the frequencies $\omega ^{(n)}$. These frequency shifts are introduced by the method used to design the graded array. Firstly, changing the geometry of a cavity from a pair of identical barriers to a pair of unequally submerged barriers increases its fundamental frequency, provided that the mean submergence of the barriers is kept constant. This effect was shown in figure 2(a). Secondly, the coupling of multiple resonators in the graded array alters the resonant frequencies in a non-transparent way.

Figure 3. Colour plots of the velocity potential $|\phi (x,z)|$ when forced by $0.1$ m waves which are incident from the left of the array of barriers. The frequencies of the incident waves shown here are the desired resonant frequencies, which were used to determine the parameters of the array. Specifically, these frequencies are (a) $\omega =1.6$ s$^{-1}$, (b) $1.5$ s$^{-1}$, (c) $1.4$ s$^{-1}$, (d) $1.3$ s$^{-1}$, (e) $1.2$ s$^{-1}$, ( f) $1.1$ s$^{-1}$, (g) $1.0$ s$^{-1}$, (h) $0.9$ s$^{-1}$ and (i) $0.8$ s$^{-1}$. The position of the strongest local energy amplification moves towards the right as the frequency decreases.

In order to analyse this phenomenon further, we consider the propagation of water waves through an infinite periodic array of uniform vertical barriers. Using the present notation, this is the case where $x^{(n)}=nl$ and $d^{(n)}=d$ for all $n\in \mathbb {Z}$, where $l >0$ and $d\in (0,H)$ are fixed. Our intention is to locally approximate the wave field within a finite graded array of vertical barriers using solutions to the infinitely periodic case, which is reasonable provided that the grading is weak (Bennetts et al. 2019). Infinite periodic arrays of uniform vertical barriers can be studied using Bloch's theorem, which motivates seeking solutions of the form

(2.17)\begin{equation} \phi(x+nl,z) = {\rm e}^{\mathrm{i}qnl}\phi(x,z). \end{equation}

In the above equation we have introduced the unknown Bloch wavenumber $q\in \mathbb {C}$. To compute the corresponding Bloch modes, (2.17) is applied to the horizontal boundary of the unit cell $\{(x,z)|x\in (-l/2,l/2)\mbox { and }z\in (-H,0)\}$, which yields

(2.18a,b)\begin{equation} \phi(l/2,z)= {\rm e}^{\mathrm{i}q l}\phi({-}l/2,z)\quad\mbox{and}\quad \frac{\partial \phi}{\partial x}(l/2,z)= {\rm e}^{\mathrm{i}q l} \frac{\partial \phi}{\partial x}({-}l/2,z). \end{equation}

These quasi-periodic boundary conditions are then solved in conjunction with (2.4a)–(2.4d), where suitable modifications are made so that $x\in (-l/2,l/2)$.

A numerical solution is obtained by solving a generalised eigenvalue problem, which incorporates the scattering matrix for a single vertical barrier (see e.g. Peter & Meylan 2009). For each frequency $\omega$, propagating Bloch modes (i.e. $q\in \mathbb {R}$) are sought in the irreducible Brillouin zone $q\in [0,{\rm \pi} /l]$. We obtain dispersion curves in the $q{-}\omega$ plane, which form so-called band diagrams. We remark that frequency intervals occupied by the dispersion curves are called passbands, and these frequencies are associated with a propagating mode which can travel through the infinite array without a change in amplitude. The slope of the dispersion curve $\textrm {d}\omega /\textrm {d}q$ determines the group velocity of this propagating mode. The frequency intervals in which there are no real-valued Bloch wavenumbers are called bandgaps. At these frequencies, waves cannot propagate without a change in amplitude.

Some examples of band diagrams with different barrier submergences and spacings are shown in figure 4. In each case, our numerical scheme identifies only one passband and one bandgap. We observe that the fundamental frequency of a pair of barriers with separation $l$ and submergence $d$ is a lower bound of the bandgap. The fundamental frequency decreases as the barriers become deeper, which coincides with the dispersion curves becoming more compressed. This mirrors a similar result found by Bennetts et al. (2018) for infinite two-dimensional arrays of C-shaped cylindrical resonators. In that setting, the resonant frequency of the isolated resonator is an upper bound of the passband, and the resonance is said to induce the bandgap (Chaplain et al. 2020). In both cases, the bandgap can be shifted by modifying the resonant substructure of the array.

Figure 4. Dispersion curves for infinite periodic arrays of uniform vertical barriers (solid lines) in which $d=(d^{(n-1)}+d^{(n)})/2$ and $l=l^{(n)}$ for (a) $n=2$, (b) $4$, (c) $6$ and (d) $8$. Here, the values of $d^{(n)}$ and $l^{(n)}$ are those of the graded array shown in figure 3. In each panel, the dotted line shows the dispersion relation (2.8), which is the dispersion curve for unobstructed waves. Notably, the bandgap begins above the fundamental frequency of a pair of barriers with separation $l$ and submergence $d$. These fundamental frequencies, which are represented by dashed lines, are (a) $\omega ^{(n)}=1.5$ s$^{-1}$, (b) $1.3$ s$^{-1}$, (c) $1.1$ s$^{-1}$ and (d) $0.9$ s$^{-1}$.

In the finite array of barriers, the potential surrounding each successive pair of barriers can be understood in terms of the corresponding Bloch problem. In particular, the dispersion curve of the infinite array of barriers with submergence $(d^{(n-1)}+d^{(n)})/2$ and spacing $l^{(n)}$ approximately determines the behaviour in the $n$th cavity for a weakly graded structure, and in this setting it is referred to as the local dispersion curve. Since the fundamental frequencies of the cavities decrease from left to right, the local dispersion curves become more compressed and the frequency interval which supports propagating modes becomes narrower. Thus low-frequency waves propagate further towards the right than high-frequency waves.

The dispersion curves can also be used to explain the local energy amplification. To demonstrate this, finite difference estimates of the group velocity of the propagating mode are tabulated in table 2 for all cavities and for a range of frequencies. Waves slow down as they propagate through the finite, graded array because the group velocity is a decreasing function of $n$, which attains a minimum before becoming undefined. This minimum typically occurs in the cavity whose fundamental frequency is the frequency of the wave. The wave can only attenuate beyond this point (which is called the turning point by Romero-García et al. (2013)), since propagation at this frequency is no longer supported by the array. A local energy amplification occurs near the turning point because the propagating mode has a relatively small group velocity, so its energy is localised for a relatively large duration. Moreover, the group velocity is a decreasing function of $\omega$ where it is defined, which means that higher-frequency waves slow down more rapidly than lower-frequency waves and amplify further towards the left, resulting in spatial separation. Thus, graded arrays of vertical barriers give rise to the rainbow-reflection effect.

Table 2. Finite difference estimates of the group velocity (m s$^{-1}$) of waves travelling through the infinite array of barriers with $d=d^{(n-1)}$ and $l=l^{(n)}$, with frequency $\omega$. Dashes (—) are used to denote the absence of a propagating Bloch mode. The group velocity is a decreasing function of both $\omega$ and $n$.

However, the local energy amplification at the turning point is not permanent, so the graded array does not give rise to the rainbow-trapping effect. Because the unit cell of the Bloch problem possesses reflective symmetry, the infinite array supports both left- and right-propagating modes at passband frequencies, which have equal and opposite group velocity. The strong coupling between these modes near the turning point causes the wave to change direction (He et al. 2012; Chaplain et al. 2020). This coupling is forbidden by rainbow-trapping devices, which can be achieved through the use of certain asymmetric substructures (Chaplain et al. 2020). Since transmission by the array of barriers is negligible, the left-travelling wave originating at the turning point reflects all of the incident energy. In the following section, we introduce a loss-inducing mechanism into the model, which will allow the localised energy to be absorbed within the structure rather than being reflected.

3. Solution of multiple absorbing cavity problem

3.1. Preliminaries

We adapt the problem considered in § 2 by introducing $N$ rigid rectangular floating bodies, or ‘pistons’, which completely occupy the horizontal gap between each adjacent pair of barriers. The pistons are restricted to move vertically in heave and are less dense than the fluid. Further, the density of each piston is $\rho _A^{(n)}/\delta ^{(n)}$, where $\rho _A^{(n)}$ is the area density of the piston with respect to the wetted surface and $\delta ^{(n)}$ is the thickness. The motion of each piston is governed by a damping term proportional to the piston's vertical velocity, with the damping coefficient denoted $\mu ^{(n)}$. A schematic of the geometry of this problem is given in figure 5.

Figure 5. Schematic of the multiple piston cavity problem. Viscous damping is symbolised by dashpots.

We assume that pistons are undergoing small, time-harmonic oscillations about their equilibria. This implies that the submergence of the bottom face of each piston is

(3.1)\begin{equation} z={{\rm Re}}(s^{(n)}\, {\rm e}^{-\mathrm{i}\omega t})-D^{(n)},\end{equation}

where $s^{(n)}$ and $D^{(n)}=\rho _A^{(n)}/\rho$ are the complex amplitude and the equilibrium submergence of the $n$th piston respectively. Here, $\rho$ is the density of the fluid. The relevant forces which govern the dynamics of the pistons are due to hydrostatic pressure, hydrodynamic pressure, gravity and viscous damping. After some algebra, the linearised equations of motion can be stated as

(3.2)\begin{equation} -\omega^{2}\rho_A^{(n)}l^{(n)}s^{(n)}-\mathrm{i}\omega\mu^{(n)} s^{(n)}+\rho g l^{(n)}s^{(n)}=\mathrm{i}\omega\rho \int_{x^{(n-1)}}^{x^{(n)}}\phi(x,-D^{(n)})\,{\rm d} x,\end{equation}

where $l^{(n)}=x^{(n)}-x^{(n-1)}$ (Wehausen 1971).

We redefine the subregions $\varOmega ^{(n)}$ to be the sets of points beneath the $n$th piston when at equilibrium, that is, $\varOmega ^{(n)} = \{(x,z)|x\in (x^{(n-1)},x^{(n)})\textrm { and }z\in (-H,-D^{(n)})\}$ for all $n\in \{1,\ldots,N\}$. The definitions of $\varOmega ^{(0)}$ and $\varOmega ^{(N+1)}$ are unchanged. The velocity potential of the fluid $\phi$ satisfies the boundary value problem given by (2.4a), (2.4b), (2.4d) and the additional conditions

(3.3a)$$\begin{gather} \partial_z\phi=\frac{\omega^{2}}{g}\phi \quad \text{for } z=0 \text{ and }x\in(-\infty,x^{(0)})\cup(x^{(N)},\infty), \end{gather}$$

(3.3b)$$\begin{gather}\partial_z\phi={-}\mathrm{i}\omega s^{(n)}\quad \mbox{where }x\in(x^{(n-1)},x^{(n)})\text{ and }z={-}D^{(n)}. \end{gather}$$

The potential is also subject to (2.5a,b), (2.6), (2.11) and (2.12).

3.2. Derivation of a system of integral equations

The expansion for the solution given in (2.7) still applies in subregions $\varOmega ^{(0)}$ and $\varOmega ^{(N+1)}$, since there is no piston at the surface of the fluid. In the remaining subregions $\varOmega ^{(n)}$, we express the general solution as the sum of a homogeneous solution $\phi _h^{(n)}$ and a particular solution $\phi _p^{(n)}$, for all $n\in \{1,\ldots,N\}$. The homogeneous solution satisfies (2.4a), (2.4b) and

(3.4)\begin{equation} \partial_z\phi_h^{(n)}(x,-D^{(n)})=0, \end{equation}

for $x\in (x^{(n-1)},x^{(n)})$. These three equations constitute a Sturm–Liouville problem in the $z$ direction. By using separation of variables and choosing an expansion centred at the midpoints of each piston $m^{(n)}=x^{(n-1)}+\frac {1}{2}l^{(n)}$ we obtain

(3.5)\begin{align} \phi_h^{(n)}(x,z)&=A_0^{(n)}+B_0^{(n)}(x-m^{(n)}) +\sum_{m=1}^{\infty} (A_m^{(n)}\cosh(\kappa_m^{(n)} (x-m^{(n)}))\nonumber\\ &\quad +B_m^{(n)}\sinh(\kappa_m^{(n)} (x-m^{(n)}))) \frac{\psi_m^{(n)}(z)}{\cosh\left(\frac{1}{2} \kappa_m^{(n)}l^{(n)}\right)}, \end{align}

where $\kappa _m^{(n)}=m{\rm \pi} /h^{(n)}$ in which $h^{(n)}=H-D^{(n)}$, and the vertical eigenfunctions have been defined as

(3.6)\begin{equation} \psi_m^{(n)}(z)=\sqrt{2}\cos(\kappa_m^{(n)}(z+H)), \end{equation}

for all $m\in \mathbb {N}$. These eigenfunctions satisfy the following orthogonality relation:

(3.7)\begin{equation} \int_{{-}H}^{{-}D^{(n)}}\psi_m^{(n)}(\xi)\psi_p^{(n)}(\xi)\, {{\rm d}} \xi = h^{(n)}\delta_{m p}, \end{equation}

for all $m,p\in \{0\}\cup \mathbb {N}$ and for all $n\in \{1,\ldots,N\}$.

The particular solution $\phi _p^{(n)}$ satisfies (2.4a), (2.4b) and (3.3b). We find a particular solution with the ansatz that $\phi _p^{(n)}$ is quadratic in both $x$ and $z$, which gives

(3.8)\begin{equation} \phi_h^{(n)}(x,z)={-}\frac{\mathrm{i}\omega s^{(n)}}{2h^{(n)}}((z+H)^{2}-(x-m^{(n)})^{2}),\end{equation}

where we have imposed symmetry about the horizontal midpoints $x=m^{(n)}$. This treatment of the inhomogeneous boundary condition is in complete analogy to the solution obtained by Yeung (1981) for the problem of an oscillating vertical circular cylinder in a fluid of finite depth.

The scattering matrix method presented in § 2.2 cannot be straightforwardly applied to the present problem, since the equations of motion of the pistons and the wave scattering problem must be solved simultaneously. Instead, we solve the system using a fully coupled system of equations. We begin by introducing the auxiliary functions

(3.9)\begin{equation} u^{(n)}(z) = \left\{\begin{array}{ll} \lim_{x\to x_j}\partial_x\phi(x,z) & \mbox{for }z<{-}d^{(n)}\\ 0 & \mbox{for }z\geq{-}d^{(n)}. \end{array}\right. \end{equation}

Applying (2.5b) as $x \to x^{(0)-}$ and $x \to x^{(N)+}$ to (2.7) yields

(3.10)$$\begin{gather} u^{(0)}(z)=\mathrm{i}k_0A_0^{(0)}\psi_0^{(0)}(z)-\sum_{m=0}^{\infty} \mathrm{i}k_mB_m^{(0)}\psi_m^{(0)}(z), \end{gather}$$

(3.11)$$\begin{gather}u^{(N)}(z)={-}\mathrm{i}k_0B_0^{(N+1)}\psi_0^{(0)}(z)+ \sum_{m=0}^{\infty} \mathrm{i}k_mA_m^{(N+1)}\psi_m^{(0)}(z). \end{gather}$$

After referring to (2.10), we obtain the following expressions for the scattered wave coefficients:

(3.12)$$\begin{gather} B_m^{(0)}=A_0^{(0)}\delta_{0m}+\frac{\mathrm{i}}{k_m H} \mathcal{P}^{(0)}(u^{(0)},\psi_m^{(0)}), \end{gather}$$

(3.13)$$\begin{gather}A_m^{(N+1)}=B_0^{(N+1)}\delta_{0m}-\frac{\mathrm{i}}{k_m H} \mathcal{P}^{(N)}(u^{(N)},\psi_m^{(0)}), \end{gather}$$

where $\mathcal {P}^{(n)}(\cdot,\cdot )$ is the real inner product over $(-H,-d^{(n)})$. Similarly, applying (2.5b) as $x \to x^{(n-1)+}$ and $x \to x^{(n)-}$ gives, after referring to (3.5) and (3.8) and invoking (3.7),

(3.14)$$\begin{gather} A_m^{(n)}=\frac{{{\rm coth}}\left(\frac{1}{2}\kappa_m^{(n)}l^{(n)}\right)}{2m{\rm \pi}} (\mathcal{P}^{(n)}(u^{(n)},\psi_m^{(n)})-\mathcal{P}^{(n-1)}(u^{(n-1)},\psi_m^{(n)})), \end{gather}$$

(3.15)$$\begin{gather}B_m^{(n)}=\frac{1}{2m{\rm \pi}}(\mathcal{P}^{(n)}(u^{(n)},\psi_m^{(n)}) +\mathcal{P}^{(n-1)}(u^{(n-1)},\psi_m^{(n)})). \end{gather}$$

Then, applying (2.5a) as $x\to x^{(n)}$ for all $n\in \{0,\ldots,N\}$ and using (3.12), (3.13), (3.14) and (3.15), we obtain the following system of integral equations:

(3.16)\begin{align} &\int_{{-}H}^{{-}d^{(0)}}(K^{(0)}(z,\xi)+K^{(1)}(z,\xi))u^{(0)}(\xi)\, {{\rm d}} \xi-\int_{{-}H}^{{-}d^{(1)}}L^{(1)}(z,\xi)u^{(1)}(\xi)\, {{\rm d}} \xi \nonumber\\ &\quad +\frac{\mathrm{i}\omega s^{(1)}}{2h^{(1)}}((z+H)^{2}-\frac{1}{4}(l^{(1)})^{2})-A_0^{(1)} +\frac{1}{2}B_0^{(1)}l^{(1)}={-}2A_0^{(0)}\psi_0^{(0)}(z), \end{align}

(3.17)\begin{align} &-\int_{{-}H}^{{-}d^{(n-1)}}L^{(n)}(z,\xi)u^{(n-1)}(\xi)\, {{\rm d}} \xi+\int_{{-}H}^{{-}d^{(n)}}(K^{(n)}(z,\xi)+K^{(n+1)}(z,\xi))u^{(n)}(\xi)\, {{\rm d}} \xi\nonumber\\ &\quad -\int_{{-}H}^{{-}d^{(n+1)}}L^{(n+1)}(z,\xi)u^{(n+1)}(\xi)\, {{\rm d}} \xi+\frac{\mathrm{i}\omega}{2}\left(\frac{s^{(n+1)}}{h^{(n+1)}} \left((z+H)^{2}-\frac{1}{4}(l^{(n+1)})^{2}\right)\right.\nonumber\\ &\quad \left.-\frac{s^{(n)}}{h^{(n)}}\left((z+H)^{2}-\frac{1}{4}(l^{(n)})^{2}\right) \right)+A_0^{(n)}+\frac{1}{2}B_0^{(n)}l^{(n)}-A_0^{(n+1)} +\frac{1}{2}B_0^{(n+1)}l^{(n+1)}=0, \end{align}

(3.18)\begin{align} &-\int_{{-}H}^{{-}d^{(N-1)}}L^{(N)}(z,\xi)u^{(N-1)}(\xi)\, {{\rm d}} \xi+\int_{{-}H}^{{-}d^{(N)}}(K^{(N)}(z,\xi)+K^{(0)}(z,\xi))u^{(N)}(\xi)\, {{\rm d}} \xi\nonumber\\ &\quad -\frac{\mathrm{i}\omega s^{(N)}}{2h^{(N)}}\left((z+H)^{2}-\frac{1}{4}(l^{(N)})^{2}\right)+A_0^{(N)} +\frac{1}{2}B_0^{(N)}l^{(N)}=2B_0^{(N+1)}\psi_0^{(0)}(z), \end{align}

where (3.17) holds for all $n\in \{1,\ldots,N-1\}$ and we have defined the integral kernels

(3.19)$$\begin{gather} K^{(0)}(z,\xi)=\sum_{m=0}^{\infty} \frac{\mathrm{i}}{k_mH} \psi_m^{(0)}(z)\psi_m^{(0)}(\xi), \end{gather}$$

(3.20)$$\begin{gather}K^{(n)}(z,\xi)=\sum_{m=1}^{\infty} \frac{{{\rm coth}}(\kappa_m^{(n)}l^{(n)})}{m {\rm \pi}}\psi_m^{(n)}(z)\psi_m^{(n)}(\xi), \end{gather}$$

(3.21)$$\begin{gather}L^{(n)}(z,\xi)=\sum_{m=1}^{\infty} \frac{\operatorname{csch}(\kappa_m^{(n)}l^{(n)})}{m {\rm \pi}}\psi_m^{(n)}(z)\psi_m^{(n)}(\xi). \end{gather}$$

3.3. Numerical solution using a Galerkin method

To find an approximate solution to the system of integral equations, we first suppose that for each $n=0,\ldots,N$, the auxiliary function $u^{(n)}$ can be expanded in terms of an appropriately chosen basis $\{v_j^{(n)}:(-H,-d^{(n)})\to \mathbb {C}|j\in \mathbb {N}\}$ over the interval $(-H,-d^{(n)})$. We then approximate each auxiliary function $u^{(n)}$ by

(3.22)\begin{equation} u^{(n)}(z) \approx \begin{cases} \displaystyle\sum_{j=1}^{M_{{aux}}} c_j^{(n)}v_j^{(n)}(z) & \mbox{for }z<{-}d^{(n)}\\ 0 & \mbox{for }z\geq{-}d^{(n)}, \end{cases} \end{equation}

where we have truncated each expansion at some sufficiently large $M_{{aux}}\in \mathbb {N}$.

After substitution of (3.22) into the system (3.16)–(3.18), we multiply each expression by $v_p^{(n)}(z)$, for some $p\in \{1,\ldots,M_{{aux}}\}$. Here, $n\in \{0,\ldots,N\}$ is the index of the barrier associated with each integral equation. For instance, (3.16) was derived by applying (2.5a) at $x^{(0)}$, so this expression is multiplied by $v_p^{(0)}(z)$ following substitution of (3.22). Integrating each of these expressions with respect to $z$ over $(-H,-d^{(n)})$ subsequently yields the following linear system of equations:

(3.23)\begin{align} &\sum_{m=0}^{\infty}\mathcal{U}_{mp}^{(0)}\mathcal{K}_{m}^{(0)}\sum_{j=1}^{M_{{aux}}} \mathcal{U}_{mj}^{(0)}c_j^{(0)}+\sum_{m=1}^{\infty}\mathcal{V}_{mp}^{(1)} \left(\mathcal{K}_{m}^{(1)}\sum_{j=1}^{M_{{aux}}} \mathcal{V}_{mj}^{(1)}c_j^{(0)}-\mathcal{L}_{m}^{(1)} \sum_{j=1}^{M_{{aux}}}\mathcal{U}_{mj}^{(1)}c_j^{(1)}\right)\nonumber\\ &\quad +\frac{\mathrm{i}\omega s^{(1)}}{2h^{(1)}}\left(\zeta_p^{(1)}-\frac{1}{4}(l^{(1)})^{2}\gamma_{p}^{(1)}\right) -A_0^{(1)}\gamma_{p}^{(1)} +\frac{1}{2}B_0^{(1)}l^{(1)}\gamma_{p}^{(1)}={-}2\mathcal{U}_{0p}^{(0)}A_0^{(0)}, \end{align}

(3.24)\begin{align} &-\sum_{m=1}^{\infty} \mathcal{U}_{mp}^{(n)}\left(\mathcal{L}_{m}^{(n)} \sum_{j=1}^{M_{{aux}}}\mathcal{V}_{mj}^{(n)}c_j^{(n-1)}- \mathcal{K}_{m}^{(n)}\sum_{j=1}^{M_{{aux}}}\mathcal{U}_{mj}^{(n)}c_j^{(n)}\right)\nonumber\\ &\quad +\sum_{m=1}^{\infty} \mathcal{V}_{mp}^{(n+1)}\left(\mathcal{K}_{m}^{(n+1)} \sum_{j=1}^{M_{{aux}}}\mathcal{V}_{mj}^{(n+1)}c_j^{(n)}-\mathcal{L}_{m}^{(n+1)} \sum_{j=1}^{M_{{aux}}}\mathcal{U}_{mj}^{(n+1)}c_j^{(n+1)}\right)\nonumber\\ &\quad +\frac{\mathrm{i}\omega}{2}\left(\frac{s^{(n+1)}}{h^{(n+1)}} \left(\zeta_p^{(n)}-\frac{1}{4}(l^{(n+1)})^{2} \gamma_{p}^{(n)}\right)- \frac{s^{(n)}}{h^{(n)}}\left(\zeta_p^{(n)}-\frac{1}{4}(l^{(n)})^{2} \gamma_{p}^{(n)}\right)\right)\nonumber\\ &\quad+A_0^{(n)} \gamma_{p}^{(n)} +\frac{1}{2}B_0^{(n)}l^{(n)} \gamma_{p}^{(n)}-A_0^{(n+1)} \gamma_{p}^{(n)} +\frac{1}{2}B_0^{(n+1)}l^{(n+1)} \gamma_{p}^{(n)}=0, \end{align}

(3.25)\begin{align} &-\sum_{m=1}^{\infty} \mathcal{U}_{mp}^{(N)}\mathcal{L}_{m}^{(N)}\sum_{j=1}^{M_{{aux}}} \mathcal{V}_{mj}^{(N)}c_j^{(N-1)}+\sum_{m=1}^{\infty} \mathcal{U}_{mp}^{(N)}\mathcal{K}_{m}^{(N)}\sum_{j=1}^{M_{{aux}}} \mathcal{U}_{mj}^{(N)}c_j^{(N)}\nonumber\\ &\quad +\sum_{m=0}^{\infty} \mathcal{V}_{mp}^{(N+1)}\mathcal{K}_{m}^{(0)}\sum_{j=1}^{M_{{aux}}} \mathcal{V}_{mj}^{(N+1)}c_j^{(N)}-\frac{\mathrm{i}\omega s^{(N)}}{2h^{(N)}}\left( \zeta_{p}^{(N)}-\frac{1}{4}(l^{(N)})^{2} \gamma_{p}^{(N)}\right)\nonumber\\ &\quad +A_0^{(N)} \gamma_{p}^{(N)} +\frac{1}{2}B_0^{(N)}l^{(N)} \gamma_{p}^{(N)}=2 B_0^{(N+1)}\mathcal{V}_{0p}^{(N+1)}, \end{align}

where (3.24) holds for all $n\in \{1,\ldots,N-1\}$ and we have defined

(3.26a––c)\begin{equation} \mathcal{K}_{m}^{(0)}=\frac{\mathrm{i}}{k_m H},\quad \mathcal{K}_{m}^{(n)}=\frac{{{\rm coth}}(\kappa_m^{(n)}l^{(n)})}{m{\rm \pi}},\quad \mathcal{L}_{m}^{(n)}=\frac{\operatorname{csch}(\kappa_m^{(n)}l^{(n)})}{m{\rm \pi}}. \end{equation}

Additionally, we have defined

(3.27a,b)$$\begin{gather} \mathcal{U}^{(0)}_{mj}=\mathcal{P}^{(0)}(\psi_m^{(0)},v_j^{(0)}), \quad \mathcal{V}^{(N+1)}_{mj}=\mathcal{P}^{(N)}(\psi_m^{(0)},v_j^{(N)}), \end{gather}$$

(3.28a,b)$$\begin{gather}\mathcal{U}^{(n)}_{mj}=\mathcal{P}^{(n)}(\psi_m^{(n)},v_j^{(n)}),\quad \mathcal{V}^{(n)}_{mj}=\mathcal{P}^{(n-1)}(\psi_m^{(n)},v_j^{(n-1)}), \end{gather}$$

(3.29a,b)$$\begin{gather}\gamma_{j}^{(n)}=\int_{{-}H}^{{-}d^{(n)}}v_j^{(n)}(\xi)\, {{\rm d}} \xi, \quad \zeta_j^{(n)}=\int_{{-}H}^{{-}d^{(n)}}(\xi+H)^{2}v_j^{(n)} (\xi)\, {{\rm d}} \xi. \end{gather}$$

The singularities at $(x^{(n)},-d^{(n)})$ motivate the following choice of basis functions:

(3.30)\begin{equation} v_j^{(n)}(\xi) = \frac{2({-}1)^{j-1}}{{\rm \pi} \sqrt{(H-d^{(n)})^{2}-(H+\xi)^{2}}}T_{2(j-1)} \left(\frac{H+\xi}{H-d^{(n)}}\right), \end{equation}

in which $T_n$ is the Chebyshev polynomial of the first kind of degree $n$. This basis was first introduced by Porter & Evans (1995) to solve the scattering of linear waves by a single surface-piercing barrier. By using (1.10.2) in Bateman (1954), the integrals in (3.27a,b) and (3.28a,b) can be evaluated as

(3.31)\begin{equation} \left.\begin{array}{c} \mathcal{U}_{mj}^{(0)}=({-}1)^{j-1}\beta_{m}^{{-}1/2} I_{2(j-1)}(k_{m}(H-d^{(0)})),\\[5pt] \mathcal{V}_{mj}^{(N+1)}=({-}1)^{j-1}\beta_{m}^{{-}1/2} I_{2(j-1)}(k_{m}(H-d^{(N)})),\\[3pt] \mathcal{U}_{mj}^{(n)}=\sqrt{2}J_{2(j-1)}(\kappa_m^{(n)}(H-d^{(n)})),\\[5pt] \mathcal{V}_{mj}^{(n)}=\sqrt{2}J_{2(j-1)}(\kappa_m^{(n)}(H-d^{(n-1)})), \end{array}\right\} \end{equation}

where $J_n$ is the Bessel function of the first kind of order $n$ and $I_n$ is the corresponding modified Bessel function. Furthermore, by applying the symmetry and orthogonality of the even Chebyshev polynomials, the integrals in (3.29a,b) can be evaluated as

(3.32)\begin{equation} \gamma_{j}^{(n)}=\delta_{1j}\quad\mbox{and}\quad \zeta_j^{(n)}=\frac{(H-d^{(n)})^{2}}{4}\begin{cases} 2 & j=1\\ -1 & j=2\\ 0 & j\geq 3. \end{cases} \end{equation}

Equations (3.23), (3.24) and (3.25) form a linear system of $(N+1)M_{{aux}}$ equations with $(N+1)M_{{aux}}+3N$ unknowns. To obtain a further $2N$ equations, we once again apply (2.5b) as $x \to x^{(n-1)+}$ and $x \to x^{(n)-}$. However, in contrast to the derivation of (3.14) and (3.15), we use the fact that $\psi _m^{(n)}$ is orthogonal to the constant function $1$ for all $m\in \mathbb {N}$. After referring to (3.29a,b), we have

(3.33)$$\begin{gather} \sum_{j=1}^{M_{{aux}}}\gamma_{j}^{(n-1)}c_j^{(n-1)}+ \frac{\mathrm{i}\omega s^{(n)}l^{(n)}}{2}-B_0^{(n)}h^{(n)}=0, \end{gather}$$

(3.34)$$\begin{gather}\sum_{j=1}^{M_{{aux}}}\gamma_{j}^{(n)}c_j^{(n)}- \frac{\mathrm{i}\omega s^{(n)}l^{(n)}}{2}-B_0^{(n)}h^{(n)}=0, \end{gather}$$

for $n\in \{1,\ldots,N\}$.

The remaining $N$ equations are derived from the equations of motion of the pistons (3.2). The integral of the hydrodynamic pressure is evaluated by invoking symmetry about ${x=m^{(n)}}$, then expressed in terms of the coefficients of the auxiliary functions by referring to (3.14), (3.22) and (3.28a,b):

(3.35)\begin{align} &\int_{x^{(n-1)}}^{x^{(n)}}\phi(x,-D^{(n)}){{\rm d} x}={-}\frac{\mathrm{i}\omega s^{(n)}}{2}h^{(n)}l^{(n)}+\frac{\mathrm{i}\omega s^{(n)}(l^{(n)})^{3}}{24~h^{(n)}}+l^{(n)}A_0^{(n)}\nonumber\\ &\quad +\frac{\sqrt{2}h^{(n)}}{{\rm \pi}^{2}} \sum_{m=1}^{\infty}\frac{({-}1)^{m}}{m^{2}}\left(\sum_{j=1}^{M_{{aux}}} \mathcal{U}_{mj}^{(n)}c_j^{(n)}-\sum_{j=1}^{M_{{aux}}} \mathcal{V}_{mj}^{(n)}c_j^{(n-1)}\right). \end{align}

After some rearrangement, substitution of (3.1) and (3.35) into (3.2) yields

(3.36)\begin{equation} \sigma^{(n)}s^{(n)}-\mathrm{i}\omega\rho l^{(n)}A_0^{(n)}-\frac{\mathrm{i}\omega\rho\sqrt{2}h^{(n)}}{{\rm \pi}^{2}} \sum_{m=1}^{\infty}\frac{({-}1)^{m}}{m^{2}}\left(\sum_{j=1}^{M_{{aux}}} \mathcal{U}_{mj}^{(n)}c_j^{(n)}-\sum_{j=1}^{M_{{aux}}} \mathcal{V}_{mj}^{(n)}c_j^{(n-1)}\right)=0, \end{equation}

where we have defined

(3.37)\begin{equation} \sigma^{(n)}=\rho_A^{(n)}l^{(n)}\omega^{2}+\rho g l^{(n)} -\frac{\omega^{2}\rho h^{(n)}l^{(n)}}{2}+\frac{\omega^{2}\rho(l^{(n)})^{3}}{24h^{(n)}} -\mathrm{i}\omega\mu^{(n)}.\end{equation}

Finally, we truncate the infinite sums in (3.23), (3.24) and (3.25) after $M_{{ker}}\sim 1000$ terms, then combine them with (3.33), (3.34) and (3.37) to obtain a matrix equation of dimension $M_{{aux}}(N+1)+3N$. This allows us to simultaneously solve for the coefficients of the auxiliary functions $c_j^{(n)}$ for all $n\in \{0,\ldots,N\}$ and $j\in \{1,\ldots,M_{{aux}}\}$, as well as the quantities $A_0^{(n)}$, $B_0^{(n)}$ and $s^{(n)}$ for all $n\in \{1,\ldots,N\}$. The size of $M_{{aux}}$ must be chosen in accordance with $M_{{ker}}\sim 1000$. By choosing $M_{{ker}}=1000$ and $M_{{aux}}=50$, we obtain four-digit accuracy in the unknown coefficients $R$, $T$ and $s^{(n)}$. The convergence properties of this class of Galerkin methods are discussed in more detail by Linton & McIver (2001).

The unknown wave amplitudes in regions $\varOmega ^{(0)}$ and $\varOmega ^{(N+1)}$ can then be recovered by substituting (3.22) into (3.12) and (3.13). Referring to (3.27a,b) then gives

(3.38a,b)\begin{equation} B_m^{(0)}=A_0^{(0)}\delta_{0m}+\mathcal{K}_{m}^{(0)}\sum_{j=1}^{M_{{aux}}} \mathcal{U}^{(0)}_{mj}c_j^{0}\quad\mbox{and}\quad A_m^{(N+1)}=B_0^{(N+1)}\delta_{0m}-\mathcal{K}_{m}^{(0)} \sum_{j=1}^{M_{{aux}}}\mathcal{V}^{(N+1)}_{mj}c_j^{N}. \end{equation}

Similarly, the remaining unknown coefficients can be recovered after additionally referring to (3.14) and (3.15):

(3.39)$$\begin{gather} A_m^{(n)}=\frac{{{\rm coth}}(\kappa_m^{(n)}l^{(n)}/2)}{2m{\rm \pi}} \left(\sum_{j=1}^{M_{{aux}}}\mathcal{U}_{mj}^{(n)}c_j^{(n)}- \sum_{j=1}^{M_{{aux}}}\mathcal{V}_{mj}^{(n)}c_j^{(n-1)}\right), \end{gather}$$

(3.40)$$\begin{gather}B_m^{(n)}=\frac{1}{2m{\rm \pi}}\left(\sum_{j=1}^{M_{{aux}}} \mathcal{U}_{mj}^{(n)}c_j^{(n)}+\sum_{j=1}^{M_{{aux}}} \mathcal{V}_{mj}^{(n)}c_j^{(n-1)}\right). \end{gather}$$

4. Energy absorption

Due to the presence of damping, the class of devices introduced in § 3 is capable of absorbing energy from the incident waves. Here, absorption is defined as $\alpha (\omega )= 1-|T|^{2}-|R|^{2}$, where the transmission and reflection coefficients $T$ and $R$ are calculated when the system is forced from the left by a plane wave of frequency $\omega$. By considering the energy flux of the incident and scattered plane waves (Linton & McIver 2001) and the time-average work done by the damping force, we can write

(4.1)\begin{equation} \alpha(\omega)=\sum_{n=1}^{N} \alpha^{(n)}(\omega), \quad\textrm{where }\alpha^{(n)}(\omega)=\frac{\omega\mu^{(n)}}{\rho k_0 H} \left|\frac{s^{(n)}}{A_0^{(0)}}\right|^{2}. \end{equation}

This second form allows us to consider the contribution of each piston to the absorption by the entire device. It is clear that the minimum value of absorption is $\alpha =0$, which occurs, for example, in the absence of any damping. Further, the maximum value of absorption is $\alpha =1$, which occurs when no energy is scattered by the device (i.e. $R=T=0$).

4.1. Absorption by a single cavity

We first consider absorption by a single cavity absorber (SCA), here defined as a single piston surrounded by a pair of barriers. When the two barriers are identical, an established theoretical result implies that the maximum value of absorption is $\alpha =0.5$ (Budal & Falnes 1975; Evans 1976; Mei 1976; Newman 1977). For any given frequency, this maximum can be attained by tuning the damping coefficient of the piston and the submergence of the barriers. This tuning is simplified by fixing $l^{(1)}=d^{(0)}/5$ and $\rho _A^{(1)}=500$ kg m$^{-2}$. The parameters which maximise absorption for the frequencies considered in § 2.3 are presented in table 3, where we use a tilde to distinguish parameters that relate to SCAs. We observe that the submergence depths are nearly identical to those in table 1, which were obtained by maximising the time-average power of the free-surface oscillations. This is because the fundamental resonance considered in § 2.3 is essentially a vertical fluid motion between the barriers, which is the fluid motion that maximises energy absorption by a heaving piston.

Table 3. The optimised parameters of a SCA which maximise absorption at $\omega ^{(n)}$, i.e. $\alpha (\omega ^{(n)})=0.5$ is attained.

We remark that larger values of $\alpha$ can be attained when $d^{(1)}>d^{(0)}$. Perfect absorption (i.e. $\alpha =1$) occurs when the wave radiated by the piston destructively interferes with the scattered wave, which is impossible for a symmetric device oscillating in one mode of motion (see Falnes (2002) for more details). However, in the limiting case where $d^{(1)}=H$, the scattered wave propagates in one direction only as no transmission is allowed. This implies that complete destructive interference of the scattered wave by the radiated wave produced by a heave-restricted piston is possible. This limiting case is analogous to the oscillating water columns positioned against a wall considered by Evans & Porter (1995), where perfect absorption at a given frequency was achieved.

4.2. Optimised absorption at $N$ frequencies

Consider now an array of absorbing cavities whose absorption spectrum has maxima at the $N=9$ prescribed frequencies considered in § 4.1. We construct this array by modifying the algorithm which was used by Jiménez et al. (2017) to design an acoustic RA. The algorithm is a sequence of optimisation problems, whereby the size of the device and the set of frequencies at which absorption is maximised both increase with each iteration. The details of the algorithm are provided in the following paragraph.

For $n$ iteratively incrementing from 1 to $N$, we maximise $\sum _{j=1}^{n} \alpha (\omega ^{(N-j+1)})$ for a device consisting of $n$ cavities, by tuning $d_n^{(j)}$ for $j\in \{0,\ldots,n\}$, and $\mu _n^{(j)}$ for $j\in \{1,\ldots,n\}$. Here, the subscript $n$ refers to the stage of the iteration. The optimisation is performed while fixing ${l_n^{(j)}=(d_n^{(j-1)}+d_n^{(j)})/10}$ for $j\in \{1,\ldots,n\}$, and with the constraints ${\tilde {d}^{(N-n)}< d_n^{(0)}<\dots < d_n^{(n)}<0.8H}$. When $n=N$, these constraints instead become ${2D^{(0)}< d_N^{(0)}<\dots < d_N^{(N)}<0.8H}$. For $n=1$ the initial guesses of the parameters are $d_n^{(0)}=d_n^{(1)}=\tilde {d}^{(N)}$ and $\mu _n^{(1)}=\tilde {\mu }^{(N)}$. At each subsequent stage $n>1$, we initialise the optimisation with ${d_n^{(0)}=\tilde {d}^{(N-n+1)}}$, ${d_n^{(j)}=d_{n-1}^{j-1}}$ for $j\in \{1,\ldots,n\}$, $\mu _n^{(1)}=\tilde {\mu }^{(N-n+1)}$ and $\mu _n^{(j)}=\mu _{n-1}^{(j-1)}$ for $j\in \{2,\ldots,n\}$. In other words, a new cavity which maximises absorption near $\omega ^{(N-n+1)}$ is added to the left of the array constructed during step $n-1$.

The absorption spectra of all of the intermediate devices obtained by this process are given in figure 6, and the parameters of the final device with $n=9$ cavities are given in table 4. As desired, we find that the peaks of the spectra of all intermediate devices coincide with the prescribed frequencies. In analysis not presented here, we find that at the $j$th peak of each spectrum, the largest contribution to absorption is by the $(n-j+1)$th piston. For the final device, we obtain $\alpha (\omega ^{(j)})>0.98$ for $j\in \{1,\ldots,N-1\}$. The constraint that $d^{(N)}<0.8H$ (which was imposed to ensure the solution converges rapidly) hinders absorption at the lowest frequency, where we obtain $\alpha (\omega ^{(N)})\approx 0.93$. In contrast, Jiménez et al. (2017) were able to obtain an acoustic device for which $\alpha =1$ at all $N$ frequencies because the viscothermal energy loss mechanism included in their model is less restrictive. We also observe that the peaks of the spectrum in figure 6(i) become increasingly sharp as $\omega$ increases. As a consequence, absorption becomes increasingly lower in the intervals between the prescribed frequencies. Our next task is to rectify this by using a broadband absorption metric to design a RA, whose absorption spectrum remains close to $1$ over the entire interval $[\omega ^{(N)},\omega ^{(1)}]$.

Figure 6. Absorption spectra of the devices obtained using the steps outlined in § 4.2, which is a modification of that used by Jiménez et al. (2017). The frequencies at which absorption is maximised at each stage of the algorithm are denoted by vertical dotted lines. As desired, the peaks of the spectra coincide with the desired peaks of absorption.

Table 4. Parameters of the device obtained using the modified algorithm of Jiménez et al. (2017). The absorption spectrum of this device is shown in figure 6(i).

4.3. Optimised broadband absorption

Using a similar metric to Jiménez et al. (2017), we quantify the broadband performance of the RA as the average value of $\alpha$ over the interval $[\omega ^{(N)},\omega ^{(1)}]$, i.e.

(4.2)\begin{equation} E=\frac{1}{\omega^{(1)}-\omega^{(N)}}\int_{\omega^{(N)}}^{\omega^{(1)}} \alpha(\omega)\, {{\rm d}} \omega. \end{equation}

It is clear that $E\in [0,1]$. Naively, we can increase $E$ by increasing the length of the cavities, since this decreases the $Q$-factors of the resonant peaks and allows them to overlap more. Broadband performance can also be increased by adding more appropriately tuned cavities to the RA (i.e. increasing $N$), since this increases the number of absorption peaks in a given interval. Both of these avenues to increasing $E$ also increase the total length of the RA given by $L=\sum _{n=1}^{N} l^{(n)}$, which is undesirable due to practical considerations. This motivates the present problem of maximising $E$ while constraining $L$, where we will fix $N=9$.

The maximisation of $E$ is conducted by simultaneously tuning the parameters $d^{(n)}$, $\mu ^{(n)}$ and $l^{(n)}$ over a bounded feasible set which is defined by inequality constraints. Specifically, we require that $L<\lambda _{0.8}/5$, where $\lambda _{0.8}\approx 86.4$ m is the wavelength of a plane wave with frequency $\omega ^{(N)}=0.8$ s$^{-1}$. We also require that $l^{(n)}$ and $\mu ^{(n)}$ are increasing, that $2D^{(1)}< d^{(0)}<\cdots < d^{(N)}<0.8H$ and that $\mu ^{(n)}<2\tilde {\mu }^{(n)}$ for all $n\in \{1,\ldots,9\}$. The resulting constrained optimisation problem is then solved using the MATLAB routine fmincon. In preliminary tests, the introduction of constraints greatly reduced the runtime of the optimisation algorithm by reducing the number of iterations required to achieve convergence. The numerical evaluation of the objective function $E$ is performed using trapezoidal quadrature. The initial guess is the device with $9$ cavities obtained in § 4.2, which has an absorption peak at $9$ discrete frequencies and for which $E\approx 0.787$. In particular, we remark that this initial guess satisfies the constraints. The RA which results from this optimisation procedure is a local maximiser of $E$, but should not be considered a global maximiser.

The parameters of the optimised RA, which attains a value of $E\approx 0.982$, are given in table 5, and its absorption spectrum is shown in figure 7. The spectrum is predictably flat over the interval $[0.8,1.6]$ s$^{-1}$, and we observe two additional absorption peaks at $\omega \approx 1.84$ s$^{-1}$ and $\omega \approx 2.07$ s$^{-1}$. By using (4.1) to partition the area under the absorption spectrum in figure 7, we find that the two additional peaks are the resonances of the third and second cavity respectively. Additionally, the first three pistons have a minor contribution to absorption over the target interval, which can be demonstrated by computing

(4.3)\begin{equation} \frac{1}{\omega^{(1)}-\omega^{(N)}}\sum_{n=1}^{3}\int_{\omega^{(N)}}^{\omega^{(1)}} \alpha^{(n)}(\omega)\, {{\rm d}} \omega<0.035. \end{equation}

However, the RA obtained by removing the first three cavities, with absorption spectrum given by the dotted line in figure 8(a), only attains a value of $E\approx 0.909$. Therefore the purpose of first three cavities cannot be understood in terms of absorption alone.

Figure 7. The absorption spectrum (solid curve) of the optimised RA. The coloured areas under the curve show the relative contribution of each piston to the total absorption by the device.

Table 5. Parameters of the device designed by maximising the broadband performance metric $E$. The absorption spectrum of this device is shown in figure 7.

Figure 8. (a) The absorption spectrum of the optimised RA where the first three cavities have been removed (dotted line) and where the damping coefficients of the first three pistons are changed to zero (solid line). The absorption spectrum of the unaltered RA also shown in figure 7 is included for comparison (dashed line). (b) The absorption spectra of devices (i) (dotted line) and (ii) (solid line) described in the text. As was the case in (a), device (ii) differs from device (i) by the inclusion of three non-absorbing cavities.

Indeed, the primary role of cavities $1$–$3$ is to intensify the local energy amplification in cavities $4$–$9$, which results in greater absorption by the pistons in these cavities. We interpret that this is because the coupling between the incident plane wave and the resonant mode in the cavity where the amplification occurs is strengthened by the inclusion of cavities $1$–$3$. To illustrate this, we consider (i) a SCA with the same parameters as cavity $4$ of the optimised RA and (ii) an array of four cavities consisting of the first four cavities of the optimised RA, but where $\mu ^{(n)}=0$ for $n\in \{1,2,3\}$. Hence both devices contain a single absorbing cavity whose parameters are identical, but device (ii) additionally contains three non-absorbing cavities positioned to the left. Cavities $5$–$9$ have been omitted from both devices in order to simplify their absorption spectra, which are contrasted in figure 8(b). Importantly, the absorption peak of device (ii) is significantly higher than that of device (i), which suggests that the resonant mode in the absorbing cavity of device (ii) is more strongly coupled to the incident plane wave. This coupling occurs sequentially via the propagating Bloch modes of the three non-absorbing cavities, through which the group velocity gradually decreases (recall from § 2.4 that the group velocity of the resonant mode is close to zero). In device (i), the coupling between the plane wave and the resonant mode must occur directly, which evidently results in a less-intense energy amplification. In figure 8(a), this experiment is also conducted without omitting cavities $5$–$9$ (see solid line), so that both of the devices being compared have six identical absorbing cavities. We observe that absorption over $[0.8,1.6]$ s$^{-1}$ is reduced much less by setting the damping coefficient of the first three pistons to zero than by removing the first three cavities. This further demonstrates the important role that cavities $1$–$3$ have in facilitating the transfer of energy from the incident plane wave to pistons $4$–$9$ in the optimised RA. We conclude that the two additional absorption peaks in figure 7 are a by-product of the geometry that cavities $2$ and $3$ must have in order to optimise this energy transfer.

5. Concluding remarks

Graded arrays of surface-piercing vertical barriers have been proposed as a model for studying the rainbow reflection of two-dimensional linear water waves. An example consisting of 10 barriers was studied, where the parameters were selected so that the fundamental frequency of each successive pair of barriers gradually decreases throughout the array. An incident plane wave causes a local energy amplification in this array in the pair of barriers with that fundamental frequency. Using band diagram calculations for the cognate infinite array problem, we demonstrated that the local energy amplification originates from the rainbow-reflection effect, since broadband wave signals slow down and become spatially separated based on frequency.

The graded array was modified to incorporate absorption by adding damped, heave-restricted pistons between each adjacent pair of barriers. We solved the resulting boundary value problem using an integral equation/Galerkin technique, giving accurate numerical results. We found parameters for devices which maximise absorption over a discrete set of frequencies by modifying the algorithm of Jiménez et al. (2017). Further tuning using a broadband absorption metric produced a device that achieves near-perfect absorption over a prescribed frequency interval. While several of the cavities in the optimised RA absorb a negligible proportion of energy over this interval, they were shown to be indirectly important for transferring energy from the incident wave to other absorbing cavities.

The study of rainbow reflection and broadband energy absorption presented here is a proof of concept that should motivate further exploration of the applications of these phenomena to water waves. For instance, time-domain simulations of broadband wave packets could be used to investigate the duration of local energy amplifications. Further, the effect of obliquely incident waves could be considered by modelling a graded array of resonant absorbers in a three-dimensional fluid domain. Future work should also seek to validate the local energy amplification predicted using linear water wave theory, which could involve the use of computational fluid dynamics software or wave-tank experiments. Finally, more realistic energy absorption models, including models of wave-to-wire energy conversion, could be explored in the context of rainbow absorption for applications to the design of wave energy converters.