Wave scattering by an array of metamaterial cylinders

Abstract In this paper, a semi-analytical model based on linear potential flow theory and an eigenfunction expansion method is developed to study wave scattering by an array of structured cylinders in water of finite depth. Each cylinder is formed by a closely spaced array of thin vertical plates, between which fluid can flow, extending through the depth. In order to consider the wave attenuation and energy dissipation in narrow gaps between the thin vertical plates, a damping mechanism is introduced at the surface of the fluid occupied by the structured cylinders. In addition to a direct calculation of the energy dissipation, an indirect method based on Kochin functions is derived with the employment of energy identities. The present model is shown to be in excellent agreement with both the published data and those obtained by using different methods. The validated model is then applied to study the effect of a pair of structured cylinders on wave focusing/blocking, scattered far-field amplitude and wave power dissipation. Results show that wave focusing/blocking can be achieved by the appropriate choice of plate alignment. The structured cylinders hold profound potential for wave power extraction.


Introduction
The interaction of water waves with impermeable vertical cylinders extending through the surface of a fluid has been an active area of study over many decades. This is partly because of its practical relevance in relation to marine structures such as the supporting columns of wind turbines, oil rigs, bridges and so on. In conjunction, the boundary-value problem that results from the mathematical description of the water-wave problem is amenable to analytic methods with particular advantage being taken of the alignment of fluid boundaries with coordinate surfaces in cylindrical polar coordinates. Consequently, it also acts as a prototype problem for many computational and experimental methods.
Under the small-amplitude (linearised) description of water waves, the scattering of incident plane waves by a single rigid vertical cylinder extending upwards through the surface from the bed of a fluid of constant depth is explicit; see MacCamy & Fuchs (1954). For cylinders extending uniformly through the depth, the dependence upon vertical coordinate is separable and, consequently, the problem is governed by the two-dimensional Helmholtz equation with the implication that the solutions have interpretations in other physical settings such as two-dimensional linearised acoustics and transverse-magnetic-polarised electromagnetics. The extension to multiple cylinders has been the subject of a number of papers (e.g. Siddorn & Eatock Taylor 2008;Zheng, Zhang & Iglesias 2018) and the theory for a finite number of arbitrarily placed cylinders in a water wave setting is described by Linton & Evans (1990), who followed and extended the original method of solution devised by Záviška (1913) and later Spring & Monkmeyer (1974) to show that the forces on the cylinders could be expressed in a particularly simple way in terms of the solution of certain infinite systems of equations. A number of interesting effects occur when waves interact with multiple vertical cylinders. For example, when then vertical axes of N ≥ 4 cylinders are equally spaced in a circular arrangement, Evans & Porter (1997) showed that large amplifications of the incident waves could occur inside the ring of cylinders. This so-called near-trapping phenomenon becomes especially dramatic as the gaps between cylinders become much smaller than the cylinder diameter, resulting in large peaks in wave forces close to certain frequencies linked to near resonance.
For long arrays of cylinders with vertical axes equally spaced along a straight line (a truncated periodic array) Maniar & Newman (1997) also discovered a near-trapping phenomenon with similar consequences on surface elevation and cylinder wave force amplification. This time, the connection was made to free oscillations that were shown to occur in the equivalent infinite periodic array (e.g. Linton & Evans 1993;Porter & Evans 1999;Thompson, Linton & Porter 2008). Notable extensions to problems involving non-circular cylinders, truncated cylinders, second-order theory and ice covered surfaces are described by Chatjigeorgiou (2011), Zheng et al. (2020b), Wolgamot, Eatock Taylor & Taylor (2015), Malenica, Eatock Taylor & Huang (1999) and Ren, Wu & Ji (2018). It is worth remarking here that if the cylinder is not uniform in the depth (i.e. it is truncated), the boundary-value problem becomes more complicated as the separable depth dependence can no longer be assumed and solutions are complicated by the need to expand over an infinite set of depth eigenfunctions (e.g. Zheng, Zhang & Iglesias 2019a;Zheng et al. 2019b).
In this paper the focus is on vertical cylinders which are no longer rigid and impermeable, but which are structured in such a way that fluid is allowed to flow inside the cylinder. The particular structure of the cylinder we choose to consider is comprised of closely spaced thin parallel array of vertical plates whose lateral edges form the outline of a cylinder when viewed from above. Thus the fluid (and waves on the surface of the fluid) can move between the plates in the direction of the plates, but there is limited motion perpendicular to the plates owing to the assumed narrowness of the gaps between adjacent plates forming the structure. The idea for the use of such a structure in the water wave context originates from Porter (2018), who used the same parallel-plate-array 'metamaterial' occupying an infinitely long rectangular domain. The terminology metamaterial is used to describe a medium which exhibits behaviour not associated with normal materials; this is manifested by some complicated form of anisotropy. A metamaterial obtains its properties from a microstructure whose length scale is much smaller than that of the underlying field variables. Properties may be obtained by direct simulation, or by deriving an effective equation governing the microstructured medium via homogenisation or a multiscale method. See for example, Mei & Vernescu (2010), Berraquero et al. (2013) and Maurel et al. (2017) for general theory and its application to structured bathymetry in water waves as devices for producing anisotropic effects in surface wave propagation. Porter (2018) shows that the plate-array metamaterial is governed by a reduced wave equation allowing waves to travel only in the direction aligned with the plate array. This anisotropy of wave propagation manifested itself both as an all-frequency perfectly transmitting negative refraction device for a particular incident wave angle or as a perfectly transmitting all-angle negative refraction device for particular frequencies. Also in Porter (2018) an outline description of the use of the metamaterial plate-array vertical cylinder was produced with some preliminary results. For example, for wave headings aligned with the plate array, the cylinder is transparent to incident waves whilst for other wave headings and frequencies the interaction is more complicated producing, in general, unsymmetric wave diffraction.
This paper develops the preliminary study of Porter (2018) on cylinders and extends the theory in two directions. First, we consider multiple cylinders and the interaction between them. Secondly, recognising that the assumed narrow fluid channels between the closely spaced plates may lead to viscous damping, we include in our model an artificial linear damping mechanism added to the free-surface dynamics. Although not directly related to the physical source of viscous damping on the vertical plate structures, artificial damping in the free-surface condition is a commonly used device whose dependence upon physical parameters, we imagine, can be parametrised via computational fluid dynamics or experimental methods. Whilst standard analytic tools (see above) can be used to consider the interaction between multiple cylinders, the addition of damping to the surface condition inside the metamaterial cylinder means we no longer enjoy a separable depth dependence and are required to deploy a full expansion in depth eigenfunctions for the velocity potential. The application of effective boundary conditions matching the flow in the exterior of the cylinders to the uni-directional flow inside the structured cylinders leads to infinite systems of equations to be solved. This procedure follows the problem statement outlined in § 2 of the paper where subsequently we derive (with algebraic details relegated to appendix B) expressions for the rate of energy dissipation due to damping and the far-field diffraction coefficient. In § 3 we describe the numerical convergence characteristics and validate the model. In § 4 we produce a set of results mainly focusing on the interaction effects between two cylinders. Particular attention is given to wave focusing, wave sheltering and energy dissipation characteristics. We draw conclusions to the work in § 5 systematically.

