In this paper, we study the interaction of water waves with a surface-piercing truncated cylindrical meta-structure consisting of two overlapping arrays of closely spaced vertical thin plates. The fluid resonance promoted in the narrow vertical channels formed by the cylindrical meta-structure is exploited by a novel design concept of the wave power converter by covering the surface of the cylinder with an array of small cuboid buoys which float in the gaps between the intersecting plate arrays. Each buoy is attached to its own spring and power takeoff damping mechanism, and the vertical displacement of individual buoys is replaced by a continuous two-dimensional function of position which follows from homogenization of the plate/fluid structure of the cylinder. Effective medium equations and boundary conditions are derived under both full depth-dependent theory and shallow-water theory, allowing semi-analytical methods to be developed to investigate the wave scattering and wave energy absorption properties of this cylindrical meta-structure. Results illustrate that the internal resonance of the cylindrical meta-structure can lead to significant wave power capture across a broad range of frequencies.

## I. INTRODUCTION

When water waves interact with a closely spaced periodic array of fixed rigid elements, they can exhibit behavior not typically observed when marine structures have smooth surfaces. For example, in a series of recent papers,^{1–5} periodic arrays of thin vertical plates protruding from the base of a fluid have been shown to produce refractive effects (including negative refraction) not possible with conventional bathymetry. This is partly due to the contrast in lengthscales between the wavelength and the spacing between elements of structure, but also due to the anisotropy built into the design of the structure allowing waves to experience different depths depending on the wave heading. Adopting terminology used across other areas of physics, this new type of offshore structure is defined as a water wave *meta-structure.*^{6} The so-called cylindrical meta-structure that is formed when the plate array extends fully through the depth and is confined within a cylindrical domain has been considered by Refs. 7 and 8. Apart from its anisotropic scattering character (e.g., incident waves propagating in directions aligned with the plate array experience no scattering), it has been shown that the plate array structure reduces the wave speed inside the cylinder leading to a resonant amplification of the elements in the system within the cylinder which produces a strong lensing effect.^{8}

In the field of fluid dynamics, meta-structures have recently been applied to power extraction from water waves. Vita *et al.*^{9} presented a two-dimensional example of the interaction of surface gravity waves with a wave energy device consisting of an array of periodic submerged harmonic oscillators. By redirecting and accelerating/decelerating the flow in inhomogeneous and anisotropic material, Pang *et al.*^{10} designed an energy harvesting device to capture the kinetic energy of low-speed water flow based on transformation hydrodynamics.^{11} In Ref. 7, a damped surface boundary condition was introduced inside the cylinder to mimic the effect of a wave energy converter (WEC) whereby it was shown that a single cylindrical meta-structure is capable of harnessing multiple times the maximum theoretical wave power of a cylindrical device of an equivalent size operating under rigid body motion (e.g., Refs. 12 and 13). Two of the current authors have been working on a project to investigate how to develop this result by replacing the surface damping condition by practical mechanical mechanisms. In Ref. 14, the damping condition was removed, and power was instead generated by an arrangement of opposing pairs of vertically buoyant hinged paddles distributed midway along each channel through the middle of the structure of the cylinder. It was demonstrated that power well in excess of the equivalent rigid body limits could be generated by a non-structured cylindrical device of the same size. Quite remarkably, it was also shown that power capture characteristics were relatively insensitive to the wave heading despite the anisotropy of the cylindrical meta-structure design.

In this paper, we consider an alternative mechanical method for generating useful power from a cylindrical meta-structure subject to incident waves. The cylindrical meta-structure is now formed by two overlapping pairs of plate arrays at right angles to one another and submerged uniformly through the surface of the fluid to different depths. The doubly periodic arrays of square section cavities that divide the surface of the fluid within the cylinder are placed with floating buoys that are attached to their own spring and damper, allowing power to be generated through the vertical motion of the buoys. This paper presents a theoretical model as a preliminary study into the use of a cylindrical meta-structure as a wave energy converter. We have therefore assumed an ideal fluid and ignored turbulent and viscous losses due to the interaction of the fluid with the structure.

This work has similarities to Garnaud and Mei,^{15,16} who also carried out a theoretical study on wave power extraction by a compact array of small floating buoys absorbing power in heave with spacings much shorter than the typical wavelength. They showed that a high efficiency of wave energy conversion can be achieved compared to the limits upon an equivalent rigid cylinder, and this can be maintained over a wide range of frequencies. The reason for this is that the independent motion of elements in the array allow energy to be captured from multiple Fourier components of the incident wave as opposed to the first two components associated with the heave and surge motion of the equivalent rigid cylinder (see Ref. 17). The present work includes the additional complexity of the plate array structure and can be seen as an extension of the previous study of Porter *et al.*^{18} who considered a single array of bottom-mounted plates extending only partially through the depth. The use of plate arrays to confine arrays of buoys to extract energy through their vertical motion was also recently proposed in Ref. 19 in a two-dimensional setting. They considered a single array of plates of increasing depth to develop multiple closely spaced fluid/mass resonances within the array and numerically demonstrated impressive broadband energy capture. Their solution was represented by a coupled system of integral equations associated with the unknown fluid velocity across the gap under each vertical plate.

