2 Preliminaries

We consider a water domain of infinite horizontal extent and finite depth $h$ , bounded below by a flat bed and above by a free surface, and containing a line array of $M$ vertical C-shaped cylinders that extend throughout the water column. A Cartesian coordinate system $\boldsymbol{X}=(x,y,z)$ defines locations in the water domain, where $\boldsymbol{x}=(x,y)$ is the horizontal coordinate and $z$ is the vertical coordinate. The vertical coordinate points directly upwards, with its origin, $z=0$ , coinciding with the equilibrium free surface, and $z=-h$ denoting the bed.

The water velocity field is defined as $\boldsymbol{u}(\boldsymbol{X},t)=\unicode[STIX]{x1D735}\text{Re}\{(gA/\text{i}\unicode[STIX]{x1D714})\unicode[STIX]{x1D719}(\boldsymbol{X})\exp (-\text{i}\unicode[STIX]{x1D714}t)\}$ , where $\unicode[STIX]{x1D719}$ is a (dimensionless) reduced complex-valued velocity potential, $A$ is the incident wave amplitude, $\unicode[STIX]{x1D714}\in \mathbb{R}$ is a prescribed angular frequency and $g\approx 9.81$  m s $^{-2}$ is the constant of gravitational acceleration. The velocity potential satisfies Laplace’s equation throughout the linearised water domain, i.e. 

where ${\mathcal{D}}=\mathbb{R}^{2}\backslash \overline{{\mathcal{C}}}$ and ${\mathcal{C}}$ is the union of the $M$ C-shaped cylinder horizontal cross-sections. It also satisfies the impermeable bed and linearised free-surface conditions,

respectively, and the Sommerfeld radiation condition in the far field, $|\boldsymbol{x}|\rightarrow \infty$ .

Without loss of generality, the array is assumed to lie along the $x$ -axis. Let the cylinders be indexed $m=1,\ldots ,M$ from left to right, and the domain occupied by the $m$ th cylinder be $\boldsymbol{X}\in {\mathcal{C}}_{m}\times (-h,0)$ , where

Here, $x_{m}=(m-1)\,W$ is the cylinder centre location along the $x$ -axis, $a_{m}$ is its radius and $\unicode[STIX]{x1D711}$ is the half-angle of its opening (identical for all cylinders), with the openings at the left-hand end of the cylinders and symmetric about the $x$ -axis, as shown in figure 1. The velocity potential satisfies a no-normal-flow condition at the cylinder surfaces, i.e. 

where $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}n\equiv \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ and $\boldsymbol{n}$ is the normal vector to the cylinder surfaces, together with a condition ensuring the correct singularity at the tips of the C-shape.

Motions are forced by a plane incident wave, with velocity potential

where $\unicode[STIX]{x1D713}$ is the incident wave direction with respect to the positive $x$ -axis ( $\unicode[STIX]{x1D713}=0$ for head-on incidence, as in figure 1), and the wavenumber $k=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}\in \mathbb{R}_{+}$ satisfies the dispersion relation $k\tanh (k\,h)=\unicode[STIX]{x1D714}^{2}/g$ , and is therefore used as a proxy for frequency. The full wave field is $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{inc}+\unicode[STIX]{x1D719}_{sca}$ , where $\unicode[STIX]{x1D719}_{sca}$ is the scattered wave field.

A modified version of the method outlined by Bennetts, Peter & Montiel (Reference Bennetts, Peter and Montiel2017) is used to solve the problem. To implement the method, the $x$ -axis containing the cylinders is divided into $M$ contiguous subintervals $(x_{m-1/2},x_{m+1/2})$ for $m=1,\ldots ,M$ , where $x_{m\pm 1/2}=x_{m}\pm W/2$ , so that the $m$ th subinterval contains cylinder $m$ only. The wave field in the $m$ th subinterval is expressed in the directional-spectrum form