Mathematical model
In this model, a number (N) of metamaterial circular cylinders conceptually deployed as an array in water of finite depth h are considered (see figure 1). A global Cartesian coordinate system Ox yz is chosen with the mean free surface coinciding with the (x, y)-plane and z measured vertically upwards. Hence the fluid bottom is at z = −h. Cylinder n with its radius denoted as R n is composed of a periodic array of infinitely thin vertical plates rotated through a clockwise angle β n relative to the Ox axis; (x n , y n , 0) denotes the horizontal position of cylinder n in the coordinate system Ox yz. Plane waves propagating at an angle β relative to the Ox axis are incident on these metamaterial cylinders. Fluid is allowed to flow in gaps between adjacent plates and waves are supported by the free surface. In addition to the global Cartesian coordinate system, the local one O n x n y n is also adopted with O n x n in parallel with the plates. The effect of these plates allows waves to propagate in the ±O n x n -direction only. Moreover, N cylindrical coordinate systems, O n r n θ n z, are chosen for the purpose of convenience of mathematical expression. Additionally, one more cylindrical coordinate system Or 0 θ 0 z is defined (not plotted in figure 1), the origin of which coincides with the coordinate system Ox yz. R n,j and α n,j denote the length and the angle, respectively, of a vector pointing from O n to O j . We assume that all amplitudes are small enough that linear theory applies and we make the usual assumptions that the fluid is inviscid, incompressible and its motion is irrotational. We denote the fluid velocity potential by Φ (x, y, z, t). It is further assumed that all motion is time harmonic with angular frequency ω. Thus, we can write where Re denotes the real part. Thus φ is the spatial velocity potential which is independent of time, i.e. t. i is the imaginary unit. The fluid domain can be divided into N interior domains, which fill the N cylinders accordingly, and an exterior domain, representing the remainder of fluid domain extending towards infinity horizontally.
The spatial velocity potential satisfies Laplace equation, 2 φ = 0 in the water, (2.2) the boundary condition at sea bed, and the boundary condition at the water surface of the exterior domain in which g denotes the acceleration due to gravity. Within the fluid in the nth cylinder, (2.2) also holds although it is confined to narrow disconnected domains bounded by thin plates aligned with the x n coordinate. Writing (2.2) in coordinates O n x n y n with rescaled in x n and y n coordinates and imposing the boundary conditions on the channel walls shows that the field within the whole of the nth cylinder is governed by an effective medium governing equation involving the reduced Laplacian It is assumed that the separation between plates is small compared to the wavelength (i.e. d p /λ 1 where d p is the distance between plates and λ is the wavelength) and also the length of the plate (i.e. d p /L p 1 where L p is the length of the plate). See also Porter (2018) and Jan & Porter (2018) who employed the same models. Equation (2.5) represents conservation of mass for an irrotational flow in which the motion perpendicular to the plates is inhibited.
Within the boundary of the cylinder and between the plates we allow for the possibility of energy dissipation and will employ the modified free-surface condition , z = 0, (2.6) (sometimes referred to as a 'damping lid' model, e.g. Kim, Koo & Hong 2014;Dinoi 2016) withν ≥ 0 within the cylinder as a means of achieving this. We identify three physical settings in which this condition applies. The first is that the surface of the fluid within the cylinder is covered with a fixed porous medium with permeability κ submerged to a small depth d. The flow through small vertical pores is assumed to be dominated by the fluid dynamic viscosity, μ, and it is appropriate to use Darcy's law (e.g. Chwang & Chan 1998) to relate the vertical fluid velocity w to the pressure gradient p z via w = −(κ/μ)( p z + ρg), where ρ represents the water density. Integrating subject to the kinematic and dynamic free-surface conditions and matching the pressure and the mass flux to an inviscid fluid described by potential flow theory beneath the porous medium readily leads to the free-surface condition where 0 < α ≤ 1 is a 'blockage coefficient' representing the fractional area of the medium occupied by pores in horizontal cross-section. The second physical setting involves the surface of the narrow channels within the cylinder being covered by floating buoys constrained to move in heave. The buoys are designed to operate as wave energy converters being connected to a power take-off mechanism with a linear damping rate c. Garnaud & Mei (2009) showed, using multiscale homogenisation theory underpinned by an assumed contrast in wavelength and buoy separation, that the effect of a compact array of buoys occupying a fraction γ ∈ (0, 1] of the area of the surface can be represented by the modified free-surface condition which coincides with (2.6) when γ = 1. Garnaud & Mei (2009) give an example of the application of this condition to an array of buoys along a rectangular channel.
The final setting arises from consideration of the viscous dissipation due to fluid interaction with the sidewalls and bottom of the narrow rectangular fluid-filled channels with a normal air-fluid free surface. Hunt (1952) and Mei, Stiassnie & Yue (2005) ( § 9, Exercise 9.2) have shown that the effect of dynamic viscosity, μ, on a plane wave of angular frequency ω propagating along a uniform channel of width d p and depth h is to shift the inviscid wavenumber from k to provided k is not close to zero and is small. The condition (2.6) can be used to generate the same effect since, ifν is small, and the velocity potential of a propagating wave along the channel is sought to in the form and allowing a connection to be made betweenν and in (2.12) and (2.9) above. The range of values ofν that we shall consider in later results may not be appropriate to all physical settings but are included to demonstrate the full range of wave interaction available under the condition (2.6).