The mathematical approach employed here adopted combines Refs. 15 and 18 using homogenization methods to replace exact governing equations and boundary conditions that apply on the microscale by effective medium equations and conditions on the macroscale. In tandem with a full depth-dependent description of the problem, we develop a shallow water approximation, again using homogenization methods, which results in a simpler and more numerically robust implement and makes the role of the physical parameters in determining wave propagation characteristics explicit. Indeed, even after using a homogenization approximation to governing equations, the solution to the full depth-dependent problem is numerically challenging with the effect of resonance in the present problem adding to the difficulties reported in Ref. 18.

The paper is arranged as follows. The modeling and assumptions are described in Sec. II of this paper. Appendixes A and C describe the homogenization approach that is used to derive the effective equations, and these equations are presented at the beginning of Secs. III and IV which describe the full depth-dependent theory and the shallow water theory, respectively. In Sec. V, we describe two independent methods for calculating the power developed by the WEC device, which are used alongside other results to validate the model in Sec. VI. In Sec. VII, we present a number of typical cases to assess the efficacy of the proposed device as a WEC including a comparison with the results of Ref. 15. Finally, a summary of the work is given in Sec. VIII.

## II. PROBLEM STATEMENT

As shown in Fig. 1, a structured cylinder of radius *a* consists of two arrays of parallel vertical thin plates that are overlapping and perpendicular to each other. The two arrays are, respectively, submerged through the free surface to depths *d _{x}* and

*d*in water of constant depth

_{y}*h*and density

*ρ*. We assume that the separation between two adjacent vertical plates,

*L*, is equal throughout both arrays and is small compared to the typical wavelength, resulting in a two-dimensional periodic array of identical open-ended vertical channels formed within the cylinder through which the fluid is allowed to flow. Floating on the surface of each of the square cross section channels within the cylinder is a cuboid buoy with sides of length

*L*and draft

*d*which is connected to its own linear damper with damping rate,

*γ*, and a linear spring with spring constant,

*σ*. The buoys are thus confined to moving in heave only (the vertical direction).

The Cartesian coordinate system, *Oxyz*, is defined with its origin, *O*, in the mean water surface, the *Ox*-axis and *Oy*-axis aligned with the two arrays of plates submerged to depths *d _{x}* and

*d*, respectively. The

_{y}*Oz*-axis coincides with the vertical axis of symmetry of the cylinder and is directed upward. A plane wave with the amplitude

*A*and angular frequency

*ω*is incident at an angle

*β*relative to the positive

*Ox*-axis. Under the action of waves, the buoys oscillate vertically, and energy is extracted via the damper. In our theory, we assume no hydrodynamical or mechanical losses. With no loss of generality, we let $ d x \u2264 d y$ since the incident wave angle is arbitrary.

Additionally, the cylindrical coordinate system, $ Or \theta z$, is employed for mathematical convenience with $ x = r \u2009 cos \u2009 \theta $ and $ y = r \u2009 sin \u2009 \theta $. The entire fluid domain can be divided into an outer region $ \Omega 1 = { r \u2265 a , 0 \u2264 \theta < 2 \pi , \u2212 h \u2264 z \u2264 0}$ and an inner region $ \Omega 2 = { r < a , 0 \u2264 \theta < 2 \pi , \u2212 h \u2264 z \u2264 \u2212 d}$. It is convenient to further divide the inner cylindrical region into three layers: $ \Omega 21 = { r < a , 0 \u2264 \theta < 2 \pi , \u2212 d x \u2264 z \u2264 \u2212 d} , \u2009 \Omega 22 = { r < a , 0 \u2264 \theta < 2 \pi , \u2212 d y \u2264 z < \u2212 d x}$, and $ \Omega 23 = { r < a , 0 \u2264 \theta < 2 \pi , \u2212 h \u2264 z < \u2212 d y}$.

*hydrodynamic*pressure (i.e., in excess of the hydrostatic background pressure, $ \u2212 \rho g z$, in which the acceleration due to the gravity

*g*acts in the negative

*z*-direction).

*L*of each square cell but, under homogenization, will eventually be regarded as a continuous function of the macroscale variables

*x*,

*y*. Thus, the linearized kinematic and dynamic conditions on the bottom of the buoy can be written as