for $0\leqslant \pm (x-x_{m})\leqslant W/2$ and $\boldsymbol{x}\notin \overline{\unicode[STIX]{x1D6FA}}_{m}$ , where $\unicode[STIX]{x1D6FA}_{m}=\{\boldsymbol{x}:(x-x_{m})^{2}+y^{2}<a_{m}\}$ , $\unicode[STIX]{x1D6E4}_{-}=\{-\unicode[STIX]{x03C0}/2+\text{i}\unicode[STIX]{x1D6FE}:\unicode[STIX]{x1D6FE}\in \mathbb{R}_{+}\}\cup \{\unicode[STIX]{x1D6FE}\in \mathbb{R}:-\unicode[STIX]{x03C0}/2\leqslant \unicode[STIX]{x1D6FE}\leqslant \unicode[STIX]{x03C0}/2\}\cup \{\unicode[STIX]{x03C0}/2-\text{i}\unicode[STIX]{x1D6FE}:\unicode[STIX]{x1D6FE}\in \mathbb{R}_{+}\}$ and $\unicode[STIX]{x1D6E4}_{+}=\unicode[STIX]{x1D6E4}_{-}+\unicode[STIX]{x03C0}$ . The amplitude functions $A_{m}^{\pm }$ and $B_{m}^{\pm }$ satisfy the scattering relations

where ${\mathcal{R}}_{m}^{\pm }$ and ${\mathcal{T}}_{m}^{\pm }$ are the operators

The scattering kernels, reflection $R_{m}^{\pm }$ and transmission $T_{m}^{\pm }$ respectively, are

and

where $D_{p,q}^{(m)}$ are the entries of the diffraction transfer matrix/operator for the $m$ th cylinder, calculated numerically using the method of Montiel et al. (Reference Montiel, Chung, Karimi and Kessissoglou2017). Using relations (2.6), the amplitude functions can be mapped from the left ( $-$ ) to the right ( $+$ ) of their subintervals via

where ${\mathcal{P}}_{m}$ is the monodromy operator for the $m$ th subinterval (in terms of the amplitude functions). Noting that $A_{m+1}^{-}(\unicode[STIX]{x1D712})=B_{m}^{+}(\unicode[STIX]{x1D712})$ and $B_{m+1}^{-}(\unicode[STIX]{x1D712})=A_{m}^{+}(\unicode[STIX]{x1D712})$ for $m=1,\ldots ,M-1$ , amplitude functions in any two subintervals, the $r$ th and $s$ th subintervals say, where $r<s$ , are thus related via

For numerical computations, the scattering kernels are sampled at uniformly distributed points along the real branches of the contours $\unicode[STIX]{x1D6E4}_{\pm }$ and truncated versions of their complex branches, thus mapping the operators ${\mathcal{R}}_{m}^{\pm }$ and ${\mathcal{T}}_{m}^{\pm }$ to so-called scattering matrices, and the monodromy operators ${\mathcal{P}}_{m}$ to the so-called transfer (or monodromy) matrices $\unicode[STIX]{x1D64B}_{m}$ ( $m=1,\ldots ,M$ ). The amplitude functions, and hence solution, are then found recursively, using the algorithm of Bennetts & Squire (Reference Bennetts and Squire2009), and setting $A_{1}^{-}(\unicode[STIX]{x1D712})$ and $A_{M}^{+}(\unicode[STIX]{x1D712})$ according to the incident wave (2.4). The Bennetts et al. (Reference Bennetts, Peter and Montiel2017)-esque method described above is chosen over the standard method based on Graf’s addition formula (Kagemoto & Yue Reference Kagemoto and Yue1986; Peter & Meylan Reference Peter and Meylan2004) to gain insights into the mechanisms underlying the large amplifications shown in figure 1, as explained in § 3.

3 Rayleigh–Bloch waves, dispersion curves and quasi-bandgaps

Figure 2 shows results that provide the above-mentioned insights into the large amplifications observed in figure 1, for which the array is defined by $M=10$ , $W=15$  m, $\unicode[STIX]{x1D711}=0.1\unicode[STIX]{x03C0}$ and $a_{m}=a_{1}(m+M-2)/(M-1)$ , where $a_{1}=3.25$  m. The largest amplification case from figure 1, with $\unicode[STIX]{x1D706}=88$  m, is chosen as an example, and figure 2(a) is a magnified version of figure 1(e). The solutions are calculated using 101 sampled points along the real branches of $\unicode[STIX]{x1D6E4}_{\pm }$ and 200 points for each complex branch, along with the truncation $\text{Im}(\unicode[STIX]{x1D712})\leqslant 4$ .