Expressions of spatial velocity potential in different domains
The standard method of eigenfunction expansions is used to solve the wave-structure interaction problem (e.g. Mei 1983).

Exterior domain
The spatial velocity potential in the exterior domain can be expressed as (e.g. Siddorn & Eatock Taylor 2008; where φ I represents the velocity potential of incident waves. The second term denotes the components contributed by the waves scattered from the N cylinders; A (n) m,l are the unknown coefficients to be determined; H m denotes the Hankel function of the first kind of order m; Z l (z) = cosh[k l (z + h)]/cosh(k l h); k 0 ∈ R + and k l ∈ iR + for l = 1, 2, 3, . . . are associated with propagating waves and evanescent waves, respectively, and they are the positive real root and the infinite positive imaginary roots of the dispersion relation for the exterior domain ω 2 = gk l tanh(k l h). (2.14) For the plane incident waves with amplitude A, angular frequency ω and wave direction β, φ I can be expressed in the coordinate systems of Ox yz and O n r n θ n z, respectively, as where J m is the Bessel function of order m (e.g. Linton & Evans 1990;. After using Graf's addition theorem for Bessel functions, φ ext can be rewritten in the cylindrical coordinate system O n r n θ n as (2.17)

Interior domain
General solutions of the reduced Laplace's equation, i.e. (2.5), inside the nth cylinder satisfying free-surface and bed boundary conditions, i.e. (2.6) and (2.4), can be expressed as . are the complex roots of the dispersion relation for the interior domains which degenerates into the dispersion relation for the exterior domain, i.e. (2.14), whenν = 0. The values of k l for l = 0, 1, 2, 3, . . . can be calculated using an analytic continuation method, starting with the corresponding roots for the case ofν = 0 (i.e. k l for l = 0, 1, 2, 3, . . .), and incrementingν to the specified value (e.g. Meylan, Bennetts & Peter 2017;Zheng et al. 2020a).
In (2.18), B n,l and C n,l are coefficients expressing the amplitude of waves propagating in each direction within the channels as a function of y n ; E n,l and F n,l are the same coefficients as a function of the angle (θ n ) at which the channel emerges at the edge of the cylinder (i.e. the angle of − − → O n P as shown in figure 2). When r n = R n , we have θ n = θ n , and which express the fact that the channels between the plates connect the cylindrical surface θ n ∈ [β n − π/2, β n + π/2] and θ n ∈ [β n + π/2, β n + 3π/2].
With consideration of the properties as given in (2.20a,b), and for the purposes of deriving a solution, we expand the functions E n,l and F n,l as where E (n) p,l and F (n) p,l are the unknown coefficients to be determined. The expression of φ (n) int as given in (2.18) can be further rewritten with the employment of (2.23)

Solution of unknown coefficients
Continuity of the field in terms of pressure and flux across the interfaces of the interior and exterior domains requires which can be used to determine the unknown coefficients A (n) m,l , E (n) p,l and F (n) p,l . The latter condition is derived from a flux balance through a small right-angled triangle with sides approximating the circular boundary of the cylinder, the perpendicular line across the entrance to a narrow channel and a channel sidewall. Detail derivation and calculation of the unknown coefficients are given in appendix A.

Wave motion
The water elevation non-dimensionalised by the incident wave amplitude can be expressed asη in which the term −νi will vanish for the exterior domain.

Far-field scattering amplitudes
In the water domain far away from an array of metamaterial cylinders, only the propagating modes exist in the scattered waves. With the asymptotic forms of H m for r 0 → ∞, H m (kr 0 ) = 2/π exp(−i(mπ/2 + π/4))(kr 0 ) −1/2 exp(ikr 0 ) for r 0 → ∞, (2.27) the scattered wave potential, i.e. the accumulative term in (2.13), can be rewritten as which can be further expressed in the global cylindrical coordinate system O 0 r 0 θ 0 z as where A S is the so-called far-field scattering amplitude that is independent of r 0 and z, and can be expressed as The energy dissipated by the N metamaterial cylinders due to damping coefficient can be calculated by (Zheng et al. 2020a) where Ω n denotes the water surface of the interior domain occupied by cylinder n, and η denotes the time-independent surface elevation. Equation (2.3) presents a straightforward way to calculate the energy dissipation by the array of metamaterial cylinders. From the view of energy identities, the energy dissipation can also be evaluated based on the spatial potentials in the exterior domain where Ω R represents an envisaged vertical cylindrical control surface with its radius denoted by r 0 = R 0 , which is large enough to enclose all the cylinders. The derivation process of (2.32) can be found in appendix B; when r 0 = R 0 → ∞, (2.32) holds as well with the control surface Ω R replaced by Ω ∞ , i.e. r 0 → ∞. It has been shown that the integral in (2.32) can be expressed in terms of Kochin functions (Falnes 2002), where H R is the Kochin function which can be expressed as follows (Falnes 2002) Therefore, the energy dissipated by the array of metamaterial cylinders can be evaluated by using an indirect method based on Kochin functions which presents a way to check the accuracy of the proposed semi-analytical model.
The energy dissipated by the cylinders can be written in non-dimensional format as where P in is the incident wave power per unit width of wave front, . (2.38)

Model validation
The effect of truncation of the infinite sums on the angular and vertical modes to finite sums over −M ≤ m ≤ M and 0 ≤ l ≤ L have been carried out and suggest that M ≥ 20 and L ≥ 5 provide sufficiently converged results for kh = 1.3. As kh becomes larger, more truncated terms of m and l may be required to obtain the converged results. Hereinafter, M = 20 and L = 5 are adopted unless otherwise specified. Forν = 0 with β n = β, i.e. when waves propagate into the cylinders with the plates aligned to the incident wave direction, the incident waves would not be affected at all. figure 3 presents the predicted wave field around a pair of metamaterial cylinders with When the metamaterial cylinders are deployed far away from each other, the wave motion at each cylinder is expected to be the same as that for an isolated single metamaterial cylinder. Figure 4 illustrates the instantaneous wave field around one of a pair of metamaterial cylinders far apart from one another. The present results are found to agree well with those for a single metamaterial cylinder in the absence of damping as investigated by Porter (2018). Due to the effect of the closely spaced array of thin vertical plates aligned with the y-axis, the scattering pattern outside the cylinder has broken the symmetry of the incident wave and is very different from that for a solid cylinder. Note that in figures 4(a) and 4(b), the cylinder redirects a 'beam' of energy rightward (i.e. in a direction perpendicular to the plate direction), which may be called a 'beaming' or FIGURE 5. Instantaneous wave field due to incident wave propagation with kh = 1.3, β = π/4 on a pair of metamaterial cylinders with R 1 /h = R 2 /h = 1.0, −x 1 /h = x 2 /h = 2.0, y 1 = y 2 = 0, β 1 = β 2 = π/2: (a) cylinders with a large damping coefficient,ν = 10 5 ; (b) cylinders with fixed solid lid at the mean water surface.
'lensing' effect, and it will be further investigated in § 4. Whilst the motion of the fluid in each channel appears to be separate from the next, there is coupling between channels on the boundary of the cylinder and this coupling seems to give rise to a slow wave with high energy propagating through the cylinder. Also, it is observed that the plates have the effect of inducing resonant-like behaviour in the fluid channels. This is particularly noticeable in figure 4(c) with kh = 1.6, where there is a large resonant amplification in the channels whose lengths are approximately half a wavelength.
Another extreme case is that whenν → ∞, the wave motion on the surface of the internal region is strictly restricted, and the wave scattering problem becomes the same for the metamaterial cylinder with a fixed solid lid at the mean water surface, the wave scattering solution of which is derived in appendix C. Comparison of the wave field for a pair of metamaterial cylinders withν = 10 5 and that of metamaterial cylinders with fixed solid lids at the mean water level is plotted in figure 5. Note that a metamaterial cylinder with a rigid lid is not the same as a vertical cylinder with a rigid cylindrical surface. In the former case, fluid is still able to flow through the cylinder. Additionally, the wave power dissipated by the metamaterial cylinders evaluated by using the direct method ((2.3)) and the indirect method ((2.36)) are presented in figure 6.
Moreover, potential flow theory based numerical simulations are carried out with the employment of a commercial boundary element method (BEM) code AQWA (ANSYS, Inc. 2011) to study wave interaction with a metamaterial circular cylinder consisting of 20 thin vertical plates (figure 7). The thickness of each plate is 0.02h and the spacing distance between the centres of adjacent plates is 0.1h. The good agreement between the semi-analytical results (figure 4a) with BEM numerical simulations (figure 7b) confirms that the homogenisation of the structured cylinder into an effective medium with effective boundary conditions is a good approximation. An obvious advantage of the semi-analytical model lies in its high computational efficiency, and, indeed, our results are much easier to compute compared to BEM numerical computations. The excellent agreement between the results shown in figures 4-7, together with figure 3, gives confidence in the present model for solving wave scattering and predicting wave dissipation by an array of circular metamaterial cylinders.

Results and discussion
In this section, the effect of the metamaterial cylinders on wave focusing/blocking and scattered far-field amplitude is investigated with the employment of the validated semi-analytical model. Additionally, wave power dissipation of the cylinders is studied, and shows to form the foundation of a wave energy device with a high 'capture width' if the artificial surface damping used in our present model were to be replaced by a mechanical energy conversion device with similar effects.
Prior to investigating performance of a pair of metamaterial cylinders, the angle responses of the scattered far-field amplitude for a single metamaterial cylinder placed at x = y = 0 with β 1 = 0, π/6, π/4, π/3, π/2 are plotted in figure 8. For the non-damping situation as shown in figure 8(a), the main peak value of the far-field scattering wave amplitude and the corresponding angle are (|A S |/A, θ 0 ) = (1.75, 0.50π), (1.21, 0.67π), (0.73, 0.76π) and (0.34, 0.85π) for β 1 = 0, π/6, π/4 and π/3, respectively, in which (θ 0 − β 1 ) ≈ 0.5π is satisfied, and moreover, the |A S |/A is vanishing at θ 0 ± 0.5π approximately. This means the cylinder bends or redirects a 'beam' of energy in a direction perpendicular to the plate direction, though there is a loss in the intensity of this beam as the angle is rotated with respect to the incident wave angle. We note that there is very little lateral scattering of wave energy either laterally or back towards the incoming wave direction. That is, the metamaterial cylinder acts rather like a transparent lens, but also one which appears to absorb wave energy laterally into the microstructure and produce and intense forward beam. The same thing still roughly happens forν = 0.1 (figure 8b). However, when the damping is too large that it works like a solid lid placed on the surface of the structured cylinder (figure 8c), the angle response of the scattered far-field amplitude is lightly dependent on β 1 , indicating that the orientation of the plates is relatively unimportant as far as the overall effect of the cylinder is on wave diffraction.

Wave focusing/blocking
When β 1 = π/2, the waves pass through the cylinders with no scattering (figures 9d and 9h). For other values of β 1 the metamaterial cylinders interacts with incident waves in a non-trivial way. Compared to the waveward surface elevation, wave motion at the leeward region and close to the cylinders is more affected by the two cylinders. When the thin plates are all aligned along the incident wave crest line (β 1 = 0, figures 9a and 9e), wave motion is suppressed at the region between the two cylinders, where two small areas ofη < 0.4 are observed. Moreover, two larger wave attenuation areas ofη < 0.4 can be found on the flanks of the pair of cylinders on the leeward side. On the other hand, wave motion is strengthened at the central leeward region, which extends to the two cylinders, forming an inverted 'Y' shape area ofη > 1.2. For the case with β 1 = −π/6 (figures 9b and 9f ), there is a wave focusing area at the central leeward region of the array as well, while the region is closer to the array, and the wave in the region is much more focused withη > 2.0. The largest wave amplitude in the computed range of the exterior region isη = 2.31, which occurs at (x/h, y/h) = (0, 1.44). Meanwhile, there is a small narrow region ofη < 0.4 observed immediately beyond each cylinder, where the smallest wave motion isη = 0.02 at (x/h, y/h) = (±1.86, 1.20). As a comparison, for the metamaterial cylinders with β 1 = π/6 as shown in figures 9(c) and 9(g), there is a much larger area ofη < 0.4 at the very leeward of the array, where waves are effectively blocked by the cylinders. At (x/h, y/h) = (±3.34, 4.86),η = 0 is obtained, meaning the incident wave can be completely blocked at specified points. The dramatic amplification and focusing effects on wave motion inside the metamaterial cylinders are observed for all the studied -4 -3 -2 -1 0 1 2 3 4 5 -5 -5 x/h -4 -3 -2 -1 0 1 2 3 4 5 -5 -5 x/h 5 4 3 2 1 0 -1 -2 -3 -4 -4 -3 -2 -1 0 1 2 3 4 5 -5 -5 x/h cases, except where β 1 = π/2. The results as given in figure 9 demonstrates that wave focusing/blocking can be achieved by a pair of metamaterial cylinders with the appropriate control to the plates alignment direction. The near-field wave motion due to the same incident waves propagating on the same metamaterial cylinders withν = 0.1 is presented in figure 10. Due to wave power dissipation of the metamaterial cylinders, the wave focusing area (η > 1.2) at the leeward region of the cylinders with β 1 = 0 (figures 9a and 9e) now mostly becomes a blocking region withη < 0.8 (figures 10a and 10e). What is more, the previous regions ofη < 0.4 now merge together, resulting in a much larger 'M' shaped region. For the case with β 1 = −π/6 (figures 10b and 10f ), asν increases from 0 to 0.1, the wave focusing region ofη > 1.2 previously located at the central leeward of the array now moves to the gap between the cylinders, and gets smaller. Whereas the wave blocking regions ofη < 0.4 grow and, as a result, they merge together into an inverted 'V' shape area. With the increase ofν from 0 to 0.1, the wave blocking region ofη < 0.4 for β 1 = π/6 breaks into two regions, and the correspondingη < 0.8 region becomes broader (figures 10c and 10g). The previous amplification and focusing effects on wave motion inside the metamaterial cylinders are now significantly weakened by the damping, except the one with β 1 = π/2. Due to the existence of damping, the incident wave is disturbed by the metamaterial cylinders with β 1 = π/2, despite very limited influence (figures 10d and 10h). Figure 11 illustrates the near-field wave motion when an extremely large dampingν = 10 5 is employed, which is equivalent to a solid lid placed on the surface of each cylinder. When a solid lid is put on the surface, it largely produces the same overall wave pattern. That is, the orientation the plates are relatively unimportant as far as the overall effect of the cylinder is on wave diffraction, which is in accordance with that obtained for the isolated cylinder (see figure 8c). The regions ofη > 1.2 are distributed at the left and right sides of the array, and in the gap between the two cylinders. Additionally, due to -4 -3 -2 -1 0 1 2 3 4 5 -5 -5 x/h -4 -3 -2 -1 0 1 2 3 4 5 -5 -5 x/h -4 -3 -2 -1 0 1 2 3 4 5 -5 -5 x/h FIGURE 12. Far-field scattering wave amplitude due to incident wave propagation with kh = 1.3, β = π/2 on a pair of metamaterial cylinders with R 1 /h = R 2 /h = 1.0, −x 1 /h = x 2 /h = 2.0, y 1 = y 2 = 0, β 2 = −β 1 : (a)ν = 0; (b)ν = 0.1; (c)ν = 10 5 .

Scattered far-field amplitude
occurring at θ 0 = 0.5π. In the range of θ 0 ∈ [π, 2.0π], |A S |/A is small regardless of the value of β 1 . Asν increases from 0 to 0.1 and 10 5 (figures 12b and 12c), the main peaks at θ 0 = 0.5π for β 1 = 0 and ±π/6 decline, whereas the |A S |/A-θ 0 curve for β 1 = π/2 rises at θ 0 = 0.5π. Meanwhile, the |A S |/A response in the range of θ 0 ∈ [π, 2.0π] gets stronger and stronger, and the peaks occurring in this range ultimately become as large as those around θ 0 = 0.5π. Figure 13 demonstrates how the energy dissipated by the metamaterial cylinders with different plate alignment directions varies with damping and incident wave direction.

Wave power dissipation
For incident waves incoming with β = π/2 (i.e. beam incidence, figure 13a), as the damping coefficientν increases from 0, η diss first increases and then decreases after reaching the maximum wave power dissipation. This is reasonable as no energy can be dissipated by the metamaterial cylinders forν = 0 and forν → ∞, and meanwhile η diss > 0 forν > 0. The corresponding optimisedν varies for the metamaterial cylinders with different value of β 1 . More specifically, the maximum wave power dissipation and the corresponding optimised damping, are (η diss ,ν) = (10.00, 0.15), (6.13, 0.25), (5.17, 0.35) and (3.13, 0.55) for β 1 = 0, −π/6, π/6, π/2, respectively. For the cylinders with any specified value of damping coefficient, the more perpendicular of the plate alignment relative to the incident wave propagation, the more energy can be dissipated. Note for the two cases with β 1 = −π/6 and β 1 = π/6, the former one performs obviously better than the latter one in terms of wave power dissipation, which might be explained from the view of wave focusing and blocking as studied in § 4.1.
It should be noted that the metamaterial cylinders may be utilised to capture wave energy if the channels are filled with buoys extracting power in heave. Correspondingly, the wave power dissipation represents useful power being consumed by the heaving buoys, i.e. the so called wave power absorption. The surface condition we have used is very similar to the ones Garnaud & Mei (2010) and Garnaud & Mei (2009) derived for arrays of small heaving buoys. For a traditional wave energy converter (WEC) consisting of an axisymmetric rigid cylinder moving in heave mode, it has a maximum capture width P diss /P in and a relative capture width kP diss /P in (i.e. η diss ) of 1/k and 1.0, respectively, which can be achieved with the motion fully optimised (Budal & Falnes 1975;Evans 1976;Newman 1976). When isolated rigid cylinders move in surge/pitch these maximum theoretical values double and in combined surge/pitch and heave, the maximum increases to three times the value for heave only. The present pair of metamaterial cylinders are found to give η diss > 2.0 over a wide range of conditions, absorbing more than two non-interacting heaving cylinders can ever get. For a wide range of incident angle and a pairs of cylinders in beam seas η diss can be significantly greater than 6 meaning that this metamaterial cylinder outperforms the theoretical maximum values for rigid cylinders operating in rigid body modes.

Conclusions
A semi-analytical model based on linear potential flow theory and the eigenfunction matching method has been developed to investigate the interaction of waves with an array of metamaterial circular cylinders consisting of a series of parallel thin plates. To consider the wave attenuation and energy dissipation at narrow gaps between the thin vertical plates, a damping mechanism is introduced at the surface of the fluid occupied by the structured cylinders. In addition to a straightforward way to calculate the energy dissipation, an indirect method is derived based on Kochin functions with the employment of energy identities.
Four case studies: a pair of metamaterial cylinders with the plates aligned to the incident wave direction; a pair of metamaterial cylinders deployed far away from each other; a pair of metamaterial cylinders with an extremely large damping adopted; and a pair of metamaterial with a specified range of damping, were carried out to validate the semi-analytical model. Additionally, potential flow theory based numerical simulations were carried out with the employment of a commercial BEM code to study wave interaction with a metamaterial circular cylinder consisting of 20 thin vertical plates. In these validation cases, the present model is in excellent agreement with both the published data and those obtained by using different methods. The validated model is then applied to investigate the influence of a pair of metamaterial cylinders on wave focusing/blocking and scattered far-field amplitude, and also their performance in wave power dissipation. The effect of a single metamaterial cylinder on the scattered far-field amplitude is studied as well. And the following conclusions may be drawn.
(i) A single cylinder acts as a 'lens' drawing in and emitting a beam of intense wave energy in a direction perpendicular to the plates forming the metamaterial cylinder. (ii) For a pair of metamaterial cylinders without any damping, wave focusing/blocking can be achieved by rotating the cylinder orientation to direct wave energy in the desired manner. The dramatic amplification and focusing effects on wave motion inside the metamaterial cylinders are observed for all the studied cases, except where the plates are aligned in the same direction of the incident wave propagation. (iii) A small amount of damping does not alter the underlying characteristics of wave scattering. For the metamaterial cylinder with an extremely large damping coefficient, it works like a solid lid placed on the surface of the structured cylinder. The orientation of the plates is relatively unimportant as far as the overall effect of the cylinder is on wave diffraction. (iv) There is an optimised damping coefficient to achieve the maximum wave power dissipation of the metamaterial cylinders. A pair of metamaterial cylinders have been shown to exceed the theoretical maximum power for traditional wave power devices moving in rigid body motion.
We explored a range of parameters which gave rise to different features of thē ν associated results. The results we have provided can be mapped into different specific physical interpretations depending on the formula used to connectν to physical parameters; the settings identified here are porous media, channels filled with small heaving buoys and viscous damping in rectangular channels. The semi-analytical model is proposed in the framework of linear potential flow theory, which does not capture viscous effects. Although an artificial linear damping mechanism is included in the model, it may not be suitable for extreme wave-structure interactions.

Appendix B. Derivation process of the energy identities
In the water domain enclosed by Ω 1 ∪ Ω 2 ∪ · · · Ω N ∪ Ω R , free water surface and the sea bed, using Green's theorem (Falnes 2002 For an array of metamaterial cylinders, in which circular solid lids are placed at the still water surface, the boundary condition at z = 0 in each interior domain is ∂φ ∂z = 0, z = 0. (C 1) The spatial velocity potential in the interior domain occupied by cylinder n can be expressed as φ (n) int (x n , y n , z) = Y 0 (z)(B n,0 ( y n )x n + C n,0 ( y n )) + ∞ l=1 Y l (z) B n,l ( y n ) exp(ik l x n ) + C n,l ( y n ) exp(−ik l x n ) = Y 0 (z)[E n,0 (θ n )r n cos(θ n − β n ) + F n,0 (θ n )] + ∞ l=1 Y l (z) E n,l (θ n ) exp(ik l r n cos(θ n − β n )) +F n,l (θ n ) exp(−ik l r n cos(θ n − β n )) , (C 2) in which and, in the same way, the functions E n,l and F n,l can be expressed by (2.21a,b).
Expression of the spatial velocity potential in the exterior domain can be found in (2.13). The same continuity conditions of the field across the interfaces of the interior and exterior domains, i.e. (2.24)-(2.25), should be satisfied as well. After inserting the expressions of the spatial velocity potentials in different domains into the continuity conditions and making use of the orthogonality properties of Z l (z), Y l (z) and exp(imθ n ), we have , (C 5) and ∞ l=0 A (n) τ,l k l H τ (k l R n )S ζ,l m,l (−1) τ H m−τ (k l R n,j )k l J τ (k l R n ) exp(i(mα j,n − τ α n,j ))S ζ,l igAk 0 ω exp(ik 0 (x n cos β + y n sin β))i τ exp(−iτβ)J τ (k 0 R n )S ζ,0 , (C 6) in which Equations (C 4) and (C 6) can be used to determine the unknown coefficients in the expressions of the spatial velocity potentials.