**n**is used to represent the normal into the fluid on fixed surfaces.

Since the separation of adjacent plates, *L*, is assumed to be much smaller than the typical wavelength, i.e., $ \epsilon = L / \lambda \u226a 1$, the method of multiple scales is applied to consider the effect of microstructure. We will also assume that $ L / d x , L / d y \u226a 1$, implying that the periodicity of the plate arrays is significantly smaller than their depth of submergence.

Appendixes A and C outline the application of the multiple scales approach in the case where the fluid depth is arbitrary [ $ \omega 2 h / g = O ( 1 )$] and when the depth is small compared to the wavelength or $ \omega 2 h / g \u226a 1$. The latter case is referred to as the shallow water limit and is considered in Appendix C. Otherwise, we say that we are fully depth-dependent for which the appropriate homogenization is performed in Appendix A.

## III. FULL DEPTH-DEPENDENT THEORY

### A. Governing equations and boundary conditions

*ω*proportional to $ \u2009 exp \u2212 i \omega t$. The superscript (0) has been dropped for convenience, and we have reverted to dimensional quantities. Therefore, the time-independent velocity potential $ \varphi = \varphi k ( x , y , z )$ in each respective region Ω

_{k}has been shown to satisfy

*z*direction, the flow component perpendicular to the plate is suppressed in $ \Omega 22$, and the flow under the cylinder is unconstrained.

*r*=

*a*, between inner and outer regions, continuity of pressure requires

*r*=

*a*results in the piecewise conditions

*r*=

*a*from within the overlapping plate region $ \u2212 d x < z < \u2212 d$.

### B. Expansions for the velocity potential

*J*,

_{n}*H*, and

_{n}*K*denote, respectively, the

_{n}*n*th order Bessel function of the first kind, Hankel function of the first kind, and modified Bessel function of the second kind. Also in the above,

*A*are undetermined coefficients,

_{jn}*k*are the roots of the dispersion equation

_{j}*k*$ ( j \u2265 1 )$ are real while $ k 0 = \u2212 i k$ and the wavenumber

_{j}*k*is real, and $ \psi j ( z )$ are vertical eigenfunctions

*z*. Its final expression can be determined by satisfying the boundary conditions (12b) and (14) and by balancing the pressure and flux on the interfaces $ z = \u2212 d x$ and $ z = \u2212 d y$, leading to

*d*, since the inertial effect of the mass of the buoy is the same as the fluid it displaces. As we proceed, we find the dependence on

*d*disappears from all calculations, as we must therefore expect.

*t*= 0, Eq. (26) is reduced to

*d*. This is since the barriers of depth

_{x}*d*are aligned with the

_{x}*x*-axis and transparent to waves traveling at an angle

*t*= 0. Likewise, for $ t = \pi / 2$, Eq. (26), less obviously becomes

*d*. Finally, if

_{y}*d*=

_{x}*d*, Eq. (26) reduces to Eq. (27) or (28) for all values of

_{y}*t*. In this latter special case, $ \kappa q ( t )$ are constant, and Eq. (23) reduces to the more familiar separation series

*τ*is real, a conventional method can be performed to find the roots $ k q ( t )$, but for the dispersion equation (26), the procedure becomes more complicated, and one robust approach is outlined in Appendix B.

### C. Reduction to a system of equations

## IV. SHALLOW WATER THEORY

### A. Governing equation and boundary conditions

*k*, in shallow water satisfies $ \nu = k 2 h$, which is the $ k h \u2192 0$ limit of the dispersion equation (19) for

*j*= 0. Also,

*τ*is defined in Eq. (13), and the two-dimensional diagonal tensor $h$ is, from Eq. (C22) in dimensional form,

*r*=

*a*are continuous. In terms of our variables, these require

*r*=

*a*for $ 0 \u2264 \theta < 2 \pi $ (see Refs. 4 and 18 for a similar application of shallow water matching conditions).

### B. Expansion, matching, and system of equations

*t*when

*d*=

_{x}*d*which is implied by geometric symmetry.

_{y}*r*=

*a*, and using the orthogonality of functions $ \u2009 exp i n \theta $ over $ 0 \u2264 \theta < 2 \pi $, we obtain

*A*between Eqs. (51) and (52) gives

_{n}## V. WAVE POWER ABSORPTION

*S*, of a cylinder of large radius

_{R}*R*extending through the depth,

*h*, gives the mean power as

*A*and $ c g = \omega / k$ if the shallow water theory is applied.

_{n}*S*through the fluid domain to the internal boundary representing the underside of the floating buoys using Green's second identity. Once this is done carefully using the different governing equations in the different subdomains and the boundary matching conditions that define the problem, we find that the mean power calculated by the near-field method is given by

_{R}*z*-dependence from the integral over