Figure 2.

(a) Log of normalised energy as in figure 1(e). (b) Corresponding energy with the scattered wave field approximated by the local Rayleigh–Bloch-wave components only. (c–g) Eigenvalue spectra of transfer matrices operating over (c) cylinder  $m=2$ , (d) 5, (e) 7, (f) 9 and (g) 10, where the closely spaced black dots approximate the continuous spectra and the red dots define the discrete spectra corresponding to Rayleigh–Bloch waves.

As shown by Thompson, Linton & Porter (Reference Thompson, Linton and Porter2008) and others for uniform line arrays of ordinary cylinders (with $\unicode[STIX]{x1D711}=0$ , i.e. full cylinders of circular cross-section), the scattered wave field along the array is dominated by so-called Rayleigh–Bloch waves, which propagate in both directions along the array and decay exponentially in the transverse direction (i.e.  $y$ -direction) away from it. As explained below, line arrays of C-shaped cylinders also support Rayleigh–Bloch waves, and the radius grading causes their properties to evolve along the array. Figure 2(b) is similar to figure 2(a), but with the scattered wave field approximated by the Rayleigh–Bloch-wave components only. The approximation is based on projecting the scattered wave field in each subinterval onto the eigenfunctions defined by the corresponding transfer matrix, and retaining only the eigenfunctions associated with local Rayleigh–Bloch waves, i.e. Rayleigh–Bloch waves associated with a particular cylinder and spacing. Notwithstanding the small discontinuities, which are an inevitable consequence of the approximation method, the approximation is highly accurate, particularly with respect to the large amplifications, thus confirming that the amplifications are due to excitation of local Rayleigh–Bloch waves.

Figure 2(c–g) shows eigenvalue spectra of the transfer matrices, $\unicode[STIX]{x1D707}\in \text{eig}\{\unicode[STIX]{x1D64B}_{m}\}$ , for $m=2$ , 5, 7, 9 and 10, in the complex plane. The closely spaced black dots along an arc of the unit circle approximate the continuous spectrum, corresponding to plane wave forcing, i.e.  $\unicode[STIX]{x1D707}=\exp \{\text{i}kW\cos \unicode[STIX]{x1D712}\}$ for $\unicode[STIX]{x1D712}\in (-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0})$ , and with eigenfunctions that display delta-function-like behaviour around the incident direction (not shown; see figure 4 by Bennetts et al. Reference Bennetts, Peter and Montiel2017). Eigenvalues in the continuous spectrum corresponding to evanescent wave forcing are omitted from the plots for clarity. The two red dots represent the discrete spectrum, corresponding to Rayleigh–Bloch waves, i.e.  $\unicode[STIX]{x1D707}=\exp \{\pm \text{i}\,\unicode[STIX]{x1D6FD}\,W\}$ , where $\unicode[STIX]{x1D6FD}$ is the Rayleigh–Bloch wavenumber, and with regular eigenfunctions (again, see figure 4 by Bennetts et al. Reference Bennetts, Peter and Montiel2017). The eigenvalue in the upper-half complex plane corresponds to a rightward-propagating Rayleigh–Bloch wave ( $+$ ), and the eigenvalue in the lower-half plane corresponds to a leftward-propagating Rayleigh–Bloch wave ( $-$ ). Increase of the sampling of the contours $\unicode[STIX]{x1D6E4}_{\pm }$ increases the number of eigenvalues in the discrete approximation of the continuous spectrum, but does not affect the eigenvalues (or eigenfunctions) defining the Rayleigh–Bloch waves, i.e. they are converged.

