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.

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.

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 dx and dy in water of constant depth 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).

FIG. 1.

(a) Sketch of the cylindrical meta-structure formed by two intersecting arrays of parallel vertical thin plates submerged to different depths; (b) vertical view; (c) side view; and (d) front view.

FIG. 1.

(a) Sketch of the cylindrical meta-structure formed by two intersecting arrays of parallel vertical thin plates submerged to different depths; (b) vertical view; (c) side view; and (d) front view.

Close modal

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 dx and dy, respectively. The 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 d y since the incident wave angle is arbitrary.

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

Linearized theory is adopted under the assumption that the incident wave amplitude is much smaller than the wavelength, i.e., the wave steepness A / λ 1. Assuming that the fluid is incompressible and inviscid, the velocity field u ( x , y , z , t ) = ( u , v , w ) is governed by the mass conservation equation
· u = 0 ,
(1)
and the linearized momentum conservation equation
u t = 1 ρ p ,
(2)
where p ( x , y , z , t ) is defined throughout as the hydrodynamic pressure (i.e., in excess of the hydrostatic background pressure, ρ g z, in which the acceleration due to the gravity g acts in the negative z-direction).
In the outer region, the linearized kinematic and dynamic conditions on the free surface of the fluid are
ζ t = w p = ρ g ζ on r > a , z = 0 ,
(3)
where ζ ( x , y , t ) is the free surface elevation. In the inner region, the vertical displacement of the buoys is represented by a two-dimensional function ξ ( x , y , t ). This is a piecewise constant on the spatial scale 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
ξ t = w ( M 2 t 2 + γ t + σ + ρ g ) ξ = p ¯ on r < a , z = d ,
(4)
where M = ρ d represents the buoy mass per unit area (by the Archimedes principle). Here, Eq. (4) applies in each cell of the structure, and p ¯ is the mean pressure within each cell acting over the base of the buoy occupying that cell. The presentation alludes to the notion developed in the paper that Eq. (4) will be applied as a continuous condition. Note that if we set M = γ = σ = d = 0, the effect of the buoy will be “switched off,” and Eq. (4) reverts to Eq. (3) with ξ ( x , y , t ) representing the free surface elevation. Moreover, the velocity field also satisfies the condition
w = 0 on z = h ,
(5)
representing no normal flow through the horizontal sea bed, and
u · n = 0 ,
(6)
on the vertical plates; n is used to represent the normal into the fluid on fixed surfaces.
If we further assume that the fluid motion is irrotational, that is,
× u = 0 ,
(7)
then there exists a velocity potential Φ ( x , y , z , t ) satisfying
u = Φ ,
(8)
implying, from Eq. (1), that Φ satisfies Laplace's equation in the fluid and, from Eq. (2), that hydrodynamic pressure in the fluid satisfies the linearized Bernoulli equation
p = ρ Φ t .
(9)

Since the separation of adjacent plates, L, is assumed to be much smaller than the typical wavelength, i.e., ε = L / λ 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 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 [ ω 2 h / g = O ( 1 )] and when the depth is small compared to the wavelength or ω 2 h / g 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.