*t*enabling us to write

*d*, in accordance with earlier remarks and with the computation of Eq. (60).

It is now possible to use either Eq. (23) or (49) in Eq. (63), and the final expression for *P _{n}* can be expressed explicitly. However, its numerical calculation is time-consuming since the expression includes a triple summation and a double integration for the full-depth theory and includes a double summation and a double integration under the shallow water theory (when $ \varphi 2$ is replaced by

*ξ*). In spite of the complexity of Eq. (63), it does provide us with an additional check on Eq. (61) to assess the accuracy of numerical calculation.

*W*, defined as the width of the incoming wavefront that contains the same amount of power as that extracted by the WEC (e.g., Ref. 13) can be used as an indicator of effectiveness of wave power extraction written as

## VI. VALIDATION

We first examine the convergence characteristics for the numerical scheme of the full depth-dependent theory given in Sec. III by varying the parameters *M* and *N* representing the number of angular and depth modes retained in the system of equations (40). A cylinder of radius $ a / h = 0.5$ with the two arrays of plates having the depths of submergence $ d x / h = 0.1$ and $ d y / h = 0.2$ subject to incident waves propagating at $ \beta = \pi / 4$ is considered. For the purpose of illustration, we choose a damping $ \gamma / ( \rho g h ) = 0.2$ and a spring $ \sigma / \rho g = 0.2$ as representative values of a real device. Figure 2 shows the variation in the capture width ratio, $ W / ( 2 a )$, calculated using the far-field method against the non-dimensional wavenumber $ \nu d y$. In Fig. 2(a), we fix *N* = 8 and so illustrate convergence with increasing *M*. As the frequency increases, more terms are typically required to adequately model the evanescent waves. In Fig. 2(b), we now fix *M* = 8 and show that convergence with increasing *N* is rapid.

Figure 3 presents the comparison of the wave power evaluated by the near-field and far-field expressions. Here, we set both truncation parameters *N* and *M* equal to 8, and two methods are shown to agree with a high degree of accuracy.

We further validate our numerical model by considering a special case where the geometry of the cylinder is unchanged but the damping rate and spring constant are set to zero (i.e., *τ* = 0) which means the surface of fluid within the cylinder is subject only to gravity. This allows us to compare the present results with those computed using the boundary element method proposed by Ref. 22 in which the scattering of multiple discrete vertical barriers with infinitesimal thickness is considered. As Huang and Porter^{23} have pointed out, the homogenization method can serve as a good approximation to large discrete arrays of plates when the separation of two adjacent plates $ L / d y < 0.5$ if $ \nu d y \u2272 0.4$. Our boundary element method was used to calculate wave interaction with 24 plates in both *x* and *y* directions. Figure 4 plots the contour of wave amplitude in the wave field. It illustrates that the homogenization method accurately captures the fluid motion in both interior and exterior regions. Specifically, compared with the boundary element method, the results from homogenization slightly overpredict the wave amplitudes within the cylinder, and the free surface elevation is not continuous at the interface of interior and exterior regions [see subplots (a) and (c)] since no pressure continuity condition is imposed on $ \u2212 d x < z < 0$ [see Eq. (15)]. On the other hand, in the boundary element method, there should be discontinuities in the surface elevation across each of the 48 plate elements in the interior domain since each plate is exactly described in the numerical computation. However, these are not evident in Fig. 4. Therefore, this also provides a justification for the homogenization approach that assumes that the quantities vary continuously in the effective medium.

Next, a comparison is made between the full depth-dependent theory and the shallow water theory. Both are homogenization approximations, and we expect good agreement for values of $ \nu h \u226a 1$. A cylinder with the same settings used for Fig. 3 is considered in Fig. 5, which plots curves of the capture width ratio $ W / ( 2 a )$ vs the non-dimensional wavenumber $ \nu d y$. The results show increasing agreement as $ \nu h \u2192 0$, as expected. The corresponding pressure distribution on the bottom of the buoy at $ \nu d y = 0.1$ (corresponding to $ \nu h = 0.5$) is given in Fig. 6 (note the compressed vertical scale). Since a long wave is considered, the hydrodynamic pressure is almost unchanged in the array.

The wave amplitude for the same cylinder given above but without internal buoys at $ \nu d y = 0.1$ is plotted in Fig. 7. Although the results from these two theories are similar, the full depth-dependent theory predicts larger wave amplification in and around the cylinder. This may be because the first-order shallow water theory does not include the inertial effects of the fluid or floating buoy in the internal domain $ \u2212 d < z < 0$ (see the neglect of $M$ in Appendix C). Indeed, with *τ* = 0, the effective shallow water governing equation Eq. (44) and effective boundary condition Eq. (47) are identical to Ref. 18 for a bottom-mounted structured cylinder with the reduced depths $ h \u2212 d y$ and $ h \u2212 d x$ replacing the plate submergence depths *D* and *d* in Ref. 18.