Figure 2(c) shows the spectrum for the transfer matrix that maps over cylinder $m=2$ , for which the Rayleigh–Bloch eigenvalues are fractionally displaced along the unit circle from the ends of the continuous spectrum. Figure 2(d–f) shows that the spectra of the transfer matrices farther along the array have identical continuous spectra, but that the Rayleigh–Bloch eigenvalues rapidly move along the unit circle away from the continuous spectrum. As the Rayleigh–Bloch eigenvalues approach $-1$ from above ( $+$ ) and below ( $-$ ), the group velocity of the Rayleigh–Bloch waves reduces, and when the eigenvalues meet at $-1$ , the Rayleigh–Bloch waves become standing waves, known as a Neumann trapped mode (Evans & Porter Reference Evans and Porter1999). Figure 2(g) shows the spectrum for the final cylinder along the array ( $m=10$ ), for which the Rayleigh–Bloch eigenvalues have jumped onto the negative branch of the real line, meaning that $\unicode[STIX]{x1D6FD}\in \unicode[STIX]{x03C0}+\text{i}\mathbb{R}_{+}$ and the Rayleigh–Bloch waves no longer propagate. (The eigenvalue associated with $-\unicode[STIX]{x1D6FD}$ is beyond the axis limits for $m=10$ .)

Therefore, the cylinder-radius grading causes the rightward-propagating Rayleigh–Bloch wave excited at the leading end of the array ( $m=1$ ) to slow down progressively along the array, until it reaches a turning point, at which the group velocity is zero. In the language of waveguide modes, at this point, the wave is cut off and ceases to propagate. The energy carried by the local Rayleigh–Bloch wave accumulates at the turning point, thereby generating large amplifications. This behaviour is broadly analogous to the effect of gradually decreasing the width of an acoustic or ocean waveguide, so that a waveguide mode is eventually cut off, with its energy reflected from the turning point, with propagation in a wedge being the archetypal example (Arnold & Felsen Reference Arnold and Felsen1983).

The location of the Rayleigh–Bloch eigenvalues in the complex plane depends on the frequency of motion, i.e. incident wavelength. In particular, as the frequency increases, the eigenvalues depart the unit circle for smaller cylinder radii, meaning that the local Rayleigh–Bloch waves propagate shorter distances along the array. Thus, as shown in figure 1, energy accumulates closer to the leading end of the array as the frequency increases, i.e. the incident wavelength decreases, so that the array spatially separates frequencies/wavelengths.

Figure 3 shows Rayleigh–Bloch dispersion curves for C-shaped cylinders in periodic line arrays. These curves are valuable for interpreting behaviours along graded arrays as, assuming that the array grading is sufficiently gradual that waves behave locally as if they are in a uniform array, one can then use the dispersion curves to estimate the turning points for different frequencies. The dispersion curves are shown in the first irreducible Brillouin zone $\unicode[STIX]{x1D6FD}\,W\in [0,\unicode[STIX]{x03C0}]$ , and curves for ordinary cylinders are shown for comparison where appropriate, along with the dispersion line (so-called light line) for the bulk medium, $k=\unicode[STIX]{x1D6FD}$ . Figure 3(a) shows dispersion curves corresponding to the array used in figures 1 and 2. Each C-shaped cylinder curve follows the light line for low frequencies $k$ /small wavenumbers $\unicode[STIX]{x1D6FD}$ , but rapidly departs from the light line as the frequency/wavenumber increases, and cuts off ( $\unicode[STIX]{x1D6FD}W=\unicode[STIX]{x03C0}$ ) at a relatively low frequency, just below the resonant frequency for the isolated cylinder, meaning that the cutoff can be tuned by altering the resonance.

The frequency interval occupied by the dispersion curve is conventionally known as a passband in periodic media. In this instance, frequencies in the passband support propagating Rayleigh–Bloch waves. For frequencies above the passband, Rayleigh–Bloch waves do not propagate, which is analogous to a so-called bandgap in periodic media, but is here referred to as a quasi-bandgap due to propagating modes in the continuous spectrum being supported by the array for these frequencies (as shown in figure 2 g). As the cylinder radius increases, the resonant frequency decreases, thus pushing the dispersion curves down, and narrowing the passbands. Therefore, a range of frequencies exists in the passband for the small cylinders at the leading end of the array that transition into quasi-bandgaps as the radius increases along the array. Local Rayleigh–Bloch waves at these frequencies meet a turning point within the array, with the associated energy trapped within the array. Doubling the radius from $a=3.25$  m to 6.5 m almost halves the passband width (reduction factor ${\approx}0.51$ ), indicating that radius grading is an effective way to capture a wide range of frequencies.