We first present the effective governing equations and boundary conditions satisfied by the leading-order velocity potential after being subjected to the homogenization theory. The detailed derivation is shown in  Appendix A. Linearity of the governing equations has been used to assume the factorization of a time-harmonic variation of angular frequency ω proportional to exp i ω t. The superscript (0) has been dropped for convenience, and we have reverted to dimensional quantities. Therefore, the time-independent velocity potential ϕ = ϕ k ( x , y , z ) in each respective region Ωk has been shown to satisfy
2 ϕ 1 = 0 in Ω 1
(10)
and
{ 2 z 2 ϕ 2 = 0 in Ω 21 , ( 2 y 2 + 2 z 2 ) ϕ 2 = 0 in Ω 22 , 2 ϕ 2 ( 2 x 2 + 2 y 2 + 2 z 2 ) ϕ 2 = 0 in Ω 23 .
(11)
Equation (11) describes that the flow at different depths behaves with different characteristics that are physically quite natural: in Ω 21, the fluid is constrained to oscillate in the z direction, the flow component perpendicular to the plate is suppressed in Ω 22, and the flow under the cylinder is unconstrained.
The combined linearized kinematic and dynamic boundary conditions can be expressed as
(12)
ϕ 1 z = ν ϕ 1 on r a , z = 0
(12a)
and
ϕ 2 z = ν 1 + τ ν d ϕ 2 on   r < a , z = d ,
(12b)
where ν = ω 2 / g and
τ = ( i ω γ + σ ) / ρ g ,
(13)
combines the effect of the spring and damper which are separated in Eq. (12b) from the effect of the mass and hydrostatic stiffness. Obviously, when τ = d = 0, the fluid has the free surface in the inner region. The velocity potentials ϕ k also satisfy the no-normal flow condition
ϕ k z = 0 on   z = h .
(14)
Moreover, at the fluid interface, r = a, between inner and outer regions, continuity of pressure requires
ϕ 1 ( a , θ , z ) = ϕ 2 ( a , θ , z ) , 0 θ < 2 π , h < z < d x ,
(15)
whilst matching fluxes across r = a results in the piecewise conditions
ϕ 1 ( a , θ , z ) r = { 0 , d x z 0 , sin θ ϕ 2 ( a , θ , z ) y , d y z < d x , ϕ 2 ( a , θ , z ) r , h z < d y , 0 θ < 2 π .
(16)
The geometric factor sin θ arises from matching the inner flux confined between parallel plates to the flux into the fluid region outside the cylinder through the circular boundary across an approximated triangular region (see Ref. 3). Note that there are no conditions relating to ϕ 2 or its derivative at r = a from within the overlapping plate region d x < z < d.
In the outer region, the velocity potential satisfying Eqs. (10), (12a), and (14) can be expressed as
ϕ 1 = ϕ inc i g A ω n = i n exp i n θ A 0 n H n ( k r ) ψ 0 ( z ) i g A ω n = i n exp i n θ j = 1 A j n K n ( k j r ) ψ j ( z ) ,
(17)
where
ϕ inc = i g A ω n = i n exp i n ( θ β ) J n ( k r ) ψ 0 ( z )
(18)
represents the incident wave, and Jn, Hn, and Kn denote, respectively, the nth 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, Ajn are undetermined coefficients, kj are the roots of the dispersion equation
ω 2 = g k j tan k j h ,
(19)
such that kj ( j 1 ) are real while k 0 = i k and the wavenumber k is real, and ψ j ( z ) are vertical eigenfunctions
ψ j ( z ) = N j 1 / 2 cos k j ( z + h )
(20)
with
N j = 1 2 ( 1 + sin 2 k j h 2 k j h ) ,
(21)
satisfying the orthogonality relation
1 h h 0 ψ j ( z ) ψ n ( z ) d z = δ j n .
(22)
In the inner region, we will utilize the general expression of the velocity potential proposed by Ref. 18 written as
ϕ 2 ( r , θ , z ) = i g A ω q = 0 π π B q ( t ) exp i κ q ( t ) r cos ( θ t ) Z q ( z , t ) d t = i g A ω q = 0 π π B q ( t ) n = i n exp i n ( θ t ) J n [ κ q ( t ) r ] Z q ( z , t ) d t .
(23)
This expression represents the potential as a sum over all possible wavenumbers κ q ( t ) integrated over all angles t [ π , π ). Since the velocity potential ϕ 2 ( r , θ , z ) is governed by Eq. (11), the vertical function Z q ( z , t ) is expressed as a piecewise function satisfying different equations in each of the three intervals of 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 = d x and z = d y, leading to
Z q ( z , t ) = { ν z + τ + 1 , d x z d , C q ( z , t ) , d y z < d x , C q ( d y , t ) cosh [ κ q ( z + h ) ] cosh [ κ q ( h d y ) ] , h z < d y ,
(24)
in which
C q ( z , t ) = ν sinh [ κ q sin t ( z + d x ) ] κ q sin t + ( τ + 1 ν d x ) cosh [ κ q sin t ( z + d x ) ] ,
(25)
and κ q = κ q ( t ) are the roots of
tanh [ κ q ( h d y ) ] = ν κ q ( τ + 1 ν d x ) sin t tanh [ κ q sin t ( d y d x ) ] κ q ( τ + 1 ν d x ) ( ν / sin t ) tanh [ κ q sin t ( d y d x ) ] .
(26)
Note that Eq. (26) does not depend upon 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.
Equation (26) is complicated, and we should pick out special cases. Thus, for t = 0, Eq. (26) is reduced to
κ q h tanh [ κ q ( h d y ) ] = ν h ( τ + 1 ν d y ) ,
(27)
which is independent of dx. This is since the barriers of depth dx are aligned with the x-axis and transparent to waves traveling at an angle t = 0. Likewise, for t = π / 2, Eq. (26), less obviously becomes
κ q h tanh [ κ q ( h d x ) ] = ν h ( τ + 1 ν d x ) ,
(28)
which is independent of dy. Finally, if dx = dy, Eq. (26) reduces to Eq. (27) or (28) for all values of t. In this latter special case, κ q ( t ) are constant, and Eq. (23) reduces to the more familiar separation series
ϕ 2 ( r , θ , z ) = i g A ω q = 0 n = i n b q n J n ( κ q r ) Z q ( z ) exp i n θ ,
(29)
in which
b q n = π π B q ( t ) exp i n t d t
(30)
represents unknown coefficients. Since Eqs. (27) and (28) have the same form as Eq. (19) if τ 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.
Now we have the expressions of the velocity potential in the inner and outer regions. After substituting Eqs. (17) and (23) into velocity matching condition (16), multiplying ψ j ( z ) exp i n θ on both sides, integrating over h z 0 and 0 θ 2 π and applying Eq. (22) and the orthogonality of functions exp i n θ, we can obtain
1 2 i q = 0 π π B q ( t ) κ q J n 1 ( κ q a ) sin t exp i ( n 1 ) t F q 0 ( 1 ) ( t ) d t + 1 2 i q = 0 π π B q ( t ) κ q J n + 1 ( κ q a ) sin t exp i ( n + 1 ) t F q 0 ( 1 ) ( t ) d t + q = 0 π π B q ( t ) κ q J n ( κ q a ) exp i n t F q 0 ( 2 ) ( t ) d t = { k J n ( k a ) exp i n β + A 0 n k H n ( k a ) , j = 0 , A j n k j K n ( k j a ) , j = 1 , 2 , ,
(31)
where standard recurrence relations for Bessel functions have been used (see Ref. 20, Sec. 9.1.27) and
F q j ( i ) ( t ) = { 1 h d y d x Z q ( z , t ) ψ j ( z ) d z , i = 1 , 1 h h d y Z q ( z , t ) ψ j ( z ) d z , i = 2 ,
(32)
which can be expressed explicitly.
Similarly, we can apply pressure matching condition (15) to the velocity potentials given in Eqs. (17) and (23), multiply on equation both sides by ψ m ( z ) exp i n θ, and integrate over the regions of validity. The orthogonality of functions exp i n θ still exists, but the orthogonality in the vertical direction no longer satisfies since the pressure is continuous only in h z d x. Thus, we have
[ J n ( k a ) exp i n β + A 0 n H n ( k a ) ] E 0 k + j = 1 A j n K n ( k j a ) E j m = q = 0 π π B q ( t ) exp i n t J n ( κ q a ) [ F q m ( 1 ) ( t ) + F q m ( 2 ) ( t ) ] d t ,
(33)
where
E j m = 1 h h d x ψ j ( z ) ψ m ( z ) d z ,
(34)
which also can be expressed explicitly.
Then, after substituting Eq. (31) into Eq. (33), we can obtain an equation system related to the unknown functions B q ( t )
q = 0 π π B q ( t ) M mnq ( t ) exp i n t d t = 2 i exp i n β E 0 m π kaH n ( k a ) ,
(35)
where
M mnq ( t ) = J n ( κ q a ) [ F q 0 ( 1 ) ( t ) + F q 0 ( 2 ) ( t ) ] + i κ q sin t 2 [ J n 1 ( κ q a ) exp i t + J n + 1 ( κ q a ) exp i t ] G mnq ( 1 ) κ q ( t ) J n ( κ q a ) G mnq ( 2 )
(36)
and
G mnq ( i ) = H n ( k a ) E 0 m k H n ( k a ) F q 0 ( i ) ( t ) + j = 1 K n ( k j a ) E j m k j K n ( k j a ) F q j ( i ) ( t ) .
(37)
Note that the series in Eq. (37) decays like O ( j 3 ) which allows us to make efficient and accurate computations.
In order to solve Eq. (35) for the unknowns B q ( t ), we take advantage of the 2 π-periodicity of M mnq ( t ) and expand it in its Fourier series as
M mnq ( t ) = 1 2 π p = M mnpq exp i ( p + n ) t exp i κ q ( t ) a ,
(38)
such that
M mnpq = π π M mnq ( t ) exp i κ q ( t ) a exp i ( p + n ) t d t ,
(39)
where the scale factor exp i κ q ( t ) a is introduced to balance the increasingly exponential behavior of the functions J n ( κ q a ) for q 1 with the aim of suppressing the influence of rounding error on numerical results. Substituting Eq. (38) into Eq. (35) and truncating the infinite system of equations results in
p = N N q = 0 M b p q M mnpq = 2 i exp i n β E 0 m π kaH n ( k a ) ,
(40)
for m = 0 , 1 , 2 , , M (angular mode truncation) and n = N , N + 1 , , N (vertical mode truncation), where
b p q = 1 2 π π π B q ( t ) exp i p t exp i κ q ( t ) a d t
(41)
with
B q ( t ) = p = N N b p q exp i p t exp i κ q ( t ) a .
(42)
In this section, we apply the shallow water theory to consider the interaction of water waves with our cylindrical meta-structure. The detailed derivation of the linear homogenized governing equation is shown in  Appendix C. Returning to dimensional variables and considering motion with a time factor of exp i ω t assumed, the time-independent surface elevation η ( x , y ) outside the cylinder satisfies
h h 2 η + ν η = 0 , in r > a ,
(43)
and the time-independent buoy elevation ξ ( x , y ) inside the cylinder (assuming constant properties of the cylinder in the vertical direction) satisfies
h · ( h h ξ ) + ν 1 + τ ξ = 0 in r < a .
(44)
The wavenumber, k, in shallow water satisfies ν = k 2 h, which is the k h 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,
h = ( h d y 0 0 h d x ) .
(45)
The physical conditions that apply at the interface of inner and outer regions are that the pressure and depth-averaged flux normal to the circular boundary r = a are continuous. In terms of our variables, these require
η = ξ
(46)
and
η r = h d x sin 2 θ d y cos 2 θ h ξ r + ( d y d x ) sin θ cos θ a h ξ θ ,
(47)
on r = a for 0 θ < 2 π (see Refs. 4 and 18 for a similar application of shallow water matching conditions).
In the outer region Ω 1, the elevation satisfying Eq. (43) can be written as
η ( r , θ ) = A n = i n exp i n θ [ J n ( k r ) exp i n β + A n H n ( k r ) ] ,
(48)
while in the inner region Ω 2, the elevation satisfying Eq. (44) can be expressed as
ξ ( r , θ ) = A π π B ( t ) exp i κ ( t ) r cos ( θ t ) d t = A π π B ( t ) n = i n exp i n ( θ t ) J n [ κ ( t ) r ] d t ,
(49)
where
κ ( t ) = ν ( 1 + τ ) ( h d x sin 2 t d y cos 2 t ) .
(50)
It can readily be confirmed that Eq. (50) is the ν 0 limit of Eq. (26), and the wavenumber is now explicit, rather than implicit and requiring roots to be numerically determined. In addition, Eq. (50) inherits the properties of Eq. (26): for example, κ ( t ) is independent of the angle t when dx = dy which is implied by geometric symmetry.
Applying the matching conditions at r = a, and using the orthogonality of functions exp i n θ over 0 θ < 2 π, we obtain
J n ( k a ) exp i n β + A n H n ( k a ) = π π B ( t ) exp i n t J n ( κ a ) d t ,
(51)
and
k J n ( k a ) exp i n β + A n k H n ( k a ) = π π B ( t ) exp i n t κ J n ( κ a ) d t + d y 2 h π π κ B ( t ) cos t exp i n t [ J n 1 ( κ a ) exp i t + J n + 1 ( κ a ) exp i t ] d t + i d x 2 h π π κ B ( t ) sin t exp i n t [ J n 1 ( κ a ) exp i t + J n + 1 ( κ a ) exp i t ] d t .
(52)
Eliminating An between Eqs. (51) and (52) gives
π π B ( t ) M n ( t ) exp i n t d t = 2 i π k a exp i n β ,
(53)
where
M n ( t ) = π π [ J n ( κ a ) H n ( k a ) κ k J n ( κ a ) H n ( k a ) ] B ( t ) exp i n t d t + 1 2 k h H n ( k a ) π π κ d y B ( t ) cos t exp i n t [ J n 1 ( κ a ) exp i t J n + 1 ( κ a ) exp i t ] d t + 1 2 i k h H n ( k a ) π π κ d x B ( t ) sin t exp i n t [ J n 1 ( κ a ) exp i t + J n + 1 ( κ a ) exp i t ] d t .
(54)
Since M n ( t ) = M n ( t + 2 π ), we expand M n ( t ) in its Fourier basis
M n ( t ) = 1 2 π p = M n p exp i ( p + n ) t ,
(55)
from which follows that
M n p = π π M n ( t ) exp i ( p + n ) t d t .
(56)
Substituting Eq. (55) into Eq. (53) and truncating the infinite system of equations, we obtain
p = N N b p M n p = 2 i π k a exp i n β ,
(57)
for n = N , N + 1 , , N, where
b p = 1 2 π π π B ( t ) exp i p t d t ,
(58)
with
B ( t ) = p = b p exp i p t .
(59)
The mean power absorbed by a wave energy device can usually be evaluated using two independent methods that make it a useful tool for checking the accuracy of the numerical scheme. The first method involves balancing the incoming and outgoing energy flux far from the device, and we call this the far-field method. Thus, integrating the time-averaged product of pressure and velocity over the surface, SR, of a cylinder of large radius R extending through the depth, h, gives the mean power as
P f = ω ρ 2 lim R Im [ S R ϕ 1 ϕ 1 * r d S ] ,
(60)
and the asterisk denotes complex conjugation. After substituting Eq. (17) or (48) into the far-field expression (60) and exploiting the theorem of stationary phase (see Ref. 21), it can be expressed as
P f = 2 ρ g A 2 c g k n = 0 { [ A 0 n exp i n β ] + | A 0 n | 2 } ,
(61)
where c g = ( ω / 2 k ) ( 1 + 2 k h / sinh 2 k h ) is the group velocity; A 0 n is replaced with An and c g = ω / k if the shallow water theory is applied.
The second method is related to the first by connecting the boundary SR 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
P n = ω ρ 2 Im [ 0 a 0 2 π ϕ 2 ϕ 2 * z r d θ d r ] z = d .
(62)
This can be interpreted as the time-averaged rate of working of the fluid pressure on the motion of the buoys. Using Eq. (12b) allows us to write
P n = ν 2 γ 2 0 a 0 2 π | ϕ ̃ 2 ( r , θ ) | 2 r d θ d r ,
(63)
where, from Eqs. (23) and (24) in d x < z < d, we have factorized the z-dependence from the integral over t enabling us to write
ϕ 2 ( r , θ , d ) = ( 1 + τ ν d ) ϕ 2 ̃ ( r , θ ) ,
(64)
and the subsequent simplification renders the final expression in Eq. (63) independent of 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 Pn 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 ϕ 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.

The capture width, 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
W = P f / n P inc ,
(65)
where P inc = 1 2 ρ g A 2 c g represents the incident wave power per unit width of the wave crest. Another important measure that determines the economical performance of WEC is the capture width ratio, W / 2 a, obtained by dividing the capture width by the characteristic dimension (in this case, the diameter) of the device.

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 β = π / 4 is considered. For the purpose of illustration, we choose a damping γ / ( ρ g h ) = 0.2 and a spring σ / ρ 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 ν 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.

FIG. 2.

The convergence of the capture width ratio W / ( 2 a ) for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4: (a) N = 8 and (b) M = 8.

FIG. 2.

The convergence of the capture width ratio W / ( 2 a ) for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4: (a) N = 8 and (b) M = 8.

Close modal

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.

FIG. 3.

Comparison of the wave power evaluated by the near-field and far-field methods for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4.

FIG. 3.

Comparison of the wave power evaluated by the near-field and far-field methods for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4.

Close modal

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 Porter23 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 ν d y 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 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.

FIG. 4.

Comparison of wave amplitude computed using full depth-dependent homogenization theory (a) and (c) and a boundary element method for a discrete plate array (b) and (d) for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1 and d y / h = 0.2 with τ = 0 at β = π / 4: (a) and (b) ν d y = 0.2; and (c) and (d) ν d y = 0.4.

FIG. 4.

Comparison of wave amplitude computed using full depth-dependent homogenization theory (a) and (c) and a boundary element method for a discrete plate array (b) and (d) for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1 and d y / h = 0.2 with τ = 0 at β = π / 4: (a) and (b) ν d y = 0.2; and (c) and (d) ν d y = 0.4.

Close modal

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 ν h 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 ν d y. The results show increasing agreement as ν h 0, as expected. The corresponding pressure distribution on the bottom of the buoy at ν d y = 0.1 (corresponding to ν 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.

FIG. 5.

Comparison of the capture width ratio, W / ( 2 a ), for full depth-dependent theory (F) and shallow water theory (S) for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4.

FIG. 5.

Comparison of the capture width ratio, W / ( 2 a ), for full depth-dependent theory (F) and shallow water theory (S) for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4.

Close modal
FIG. 6.

Comparison of pressure distribution on the bottom of the buoy between full depth-dependent theory (a) and shallow water theory (b) for a cylindrical meta-structure of a / h = 2.0 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at ν d y = 0.1 and β = π / 4.

FIG. 6.

Comparison of pressure distribution on the bottom of the buoy between full depth-dependent theory (a) and shallow water theory (b) for a cylindrical meta-structure of a / h = 2.0 , d x / h = 0.1, and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at ν d y = 0.1 and β = π / 4.

Close modal

The wave amplitude for the same cylinder given above but without internal buoys at ν 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 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 d y and h d x replacing the plate submergence depths D and d in Ref. 18.

FIG. 7.

Comparison of wave amplitude between full depth-dependent theory (a) and shallow water theory (b) for a cylindrical meta-structure of a / h = 2.0 , d x / h = 0.1, and d y / h = 0.2 with τ = 0 at ν d y = 0.1 and β = π / 4.

FIG. 7.

Comparison of wave amplitude between full depth-dependent theory (a) and shallow water theory (b) for a cylindrical meta-structure of a / h = 2.0 , d x / h = 0.1, and d y / h = 0.2 with τ = 0 at ν d y = 0.1 and β = π / 4.

Close modal

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 ν 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 dx for an incident wave angle β = π / 4 when ν d y = 0.4 and τ = 0. From  Appendix B, we have known that as dx approaches dy, 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.

FIG. 8.

Wave amplitude for a cylindrical meta-structure of a / h = 0.5 and d y / h = 0.2 with τ = 0 at β = π / 4 and ν d y = 0.4: (a) d x / h = 0.0; (b) d x / h = 0.1; and (c) d x / h = 0.2.

FIG. 8.

Wave amplitude for a cylindrical meta-structure of a / h = 0.5 and d y / h = 0.2 with τ = 0 at β = π / 4 and ν d y = 0.4: (a) d x / h = 0.0; (b) d x / h = 0.1; and (c) d x / h = 0.2.

Close modal

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 γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2, dx 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 dx. As the incident wave angle β increases, the incident wave direction is gradually aligned with the array of plates in the y direction. As this happens, the value of dx has a greater influence, and more power is extracted at low frequencies as the value of dx is increased, although interestingly at higher frequencies results again suggest a lower value of dx is optimal. Note that when dx = dy, the wave power extraction is independent of the incident wave direction since the description of the structured cylinder under homogenization is axisymmetric.

FIG. 9.

The variation of wave power generated by a cylindrical meta-structure of a / h = 0.5 and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2: (a) β = 0; (b) β = π / 4; and (c) β = π / 2. (F: full depth-dependent theory; S: shallow water theory.)

FIG. 9.

The variation of wave power generated by a cylindrical meta-structure of a / h = 0.5 and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2: (a) β = 0; (b) β = π / 4; and (c) β = π / 2. (F: full depth-dependent theory; S: shallow water theory.)

Close modal

The results of the shallow water theory in the range of v d y [ 0 , 0.5 ] are also plotted in Fig. 9 for comparison. We can see that the wave power is little affected by dx. When d x / h = 0 and β = π / 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 β = π / 4 and ν 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 dx = 0, the contours are roughly aligned with the remaining single array of plates. As dx approaches dy, the contours become increasingly parallel to the direction of the incident wave crest.

FIG. 10.

The vertical displacement of buoys for a cylindrical meta-structure of a / h = 0.5 and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4 and ν d y = 0.4: (a) d x / h = 0.0; (b) d x / h = 0.1; and (c) d x / h = 0.2.

FIG. 10.

The vertical displacement of buoys for a cylindrical meta-structure of a / h = 0.5 and d y / h = 0.2 with γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4 and ν d y = 0.4: (a) d x / h = 0.0; (b) d x / h = 0.1; and (c) d x / h = 0.2.

Close modal

Then, the curves of wave power extracted from two cylindrical meta-structures of a / h = 0.5 with different plate submergence at β = π / 4 are shown in Fig. 11. The settings of spring and damping are the same above, i.e., γ / ( ρ g h ) = 0.2 and σ / ρ 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 ν c h = 3.0. Since the inner free surface is covered by the buoys, the resonance is suppressed by the damping mechanism.

FIG. 11.

The variation of wave power generated by cylindrical meta-structures of a / h = 0.5 with different plate submergence for γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4.

FIG. 11.

The variation of wave power generated by cylindrical meta-structures of a / h = 0.5 with different plate submergence for γ / ( ρ g h ) = 0.2 and σ / ρ g = 0.2 at β = π / 4.

Close modal

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 , d x / h = 0.1, and d y / h = 0.2 at β = π / 4 with different damping rates and spring constants. In Fig. 12(a), where γ / ( ρ 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 σ / ρ g = 0.2, we see that different damping parameters work best at different frequencies.

FIG. 12.

Effect of damping rate γ and spring constant σ on the efficiency of wave power extraction for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 at β = π / 4 with (a) γ / ( ρ g h ) = 0.4 and (b) σ / ρ g = 0.2.

FIG. 12.

Effect of damping rate γ and spring constant σ on the efficiency of wave power extraction for a cylindrical meta-structure of a / h = 0.5 , d x / h = 0.1, and d y / h = 0.2 at β = π / 4 with (a) γ / ( ρ g h ) = 0.4 and (b) σ / ρ g = 0.2.

Close modal

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 γ / ( ρ g h ) = 0.5 and σ / ρ g = 0.0 at β = π / 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 1.0 when a / h = 1 or k h 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.

FIG. 13.

The variation of wave power extraction for a cylindrical meta-structure of d x / h = 0.1 , d y / h = 0.2 with γ / ( ρ g h ) = 0.5 and σ / ρ g = 0.0 at β = π / 4 against non-dimensional wavenumber νh: (a) non-dimensional capture width kW; and (b) capture width ratio W / ( 2 a ) [comparison between the present model (P) and Ref. 15 (G&M)].

FIG. 13.

The variation of wave power extraction for a cylindrical meta-structure of d x / h = 0.1 , d y / h = 0.2 with γ / ( ρ g h ) = 0.5 and σ / ρ g = 0.0 at β = π / 4 against non-dimensional wavenumber νh: (a) non-dimensional capture width kW; and (b) capture width ratio W / ( 2 a ) [comparison between the present model (P) and Ref. 15 (G&M)].

Close modal

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.

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).

The authors have no conflicts to disclose.

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).

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

In this appendix, we will present a detailed derivation of effective governing equations and boundary conditions using a homogenization approach, without making any assumptions about the depth of the fluid. Since homogenization is a spatial operation, we could retain full time dependence in the derivation below. However, we only make single-frequency computations and are in the interest of simplifying the presentation of the method, so we first let Φ ( x , y , z , t ) = { ϕ ( x , y , z ) exp i ω t }. Similar factorizations of time can be applied to all other time-dependent functions such that the full problem outlined in Sec. II may be posed entirely in terms of ϕ with
2 ϕ = 0 in the fluid ,
(A1)
the combined kinematic and dynamic conditions
ϕ z ν ϕ = 0 on z = 0 , r > a
(A2)
and
ϕ z = ν 1 + τ ν d ϕ ¯ on z = d , r < a ,
(A3)
where ν = ω 2 / g and τ = ( i ω γ + σ ) / ρ g,
ϕ z = 0 on z = h ,
(A4)
and n · ϕ = 0 on all vertical sides of the structured plate array. Our depth assumption translates to ν h = O ( 1 ). In Eq. (A3), ϕ ¯ represents the spatial average of ϕ across a single cell, following the definition of p ¯ introduced in Eq. (4).
In Ω 21, since there are two intersecting arrays of closely spaced parallel plates, we introduce the multiple scales
x = L X + λ x , y = L Y + λ y , z = h z ,
(A5)
where ( x , y , z ) are the macroscopic field variables while (X, Y) are microscopic field variables based on a single cell such that 0 X 1 and 0 Y 1. Simultaneously, we shall define ε = L / λ and expand variables as
ϕ ( x , y , z ) = ϕ ( 0 ) ( x , y , z ; X , Y ) + ε ϕ ( 1 ) ( x , y , z ; X , Y ) + ε 2 ϕ ( 2 ) ( x , y , z ; X , Y ) + .
(A6)
From Eq. (A1), we have
[ ( 2 X 2 + 2 Y 2 ) + 2 ε ( 2 x X + 2 y Y ) + ε 2 ( 2 x 2 + 2 y 2 + λ 2 h 2 2 z 2 ) ] ( ϕ ( 0 ) + ε ϕ ( 1 ) + ) = 0 ,
(A7)
and λ / h is assumed to be O(1). From Eq. (A3), we have
z ( ϕ ( 0 ) + ε ϕ ( 1 ) + ) = ν h 1 + τ ν d 0 1 0 1 ( ϕ ( 0 ) + ε ϕ ( 1 ) + ) d X d Y ,
(A8)
on z = d / h. The velocity potential also satisfies the no-normal flow condition on either side of each plate, which implies that
( X + ε x ) ( ϕ ( 0 ) + ε ϕ ( 1 ) + ) = 0 on X = 0 , 1 , ( Y + ε y ) ( ϕ ( 0 ) + ε ϕ ( 1 ) + ) = 0 on Y = 0 , 1.
(A9)
We can easily confirm that the zero-order velocity potential is independent of X and Y [i.e., ϕ ( 0 ) = ϕ ( 0 ) ( x , y , z )] due to the fact that
Ω | H ϕ ( 0 ) | 2 d X d Y = Ω ϕ ( 0 ) N · H ϕ ( 0 ) d S Ω H 2 ϕ ( 0 ) d X d Y = 0 ,
(A10)
after use of Eqs. (A7) and (A9), where H = ( / X , / Y ) , Ω = { 0 X 1 , 0 Y 1 } , Ω is the boundary of Ω and N is the unit normal in terms of (X, Y) variables out of Ω.
At O ( ε ), we have from Eq. (A7)
( 2 X 2 + 2 Y 2 ) ϕ ( 1 ) = 2 ( 2 x X + 2 y Y ) ϕ ( 0 ) = 0 ,
(A11)
again and
X ϕ ( 1 ) = x ϕ ( 0 ) on X = 0 , 1 , Y ϕ ( 1 ) = y ϕ ( 0 ) on Y = 0 , 1.
(A12)
After applying the method of separation of variables, we can obtain the expression of the first-order velocity potential
ϕ ( 1 ) = ϕ ( 0 ) x X ϕ ( 0 ) y Y + f ( x , y , z ) ,
(A13)
where f ( x , y , z ) is an arbitrary function of the macroscopic coordinates.
Finally, we consider the second-order problem O ( ε 2 ) and Eq. (A7) gives us
( 2 X 2 + 2 Y 2 ) ϕ ( 2 ) = 2 ( 2 x X + 2 y Y ) ϕ ( 1 ) ( 2 x 2 + 2 y 2 + λ 2 h 2 2 z 2 ) ϕ ( 0 ) = ( 2 x 2 + 2 y 2 λ 2 h 2 2 z 2 ) ϕ ( 0 ) ,
(A14)
where Eq. (A13) has been used, and
X ϕ ( 2 ) = x ϕ ( 1 ) on X = 0 , 1 , Y ϕ ( 2 ) = y ϕ ( 1 ) on Y = 0 , 1.
(A15)
Integrating Eq. (A14) over 0 X 1 and 0 Y 1 and applying the boundary condition (A15) with Eq. (A13), we finally obtain
2 z 2 ϕ ( 0 ) = 0 , d x / h < z < d / h ,
(A16)
as the governing equation for the leading-order problem subject to the boundary condition
ϕ ( 0 ) z = ν h 1 + τ ν d ϕ ( 0 ) on z = d / h ,
(A17)
from Eq. (A8) at O ( ε 0 ).
In Ω 22, a similar procedure can be performed. Since the fluid is only separated by the vertical plates aligned in the y direction, we rescale with
x = L X + λ x , y = λ y , z = h z
(A18)
with the understanding that the 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
( 2 y 2 + λ 2 h 2 2 z 2 ) ϕ ( 0 ) = 0 ,
(A19)
again for ϕ ( 0 ) = ϕ ( 0 ) ( x , y , z ).

For the leading-order governing equation in Ω 23 or outer region Ω 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 = d x and z = d y on the interface between Ω 21 and Ω 22 and between Ω 22 and Ω 23 which involves rescaling z on the lengthscale L via z = d x , y + L Z reveals that, at leading order, ϕ z ( 0 ) and ϕ ( 0 ) are continuous across z = d x , y. Readers are directed to Ref. 5 for details of the application of this process in a similar problem.

In the numerical process of finding the roots of Eq. (26), we first solve the equation with Im [ τ ] = 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 κ q ( t ) is a root of Eq. (26), then so is κ q ( t ). Since t is integrated over all angles, π t < π, the real positive root κ 0 ( t ) representing the propagating wave and an infinite number of pure imaginary roots κ q ( t ) = i κ ̂ q ( t ) on the positive imaginary axis representing the evanescent waves need to be considered. Furthermore, it can be observed from Eq. (26) that κ q ( t ) = κ q ( t ) and κ q ( π t ) = κ q ( t ) which confine us to find the roots κ q ( t ) only in the range of t [ 0 , π / 2 ] and helps reduce the numerical effort.

From Eq. (27), we can see that the real root κ 0 ( t ) will tend to infinity as ν d y 1 + [ τ ] and only exists for 0 t π / 2 when ν d y < 1 + [ τ ], which means that waves do not propagate within the inner region when ν > ν c and ν c d y = 1 + [ τ ] 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 [ 0 , π / 2 ] when Im [ τ ] = 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 κ q ( t ) move off the real or imaginary axis into the complex plane. We can take the root at Im [ τ ] = 0 as the initial value [denoted as κ 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 κ q 0 ( t ) as an initial value and obtain a new root κ q 1 ( t ) after iteration; (ii) if | ( κ q 1 ( t ) κ q 0 ( t ) ) / κ q 0 ( t ) | is less than a threshold Δ κ [which can guarantee that for a certain q the root does not jump to other branches (see Fig. 15)], κ 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 [ τ ] ] to get κ q 2 ( t ); and (iii) if we still have | ( κ q 2 ( t ) κ q 0 ( t ) ) / κ q 0 ( t ) | > Δ κ, the interval is halved again until the termination condition is satisfied—otherwise, we will take κ q 2 ( t ) as an initial value and repeat procedure (i) on the interval [ Im [ τ ] / 2 , Im [ τ ] ].

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 κ 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 + [ τ ], we do not have the corresponding initial value to find κ 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 κ 0 ( t ) for t = π / 4 and t = π / 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 dx = dy, all the curves of κ 0 ( 0 ) for three cases are the same as the black dot dash line in Fig. 14(a) or 14(b). For a certain t, the real root, κ 0 ( t ), increases with the frequency until ν d x = 1 since the real root exists only when ν d x < 1 as mentioned above. In addition, when dx tends to dy, κ 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.

FIG. 14.

The variation of the real root κ 0 ( t ) of the dispersion equation (26) with τ = 0 and d y / h = 0.2: (a) t = π / 4 and (b) t = π / 2.

FIG. 14.

The variation of the real root κ 0 ( t ) of the dispersion equation (26) with τ = 0 and d y / h = 0.2: (a) t = π / 4 and (b) t = π / 2.

Close modal

Furthermore, if τ has an imaginary part, κ q ( t ) will go into the complex plane. In Fig. 15, we let σ = 0 and γ / ( ρ g h ) = 0.2 and present the locations of the first four root κ q ( π / 4 ) ( q = 0 , 1 , 2 , 3 ) at each step as Im [ τ ] increases from zero. κ q ( π / 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.

FIG. 15.

The paths of the first four root κ q ( π / 4 ) ( q = 0 , 1 , 2 , 3 ) as Im [ τ ] increases from zero for ν d y = 2.5 , d y / h = 0.2, σ = 0, and γ / ( ρ g h ) = 0.2.

FIG. 15.

The paths of the first four root κ q ( π / 4 ) ( q = 0 , 1 , 2 , 3 ) as Im [ τ ] increases from zero for ν d y = 2.5 , d y / h = 0.2, σ = 0, and γ / ( ρ g h ) = 0.2.

Close modal

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., μ = h / λ 1. Additionally, the spacing of plates should be sufficiently small compared to the water depth, i.e., L / h 1, and we specifically require ε = O ( μ 2 ).

Although we treat various quantities including the depths h, dx, and dy as fixed in the derivation below, a more general derivation (e.g., see Ref. 4) allows these (in addition to mechanical parameters M, σ, 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.

We first nondimensionalize the variables in the governing equations and boundary conditions given in Sec. II as
( h , d x , d y ) = h ( 1 , d x , d y ) , ξ = A ξ , p = ρ gAp , u h = ( A / h ) g h u h , w = ( A / λ ) g h w , t = ( λ / g h ) t , M = ( ρ λ 2 / h ) M , γ = ( ρ λ g / h ) γ , σ = ρ g σ ,
(C1)
where 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 = ρ d and d < h, then M = O ( μ 2 ) = O ( ϵ ), and we therefore write M = ϵ M .
In Ω 21, we introduce multiple scales in the horizontal coordinates scaled by
x = L X + λ x , y = L Y + λ y , z = h z ,
(C2)
and expand variables in the parameter ε = L / λ 1 with
{ u , p } = { u , p } ( 0 ) ( x , y , z , t ; X , Y ) + ε { u , p } ( 1 ) ( x , y , z , t ; X , Y ) + ,
(C3)
and since the buoy elevation is constant in each cell,
ξ = ξ ( 0 ) ( x , y , z , t ) + ε ξ ( 1 ) ( x , y , z , t ) + .
(C4)
For clarity, we drop the primes on the both dependent and independent variables hereafter. The mass conservation equation (1) can be written as
H · ( u h ( 0 ) + ε u h ( 1 ) + ) + ε h · ( u h ( 0 ) + ε u h ( 1 ) + ) + ε z ( w ( 0 ) + ε w ( 1 ) + ) = 0 ,
(C5)
where we have use the notation H = ( / X , / Y ) and h = ( / x , / y ). The momentum equation (2) is expressed as
t ( u h ( 0 ) + ε u h ( 0 ) + ) = ( 1 ε H + h ) ( p ( 0 ) + ε p ( 1 ) + ) , ε t ( w ( 0 ) + ε w ( 1 ) + ) = z ( p ( 0 ) + ε p ( 1 ) + ) ,
(C6)
where ε = O ( μ 2 ) has been applied. The corresponding boundary conditions (linearized) on z = d / h are
t ( ξ ( 0 ) + ε ξ ( 1 ) + ) = w ( 0 ) + ε w ( 1 ) +
(C7)
and
( ϵ M 2 t 2 + γ t + σ + 1 ) ( ξ ( 0 ) + ε ξ ( 1 ) + ) = 0 1 0 1 ( p ( 0 ) + ε p ( 1 ) + ) d X d Y
(C8)
with
[ ( u h , w ) ( 0 ) + ε ( u h , w ) ( 1 ) + ] · N = 0 ,
(C9)
on 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
H w ( 0 ) = 0 , and u ( 0 ) Y v ( 0 ) X = 0.
(C10)
We continue by considering terms at O ( ε 0 ) and making use of Eq. (C10) to show that w ( 0 ) = w ( 0 ) ( x , y , z , t ) and by using Eq. (C5) with Eqs. (C9) and (C10) to deduce that
u ( 0 ) = v ( 0 ) = 0.
(C11)
Next, the terms at the order O ( ε ) in Eq. (C5) integrated over 0 < X < 1 and 0 < Y < 1 and applying Eq. (C9) result in
w ( 0 ) z = 0 ,
(C12)
so that now w ( 0 ) = w ( 0 ) ( x , y , t ). From the kinematic boundary condition in Eq. (C7), we have
w ( 0 ) ( x , y , t ) = ξ ( 0 ) t .
(C13)
Taking the leading order terms in Eq. (C6), we see that p ( 0 ) = p ( 0 ) ( x , y , t ). Using this in the dynamic boundary condition in Eq. (C8), we can obtain
p ( 0 ) ( x , y , t ) = ( γ t + σ + 1 ) ξ ( 0 ) ( x , y , t ) .
(C14)
In the region Ω 22 occupied by plates parallel to the y-axis, we can apply the same treatment to the governing equations and boundary conditions, after replacing the scaling in Eq. (C2) by
x = L X + λ x , y = λ y , z = h z ,
(C15)
and correspondingly drop the dependence on Y from the dependent variables. We can deduce that u ( 0 ) is independent of the microscopic coordinate X. In addition, we also have
u ( 1 ) X + v ( 0 ) y + w ( 0 ) z = 0
(C16)
and
p ( 0 ) y + v ( 0 ) t = 0.
(C17)
From the vertical component of momentum, p ( 0 ) is independent of 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 d y / h < z < d x / h, we have
d y d x h v ( 0 ) y = w ( 0 ) ( x , y , d y / h , t ) w ( 0 ) ( x , y , d x / h , t ) ,
(C18)
where w ( 0 ) is independent of X as a consequence of the leading order contribution from the irrotationality condition Eq. (7) under the multiple scales expansion.
Analogously, in Ω 23 which is the lower fluid region and free from barriers, we do not require multiple scales and now the leading order contributions to the momentum equation provide us with
h p ( 0 ) + t u h ( 0 ) = 0 and p ( 0 ) z = 0 ,
(C19)
such that p ( 0 ) = p ( 0 ) ( x , y , t ) and hence also u ( 0 ) = u ( 0 ) ( x , y , t ) and v ( 0 ) = v ( 0 ) ( x , y , t ). Thus, if we integrate the mass conservation equation over 1 < z < d y / h, we obtain
w ( 0 ) ( x , y , d y / h , t ) = ( 1 d y h ) ( u ( 0 ) x + v ( 0 ) y ) ,
(C20)
where the boundary condition of no-normal flow on the sea bed has been applied.
Matching conditions across the interfaces z = d x / h and z = d y / h can be considered carefully using methods described in Ref. 4 and result in the requirement that the leading order pressure and normal component of velocity should be continuous across both interfaces. This information first allows us to connect Eqs. (C13), (C18), and (C20) to result in
ξ ( 0 ) t = h · ( h u h ( 0 ) ) ,
(C21)
where
h = ( 1 d y / h 0 0 1 d x / h )
(C22)
is a two-dimensional tensor. Second, connecting p ( 0 ) from Eqs. (C14) and (C19) gives
t u h ( 0 ) = h ( γ t + σ + 1 ) ξ ( 0 ) .
(C23)
Combining Eq. (C21) with Eq. (C23) and eliminating the leading-order horizontal velocity u h ( 0 ), we can get
2 t 2 ξ ( 0 ) = h · [ h h ( γ t + σ + 1 ) ξ ( 0 ) ] .
(C24)

The governing equation in Ω 1 can be obtained by letting γ = σ = d x = d y = 0 in Eq. (C24), which reduces to the standard shallow water equation over constant depth h. Despite setting h, γ, σ, dx, and dy 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.

1.
C. P.
Berraquero
,
A.
Maurel
,
P.
Petitjeans
, and
V.
Pagneux
, “
Experimental realization of a water-wave metamaterial shifter
,”
Phys. Rev. E
88
,
051002
(
2013
).
2.
A.
Maurel
,
J.-J.
Marigo
,
P.
Cobelli
,
P.
Petitjeans
, and
V.
Pagneux
, “
Revisiting the anisotropy of metamaterials for water waves
,”
Phys. Rev. B
96
,
134310
(
2017
).
3.
R.
Porter
, “
Plate arrays as a perfectly-transmitting negative-refraction metamaterial
,”
Wave Motion
100
,
102673
(
2021
).
4.
C.
Marangos
and
R.
Porter
, “
Shallow water theory for structured bathymetry
,”
Proc. R. Soc. A
477
,
20210421
(
2021
).
5.
R.
Porter
and
C.
Marangos
, “
Water wave scattering by a structured ridge on the sea bed
,”
Ocean Eng.
256
,
111451
(
2022
).
6.
V.
Gupta
,
S.
Adhikari
, and
B.
Bhattacharya
, “
Exploring the dynamics of hourglass shaped lattice metastructures
,”
Sci. Rep.
10
,
20943
(
2020
).
7.
S.
Zheng
,
R.
Porter
, and
D.
Greaves
, “
Wave scattering by an array of metamaterial cylinders
,”
J. Fluid Mech.
903
,
A50
(
2020
).
8.
H. J.
Putley
,
S.
Guenneau
,
R.
Porter
, and
R. V.
Craster
, “
A tuneable electromagnetic metagrating
,”
Proc. R. Soc. A
478
,
20220454
(
2023
).
9.
F. D.
Vita
,
F. D.
Lillo
,
F.
Bosia
, and
M.
Onorato
, “
Attenuating surface gravity waves with mechanical metamaterials
,”
Phys. Fluids
33
,
047113
(
2021
).
10.
H.
Pang
,
A.
Feng
,
Y.
You
, and
K.
Chen
, “
Design of novel energy harvesting device based on water flow manipulation
,”
Phys. Fluids
34
,
093609
(
2022
).
11.
H.
Pang
,
A.
Feng
,
Y.
You
, and
K.
Chen
, “
Hydrodynamic manipulation cloak for redirecting fluid flow
,”
Phys. Fluids
34
,
053603
(
2022
).
12.
D. V.
Evans
, “
A theory for wave-power absorption by oscillating bodies
,”
J. Fluid Mech.
77
,
1
25
(
1976
).
13.
K.
Budar
and
J.
Falnes
, “
A resonant point absorber of ocean-wave power
,”
Nature
256
,
478
479
(
1975
).
14.
J.
Huang
and
R.
Porter
, “
Wave power absorption by a metamaterial cylinder with internal paddle power take-off system
,”
Appl. Ocean Res.
128
,
103315
(
2022
).
15.
X.
Garnaud
and
C. C.
Mei
, “
Wave-power extraction by a compact array of buoys
,”
J. Fluid Mech.
635
,
389
413
(
2009
).
16.
X.
Garnaud
and
C. C.
Mei
, “
Comparison of wave power extraction by a compact array of small buoys and by a large buoy
,”
IET Renewable Power Gener.
4
,
519
530
(
2010
).
17.
R.
Porter
,
S.
Zheng
, and
D.
Greaves
, “
Extending limits for wave power absorption by axisymmetric devices
,”
J. Fluid Mech.
924
,
A39
(
2021
).
18.
R.
Porter
,
S.
Zheng
, and
H.
Liang
, “
Scattering of surface waves by a vertical truncated structured cylinder
,”
Proc. R. Soc. A
478
,
20210824
(
2022
).
19.
B.
Wilks
,
F.
Montiel
, and
S.
Wakes
, “
Rainbow reflection and broadband energy absorption of water waves by graded arrays of vertical barriers
,”
J. Fluid Mech.
941
,
A26
(
2022
).
20.
M.
Abramowitz
and
I. A.
Stegun
,
Handbook of Mathematical Functions
(
Dover
,
New York
,
1965
).
21.
C. C.
Mei
,
The Applied Dynamics of Ocean Surface Waves
(
World Scientific
,
1989
).
22.
B.
Teng
and
J.
Huang
, “
Examination of extraordinary transmission of waves propagation through gaps of vertical thin barriers in channels by a hypersingular boundary element method
,”
China Ocean Eng.
33
,
509
521
(
2019
).
23.
J.
Huang
and
R.
Porter
, “
Water wave propagation through arrays of closely spaced surface-piercing vertical barriers
,”
J. Fluid Mech.
960
,
A20
(
2023
).
24.
A. F. d. O.
Falcão
, “
Wave energy utilization: A review of the technologies
,”
Renewable Sustainable Energy Rev.
14
,
899
918
(
2010
).
Published open access through an agreement withJISC Collections