## VII. RESULTS

In this section, we focus on the operation of the cylindrical meta-structure as a WEC device and consider a range of parameters that are indicative of the typical performance of what we imagine to be a realistic device and wave conditions. In particular, the choice $ d y / h = 0.2$ and $ a / h = 0.5$ is suggestive of a device of 10 m in radius having plates submerged to a depth 4 m. The maximum upper value of $ \nu d y = 1$ used in the plots is then suggestive of a minimum wavelength of approximately 25 m.

First, Fig. 8 plots the contour of wave scattering of three cylindrical meta-structures of $ a / h = 0.5$ and $ d y / h = 0.2$ with different values of *d _{x}* for an incident wave angle $ \beta = \pi / 4$ when $ \nu d y = 0.4$ and

*τ*= 0. From Appendix B, we have known that as

*d*approaches

_{x}*d*, the wavelength in the structure will decrease, and the oscillation within the cylinder becomes increasingly strong. Thus, Fig. 8(c) shows a large wave amplitude in the undamped cylindrical meta-structure, which is close to three times the incident wave amplitude.

_{y}Curves of power, represented by the capture width ratio $ W / 2 a$, generated by the same three cylindrical meta-structures used in Fig. 8 but with constrained buoys at three different incident wave angles *β* are shown in Fig. 9. For the particular spring and damper parameters $ \gamma / ( \rho g h ) = 0.2$ and $ \sigma / \rho g = 0.2$, *d _{x}* has limited influence on the generated power when the incident wave angle is perpendicular to the deeper array of plates at low frequencies while for higher frequencies power is increased by reducing

*d*. As the incident wave angle

_{x}*β*increases, the incident wave direction is gradually aligned with the array of plates in the

*y*direction. As this happens, the value of

*d*has a greater influence, and more power is extracted at low frequencies as the value of

_{x}*d*is increased, although interestingly at higher frequencies results again suggest a lower value of

_{x}*d*is optimal. Note that when

_{x}*d*=

_{x}*d*, the wave power extraction is independent of the incident wave direction since the description of the structured cylinder under homogenization is axisymmetric.

_{y}The results of the shallow water theory in the range of $ v d y \u2208 [ 0 , 0.5 ]$ are also plotted in Fig. 9 for comparison. We can see that the wave power is little affected by *d _{x}*. When $ d x / h = 0$ and $ \beta = \pi / 2$, the results from the shallow water theory present a good prediction on wave power extraction since the incident wave angle is parallel to the plates such that the water depth becomes a less important factor.

Figure 10 shows the corresponding distribution of buoy vertical displacement based on the full depth-dependent theory at $ \beta = \pi / 4$ and $ \nu d y = 0.4$. As expected, amplitudes of motion on the waveward side are larger than those at the leeward side, but these remain at the same order as the incident wave amplitude. When *d _{x}* = 0, the contours are roughly aligned with the remaining single array of plates. As

*d*approaches

_{x}*d*, the contours become increasingly parallel to the direction of the incident wave crest.

_{y}Then, the curves of wave power extracted from two cylindrical meta-structures of $ a / h = 0.5$ with different plate submergence at $ \beta = \pi / 4$ are shown in Fig. 11. The settings of spring and damping are the same above, i.e., $ \gamma / ( \rho g h ) = 0.2$ and $ \sigma / \rho g = 0.2$. It shows that when the plates submerge to a deeper depth, more wave power is generated at low frequencies, but the situation is reversed at high frequencies. In addition, for the case of $ d h / h = 2 d x / h = 0.4$, the critical frequency occurs at $ \nu c h = 3.0$. Since the inner free surface is covered by the buoys, the resonance is suppressed by the damping mechanism.

Next, we consider the effect of damping rate *γ* and spring constant *σ* on the efficiency of wave power extraction. Figure 12 shows the capture width ratio $ W / ( 2 a )$ for a cylindrical meta-structure of $ a / h = 0.5 , \u2009 d x / h = 0.1$, and $ d y / h = 0.2$ at $ \beta = \pi / 4$ with different damping rates and spring constants. In Fig. 12(a), where $ \gamma / ( \rho g h ) = 0.4$ is fixed, the capture width ratio is shown to increase with decreasing spring constant, suggesting that fixing additional springs to the buoys may be unnecessary for the proposed scheme. In Fig. 12(b), where the spring constant is fixed at $ \sigma / \rho g = 0.2$, we see that different damping parameters work best at different frequencies.