Figure 3.

Rayleigh–Bloch dispersion curves for C-shaped cylinders (solid blue curves) and ordinary cylinders (solid grey), with the light line (broken black) given for reference. (a) Varying cylinder radius: $a=3.25$  m (i.e. cylinder  $m=1$ ; dark blue/dark grey lines), $a=4.7$  m ( $m=5$ ; mid-blue/mid-grey lines) and $a=6.5$  m ( $m=10$ ; light blue/light grey lines). (b) Varying cylinder spacing for $a=4.7$  m:  $W=15$  m (dark blue/dark grey lines), $W=18.75$  m (mid-blue/mid-grey lines) and $a=22.5$  m (light blue/light grey lines). (c) Varying opening half-angle: $\unicode[STIX]{x1D711}=\unicode[STIX]{x03C0}/20$ (dark blue line), $\unicode[STIX]{x1D711}=\unicode[STIX]{x03C0}/10$ (mid-blue line) and $\unicode[STIX]{x1D711}=\unicode[STIX]{x03C0}/5$ (light blue line).

For the ordinary cylinders, the dispersion curves follow the light line for a wide range of frequencies/wavenumbers in the Brillouin zone, before departing the light line and cutting off just below $k\,W=\unicode[STIX]{x03C0}$ . Therefore, only relatively high frequencies/short wavelengths fail to excite propagating Rayleigh–Bloch waves along the array, and thus can be captured. Moreover, doubling the radius of the ordinary cylinders from $a=3.25$  m to 6.5 m reduces the passband width by a factor of ${\approx}0.94$ only, and increasing the cylinder from 4.7 m to 6.5 m reduces the width by less than 1 % (curves are virtually indistinguishable). It follows that grading the radius of ordinary cylinders captures only a narrow range of (high) frequencies and is less effective than using the C-shaped resonators.

Figure 3(b) shows the effect of changing the cylinder spacing, $W$ , on the dispersion curves, for a cylinder radius of $a=4.7$  m, i.e. cylinder $m=5$ in figures 1 and 2, noting that different radii produce cognate results. For the C-shaped cylinders, the passband width is insensitive to changes in the array spacing (increasing from $W=15$  m to 22.5 m reduces the width by ${\approx}4\,\%$ only), as the upper limit of the passband is controlled by the resonant frequency for an isolated cylinder, which is independent of the spacing. Therefore, grading the spacing of identical C-shaped cylinders is an ineffective approach to capture a wide range of frequencies. In contrast, quasi-bandgaps for ordinary cylinders are generated by Bragg scattering mechanisms, so that passband widths can be controlled by the array spacing, noting that this approach is used in chirped sonic arrays (Romero-García et al. Reference Romero-García, Picó, Cebrecos, Sánchez-Morcillo and Staliunas2013; Cebrecos et al. Reference Cebrecos, Picó, Sánchez-Morcillo, Staliunas, Romero-García and Garcia-Raffi2014). Increase of the spacing from $W=15$  m to 22.5 m reduces the passband width by a factor of ${\approx}0.71$ , providing a wide frequency capture range, albeit for relatively high frequencies.

Figure 3(c) shows the effect of changing the cylinder opening angle on the dispersion curves, which is relevant for the C-shaped cylinders only. Decrease of the half-angle from $\unicode[STIX]{x1D711}=\unicode[STIX]{x03C0}/5$ to $\unicode[STIX]{x03C0}/20$ decreases the passband width by an appreciable factor of ${\approx}0.77$ , although the reduction is small in comparison with the proportional changes in passband width given by varying the cylinder radius. Therefore, grading the opening angle along the array provides an alternative method to generate turning points along the array, and could potentially be combined with radius grading to enhance the frequency capture width (not pursued here).