Finally, we present the curves of wave power extraction, represented by the capture width *kW* in Fig. 13(a) and the capture width ratio $ W / 2 a$ in Fig. 13(b), against non-dimensional wavenumber *νh* for a cylindrical meta-structure with plate depths set at $ d x / h = 0.1$ and $ d y / h = 0.2$ with $ \gamma / ( \rho g h ) = 0.5$ and $ \sigma / \rho g = 0.0$ at $ \beta = \pi / 4$. For sufficiently low frequencies, the non-dimensional capture width, *kW*, and capture width ratio, $ W / ( 2 a )$, increase with the radius of the cylindrical meta-structure. Reference 12 was one of a number of papers who simultaneously proved a theoretical limit of *kW* = 1 for an axisymmetric device of any size absorbing in heave. The results therefore demonstrate that our proposed WEC can exceed this value across a wide range of frequencies: for example, $ k h \u2273 1.0$ when $ a / h = 1$ or $ k h \u2273 0.6$ when $ a / h = 2.0$. In Fig. 13(b), we have superimposed the results of Ref. 15 who considered a power absorption from a compact array of floating buoys operating in heave. Although the designs are clearly not the same, comparisons can be made since both consider the wave power absorption from a array of small buoys arranged in a circular array. To ensure the comparison is fair, the same damping ratio and spring constant have been chosen as in Ref. 15. The results illustrate that the present WEC appears to exploit the resonance promoted by the structured plate array, which is especially strong at low frequencies.

## VIII. CONCLUSION

In this paper, we have considered a cylindrical wave energy device consisting of two intersecting arrays of identical thin vertical plates that are immersed through the surface at right angles to one another and submerged to different depths. Floating buoys connected to springs and dampers are placed within each of the narrow vertical channels formed by the overlapping plate array and operate as wave energy absorbers. This WEC design falls into the minor “Many-body systems” classification of WEC according to the authoratitive review of Ref. 24. Although few concepts currently fall within this category the present work in conjunction with the contributions of, for example, Refs. 15 and 19 highlight the potential that such devices offer, which is far in excess of a cylinder of the same size operating in rigid body motion.

The device could be operated in shallow or deep water, and we have developed approximate homogenization methods to study its operation in both settings. The advantages of using shallow water theory are that numerical results are more robust and efficient to compute and the role of physical parameters, such as the depths of submergence of the plates, is easier to identify. Additionally, it is possible to use shallow water equations to explore the use of plate arrays of slowly varying depth within the cylinder, something we have not pursued here. The results have shown shallow water theory to be useful but limited when compared to the more accurate, complicated, and numerically challenging fully depth dependent homogenization model. This model has been validated in different ways including comparison with results from boundary element method computations made for an exact description of the plate arrays.

Numerical results for the operation of the cylinder as a WEC have been used to assess its potential by considering the influence of the key geometric and wave parameters on the operation of the device. The general conclusions are as follows: (i) the plate arrays are useful in enhancing the capture of energy from low frequency waves; (ii) tuned springs are not required for optimal power absorption; (iii) larger diameter cylinders both increase the capture width ratio and move its peak to lower wave frequencies; and (iv) there is some directional dependence to the device, but it appears to be relatively weak at low frequencies.

The modeling has been performed under the assumption of no hydrodynamical or mechanical losses, and it is not clear if viscosity and/or turbulence shed from the edge of the thin plates will have a significant effect on our predictions. Experimental testing of this device is planned and will allow us to assess these issues.

## ACKNOWLEDGMENTS

J.H. and R.P. gratefully acknowledge the EPSRC for supporting this work through Grant No. EP/V04740X/1. S.Z. gratefully acknowledges the State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University for supporting part of this work through the Open Research Fund Program (Grant No. HESS-1902).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jin Huang:** Data curation (equal); Formal analysis (equal); Methodology (equal); Validation (equal); Writing – original draft (equal). **Richard Porter:** Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Siming Zheng:** Methodology (equal); Validation (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: DERIVATION OF GOVERNING EQUATIONS AND BOUNDARY CONDITIONS UNDER FULL DEPTH-DEPENDENT THEORY

*X*,

*Y*) are microscopic field variables based on a single cell such that $ 0 \u2264 X \u2264 1$ and $ 0 \u2264 Y \u2264 1$. Simultaneously, we shall define $ \epsilon = L / \lambda $ and expand variables as

*O*(1). From Eq. (A3), we have

*X*and

*Y*[i.e., $ \varphi ( 0 ) = \varphi ( 0 ) ( x \u2032 , y \u2032 , z \u2032 )$] due to the fact that

**N**is the unit normal in terms of (

*X*,

*Y*) variables out of Ω.

*y*direction, we rescale with

*Y*dependence is dropped from the corresponding expansion in Eq. (A6) and we can follow the same process to give the leading order governing equation

For the leading-order governing equation in $ \Omega 23$ or outer region $ \Omega 1$, we have no need to introduce a microscale, so the leading order governing equation remains as Eq. (A1) and the bottom condition Eq. (A4) also applies without change.

A careful analysis of the two intermediate region close to $ z = \u2212 d x$ and $ z = \u2212 d y$ on the interface between $ \Omega 21$ and $ \Omega 22$ and between $ \Omega 22$ and $ \Omega 23$ which involves rescaling *z* on the lengthscale *L* via $ z = \u2212 d x , y + L Z$ reveals that, at leading order, $ \varphi z ( 0 )$ and $ \varphi ( 0 )$ are continuous across $ z = \u2212 d x , y$. Readers are directed to Ref. 5 for details of the application of this process in a similar problem.

### APPENDIX B: NUMERICAL PROCEDURE FOR SOLVING DISPERSION EQUATION (26)

In the numerical process of finding the roots of Eq. (26), we first solve the equation with $ Im [ \tau ] = 0$, and all roots now locate on the real or imaginary axis. For a general value of *t*, there are two sequences of discrete roots satisfying Eq. (26). That is, if $ \kappa q ( t )$ is a root of Eq. (26), then so is $ \u2212 \kappa q ( t )$. Since *t* is integrated over all angles, $ \u2212 \pi \u2264 t < \pi $, the real positive root $ \kappa 0 ( t )$ representing the propagating wave and an infinite number of pure imaginary roots $ \kappa q ( t ) = i \kappa \u0302 q ( t )$ on the positive imaginary axis representing the evanescent waves need to be considered. Furthermore, it can be observed from Eq. (26) that $ \kappa q ( \u2212 t ) = \kappa q ( t )$ and $ \kappa q ( \pi \u2212 t ) = \kappa q ( t )$ which confine us to find the roots $ \kappa q ( t )$ only in the range of $ t \u2208 [ 0 , \pi / 2 ]$ and helps reduce the numerical effort.

From Eq. (27), we can see that the real root $ \kappa 0 ( t )$ will tend to infinity as $ \nu d y \u2192 1 + \u211c [ \tau ]$ and only exists for $ 0 \u2264 t \u2264 \pi / 2$ when $ \nu d y < 1 + \u211c [ \tau ]$, which means that waves do not propagate within the inner region when $ \nu > \nu c$ and $ \nu c d y = 1 + \u211c [ \tau ]$ defined as the critical frequency. This issue has been the subject of a separate paper by the authors (see Ref. 23).

After determining the positions of roots for each $ t \u2208 [ 0 , \pi / 2 ]$ when $ Im [ \tau ] = 0$, we consider *τ* to be a complex number: the buoys are now resisted by the damper. As soon as *τ* takes on an imaginary component, all the roots $ \kappa q ( t )$ move off the real or imaginary axis into the complex plane. We can take the root at $ Im [ \tau ] = 0$ as the initial value [denoted as $ \kappa q 0 ( t )$], and obtain the final result by taking the imaginary part of *τ* into account. For efficient calculation, a self-adaptive method is adopted as follows: (i) We take $ \kappa q 0 ( t )$ as an initial value and obtain a new root $ \kappa q 1 ( t )$ after iteration; (ii) if $ | ( \kappa q 1 ( t ) \u2212 \kappa q 0 ( t ) ) / \kappa q 0 ( t ) |$ is less than a threshold $ \Delta \kappa $ [which can guarantee that for a certain *q* the root does not jump to other branches (see Fig. 15)], $ \kappa q 1 ( t )$ is the root of Eq. (26) and the calculation stops—otherwise, we solve Eq. (26) with the imaginary part of *τ* taking the value on the midpoint of the interval $ [ 0 , Im [ \tau ] ]$ to get $ \kappa q 2 ( t )$; and (iii) if we still have $ | ( \kappa q 2 ( t ) \u2212 \kappa q 0 ( t ) ) / \kappa q 0 ( t ) | > \Delta \kappa $, the interval is halved again until the termination condition is satisfied—otherwise, we will take $ \kappa q 2 ( t )$ as an initial value and repeat procedure (i) on the interval $ [ Im [ \tau ] / 2 , Im [ \tau ] ]$.

It should be noted that when frequency approaches the critical frequency *ν _{c}*, under certain conditions (for example,

*t*= 0), much computational effort is required to determine $ \kappa 0 ( t )$ since the real initial value for the

*q*= 0 root tends to infinity. In addition, since the real initial root no longer exists when $ v d y > 1 + \u211c [ \tau ]$, we do not have the corresponding initial value to find $ \kappa 0 ( t )$ in the complex plane. In order to overcome these two issues, the roots of Eq. (26) at the relatively low frequency can be taken as the initial values, and obtain the final result also by a similar self-adaptive method with increasing the frequency. Generally, the above numerical method is robust.