4 Amplification spectra

Figure 4 quantifies the overall energy amplification produced by the graded array over ranges of incident wavelengths and directions, using metrics analogous to the Q-factor familiar in assessing energy gains (or losses) given by arrays of interacting wave-energy converters (Falnes Reference Falnes1980). Figure 4(a) shows ${\mathcal{Q}}_{arr}={\mathcal{E}}/{\mathcal{E}}_{inc}$ on a logarithmic scale, where

are respectively the scaled energy contained with the C-shaped cylinders along the array and the energy of the incident field over the same area. Therefore, ${\mathcal{Q}}_{arr}$ quantifies the overall amplification given by the array.

For head-on incidence, $\unicode[STIX]{x1D713}=0$ , the array amplifies the incident energy by over an order of magnitude for wavelengths $\unicode[STIX]{x1D706}\in (60~\text{m},98~\text{m})$ , with maximum amplification ${\mathcal{Q}}_{arr}\approx 10^{1.53}\approx 33.6$ . For wavelengths $\unicode[STIX]{x1D706}<60$  m, the amplification reduces as the wavelength decreases, due to the associated frequencies lying in quasi-bandgaps for the cylinders at the leading end of the array, so that only small quantities of wave energy penetrate the array. For wavelengths shorter than those shown, higher-order resonances alter this simple trend. For wavelengths $\unicode[STIX]{x1D706}>98$  m, the amplification asymptotes towards unity as the influence of the cylinder on the waves reduces. The behaviour is similar for non-head-on incidence, with some reduction in amplification, as Rayleigh–Bloch waves are not as strongly excited due to the loss of symmetry in the incident field with respect to the axis of the array. Amplifications of over an order of magnitude exist up to $\unicode[STIX]{x1D719}\approx 0.18\unicode[STIX]{x03C0}$ , and the wavelength interval for which the amplification is over an order of magnitude is at least 25 m long up to $\unicode[STIX]{x1D719}=0.15\unicode[STIX]{x03C0}$ .

Figure 4.

Logarithmic Q-factors as functions of incident wavelength and direction: (a)  ${\mathcal{Q}}_{arr}={\mathcal{E}}/{\mathcal{E}}_{inc}$ , where ${\mathcal{E}}$ is the integrated energy within the C-shaped cylinders along the array and ${\mathcal{E}}_{inc}$ is the incident energy over the same domain; (b)  ${\mathcal{Q}}_{grd}={\mathcal{E}}/{\mathcal{E}}_{0}$ , where ${\mathcal{E}}_{0}$ is the integrated energy within equivalent isolated C-shaped cylinders.

Figure 4(b) shows ${\mathcal{Q}}_{grd}={\mathcal{E}}/{\mathcal{E}}_{0}$ on a logarithmic scale, where

with $\unicode[STIX]{x1D719}_{m}$ the velocity potential for the $m$ th cylinder in isolation, i.e. with no surrounding cylinders. Therefore, ${\mathcal{Q}}_{grd}$ quantifies the overall energy amplification given by the radius grading, independent of the amplification due to the cylinder resonances. As indicated by figures 1 and 2, for head-on incidence, the grading is most effective for wavelengths that excite resonances in the cylinders towards the trailing end of the array. The maximum overall amplification due to grading is ${\mathcal{Q}}_{grd}\approx 10^{0.77}\approx 5.87$ for $\unicode[STIX]{x1D706}\approx 92$  m. The amplification is positive for $\unicode[STIX]{x1D706}>50$  m, with negative amplifications associated with quasi-bandgaps for shorter waves, as noted above. Similarly, the strength of amplification due to grading slowly decreases as the incident wave direction moves away from head-on incidence, with the grading more than trebling the overall amplification, i.e.  ${\mathcal{Q}}_{grd}>3$ , for wavelengths around $\unicode[STIX]{x1D706}=92$  m up to $\unicode[STIX]{x1D713}\approx 0.15\unicode[STIX]{x03C0}$ .