In order to illustrate the above procedure, a case is examined. First, *τ* is taken to zero, and Fig. 14 presents the variation in the real root $ \kappa 0 ( t )$ for $ t = \pi / 4$ and $ t = \pi / 2$ against frequency *ν* with a fixed value of $ d y / h = 0.2$ and three different values of $ d x / h$. It should be noted that since Eq. (26) will be reduced to Eq. (27) when *t* = 0 or *d _{x}* =

*d*, all the curves of $ \kappa 0 ( 0 )$ for three cases are the same as the black dot dash line in Fig. 14(a) or 14(b). For a certain

_{y}*t*, the real root, $ \kappa 0 ( t )$, increases with the frequency until $ \nu d x = 1$ since the real root exists only when $ \nu d x < 1$ as mentioned above. In addition, when

*d*tends to

_{x}*d*, $ \kappa 0 ( t )$ quickly become a large value at high frequencies. It indicates that a very small wavelength appears and results in a large resonant amplification in the cylindrical meta-structure, which may violate the underlying assumption that the plate spacing is much smaller than the wavelength.

_{y}Furthermore, if *τ* has an imaginary part, $ \kappa q ( t )$ will go into the complex plane. In Fig. 15, we let *σ* = 0 and $ \gamma / ( \rho g h ) = 0.2$ and present the locations of the first four root $ \kappa q ( \pi / 4 ) ( q = 0 , 1 , 2 , 3 )$ at each step as $ Im [ \tau ]$ increases from zero. $ \kappa q ( \pi / 4 )$ starts from the real or imaginary axis and moves to a certain point in the complex plane. It also can be seen that since a self-adaptive method is performed, the final result can be determined after only a few steps which set the stage for efficient computation within the whole code.

### APPENDIX C: DERIVATION OF GOVERNING EQUATION FOR THE SHALLOW WATER THEORY

In this appendix, we also apply the homogenization method to consider the interaction of water waves with the cylindrical meta-structure under a shallow water approximation. In addition to the assumptions given in Sec. II, we should further assume that the water depth *h* is small compared to the wavelength *λ*, i.e., $ \mu = h / \lambda \u226a 1$. Additionally, the spacing of plates should be sufficiently small compared to the water depth, i.e., $ L / h \u226a 1$, and we specifically require $ \epsilon = O ( \mu 2 )$.

Although we treat various quantities including the depths *h*, *d _{x}*, and

*d*as fixed in the derivation below, a more general derivation (e.g., see Ref. 4) allows these (in addition to mechanical parameters $M$,

_{y}*σ*, and

*γ*) to be functions of the macroscopic variables

*x*,

*y*. The shallow water theory is thus more versatile than the full-depth theory, and for this reason, we have chosen to retain explicit time dependence, rather than factorizing a single frequency component as we did in Appendix A.

*A*is a wave amplitude that has been assumed to be sufficiently small compared to other lengthscales to allow the neglect of non-linear terms in the governing equations and boundary conditions presented in Sec. II. We have also written $ u = ( u h , w )$ such that $ u h = ( u , v )$. We note that since $ M = \rho d$ and

*d*<

*h*, then $ M \u2032 = O ( \mu 2 ) = O ( \u03f5 )$, and we therefore write $ M \u2032 = \u03f5 M \u2033$.

*X*= 0, 1, $ 0 < Y < 1$, and on

*Y*= 0, 1 for $ 0 < X < 1$, where

**N**is the unit normal in (

*X*,

*Y*) variables on those four sides. Additionally, we use Eq. (7) which, at leading order, gives the conditions

*y*-axis, we can apply the same treatment to the governing equations and boundary conditions, after replacing the scaling in Eq. (C2) by

*Y*from the dependent variables. We can deduce that $ u ( 0 )$ is independent of the microscopic coordinate

*X*. In addition, we also have

*z*and so from Eq. (C17), we can see $ v ( 0 ) = v ( 0 ) ( x , y , t )$. Thus, after integrating Eq. (C16) across $ 0 < X < 1$ and over $ \u2212 d y / h < z < \u2212 d x / h$, we have

*X*as a consequence of the leading order contribution from the irrotationality condition Eq. (7) under the multiple scales expansion.

The governing equation in $ \Omega 1$ can be obtained by letting $ \gamma = \sigma = d x = d y = 0$ in Eq. (C24), which reduces to the standard shallow water equation over constant depth *h*. Despite setting *h*, *γ*, *σ*, *d _{x}*, and

*d*to constant values in the present study, the same governing equation can be shown to hold when these are allowed to vary on the scale of the wavelength.

_{y}## REFERENCES

*Handbook of Mathematical Functions*