Ionut Trestian

Kévin Huguenin
email: kevin.huguenin@epfl.ch

Ling Su

Aleksandar Kuzmanovic

Understanding Human Movement Semantics: A Point of Interest Based Approach *

Keywords: H.2.8 [Database Management]: Database Applications-Data mining; H.1.2 [Information Systems]: User/Machine Systems-Human information processing Human Factors, Measurement GPS trajectories analysis, Mobility, Point of interest

The recent availability of human mobility traces has driven a new wave of research on human movement with straightforward applications in wireless/cellular network. In this paper we revisit the human mobility problem with new assumptions. We believe that human movement is not independent of the surrounding locations, i.e. the points of interest that they visit; most of the time people travel with specific goals in mind, visit specific points of interest, and frequently revisit favorite places. Using GPS mobility traces of a large number of users located across two distinct geographical locations we study the correlation between people's trajectories and the differently spread points of interest nearby.

INTRODUCTION

Human movement has been studied extensively in the past several years. This is mainly due to the rather obvious implications that results might have in fields such as urban planning, disease prevention, mobile advertising, mobile infrastructure placement, etc. As such, researchers have tried to understand statistical properties, and also proposed mathematical models that try to capture particular aspects of human movement [START_REF] Gonzalez | Understanding Individual Human Mobility Patterns[END_REF], for example by relating to the movements of banknotes [START_REF] Brockmann | The Scaling Laws of Human Travel[END_REF] or subatomic particles. But the matter of fact is that, contrary to banknotes and subatomic particles, human mobility is driven by concrete goals such as buying groceries, dining at a restaurant, watching a movie, etc. And although it is indeed true that models (e.g. Lévy flights) manage to capture important aspects of movement (e.g., the distribution of trip sizes), they fail to capture the particular semantics incorporated in human movement.

Furthermore, as we noted above, it is sensible to believe that most human movement is indeed driven by certain purposes. We can therefore only note that specific and important characteristics of human movement that should get captured in the aforementioned models are linked to two important factors. The first factor consists of people's daily routines (such as visiting a coffee shop in the morning, going to work, and/or shopping or going to a restaurant in the evening). A second factor is the actual location of the different points of interest that humans visit. The combination of these two factors could explain the characteristics that older studies on human movement have observed.

In this paper we aim to answer the following research questions. What are the actual semantics of human movement? What are the locations that people visit, with what frequency, and at what time scales? What is the influence of the spread of these locations on human movement?

DATASETS

In an effort to make our study more generic we use two datasets containing user trajectories that were collected in different regions of the globe. The first dataset we use is a GPS trajectory dataset collected by Microsoft Research Asia in their GeoLife project [START_REF]GeoLife GPS Trajectories[END_REF]. It contains the data of 155 users that were monitored over a period of over two years. The second dataset we are using is a dataset collected in the United States from the backbone of a mobile provider network. It contains the trajectories of 4,429 users that fielded applications that report their GPS locations to the provider. The above datasets are enhanced with the GPS coordinates of points of interest from China and the United States. Statistics for some of the categories we extracted can be seen in Table 1. We show the points of interest for China and the states of Illinois, Indiana, and Michigan (where most of our American users are).

Our analysis requires a proper formalization of the notion of interactions between users and possibly grouped points of interest. The first stage consists in turning 1D trajectories and positions into 2D areas in such a way that an interaction can be thought of as intersecting areas over a time interval. User positions and points of interest are therefore replaced by discs of radius r and rp respectively.

(a) Merging businesses. To capture the fact that businesses are often grouped, e.g. in shopping malls, we use a region-based segmentation approach defining meta points of interest for points of interest of the same type. When two (meta) points of interest are close enough, i.e. the intersection between their respective discs is greater than a given threshold, they are merged into a single meta point of interest represented by the minimum enclosing disc (see Fig. 1a).

We define a model of interaction inspired from physics: In our analysis we used the area of intersection as a power function. The energy received by a point of interest is obtained by summing the energy of its interactions with all the segments of all users' trajectories. Such an energy-based approach captures both the spatial and temporal aspects of interactions. In case of an interaction with a meta point of interest, the energy of the interaction is evenly split among the points of interest composing the meta point of interest.

EVALUATION

By using the methodology we defined above, we evaluate user interaction with points of interest to answer the question about which point of interest types are more popular. In general, people do not interact with different types of points of interest in the same way. In particular, people visit certain places frequently depending on their function: grocery shopping, dining at a restaurant, watching a game in a stadium, etc.

Fig. 2 shows the results. The histograms are scaled to the highest energy value observed for each dataset, in both cases shops. We purposely do not show values for companies because values for interactions that users experience with the companies they work at dwarf other values.

A few findings from the figure are as follows. The most important is that the results are surprisingly consistent across the two datasets even though they were collected in different parts of the globe. Further, restaurants and shops lead the energy-of-interaction race in both figures. They hold the most energy out of all the analyzed point of interest types. This is not surprising as they are part of people's daily activities. In addition, although many more users are seen at stadiums for the duration of the trace, users do not go to stadiums frequently. In fact libraries and stadiums compare in terms of energy of interaction most probably because the users that do visit libraries do so more frequently than the users that visit stadiums.

Further insights can be found in [START_REF] Trestian | Understanding Human Movement Semantics: A Point of Interest Based Approach[END_REF], namely (i) time of day effects for point of interest types, and (ii) the number of users seen at specific businesses that carry publicly available WiFi access points (e.g., McDonald, Starbucks) that are potential targets for data offloading for mobile users.

CONCLUSIONS

In this paper we have proposed a new approach to analyzing the mobility of humans. Noting the lack of expressiveness incorporated in to existing studies, we proposed and performed the first-of-its-kind joint analysis of human mobility correlated with the surrounding environment, i.e. the points of interest that they visit. Our key finding is in demonstrating user affinity (with different visiting frequency depending on point of interest type) towards specific points of interest and specific businesses. We believe our results show noteworthy promise for further research in this area, clearing the way for future advances in understanding basic human behavior and impacting problems related to mobile transfer scheduling algorithms, mobile infrastructure placement, mobile social networking, mobile advertising, etc.

Figure 1 :

 1 Figure 1: Methodology of our analysis.

  given two successive snapshots (ti, xi) and (ti+1, xi+1) of a user and a point of interest p at position xp, the power of the interaction at time t is defined as a function p(•) of the distance between the user's position x(t) and the position xp of the point of interest. The energy E of the interaction during the time interval [ti, ti+1] is obtained through integration: E((ti, xi), (ti+1, xi+1), xp) = t i+1 t i p(d(x(t), xp)) dt .

Figure 2 :

 2 Figure 2: Energy of interaction for users interacting with the specific point of interest types.

Table 1 : Point of interest statistics.

 1

	Dataset	Illinois Indiana Michigan China
	Restaurants	19,400	9,231	14,671	116,095
	Shops	56,768	25,454	46,199	267,541
	Govt. Offices	4,645	2,890	5,185	137,837
	Hospitals, Clinics	748	466	963	54,545
	Libraries	692	372	585	40,361
	Stadiums	71	24	61	17,868
	Population [M]	12.91	6.42	9.96	1,331
	Size [1000 sq. mi]	57.9	36.4	96.7	3,705


G C † A N D C H I A N

Wave-power extraction by a compact array of buoys

The majority of existing single-unit devices for extracting power from sea-waves relies on resonance at the peak frequency of the incident wave spectrum. Such designs usually call for structural dimensions not too small compared to a typical wavelength and yield high efficiency only within a limited frequency band. A recent innovation in Norway departs from this norm by gathering many small buoys in a compact array. Each buoy is too small to be resonated in typical sea conditions. In this article a theoretical study is performed to evaluate this new design. Within the framework of linearisation, we consider a periodic array of small buoys with similarly small separation compared to the typical wave length. The method of homogenization (multiple scales) is used to derive the equations governing the macro-scale behaviour of the entire array. These equations are then applied to energy extraction by an infinite strip of buoys, and by a circular array. In the latter case advantages are found when compared to a single buoy of equal volume.

Introduction

The prevailing ideas of wave-power extraction are based on matching the impedance of the extracting device to the characteristics of the incident wave. In particular for a single unit of an oscillating body, or for an oscillating water column, the device should be resonated at the peak frequency of the incoming wave and the extraction rate should equal that of the radiation damping. High efficiency is attainable in a limited frequency bandwidth around resonance. Ideas have been proposed to broaden the bandwidth by the method of phase control [START_REF] Budal | Power from Sea Waves, chap. Interacting point absorbers with controlled motion[END_REF]) or by combining several devices of different impedances into one. Usually such devices must be sufficiently large to operate near the peak of the sea spectrum, and very small bodies can be resonated only at frequencies above the usual range of the energetic sea.

Recently Fred Olsen and ABB Power Systems Inc. (http://www02.abb.com) in Norway have designed a system called FO3 which consists of a rig with many small floating cylinders hanging underneath it. Energy is absorbed from the waves as they set the cylinders into vertical motion which then activates a hydraulic system driving a generator to produce electricity. Currently being tested is a 1:3-scale research model which measures 12 by 12 meters and is 8 meters high. It is estimated that the full-scale model can produce 2.52 MW from 6-meter high waves with a period of 9 seconds, comparable to the capacity of a wind turbine. Eventually a large array of many rigs can be installed over a large sea surface area and produce much more electricity.

In this paper we report a theoretical evaluation of this novel concept by examining a compact array of small buoys with spacings much shorter that the typical wavelength. Based on linearised theory of small amplitude waves we first employ the method of homogenization (i.e. multiple scales) to derive effective equations governing the dynamics on the macro-scale of the wave length. We show that for buoys of dimensions and spacing small compared to the water depth and wavelength, their presence and motion are manifest in a modified free surface condition on the wavelength scale. Explicit results are obtained for a long array of finite width attacked by normal incident sea, as well as a circular array of large radius. The dynamics and the energy efficiency are then compared with those of single buoys.

Specifically we shall consider a square array of small and identical buoys floating on the surface of the sea of constant mean depth h * , as shown in figure 1. Each buoy is a vertical cylinder of circular cross section of radius a * and draft H * , spaced at the distance d * from centre to centre. Assuming monochromatic waves of frequency ω * , the wave number k * of the incident waves is given by the real root of the dispersion relation

ω * 2 = gk * tanh(k * h * ) (1.1)
The incoming wave length and the sea depth are assumed to be comparable but both are much greater than the buoy radius a * , the draft H * and the separation distance d * . i.e.,

a * h * ≡ µ ≪ 1, O(a * ) = O(H * ) = O(d * ), k * h * = O(1) (1.2)
Wave energy is extracted from the heaving oscillation of each buoy through an absorbing device anchored to the seabed or attached to a fixed supporting structure.

Linearised governing equations

We employ the following symbols for physical domains: Ω F is the fluid domain, S F is the free surface, S W is the lateral surface of the buoys, and S B is the bottom surface of the buoys. Let us denote all physical variables with asterisks. Assuming irrotational flow and infinitesimal waves, the velocity potential in water is governed by Laplace's equation

∆ * Φ * = ∂ 2 Φ * ∂x * 2 + ∂ 2 Φ * ∂y * 2 + ∂ 2 Φ * ∂z * 2 = 0, x * ∈ Ω F (2.1)
The total pressure inside water is given by Bernoulli's equation:

p * = -ρ ∂Φ * ∂t * -ρgz * (2.2)
On the free surface z * = η * (x * , y * , t * ), the kinematic boundary condition is

∂Φ * ∂z * = ∂η * ∂t * , x * ∈ S F (2.3)
and the dynamic boundary condition is,

gη * + ∂Φ * ∂t * = 0, x * ∈ S F (2.4)
As the sea-surface pressure is assumed to be constant. On the sea bed the vertical velocity vanishes, which gives

∂Φ * ∂z * = 0, z * = -h * (2.5)
On the side wall of the buoy, there is no normal velocity:

∂Φ * ∂r * = 0, x * ∈ S W (2.6)
where r * is the local radial coordinate from the axis of a cylindrical buoy. We assume that all buoys are installed on a large stationary frame or platform which is held fixed above the sea surface. On the flat bottom of the buoy the kinematic condition is

∂Φ * ∂z * = ∂ζ * ∂t * , x * ∈ S B (2.7)
where ζ * (t * ) is the unknown vertical displacement of the buoy. Modelling the energy extraction device as a linear load force

-λ * ∂ζ * ∂t * (2.8)
on a moving buoy with a constant coefficient λ * , the conservation law of vertical momentum of the buoy serves as the dynamic condition

M * ∂ 2 ζ * ∂t * 2 + λ * ∂ζ * ∂t * + πa * 2 ρgζ * = -ρ SB ∂Φ * ∂t * dS * (2.9)
where M * = ρπa * 2 H * is the buoy mass and H * its draft by Archimedis principle.

Let us introduce normalized variables as follows,

x * i = a * x ′ i , t * = t ′ h * g , Φ * = A * gh * Φ, η * = A * η, ζ * = A * ζ (2.10)
where A * is the amplitude of the incoming wave. Let us rewrite the governing equations.

Note that the length scale is the small radius of the buoy (the micro-scale). In normalized form, Laplace's equation (2.1) is unchanged. The free surface condition (2.3) becomes

µ ∂η ∂t ′ = ∂Φ ∂z ′ , x ′ ∈ S F (2.11)
where

µ = a * h * ≪ 1 (2.12)
is the key parameter in this study. Equation (2.4) gives,

η + ∂Φ ∂t ′ = 0, x ′ ∈ S F (2.13)
Equations (2.11) and (2.13) can be combined into

∂Φ ∂z ′ + µ ∂ 2 Φ ∂t ′ 2 = 0, x ′ ∈ S F (2.14)
and the condition on the seabed now reads

∂Φ ∂z ′ = 0, z ′ = - h * a * = - 1 µ (2.15)
As 1/µ ≫ 1, this micro-scale boundary condition is effectively applied at z ′ → -∞. On the buoy, the kinematic conditions are

∂Φ ∂r ′ = 0 ′ , x ′ ∈ S W (2.16) and ∂Φ ∂z ′ = µ ∂ζ ∂t ′ , x ′ ∈ S B (2.17)
The dynamic condition (2.9) now reads

a * h * H * a * ∂ 2 ζ ∂t ′2 + λ * g/h * ρ * gπa * 2 ∂ζ ∂t ′ + ζ = - SB ∂Φ ∂t ′ dS ′ π (2.18) Defining H ′ = H * a * = O(1), λ = λ * g/h * ρ * gπa * 2 = O(1) (2.19)
we change the same dynamic condition to dimensionless form,

µH ′ ∂ 2 ζ ∂t ′ 2 + λ ∂ζ ∂t ′ + ζ = - 1 π SB ∂Φ ∂t ′ dS ′ (2.20)
which can be combined with the kinematic condition (2.17) to give,

µH ′ ∂ 2 ∂t ′2 + λ ∂ ∂t ′ + 1 ∂Φ ∂z ′ = - µ π SB ∂ 2 Φ ∂t ′ 2 dS ′ (2.21)

Multiple-scale approximation

Our main objective is to consider the collective effects of many small buoys on the dynamics over a much larger region of dimensions comparable to the sea depth or to the wavelength. In view of the contrast of scales we seek an asymptotic approximation by the method of multiple scales, and define the slow (macro-scale) coordinates without primes by

x = µx ′ (3.1)
Let us denote by ∇ ′ and ∆ ′ the gradient and Laplacian on the micro-scale and ∇ and ∆ the corresponding operators on the macro-scale. We next introduce the expansions

Φ = e -iωt ′ φ 0 (x ′ , x) + µφ 1 (x ′ , x) + µ 2 φ 2 (x ′ , x) + . . . (3.2) η = e -iωt ′ η 0 (x ′ , y ′ , x, y) + µη 1 (x ′ , y ′ , x, y) + µ 2 η 2 (x ′ , y ′ , x, y) + . . . (3.3) ζ = e -iωt ′ ζ 0 (x ′ , y ′ , x, y) + µη 1 (x ′ , y ′ , x, y) + µ 2 ζ 2 (x ′ , y ′ , x, y) + . . . (3.4)
where ω is the dimensionless frequency normalized according to

ω = ω * h * g (3.5)
Referring to the dimensionless governing equations in §1, we get from Laplace's equation,

∆ ′ + 2µ∇ ′ • ∇ + µ 2 ∆ (φ 0 + µφ 1 + ...) = 0, x ′ ∈ Ω F
(3.6) The combined free surface condition becomes

∂ ∂z ′ + µ ∂ ∂z -µω 2 (φ 0 + µφ 1 + ...) = 0, x ′ ∈ S F (3.7)
while the kinematic condition is

∂ ∂z ′ + µ ∂ ∂z + • • • (φ 0 + µφ 1 + µ 2 φ 2 + • • • ) = -iµω(η 0 + µη 1 + • • • ), x ′ ∈ S F (3.8)
On the side wall of the buoy we have

∂ ∂r ′ + µ ∂ ∂r (φ 0 + µφ 1 + ...) = 0, x ′ ∈ S W (3.9)
and on the sea bed

∂ ∂z ′ + µ ∂ ∂z (φ 0 + µφ 1 + ...) = 0, z ′ = - 1 µ (3.10)
At the bottom of the buoy, the kinematic condition (2.17) gives

∂ ∂z ′ + µ ∂ ∂z (φ 0 + µφ 1 + φ 2 + • • • ) = -iµω(ζ 0 + µζ 1 + • • • ), z ′ = -H ′ (3.11)
while the dynamic condition gives

-µω 2 H ′ -iλω + 1 ζ 0 + µζ 1 + µ 2 ζ 2 + • • • = iω π SB (φ 0 + µφ 1 + µ 2 φ 2 + • • • ) dS ′ (3.12)
From the combined buoy condition, we get

-µω 2 H ′ -iλω + 1 ∂φ 0 ∂z ′ + µ ∂φ 0 ∂z + ∂φ 1 ∂z ′ + µ 2 ∂φ 1 ∂z + ∂φ 2 ∂z ′ + • • • -µ ω 2 π SB (φ 0 + µφ 1 + µ 2 φ 2 + • • • ) dS ′ = 0 (3.13)
Let us also expand the factor

1 1 -iλω -µH ′ ω 2 = j=0 µ j F j (ω) (3.14)
where

F 0 (ω) = 1 1 -iλω , F 1 (ω) = H ′ ω 2 (1 -iλω) 2 , etc. (3.15)
By separating the orders, a series of micro-scale boundary value problems are then obtained at the orders O(1), O(µ) and O(µ 2 ).

Leading order (O(1))

The governing conditions are homogeneous

∆ ′ φ 0 = 0, x ′ ∈ Ω (3.16a) ∂φ 0 ∂n ′ = 0, x ′ ∈ S F ∪ S W ∪ S B ∪ S b (3.16b)
where S b denotes the sea bed at z ′ = -µ -1 ≪ -1. Let us define a unit cell of the array as shown in figure 1. Because there are a large number of periods in the array, we impose the condition that on the micro-scale, the solution is periodic, i.e.

φ 0 (x ′ , y ′ , z ′ , x) = φ 0 (x ′ + d ′ , y ′ , z ′ , x) (3.17a) φ 0 (x ′ , y ′ , z ′ , x) = φ 0 (x ′ , y ′ + d ′ , z ′ , x) (3.17b) with d ′ ≡ d * /a * (3.18)
being the centre-to-centre distance between adjacent buoys. The leading-order solution is clearly independent of the micro-scale,

φ 0 = φ 0 (x) (3.19)
and the dependence on the macro-scale is yet to be found. It follows from (2.4) that

η 0 = iω φ 0 | z=0 (3.20)
independently of the presence of the buoys. In the buoy-area (3.12) gives the buoy displacement

ζ 0 = iωF 0 (ω) φ 0 | z=0 , x ∈ SB . (3.21)
Both η 0 and ζ 0 are independent of the micro-scale coordinates, and they are related by

ζ 0 = F 0 η 0 x ′ ∈ SB . (3.22)
inside the buoy area.

First order (O(µ))

Using (3.19), we get from (3.6) that

∆ ′ φ 1 = 0, x ′ ∈ Ω F (3.23a)
and from (3.7) that

∂φ 1 ∂z ′ = - ∂φ 0 ∂z -ω 2 φ 0 , x ′ ∈ S F (3.23b) Equation (3.13) becomes ∂φ 1 ∂z ′ = - ∂φ 0 ∂z -ω 2 F 0 φ 0 , x ′ ∈ S B (3.23c)
We also have on the side wall of the buoy,

∂φ 1 ∂r ′ = - ∂φ 0 ∂r = -n i ∂φ 0 ∂x i , x ′ ∈ S W (3.23d )
where n = (n 1 , n 2 ) = (cos θ, sin θ) denotes the unit vector normal to the side wall, and

∂φ 1 ∂z ′ = - ∂φ 0 ∂z , z ′ = - 1 µ (3.23e)
on the seabed. In addition we require micro-scale periodicity on the cell boundaries.

Once φ 1 is found, ζ 1 follows from (3.12). As it is usual in homogenization analysis, the macro-scale physics at the leading order is found by requiring the solvability of the inhomogeneous micro-scale problem at a higher order. The micro-scale cell problem for φ 1 is inhomogeneous. By applying Gauss' theorem (or, equivalently applying Green's formula to φ 0 and φ 1 over a unit cell) φ 1 over the cell volume, we get

∂Ω ∂φ 1 ∂n ′ dS ′ = 0 (3.24)
where ∂Ω is the boundary of the cell. This is just the condition of solvability for the inhomogeneous problem of φ 1 . Since

SW ∂φ 1 ∂r ′ dS ′ = -µ SW ∂φ 0 ∂r dS ′ = -µ∇φ 0 • SW e r dS ′ = 0 (3.25) we must have SB ∂φ 1 ∂z ′ dS ′ = - SF ∂φ 1
∂z ′ dS ′ which gives at the leading order:

(1 -f ) ∂φ 0 ∂z -ω 2 φ 0 + f ∂φ 0 ∂z -ω 2 F 0 φ 0 = 0, z = 0 (3.26)
where for circular buoys

f ≡ πa * 2 d * 2 = π d ′2 , with 0 < f < π 4 (3.27)
is the area fraction of solid, or the packing ratio. Hence we have

∂φ 0 ∂z -ω 2 [1 + f (F 0 -1)] φ 0 = 0, z = 0, x ∈ SB . (3.28)
This is a key result of our approximation and gives the macro-scale boundary condition over the part of the mean sea surface covered by buoys. In the open water with no buoy, f = 0, (3.28) reduces to the familiar condition on the free surface:

∂φ 0 ∂z -ω 2 φ 0 = 0, x ∈ S F (3.29)
Because of (3.28), (3.23b) and (3.23c) can be rewritten as

∂φ 1 ∂z ′ = - ∂φ 0 ∂z 1 - 1 1 -f (F 0 -1) ≡ -β ∂φ 0 ∂z , x ∈ S F (3.30) and ∂φ 1 ∂z ′ = - ∂φ 0 ∂z 1 - F 0 1 + f (F 0 -1) ≡ -β ′ ∂φ 0 ∂z , x ∈ S B (3.31)
which define β and β ′ . In view of the forms of the boundary conditions, the solution of the micro-scale problem for φ 1 in a unit cell can be sought in the form

φ 1 (x ′ , x) = - 3 j=1 N j (x ′ ) ∂φ 0 ∂x j (3.32)
Then the horizontal components N 1 , N 2 are governed by the following boundary-value problems in the unit cell,

∆ ′ N j = 0, x ′ ∈Ω F (3.33a) ∂N j ∂z ′ = 0, x ′ ∈S F ∪ S B ∪ S b (3.33b) ∂N j ∂r ′ = n j , x ′ ∈S W (3.33c)
where the outerward normal to

S W is n = (n 1 , n 2 , 0). The vertical component N 3 is governed instead by ∆ ′ N 3 = 0, x ′ ∈Ω F (3.34a) ∂N 3 ∂z ′ = β, x ′ ∈S F (3.34b) ∂N 3 ∂z ′ = 0, x ′ ∈S b (3.34c) ∂N 3 ∂z ′ = β ′ , x ′ ∈S B (3.34d ) ∂N 3 ∂r ′ = 0, x ′ ∈S W (3.34e)
The solutions are made unique by adding the constraint

ΩF N j (x) dV ′ = 0, j = 1, 2 (3.35)
and

N 3 = 0 at a point x ′ = x ′ b = (0, 0, -µ -1
) on the seabed. Being periodic in (x ′ , y ′ ), the harmonic functions N j (x ′ ) are expected to diminish exponentially in z ′ . For confirmation we have performed a numerical simulation using Finite Elements. The results, given in figure 14 in Appendix A, show indeed that for a sufficiently slender cell the solutions N i are highly localized near the buoy. In view of (3.32), a consequence is that,

∂φ 1 ∂z ′ → 0 as z ′ → -∞ (3.36)
which in turn implies

∂φ 0 ∂z = 0, z = -1 (3.37)
because of (3.23e). This provides the seabed boundary condition for the macro-scale problem.

Second order (O(µ 2 )) and the macro-scale problem

At the second order the micro-scale problem for φ 2 is again inhomogeneous,

∆ ′ φ 2 = -2∇ ′ • ∇φ 1 -∆φ 0 , x ′ ∈ Ω F (3.38a) ∂φ 2 ∂z ′ = - ∂φ 1 ∂z + ω 2 φ 1 , x ′ ∈ S F (3.38b) ∂φ 2 ∂z ′ = - ∂φ 1 ∂z -ω 2 F 0 φ 1 -ω 2 F 1 φ 0 , x ′ ∈ S B (3.38c) ∂φ 2 ∂z ′ = - ∂φ 1 ∂z , z ′ = -µ -1 (3.38d ) ∂φ 2 ∂r ′ = - ∂φ 1 ∂r , x ′ ∈ S W (3.38e) As |N i | → 0 for z ′ → -µ -1 , (3.38d) reduces to : ∂φ 2 ∂z ′ = 0, z ′ = -µ -1 (3.39)
We now apply Green's formula for φ 0 and φ 2 in the unit cell and invoke their governing conditions on the micro-scale to get

ΩF (∆φ 0 + 2∇ ′ • ∇φ 1 ) dV ′ = SF ∂φ 1 ∂z -ω 2 φ 1 dS ′ + SB ∂φ 1 ∂z -ω 2 F 0 φ 1 -ω 2 F 1 φ 0 dS ′ - SW ∂φ 1 ∂r dS ′ (3.40)
Using the fact that φ 1 vanishes with N i outside the vertical distance of O(1) from z = 0, and that the cell volume

|Ω F | = O(1/µ
) is much greater than unity, we conclude that:

ΩF ∆φ 0 dV ′ = 0 Because φ 0 (x, t) is independent of x ′ , we conclude that ∆φ 0 = 0, -1 < z < 0 (3.41)
Thus φ 0 is harmonic on the macro-scale. In summary, in the region with buoys, the macro-scale variation of φ 0 (x) is governed by (3.41) in the fluid region, subject to the boundary condition (3.28) on z = 0 in the buoy-covered area, and (3.37) on the seabed. In the open water without buoys, condition (3.28) must be replaced by (3.29), while (3.41) and (3.37) still apply. Note that due to the small draft H ′ , buoy inertia, hence resonance, is unimportant.

The homogenization analysis for finding the macro-scale behaviour can in principle be extended to periodic buoys of any shape. Once the macro-scale is completely determined, one can also derive the micro-scale fluctuations by solving the cell problems for the vector N (x ′ ). Then φ 1 (x, x ′ ) can be found according to (3.32) and used to calculate wave forces on each buoy hence the individual apparent mass and radiation damping matrices. Such effort is needed for design, but is omitted here.

We shall now apply these macro-scale equations to examine wave power extraction from one-and two-dimensional arrays in response to a plane incident wave train arriving from x ∼ -∞.

Vertical eigenfunctions

As it is well known, the general solution in the open water region where f = 0 can be expressed as a series of the form

φ 0 (x) = ∞ n=0 ψ n (x, y)f n (z) (4.1)
where

f 0 = c 0 cosh(k 0 (z + 1)), f n = c n cos(κ n (z + 1)) (4.2)
are real orthogonal eigenfunctions in -1 < z < 0, and (k 0 , k 1 , . . . ) are the eigenvalues of the dispersion relation,

ω 2 = k n tanh(k n ), n = 0, 1, . . . (4.
3) In particular k 0 is the positive real root and k n ≡ iκ n is the n-th imaginary root i.e.,

ω 2 = k 0 tanh(k 0 ), ω 2 = -κ n tan(κ n ), n = 1, 2, 3, . . . (4.4)
With the choice of

c 0 = 2 1 + ω -2 sinh 2 k 0 , c n = 2 1 -ω -2 sin 2 (κ n ) (4.5)
the vertical eigenfunctions are orthonormal,

f n |f m ≡ 0 -1 f n (z)f m (z) dz = δ nm (4.6)
Furthermore the horizontal factors ψ n must satisfy Helmholtz equations in the horizontal plane

∂ 2 ∂x 2 + ∂ 2 ∂y 2 + k 2 0 ψ 0 = 0, (4.7a) ∂ 2 ∂x 2 + ∂ 2 ∂y 2 -κ 2 n ψ n = 0, n = 1, 2, 3, ... (4.7b)
In the region of wave absorbing buoys we also assume

φ 0 (x) = ∞ n=0 Ψ n (x, y)F n (z) (4.8)
It can be shown that the eigenfunctions {F n }, n = 0, 1, 2, . . . are the solutions of the boundary value problem

F ′′ n (z) -K 2 n F n (z) = 0 -1 < z < 0 F ′ n -σ 2 F n = 0 z = 0 F ′ n = 0, z = -1, (4.9)
where σ is defined by

σ 2 ≡ ω 2 [f F 0 (ω) + (1 -f )]
(4.10) and is complex due to energy extraction. Therefore the eigenfunctions F n are complex

F n = C n cosh K n (z + 1) (4.11)
The eigenvalue K n is the n-th complex root of the relation This type of dispersion relation with complex σ 2 arises also for waves through a porous media and has been studied by [START_REF] Dalrymple | Reflection and transmission from porous structures under oblique wave attack[END_REF] and [START_REF] Mciver | The dispersion relation and eigenvalue expansion for water waves in porous structures[END_REF].

σ 2 = K n tanh K n (4.12) -1 -0.8 -0.6 -0.4 -0.2 0 -1.5 -1 -0.
It is straightforward to show that the set {F n } is orthogonal. By choosing the coefficients {C n } to be

C n = 2 σ -2 sinh 2 (K n ) + 1 (4.13)
the eigenfunctions {F n } are also orthonormal,

F n |F m ≡ 0 -1 F n (z)F m (z)dz = δ nm (4.14)
Since K n is complex, the square root above is defined such that if the complex radical is z = re iθ , its phase is limited to the range -π < θ ≤ π. [START_REF] Dalrymple | Reflection and transmission from porous structures under oblique wave attack[END_REF]; [START_REF] Mciver | The dispersion relation and eigenvalue expansion for water waves in porous structures[END_REF] showed that this set of functions is a complete basis provided the eigenvalues K n are distinct, which is in general the case. For a given frequency ω, packing ratio f and damping rate λ, σ is first defined. K n and F n are found numerically. Before employing an usual iterative algorithm to solve the complex transcendental equation, a good initial guess of the solution is needed. For this purpose we solved the eigenvalue problem governed by (4.9) by the numerical method of finite elements with a regular mesh and third-order Laplace elements. The resulting K n 's are used as a initial guesses for further iteration of (4.12). Sample F n 's are shown in figure 2. Sample eigenvalues are given in Table 1.

Note that for f = 0, k is purely real and k n = iκ n , n = 1, 2, 3... are purely imaginary. For f ≪ 1, K 0 is almost real and K n , n = 1, 2, 3, ... are almost imaginary. Perturbation solutions of (4.12) have been used to confirm the values in Table 1 where f = 0.2.

We now apply these results to examine two simple arrays.

A long array of energy-absorbing buoys

Referring to figure 3, let us first consider a long array of width L with its edges parallel to the crests of incoming plane waves. Assuming an incoming wave of unit amplitude, the velocity potential in the open water on the incidence side (zone I) is where k 0 is real and k n = iκ n , n = 1, 2, 3, ... are imaginary roots of the dispersion relation.

φ I (x, z) = -i ωf 0 (0) e ik0x f 0 (z) + ∞ n=0 R n e -iknx f n (z) , -∞ < x < 0 (5.1) Kn ω = 0.5 ω = 1 ω = 2 K1 0.
In zone II of the buoys, the potential is

φ II (x, z) = -i ωf 0 (0) ∞ n=0 B n e iKnx + B ′ n e -iKnx F n (z), 0 < x < L (5.2)
and in the open water on the transmission side (zone III) we have

φ III (x, z) = -i ωf 0 (0) ∞ n=0 T n e iknx f n (z), L < x < ∞ (5.3)
The eigenvalues (k n , K n ) and eigenfunctions (f n , F n ) have been defined in §4. Let us introduce

U (z) = ∂φ 0 ∂x (0, z), U ′ (z) = ∂φ 0 ∂x (L, z) (5.4)
as the horizontal velocities at x = 0 and x = L respectively. Requiring flux continuity and using the orthogonality of eigenfunctions, we find

R 0 = 1 - U |f 0 ik 0 R n = - U |f n ik n (5.5a) B n = - (U ′ -e -iKnL U )|F n 2K n sin(K n L) B ′ n = - (U ′ -e iKnL U )|F n 2K n sin(K n L) (5.5b) T n = U ′ |f n ik n e iknL
(5.5c)

We further require continuity of pressure (i.e., of potentials) at x = 0

f 0 (z) + 1 - U |f 0 ik 0 f 0 (z) - n≥1 U |f n ik n f n (z) = - n≥0 (U ′ -cos(K n L)U )|F n K n sin(K n L) F n (z) (5.6) and at x = L - n≥0 (U ′ cos(K n L) -U )|F n K n sin(K n L) F n (z) = n≥0 U ′ |f n ik n f n (z) (5.7)
These are a pair of integral equations for U (z) and U ′ (z) in -1 < z < 0. Let their solutions be represented by the following orthonormal expansions

U = m U m F m , U ′ = m U ′ m F m , -1 < z < 0 (5.8)
with unknown coefficients, and let

f n = m M nm F m where f n |F m = M nm
(5.9)

The matrix elements M nm can be obtained explicitly,

M nm = ω 2 f (1 -F) c n cosh(k n )C m cosh(K m ) (k 2 n -K 2 m )
(5.10) Equations (5.6) and (5.7) become 2f 0 (z)

n,q

M nq U q 1 ik n f n (z) = n - U ′ n K n sin(K n L) + U n K n tan(K n L) F n (z) (5.11) n,q M nq U ′ q 1 ik n f n (z) = n - U ′ n K n tan(K n L) + U n K n sin(K n L)
F n (z) (5.12)

By taking the scalar product with F p for p = 0, 1, 2, 3, ... in turn we obtain from (5.11) and (5.12)

2M 0p - n,q M nq U q 1 ik n M np = - 1 K p sin(K p L) p U ′ n + 1 K p tan(K p L) U p (5.13)
and n,q

M nq U ′ q 1 ik n M np = - 1 K p tan(K p L) U ′ p + 1 K p sin(K p L) U p (5.14)
The expansion coefficients U n , U ′ q are solved numerically after truncation. Afterwards we get the buoy displacement ζ 0 from the expression of φ II . The transmission and reflection coefficients follow from (5.5c) and (5.5a) : The dimensionless power-extraction efficiency is

T ≡ T 0 = M 0q U ′ q ik 0 e ik0L (
E = 1 -|T | 2 -|R| 2 .
(5.17)

Figure 4 shows the amplitude and phase of the free surface elevation inside the buoy region for array width of L = 1 and L = 5 according the macro-scale normalization. The buoy displacement is simply proportional to that of the free surface displacement in the same region by the complex reduction factor F 0 whose magnitude

|F 0 | = 1 1 + (λω) 2
(5.18) is smaller for higher extraction rate and frequency. Note first that there is no resonance. b) is due to interference by strong reflection, similar to the case of a finite shelf (cf. [START_REF] Mei | Theory and application of ocean surface waves[END_REF], p. 149). In the transmission coefficient, this oscillatory behaviour is less prominent due to energy extraction. We have indeed checked that in the limit of extremely strong load force, λ ≫ 1, the buoys no longer move. In this case reflection is the strongest and the oscillatory variation in T is recovered.

While it is not surprising that a larger L gives a higher efficiency, as shown in figure 7(b), it is nevertheless interesting that the gain of energy extraction with a wider array is more significant at low frequency. In practical situations k 0 = k * 0 h * will likely be between 0 and 3. Our predictions can help the designer to choose the proper width by considering both efficiency and construction economy.

In general scattering is significant, hence the maximum efficiency of energy extraction is somewhat lower than that a large beam-sea device such as a Salter's duck (see [START_REF] Mynett | Characteristics of Salter's cam for extracting energy from ocean waves[END_REF]).

A circular array

The solution

Now let many buoys be gathered inside a circular area of radius R. First, it is well known that the incident plane wave in the direction of x can be expanded as a sum of partial waves (see e.g. [START_REF] Abramowitz | Handbook of mathematical functions[END_REF])

φ i (x) = -i ωf 0 (0) f 0 (z)e ik0x = -i ωf 0 (0) f 0 (z) ∞ m=0 ε m i m J m (k 0 r) cos(mθ)
where ε 0 = 1 and ε n = 2 for n = 1, 2, 3, . . . are the Jacobi symbols. Let us express the total solution as

φ = ∞ m=0 φm (r, z) cos(mθ) (6.1)
In the open water, the m-th mode potential φm can be written as φm = -i ωf 0 (0)

ε m i m J m (k 0 r)f 0 (z) + ∞ n=0 a n,m ψ n,m (r)f n (z) , r > R (6.2) with ψ n,m (r) = H (1) m (k 0 r) for n = 0 K m (K n r) = H (1) m (iK n r) for n = 1, 2, . . . (6.3)
where H

(1) m is the first Hankel function of order m. The first term in (6.2) corresponds to the incident wave and the series to the scattered/radiated waves. In the circular region of buoys, 0 < r < R, we can expand the potential as:

φm = -i ωf 0 (0) ∞ n=0 b n,m Ψ n,m (r)F n (z), 0 < r < R, (6.4) with Ψ n,m (r) = J m (K n r)
where (f n , k n ) and (F n , K n ) are the same as before. Let us denote the common radial flux along r = R by

U m (θ, z) = ∂ φm ∂r , r = R. (6.5)
The expansion coefficients are found in terms of U m by orthogonality:

a 0,m = U m |f 0 -ε m i m k 0 J ′ m (k 0 R) ψ ′ 0,m (R) (6.6) a n,m = U m |f n ψ ′ n,m (R) (6.7) b n,m = U m |F n Ψ ′ n,m (R) (6.8)
which ensures the continuity of radial flux. Continuity of pressure (i.e. potential) at r = R requires that

ε m i m J m (k 0 R) - k 0 J ′ m (k 0 R) ψ ′ 0,m (R) ψ 0,m (R) f 0 (z) + ∞ n=0 (U m |f n ) ψ ′ n,m (R) ψ n,m (R)f n (z) = ∞ n=0 U m |φ n Ψ ′ n,m (R) Ψ n,m (R)F n (z) (6.9)
Introducing the expansions

f i = j M ij F j , U m = j U j,m F j (6.10) We get J m (k 0 R) - k 0 J ′ m (k 0 R) ψ ′ 0,m (R) ψ 0,m (R) j M 0j F j (z) + i,j,k ψ k,m (R) ψ ′ k,m (R) M ki U i,m M kj F j (z) = n Ψ n,m (R) Ψ ′ n,m (R)
U n,m F n (z) (6.11) for m = 0, 1, 2, . . . . By taking the scalar product with F p , we finally obtain for any value of m

j k M kp ψ k,m (R) ψ ′ k,m (R) M kj - Ψ p,m (R) Ψ ′ p,m (R) δ p,j U j,m = -ε m i m J m (k 0 R) - k 0 J ′ m (k 0 R) ψ ′ 0,m (R) ψ 0m (R) M 0p (6.12)
which is a matrix equation for the unknown vector U j,m for every m. Numerical computations can be carried out after truncation of the series. After solving for U j,m , the velocity U m hence φ m are found. Combining (6.2) and (3.21), we get the displacement of the buoys

ζ 0 (r, θ) = F 0 (ω) ∞ m=0 ∞ n=0 b n,m Ψ n,m (r) F n (0) f 0 (0) cos(mθ) (6.13)
Again it is proportional to the free surface displacement in the same area according to (3.22). Hence we only show in figure 9 the free surface displacement in and outside the buoy area, for two arrays of radii R = 1 and R = 5. For the smaller array the displacement is relatively uniform and less than 1. For the larger array, the displacement is significantly reduced on the leeward side.

Energy absorption

One can evaluate the extracted energy by calculating the total energy flux into a large circular cylindrical surface of radius r ≫ R. In physical variables the power output is

P * = 2π 0 0 -h * Re(iρω * φ * ) ∂Re(φ * ) ∂r * r * dz * dθ = ρ g h * A * 2 gh * h * 2π 0 0 -1 1 2 Re iφ ∂φ † ∂r r dz dθ (6.14)
where the overline denotes time averaging over a period and dagger indicates complex conjugate. Use has been made of the normalization defined in (2.10). The implied normalization for power output is

P * = P ρ g h * A * 2 gh * 2 (6.15)
In contrast, the power flux per unit length of the incoming wave crest is,

1 2k * 0 ρgA * 2 C * g = 1 2 h * ρgA * 2 g h * h * 1 k 0 dω dk 0 = ρgA * 2 g h * h * 2 C g 2k 0
where C g = C * g / √ gh * is the dimensionless group velocity of the incoming plane wave. As in early theories the capture width W * can be defined as the ratio of the absorption rate to the influx rate of wave power within unit length of the incoming wave front. As a measure of effectiveness, k * 0 W * represents the fraction of a wavelength where the incoming power is depleted,

k * 0 W * = k 0 W = P * 1 2k * 0 ρgA * 2 C * g = 2k 0 P C g (6.16)
Using the asymptotic expansions of Bessel functions for large k 0 r, we get from (6.2), φ m (r) ≈ A m 2 πk 0 r e i(k0r-π/4) -i ωf 0 (0) f 0 (z) for k 0 r ≫ 1 (6.17) so that φ(r) ≈ -i ωf 0 (0) e ik0r cos(θ) + m A m 2 πk 0 r e i(k0r-π/4) cos(mθ) f 0 (z) (6.18)

where the modal amplitudes A m can be computed from the solution using the asymptotic expression of the Hankel functions:

A m = a 0,m i -m
Using the method of stationary phase it can be shown that :

P = 2 1 ωf 0 (0) 2   |A 0 | 2 + 1 2 m≥1 |A m | 2 + Re   m≥0 A m     (6.19)
Details are similar to that in [START_REF] Mei | Theory and application of ocean surface waves[END_REF], p.381, and omitted. The capture width is therefore:

k * 0 W * = k 0 W = 4k 0 ωC g f 0 (0) 2   |A 0 | 2 + 1 2 m≥1 |A m | 2 + Re   m≥0 A m    
Finally, using the expression for f 0 and the dispersion relation, we find

C g f 0 (0) 2 = k 0 ω hence k 0 W = 4   |A 0 | 2 + 1 2 m≥1 |A m | 2 + Re( m≥0 A m )   (6.20)
The same result can also be derived by calculating the rate of work done on the heaving buoys. Another measure of effectiveness is the ratio of the extraction rate to the influx rate across the entire diameter of the array, W/2R which is expected to be less than unity.

For evaluating the merits of the compact array let us first recall some results known for a single buoy : (i) The optimal k 0 W is 1 at best for a heaving buoy of any size. If all three degrees of freedom are used to extract energy then maximum k 0 W = 3 [START_REF] Newman | Absorption of wave energy by elongated bodies[END_REF], [START_REF] Falnes | Ocean Waves and Oscillating Systems[END_REF], [START_REF] Mei | Theory and application of ocean surface waves[END_REF]). (ii) The peak value of k 0 W occurs at k 0 a b = O(1). In other words, the peak occurs at higher k 0 = k * h * for smaller a b = a * b /h * . (iii) The curve of k 0 W versus k 0 = k * h * has a broader peak for a smaller a b . Properties (ii) and (iii) are based on numerical computations via the eigenfunction expansion method of [START_REF] Black | Radiation and scattering of water waves by rigid bodies[END_REF] and are confirmed by approximate reasoning in Appendix B.

In light of these let us present the results for a circular array of buoys. Figure 10 shows the dependence of the two measures of effectiveness on the extraction rate λ. For two different array radii R, the greatest k 0 W and W/2R are achieved at around the same extraction rate of λ = 0.5. The optimal rate of extraction depends slightly on the frequency/wave-number of the incoming wave.

Figure 11 shows that for a fixed packing ratio and damping rate, the capture width kW and efficiency W/2R naturally increase with the radius of the array. More important, the bandwidth of both quantities is very large for all array sizes.

Figure 12 shows that the capture width increases monotonically with the packing ratio f , and with the incoming wave frequency. Recall that for circular buoys in a square array the maximum packing ratio is f ≤ π/4 ≈ 0.8.

Finally let us compare a large buoy whose radius and draft are equal, with a buoy array of the same total displaced volume πf R 2 H where H ≡ H * /h * and is taken to be H = 0.1 for illustration †. Then the radius and draft of the large buoy are both a b = (f R 2 H) 1/3 . Figure 13 compares the capture widths over a wide range of frequencies. The solid curves gives the capture width for an array for different radii R, with fixed f = 0.2 and λ = 0.5. The dashed curves represent the capture width for a single-buoy absorber of radius a b . In the range of 0 < k 0 (= k * 0 h * ) < 6 the maximum k 0 W is at most unity for a single heaving buoy, and can be 3 if roll and sway can also be resonated. Note however that the band width of a single buoy is always much narrower. Thus the circular buoy array is potentially more advantageous from the technical viewpoint of efficiency.

Conclusions

Stimulated by a recent invention in Norway, we have developed a theory for the hydrodynamics and power-extraction efficiency of a compact array of small buoys. The typical wave length is assumed to be comparable to the overall radius of the array but much greater than the dimensions of individual buoys. For a periodic array the two-scale method of homogenization leads to an effective equation governing the spatial average. The energy-absorbing efficiency is studied for a long strip of buoys and for a circular array. The latter geometry is shown to be potentially advantageous, having good efficiency over a broad range of frequencies, unlike that of one large buoy. The theory can be readily modified for wave interaction with broken ice floes on the sea surface, if the ice floes are idealized as identical floating bodies in a periodic array.

Finally, we stress that the homogenization theory employed here is effective only when two sharply different scales exist. When both the buoy dimension and the spacing are not small compared to the wavelength, direct numerical methods are available but require greater computational effort. For simple geometries such as vertical circular cylinders, formally exact theories have been reported by [START_REF] Linton | The interaction of waves with arrays of vertical cylinders[END_REF]; [START_REF] Manihar | Wave diffraction by a long array of cylinders[END_REF]; [START_REF] Linton | The scattering of water waves by an array of circular cylinders in a channel[END_REF]; [START_REF] Chamberlain | Water wave scattering by finite arrays of circular structures[END_REF]; [START_REF] Linton | The radiation and scattering of surface waves by a vertical circular cylinder in a channel[END_REF]; [START_REF] Linton | The scattering of water waves by an array of circular cylinders in a channel[END_REF] for an infinite or semi-infinite line of fixed vertical cylinders with finite radius. These methods still call for significant numerical work. Approximate theories on the interaction of water waves with many floating objects have been given by [START_REF] Falnes | Radiation impedance matrix and optimum power absorption for interacting oscillators in surface waves[END_REF]; [START_REF] Falnes | Wave-power absorption by parallel rows of interacting oscillating bodies[END_REF]; [START_REF] Falnes | Wave-power absorption by an array of attenuators oscillating with unconstrained amplitudes[END_REF][START_REF] De | Wave-power absorption by a periodic linear array of oscillating water columns[END_REF] for large separation and weak hydrodynamic interactions without Bragg scattering. For Bragg scattering by an array of very slender vertical cylinders, the present approximation leads to explicit analytical results and the accuracy has been numerically confirmed by [START_REF] Li | Bragg scattering by a line array of small cylinders in a waveguide. part 1. linear aspects[END_REF] using the method of finite elements.

We thank Drs. Yuming Liu and Dick Yue of MIT for illuminating discussions. This (1) at resonance. In the plot of k 0 W vs k 0 , resonance is at higher value of k 0 if a b is smaller, as shown in Figure 13. The dimensionless capture width can be shown to be

k 0 W = k 0 C g λ g ω 2 |F D z | 2 ω 2 (λ zz + λ g ) 2 + (πa 2 b -ω 2 (πa 2 b H + µ zz )) 2 (B 3)
where F D represents the diffraction force, λ zz the radiation damping coefficient (normalized according to λ zz = λ * zz /ρg 1/2 h * 5/2 ) and µ zz the hydrodynamic mass. Using the fact that the capture width is at its peak value of unity when resonance occurs and λ g = λ zz , (B 3) can be approximated by

k 0+ -k 0-∝ a -1 b (B 8)
Consequently the peak width is larger for a smaller buoy, as shown in Figure 13.

Figure 1 .

 1 Figure 1. Geometry of the array of buoys.

  Real part of eigenfunction Fn (b) Imaginary part of eigenfunction Fn with n labelled next to the curves.with n labeled next to the curves.

Figure 2 .

 2 Figure2. First few vertical eigenfunctions in the buoy domain according to (4.9). For ω = 1, λ = 1 and f = 0.2

Figure 4 .

 4 Figure 4. Free surface elevation inside an array of length L = 1 (a) and L = 5 (b). Solid curve: Amplitude. Dashed curve: Phase difference from the undisturbed plane wave, in radians. The parameters are f = 0.2, λ = 0.5 and k0 = 1.

  For a fixed width L, the reflection coefficient R increases with the extraction rate λ, as shown in figure 5. Both the transmission coefficient T and the extraction efficiency E reach optimal? maximum values for some intermediate extraction rate around λ = 0.5 as shown in figure 7(a). The precise optimal value is around 0.5 and can be determined numerically. For a fixed extraction rate, the effects of array width L on the transmission and reflection coefficients are shown in figure 6(b). The corresponding extraction efficiency is shown in figure 7(b). The oscillatory variation of the reflection coefficient shown in figures 5(b) and 6(

Figure 5 .

 5 Figure 5. Transmission and reflection coefficients for an array of buoys with various extraction rates λ, as indicated by numbers next to each curve. The packing ratio is f = 0.2 and L = 1.

Figure 6 .

 6 Figure 6. Transmission and reflection coefficients for a buoy array with various array width L, as indicated next to each curve. The packing ratio is f = 0.2 and λ = 0.5

Figure 7 .Figure 8 .

 78 Figure 7. Variation of the extraction efficiency of a long array of finite width with (a) the extraction rate λ for a given array width and (b) the array width L for a given λ. Values of the varying parameter are indicated by numbers next to each curve. The packing ratio is f = 0.2.

Figure 9 .

 9 Figure 9. Free surface elevation in the neighbourhood of a circular array of buoys. The circumference of the array is represented by the bold circle. λ = 0.5, f = 0.2 and k0 = 1. Waves are incident from the left. Lighter colour corresponds to larger displacement. The phase difference indicated by varying shades of gray is measured in radians.

Figure 10 .

 10 Figure 10. Dependence of effectiveness on the extraction rate λ whose values are indicated next to the curves. The packing ratio is f = 0.2.

Figure 11 .

 11 Figure 11. Dependence of the effectiveness on the array radius R whose values are indicated next to the curves. (λ = 0.5 and R = 1).

  research has been supported by a grant from MASDAR Institute of Science and Technology in the program of MIT-Abu Dhabi Alliance. Partial support has been received from US-Israel Bi-National Science Foundation and from an ignition grant by MIT Earth System Initiative. Ratio of capture width to array radius.

Figure 12 .Figure 13 .

 1213 Figure 12. Dependence of effectiveness on the packing ratio f whose values are indicated next to the curves. (λ = 0.5 and R = 1)

Figure 14 .

 14 Figure 14. Finite element solutions of (3.33) and (3.34) for a tall cell with H ′ = 1, d ′ = 4 and different µ. Note the exponential decay with depth.

ga 3 b

 3 where 1 -αHω 2 ≈ 0 and α = O(1) is a constant. The values of ω 2 = ω 2 ± when k 0 W = 1/2 on both sides of the peak are found to be ω 2 = O(k 0 ). For small k 0 a b it is known that[START_REF] Mei | Theory and application of ocean surface waves[END_REF] Numerical computations show that this order of magnitude is still valid k 0 a b = O(1), henceλ zz = O k 0 a b C gSince at the peak k 0 a b = O(1), which implies ω = O(a the peak width of the k 0 W vs. k 0 curve is

Table 1 .

 1 First ten eigenvalues of (4.12) for λ = 1 and f = 0.2.

	5107 + 0.0230i 1.1165 + 0.0835i 3.3669 + 0.3159i
	K2 0.0067 + 3.0634i 0.0357 + 2.8342i 0.0545 + 2.1332i
	K3 0.0032 + 6.2448i 0.0163 + 6.1376i 0.0449 + 5.7538i
	K4 0.0021 + 9.3992i 0.0107 + 9.3286i 0.0322 + 9.0697i
	K5 0.0016 +12.5472i 0.0080 +12.4945i 0.0247 +12.2996i
	K6 0.0013 +15.6927i 0.0064 +15.6505i 0.0200 +15.4944i
	K7 0.0011 +18.8368i 0.0053 +18.8017i 0.0168 +18.6715i
	K8 0.0009 +21.9802i 0.0046 +21.9502i 0.0144 +21.8385i
	K9 0.0008 +25.1232i 0.0040 +25.0969i 0.0126 +24.9991i
	K10 0.0007 +28.2658i 0.0035 +28.2425i 0.0113 +28.1555i

L I II III Figure 3. Cross section of an infinitely long array.

† Present address : Laboratoire d'Hydrodynamique -Ecole Polytechnique, 91128 Palaiseau, France

† The dimensionless draft H of small buoys does not influence the energy extraction, but a value is chosen here to define the total volume for the array.

Appendix A. Numerical confirmation of the localization of N j (x)

The distributions of the N j in a cell are computed by the Finite Element Method, and plotted in Figure 14. For N 3 the constraint N 3 (0, 0, -1) = 0 was imposed for uniqueness. It can be seen that all solutions diminish rapidly with the depth.

Appendix B. Order estimate for a single buoy


Xavier Garnaud

Chiang C Mei
email: ‡ccmei@mit.edu

Bragg scattering and wave-power extraction by an array of small buoys

Keywords: Periodic buoy array, Multiple scattering and radiation, Bragg resonance, Wave power extraction, Homogenisation

Future designs of systems for power extraction from ocean waves will likely involve a periodic array of absorbing units. We report an asymptotic theory of scattering and radiation by a linear array of heaving buoys in a channel and attached to power-takeoff devices. The spacing between buoys is assumed to be comparable to the incident wavelength and sea depth but much greater than the buoy size. The effects of extraction rate on the buoy motion, transmission and reflection coefficients for a range of frequencies in and outside the band gap are studied. It is found that strong reflection for frequencies inside the band gap of Bragg resonance reduces the extraction efficiency significantly. For comparison an alternate theory for the efficiency away from the band gap is derived by using Froude-Krylov approximation. The predictions confirms and complements the asymptotic theory.

Introduction

Extensive theoretical studies have been devoted to the potential of power extraction from sea waves by an isolated unit such as a buoy, a raft or an oscillating water column (see reviews by [START_REF] Newman | Absorption of wave energy by elongated bodies[END_REF]; [START_REF] Falnes | Ocean Waves and Oscillating Systems[END_REF]; [START_REF] Mei | Theory and application of ocean surface waves[END_REF]; [START_REF] Cruz | Ocean Wave Energy[END_REF]; [START_REF] Mccormick | Ocean Wave Energy Conversion[END_REF]). To achieve power output comparable to a conventional power plant or a wind-turbine farm, a large array of absorbing units is necessary. The possible effects of hydrodynamic interactions among units in any geometrical deployment are therefore of design interest.

Several methods have been developed to compute the scattering and radiation by a finite number of stationary or floating bodies. In particular [START_REF] Falnes | Radiation impedance matrix and optimum power absorption for interacting oscillators in surface waves[END_REF]; [START_REF] Falnes | Wave-power absorption by parallel rows of interacting oscillating bodies[END_REF]; [START_REF] Falnes | Wave-power absorption by an array of attenuators oscillating with unconstrained amplitudes[END_REF] have examined the case of large separation where hydrodynamic interactions between bodies are weak. [START_REF] Kagemoto | Interactions among multiple three-dimensional bodies in water waves: an exact algebraic method[END_REF] have used eigenfunction expansions and addition theorems of Bessel functions to derive a numerical method to treat scattering by a few fixed vertical cylinders of circular cross sections. Infinite and semi-infinite lines of vertical cylinders have been treated semi-numerically by [START_REF] Linton | The radiation and scattering of surface waves by a vertical circular cylinder in a channel[END_REF]; Linton and Mclver (1996) using multipole expansions. Numerical studies of a finite number of cylinders have been reported by Linton and Mclver (1996); [START_REF] Chamberlain | Water wave scattering by finite arrays of circular structures[END_REF]; [START_REF] Peter | Water-wave scattering by a periodic array of arbitrary bodies[END_REF]; [START_REF] Peter | Water-wave scattering by a semi-infinite periodic array of arbitrary bodies[END_REF]; [START_REF] Mciver | Wave interaction with arrays of structures[END_REF]; [START_REF] Siddorn | Diffraction and independent radiation by an array of floating cylinders[END_REF] and [START_REF] Mavrakos | Comparison of methods for computing hydrodynamic characteristics of arrays of wave power devices[END_REF] who also considered wave power extraction without giving numerical results. These methods can in principle deal with several cylinders of different sizes, but become computationally intensive for a large number of cylinders.

It is well known that, if resonated, a single buoy with one degree of freedom (e.g., heave) can absorb all the energy within the length of 1/k of the incoming wave front. If roll is also allowed and optimized then the length can be doubled. Thus for high efficiency the ideal spacing d between adjacent buoys is 1/k, i.e., kd = 1. For environmental and navigational considerations, future wave-power farms will likely involve two-dimensional arrays with much larger spacing. Now it is known in general physics that when kd = nπ where n is an integer, waves are strongly scattered due to Bragg resonance. It is therefore of interest to predict the effect of Bragg resonance on energy absorption and buoy response.

In this article, we shall first extend the multiple-scale analysis of Li and Mei (2007a) for a periodic array of fixed slender piles to small movable buoys which scatter, radiate and absorb energy by its heave motion only. To simulate their potential for power absorption we assume that each buoy is attached to a linear device which converts the mechanical energy of the buoy to electricity. No phase control of the power-takeoff system is assumed. Analytical results for the scattering coefficients, energy absorption rate and buoy motion will be derived and discussed.

Scales and normalization

Consider a linear array of small buoys in a long channel of constant width d and mean depth h, as shown in figure 1. Simple harmonic waves arrive from one end of the channel x ∼ -∞. In the framework of potential theory for inviscid and incompressible fluid, the mathematical problem is equivalent to an infinitely long strip of buoys in an rectangular lattice, attacked by a plane wave with crests parallel to the edge of the strip.

The geometry has three sharply different length scales: the small radius a of the buoy, the water depth h, and the large horizontal extent of the array. These scales will be respectively referred to as the micro-scale for the near field, the meso-scale for the intermediate field, and the macro-scale for the far field. In all three fields the vertical displacements of the free surface (η) and of the buoys (ζ) are characterized by the incident wave amplitude A. Let us first introduce the following dimensionless variables, distinguished by primes, for the intermediate field,

Φ = A ghΦ ′ , p = ρgA p ′ , (η, ζ) = A(η ′ , ζ ′ ), (2.1) (x, y, z, d) = h(x ′ , y ′ , z ′ , d ′ ), t = t ′ h g . (2.2)
We shall assume that the wavelength, mean depth and buoy separation are of the same order of magnitude, so that

k ′ = kh = O(1), d ′ = d h = O(1), (2.3)
are both of order unity but the buoy radius and drafts are small:

a h ≡ µ ≪ 1, and H h = O(µ) ≪ 1. (2.4)
In the near field of a buoy the physics is controlled by the much smaller radius a, hence it is proper to employ the micro-scale coordinates, distinguished by bars and defined by

(x, ȳ, z) = 1 a (x, y, z) = 1 µ (x ′ , y ′ , z ′ ). (2.5)
It is known that the scattered wave from a small cylinder is smaller than the incident wave by a factor of order (ka) 2 , as in the case sound (see, e.g., [START_REF] Mei | Theory and application of ocean surface waves[END_REF]). It will also be shown in the next section that a small heaving buoy extracts a fraction (ka) 2 of the incoming energy. Therefore, significant scattering and radiation effects can be expected to accumulate in an array with O(1/µ 2 ) buoys, or after a distance of h/µ 2 . We therefore define the following macro-scale coordinates for the far-field to describe the evolution of wave envelopes:

(X, Y, T ) = µ 2 (x ′ , y ′ , t ′ ).
(2.6)

Away from Bragg resonance

As a reference for later comparison, we first study the performance of an array of buoys well separated from each other without Bragg resonance. In this case, interaction between buoys can be neglected as a first approximation. Consider one small buoy in a plane incident wave of frequency ω and amplitude A. The potential of the incoming wave is

Φ I = φ I e -iωt with φ I = Ag iω cosh k(z + h)) cosh(kh) e ikx , (3.1)
where ω and k are related by the dispersion relation

ω 2 = gk tanh(kh). (3.2)
Because of the small size of buoys, the scattered and radiated waves are negligible.

Froude-Krylov approximation can be applied so that the hydrodynamic pressure on

Article submitted to Royal Society each buoy is dominated by the undisturbed incoming wave (see [START_REF] Newman | Absorption of wave energy by elongated bodies[END_REF]).

The vertical exciting force on the buoy is therefore iρω

SB φ I (0, 0, 0) dS = ρgAπa 2 . (3.3)
Let us assume that an energy extraction device is attached to each buoy and exerts a load force iωλ g ζ where λ g denotes the extraction rate. Since the added buoyancy force due to heave is -πa 2 ρgζ, Newton's law gives

-M ω 2 ζ = πρga 2 A + iωλ g ζ -πρga 2 ζ, (3.4)
where M = ρπa 2 H by Archimedis principle. It follows that

ζ A = 1 1 - iωλg ρgπa 2 -ω 2 H g = 1 1 - iωλg ρgπa 2 + O(µ). (3.5)
Use has been made of the fact that ω 2 H/g ∼ kH = O(µ) ≪ 1. Thus the inertia of a small buoy is relatively unimportant. As a consequence, the draft of the buoy H is much less relevant than its lateral dimension a. The time-averaged rate of energy extraction by a single buoy is given by

λ g ∂ ∂t Re (ζe -iωt ) 2 = 1 2 ω 2 λ g |ζ| 2 .
(3.6)

The fraction of power extracted by one buoy from the incident energy flux across a channel width d is then

E = E E inc = 1 2 ω 2 λ g |ζ 2 | 1 2 ρgA 2 C g d , (3.7)
where E inc denotes the rate of energy influx over a width d and

C g = ω 2k 1 + 2kh sinh(2kh) (3.8)
is the group velocity of the incident wave.

Let the normalised frequency and extraction rate be denoted by

ω ′ ≡ ω g h , λ ′ g ≡ λ g πρa 2 √ gh ,
then in dimensionless form

ω ′ 2 = k ′ tanh(k ′ ), C ′ g = ω ′ 2k ′ 1 + 2k ′ sinh 2k ′ , and ζ ′ = 1 1 -iω ′ λ ′ g + O(µ) (3.9)
Consider now a large number of buoys along the centreline of a channel and spaced at the same distance d. On the macro-scale of the total array, the mathematical effect of many point absorbers can be replaced by a continuous distribution. The fraction of the incoming energy absorbed within a unit distance of the macroscale coordinates must be

Article submitted to Royal Society

h d πa 2 h 2 λ ′ 2 g ω ′2 | ζ| 2 d ′ C ′ g = µ 2 πā d ′ 2 C ′ g λ ′ g ω ′2 1 + ω ′2 λ ′ 2 g ≡ µ 2 D.
(3.10)

In terms of the macro coordinate X = µ 2 x ′ , the spatial rate of change of E(X) is

dE dX = -D E, i.e. E(X) = e -DX . (3.11)
Thus the fraction of energy remaining at the end of the array X = L (i.e. x ′ = L/µ 2 ) is E(L) = e -DL and the extraction efficiency is

E 0 ≡ 1 -E(L) = 1 -e -DL .
(3.12)

Clearly the efficiency depends on the number of buoys N = L/(µ 2 d ′ ), the frequency of the incoming wave through ω ′ and C ′ g , and the extraction rate λ ′ g . For sufficiently large L, E 0 approaches unity. The dependence of efficiency on k ′ and the extraction rate is plotted in figure 2(a) for a fixed L = 1. For a given L, the optimal extraction rate for maximum efficiency is found from

dE 0 dλ ′ g = 0, or equivalently dD dλ ′ g = 0,
which gives the optimal extraction rate

λ ′ g,opt = 1 ω ′ . (3.13)
Thus the optimum extraction rate should be higher for longer waves. Figure 2(b) shows the efficiency for a few L/d ′ , when the extraction rate is optimized for every k ′ . One sees that at high frequencies, a long enough array can extract all the incoming energy, and a larger array is better for low-frequency waves.

We now turn to the physics of Bragg resonance.

Scattering by an array of fixed buoys (a) Linearised dimensionless equations

As it is standard in linearised theories, the problem of wave-body interaction is equivalent to the sum of two hydrodynamic problems : scattering by stationary bodies and radiation by body motion. The two problems are coupled by the dynamics of the floating bodies. Let us first study the diffraction of a wave-train by fixed buoys. The following symbols for different parts of the physical domains are employed: Ω F for the fluid domain, S F for the free surface, S W for the lateral surface of the buoys, and S B for the bottom surface of the buoys. The seabed is the entire horizontal plane z ′ = -1. In terms of the meso-scale coordinates, the dimensionless governing equations for the scattering potential Φ ′ is

∆ ′ Φ ≡ ∂ 2 Φ ′ ∂x ′ 2 + ∂ 2 Φ ′ ∂y ′ 2 + ∂ 2 Φ ′ ∂z ′ 2 = 0, x ′ ∈ Ω F , (4.1a) ∂ ∂z ′ + ∂ 2 ∂t ′2 Φ ′ = 0, x ′ ∈ S F , (4.1b) ∂Φ ′ ∂z ′ = 0, x ′ ∈ S B , (4.1c
)

∂Φ ′ ∂r ′ = 0, x ′ ∈ S W , (4.1d) ∂Φ ′ ∂z ′ = 0, z ′ = -1, (4.1e) ∂Φ ′ ∂y ′ = 0, y ′ = ± d ′ 2 . (4.1f)
We shall now find the law for the slow evolution of the envelope from the fast variations between buoys, by the method of multiple scales (homogenization).

(b) Envelope equations by multiple-scale analysis Substituting into (4.1) the following expansion

Φ ′ S = φ 0 + µ 2 φ 1 + • • • e -iω ′ t ′ , (4.2)
with φ n , n = 0, 1, 2, • • • being dimensionless functions of (x ′ , y ′ , z ′ ; X; T ), we get from (4.1a) and (4.1b) the following equations

∆ ′ + 2µ 2 ∂ 2 ∂X∂x ′ + µ 4 ∂ 2 ∂X 2 φ 0 + µ 2 φ 1 + • • • = 0, x ′ ∈ Ω F , (4.3a) ∂ ∂z ′ -ω ′2 -2µ 2 iω ′ ∂ ∂T -µ 4 ∂ 2 ∂T 2 φ 0 + µ 2 φ 1 + • • • = 0, x ′ ∈ S F . (4.3b)
From successive orders of these and the remaining equations in (4.1) we obtain boundary-value problems for the wave-scale variations in a unit cell. Because each cell is one period in vary large array, we invoke Bloch's theorem [START_REF] Ashcroft | Solid State Physics[END_REF]) which states that the solution φ n should be of the form e ±ikx ′ f (x ′ ) where f (x ′ ) is periodic in x ′ with the period d ′ . Since we shall focus on the state of Bragg resonance with kd ′ = π, Bloch's theorem implies

φ n - d ′ 2 , y ′ = φ n d ′ 2 , y ′ , and ∂φ n ∂x ′ - d ′ 2 , y ′ = - ∂φ ∂x ′ d ′ 2 , y ′ . (4.4) (i) Leading order
At the leading order O(µ 0 ), the governing equations are homogeneous:

∆ ′ φ 0 = 0, x ′ ∈ Ω F , (4.5a) ∂ ∂z ′ -ω ′2 φ 0 = 0, x ′ ∈ S F , (4.5b) ∂φ 0 ∂z ′ = 0, z ′ = -1, (4.5c) ∂φ 0 ∂y ′ = 0, y ′ = ±d ′ /2. (4.5d)
As reasoned in Li and Mei (2007a,b), since µ = a/h ≪ 1, H ′ ≪ 1 and d ′ = O(1), the areas of the buoy surfaces S B and S W are of order O(µ 2 ). Hence the buoys have negligible effects on the waves at the leading order, i.e., the boundary conditions on S B and S W are ineffective until at higher orders. The formal solution is the sum of left-and right-going plane waves in free space

φ 0 = α + (X, T )Z(z ′ )e ik ′ x ′ + α -(X, T )Z(z ′ )e -ik ′ x ′ , (4.6)
where

Z(z ′ ) = 1 iω ′ cosh(k ′ (z ′ + 1)) cosh(k ′ ) . (4.7)
The no-flux boundary conditions at y = ±d/2 and Bloch conditions (4.4) are trivially satisfied. The envelope functions α ± (X, T ) are yet to be found.

(ii) First order From (4.3) the first order O(µ) meso-scale problem in the unit cell is inhomogeneous, and governed by

∆ ′ φ 1 = -2 ∂ 2 φ 0 ∂X∂x ′ , x ′ ∈ Ω F , (4.8a) ∂ ∂z ′ -ω ′2 φ 1 = 2iω ′ ∂φ 0 ∂T , x ′ ∈ S F , (4.8b
)

∂φ 1 ∂z ′ = - 1 µ 2 ∂φ 0 ∂z ′ , x ′ ∈ S B , (4.8c
)

∂φ 1 ∂r ′ = - 1 µ 2 ∂φ 0 ∂r ′ , x ′ ∈ S W , (4.8d
)

∂φ 1 ∂z ′ = 0, z ′ = -1, (4.8e) ∂φ 1 ∂y ′ = 0, y ′ = ± d ′ 2 , (4.8f)
together with the Bloch conditions (4.4). Despite the large factor µ -2 , the integrated effects of the boundary value in (4.8c) and (4.8d) are of order 1 since the area of the buoy surface is of order µ 2 . An important estimate of φ 1 near and on the buoy may be deduced. In terms of the meso-scale coordinates, we have, in the neighborhood of the buoy,

∂φ 1 ∂n ′ ∼ φ 1 µ ≫ φ 1 .
Now (4.8c) and (4.8d) imply that

∂φ 1 ∂n ′ ∼ 1 µ 2 , (4.9) hence φ 1 ∼ 1 µ (4.10)
in the neighbourhood of the buoy. We next apply Green's formula to φ 1 and the homogeneous solutions

ψ ∓ ≡ Z(z ′ ) exp(±ik ′ x ′ ) (4.11) over the volume of cell Ω F : ΩF φ 1 ∆ ′ ψ ∓ -ψ ∓ ∆ ′ φ 1 dV ′ = ∂ΩF φ 1 ∂ψ ∓ ∂n ′ -ψ ∓ ∂φ 1 ∂n ′ dS ′ , (4.12)
which can be rewritten after using the boundary conditions in (4.8) as

2 ΩF Z(z ′ )e ∓ik ′ x ′ ∂ 2 φ 0 ∂X∂x ′ dV ′ = -2iω ′ SF ψ ∓ ∂φ 0 ∂T dS ′ + SB ∪SW 1 µ 2 ψ ∓ ∂φ 0 ∂n ′ + φ 1 ∂ψ ∓ ∂n ′ dS ′ . (4.13)
Article submitted to Royal Society This is the solvability condition for the inhomogeneous problem of φ 1 and gives a constraint on φ 0 . To the leading order, the left hand side of (4.13) can be simplified to

2 ΩF = ±2ik ′ d ′2 ∂α ± ∂X 0 -1 Z(z ′ ) 2 dz ′ + O(µ 3 ).
The integral over the free surface on the right hand side gives

SF = - 2d ′2 iω ′ ∂α ± ∂T + O(µ 2 ).
In view of (4.10), the second term in the surface integral over the buoy can be neglected with an error of O(µ). The remaining surface integral on the buoy is

SB ∪SW ≈ 1 µ 2 SB ∪SW ψ ∓ ∂ [α + ψ + + α -ψ -] ∂n ′ dS ′ = - 1 µ 2 SW ik ′ ψ ∓ α + ψ + -α -ψ -e r .e x dS ′ + 1 µ 2 SB ψ ∓ dZ(z ′ ) dz ′ α + e ik ′ x ′ + α -e -ik ′ x ′ dS ′ .
Let us now use the fact that the buoy draft is small

H ′ ≡ H/h = O(µ) so that Z ≈ 1 iω ′ and dZ dz ′ ≈ -iω ′ on S W ∪ S B (4.14)
from (4.7). In terms of the micro-scale variables r = r/a and z = z/a so that dS ′ /µ 2 = d S, we have

SB∪SW ≈ -(α + -α -) 2π 0 0 - H ik ′ cos(θ) dz dθ -(α + + α -) 2π 0 1 0 r dr dθ ≈ -π(α + + α -).
Hence (4.13) becomes

2d ′2 iω ′ ∂α ± ∂T ± 2ik ′ d ′2 ∂α ± ∂X 0 -1 Z(z) 2 dz = -π α + + α -,
which can be simplified to the coupled-mode equations on the macro-scale

∂α ± ∂T ± C ′ g ∂α ± ∂X = -iΩ 0 α + + α - (4.15)
where Ω 0 is the dimensionless coupling constant Equation (4.15) can be rewritten in physical coordinates as

Ω 0 = πω ′ 2d ′2 . (4.16) Article submitted to Royal Society -2 0 2 4 6 -4 -3 -2 -1 0 1 2 3 4 Ω/Ω0 C ′ g K S Ω 0
∂α ± ∂t ± C g ∂α ± ∂x = -i ωa 2 2d 2 α + + α - (4.17)
These equations govern the envelopes of the forward and backward waves in an array of fixed small buoys, and are similar to those for wave propagation in an array of slender piles or over a periodic seabed (see Li and Mei (2007a,b); [START_REF] Naciri | Bragg scattering of water waves by a doubly periodic seabed[END_REF]). The difference with Li and Mei (2007a,b) is due to the fact that in the present problem the scattering effect comes from the bottom of the buoy and not from the lateral wall.

(c) Macro-scale dispersion relation and band gap

As the first application of (4.15) let us consider periodically modulated wave envelopes in an infinite array, α ± ∝ e i(±KSX-ΩT ) , (4.18)

where Ω corresponds to a frequency detuning of µ 2 Ω. The homogeneous equations (4.15) have non-trivial solutions only if

C ′2 g K 2 S = Ω(Ω -2Ω 0 ), (4.19)
which is plotted in figure 4. For real Ω three regions can be distinguished. In either Ω ≤ 0 or Ω ≥ 2Ω 0 , K S is real, hence waves propagate. However, inside the band gap defined by 0 ≤ Ω ≤ 2Ω 0 , K S is imaginary, hence propagation is forbidden.

The wave train can only decay in distance. The spatial rate of amplitude decay is proportional to Im(K S ), which is maximum at Ω = Ω 0 . The presence of the band gap and its limits have been confirmed by numerical computation of the band structure but without the small-buoy assumption, using a standard numerical approach of solid-state physics [START_REF] Garnaud | Wave-power absorption by buoys Master of Science Thesis[END_REF]).

(d ) Scattering by an array of finite width

Let there be a finite array of fixed buoys in 0 < X < L, where L corresponds to the physical width (L/µ 2 )h. An incoming wave slightly detuned from Bragg

Article submitted to Royal Society resonance φ in = Z(z ′ )e ik ′ x ′ e i(K0X-ΩT ) (4.20) arrives from X ∼ -∞ with K 0 ≡ C ′ g Ω.
Let the scattering potential be

φ 0 = e ik ′ x ′ e iK0X + C R e -ik ′ x ′ e -iK0X Z(z ′ )e -iΩT , X < 0, (4.21)
on the incidence side and

φ 0 = C T e ik ′ x ′ e iK0(X-L) Z(z ′ )e -iΩT , X > L, (4.22)
on the transmission side. C R and C T are respectively the complex reflection and transmission coefficients. The potential inside the array is

φ 0 = Z(z ′ ) α + (X)e ik ′ x ′ + α -(X)e -ik ′ x ′ e -iΩT , (4.23)
where α ± are defined by

α ± = α ± e -iΩT . (4.24)
By differentiating (4.15), we find

d 2 α ± dx 2 + K 2 S α ± = 0, (4.25)
with

K S ≡ Ω S C ′ g and Ω S ≡ Ω(Ω -2Ω 0 ). (4.26)
Requiring continuity of the leading-order pressure and horizontal velocity at the edges of the array, we must have

φ 0 | X=0 -= φ 0 | X=0 + , φ 0 | X=L -= φ 0 | X=L + , (4.27a) ∂φ 0 ∂x ′ X=0 - = ∂φ 0 ∂x ′ X=0 + , ∂φ 0 ∂x ′ X=L - = ∂φ 0 ∂x ′ X=L + , (4.27b)
where L + means slightly greater than L and L -slightly less that L. Using (4.21) and (4.23) in (4.27) we obtain

1 + C R = α + (0) + α -(0), C T = α + (L) + α -(L), (4.28a) 1 -C R = α + (0) -α -(0), C T = α + (L) -α -(L). (4.28b) It follows that α -(L) = 0, α + (0) = 1. (4.29a)
The general solutions of (4.25) are of the form

α ± = C ± 1 cos(K S X) + C ± 2 sin(K S X).
Article submitted to Royal Society The six coefficients C R , C T , C ± 1 and C ± 2 can be found from the four conditions (4.28a), (4.28b) and the coupled equations (4.15). Finally the amplitudes inside the array are given by

α + (X) = i(Ω -Ω 0 ) sin(K S (X -L))) + Ω S cos(K S (X -L)) -i(Ω -Ω 0 ) sin(K S L) + Ω S cos(K S L) (4.30a) and α -(X) = iΩ 0 sin(K S (X -L)) -i(Ω -Ω 0 ) sin(K S L) + Ω S cos(K S L) (4.30b)
In the open sea, the transmission and reflection coefficients are

C T = α + (L) = Ω S -i(Ω -Ω 0 ) sin(K S L) + Ω S cos(K S L) , (4.31a)
and

C R = α -(0) = -iΩ 0 sin(K S L) -i(Ω -Ω 0 ) sin(K S L) + Ω S cos(K S L) . (4.31b)
The numerical results are qualitatively similar to those for vertical piles. As shown in figure 5, the scattering coefficients depend strongly on the detuning frequency Ω. Outside the band gap, Ω/Ω 0 < 0 or Ω/Ω 0 > 2, the scattering effects of the array are weak, as C R is small and C T is close to unity. Both are oscillatory in Ω/Ω 0 . Inside the band gap, 0 ≤ Ω/Ω 0 ≤ 2, propagation is inhibited. The reflection coefficient is close to 1 and the transmission coefficient close to 0. Inside the array the free surface profile is oscillatory in space if Ω is outside the gap and exponentially attenuating if inside, as shown in figure 6.

Next let us examine the effects of buoy motion on waves.

Envelope of radiated waves

In this section, we allow the small buoys to heave either freely or partially constrained by energy-absorbing devices. Since k ′ d ′ = π, if one buoy goes up then its immediate neighbours must go down, and vice versa. Across the wide array of many buoys, the amplitude of waves and buoy displacements will be a slowly

Article submitted to Royal Society varying function of space. The normalised heave amplitude of the m th buoy can therefore be expressed as

ζ ′(m) = (-1) m ζ 0 (X, T ) + µ 2 ζ 1 (X, T ) + • • • e -iω ′ t ′ .
(5.1)

The buoy displacement ζ 0 (X, T ) is yet unknown and will be found later by additional account of the buoy dynamics. Let the centre of the m th buoy be located at x ′ m = md ′ . For mathematical convenience, we express the term (-1) m as

(-1) m = e ik ′ x ′ m = e imπ (5.2)
in the kinematic boundary conditions on the buoys. Denoting the radiation potential by Φ R and expanding

Φ ′ R = Φ 0 + µ 2 Φ 1 + • • • e -iω ′ t ′ , (5.3)
the governing equations become W are respectively the bottom and lateral boundaries of the m th buoy. As shown in figure 3, we require the radiation potential to be anti-periodic in x ′ with the period d ′ ,

∆ ′ + 2µ 2 ∂ 2 ∂X∂x ′ + µ 4 ∆ Φ 0 + µ 2 Φ 1 + • • • = 0, x ′ ∈ Ω F , (5.4a) ∂ ∂z ′ -ω ′ 2 -2µ 2 iω ′ ∂ ∂T -µ 4 ∂ 2 ∂T 2 Φ 0 + µ 2 Φ 1 + • • • = 0, x ′ ∈ S F , (5.4b) ∂ Φ 0 + µ 2 Φ 1 + • • • ∂z ′ +iω ′ e ik ′ x ′ m ζ 0 + µ 2 ζ 1 + • • • = 0, x ′ ∈ S (m) B , (5.4c) ∂ Φ 0 + µ 2 Φ 1 + • • • ∂z ′ = 0, z ′ = -1, (5.4d) ∂ ∂r ′ + µ 2 ∂ ∂R Φ 0 + µ 2 Φ 1 + • • • = 0, x ′ ∈ S (
Φ i (x ′ + d ′ , y ′ , z ′ ) = -Φ i (x ′ , y ′ , z ′ ), ∂Φ i ∂x ′ (x ′ + d ′ , y ′ , z ′ ) = - ∂Φ i ∂x ′ (x ′ , y ′ , z ′ ). (5.5)
Article submitted to Royal Society for i = 1, 2, • • • . From (5.4), it is evident that the zeroth order radiation potential is of the form

Φ 0 = β + (X, T )e ik ′ x ′ + β -(X, T )e -ik ′ x ′ Z(z ′ ), (5.6)
where the long-scale functions β ± represent the unknown amplitudes of the propagating waves. The first-order problem is governed by

∆ ′ Φ 1 = -2 ∂ 2 Φ 0 ∂X∂x ′ , x ′ ∈ Ω F , (5.7a) ∂Φ 1 ∂z ′ -ω ′2 Φ 1 = 2iω ′ ∂Φ 0 ∂T , x ′ ∈ S F , (5.7b
)

∂Φ 1 ∂z ′ = -iω ′ ζ 0 e ik ′ x ′ n µ 2 - 1 µ 2 ∂Φ 0 ∂z ′ , x ′ ∈ S (m) B , (5.7c
)

∂Φ 1 ∂z ′ = 0, z ′ = -1, (5.7d) ∂Φ 1 ∂r ′ = - 1 µ 2 ∂Φ 0 ∂r ′ , x ′ ∈ S (m)
W .

(5.7e)

Note that as the buoy radii and drafts are of order µ, we have

x ′ = x ′ m + O(µ) for x ′ ∈ S (m) W ∪ S (m)
B . The problem (5.7) is similar to the scattering problem (4.8) except for the additional term proportional to the buoy displacement, which is also anti-periodic with a period d ′ .

As in the scattering problem, we derive the solvability condition for Φ 1 by applying Green's formula to Φ 1 and ψ

± = Z(z ′ ) exp(±ik ′ x ′ ) over a cell Ω F .
It is easily checked that the only change is in the surface integral over the buoys

SB ∪SW ≈ 1 µ 2 SB ∪SW ψ ∓ ∂ [β + ψ + + β -ψ -] ∂n ′ dS ′ - 1 µ 2 SB ψ ∓ iω ′ ζ 0 e ik ′ x ′ m dS ′ = - 1 µ 2 SW ik ′ ψ ∓ β + ψ + -β -ψ -e r .e x dS ′ + 1 µ 2 SB ψ ∓ dZ(z ′ ) dz ′ β + e ik ′ x ′ m + β -e -ik ′ x ′ m dS ′ - 1 µ 2 SB ψ ∓ iω ′ ζ 0 e ik ′ x ′ m dS ′ . (5.8)
where

S B ≡ ∪ m S (m)
B . Using (4.14), we get

SB ∪SW ≈ -(β + -β -) 2π 0 0 - H ik ′ cos(θ) dz dθ -(β + + β --ζ 0 ) 2π 0 1 0 r dr dθ ≈ -π(β + + β --ζ 0 ).
(5.9)

Article submitted to Royal Society

This finally gives us the envelope equations

∂β ± ∂T ± C ′ g ∂β ± ∂X = -iΩ 0 β + + β --ζ 0 (5.10)
where Ω 0 was defined in (4.16). In physical terms, (5.10) reads

∂β ± ∂t ± C g ∂β ± ∂x = -i ωa 2 2d 2 β + + β -- ζ 0 A .
(5.11)

The pair of equations in (5.10) expresses the coupling between the amplitude of the right-and left-going waves and the unknown buoy motion.

We must now examine the buoy displacement induced by the waves, in order to relate ζ 0 to α ± and β ± in (5.10).

Buoy dynamics

The forcing term ζ ′ in the long-scale equation for the radiation problem are related to the scattering and radiation potentials, given respectively by (4.6) and (5.6). By Froude-Krylov approximation, the leading order vertical forces on a buoy by the scattering and radiation potentials are given respectively by

f ′ m = iω ′ S (m) B φ 0 dS = iπµ 2 ω ′ (-1) m iω ′ α + + α -, (6.1)
F ′ m = iω ′ S (m) B Φ 0 dS = iπµ 2 ω ′ (-1) m iω ′ β + + β -, (6.2)
where the forces are normalised according to

(f m , F m ) = ρgh 2 A(f ′ m , F ′ m ). (6.3)
Use is made of the fact that the neighboring buoys move in opposite phases.

Let us assume that the energy extraction device exerts a force

λ ′ g ∂ζ ′(m) ∂t ′ (6.4)
on the m th buoy. Applying Newton's law to the m th buoy, we get

M ω 2 ζ ′(m) = f m + F m + iωλ ′ g ζ ′(m) -πρa 2 ζ ′(m) , (6.5)
which gives us in dimensionless form

-ω ′2 µ 3 π Hζ 0 (-1) m = πµ 2 α + + α -(-1) m + πµ 2 β + + β -(-1) m + iπω ′ µ 2 ζ 0 (-1) m λ ′ g -πµ 2 ζ 0 (-1) m . (6.6)
Again the mass of the buoy can be ignored so that

ζ 0 = G α + + α -+ β + + β -, (6.7)
Article submitted to Royal Society

where for brevity we denote

G ≡ 1 1 -iω ′ λ ′ g + O(µ), (6.8)
which is the same as (3.9). Using (6.7) we can rewrite (5.10) as:

∂β ± ∂T ± C ′ g ∂β ± ∂X = -iΩ 0 (1 -G)(β + + β -) -G(α + + α -) (6.9)
which couple the radiation and scattering components. The scattering amplitudes α ± are already found in the previous section and serve as forcing terms here.

We can now study the envelopes due to waves interacting with a finite array of energy-extracting buoys.

Radiation by a finite array of energy-extracting buoys (a) Solution for radiation amplitudes

Consider again an array of finite but large width (L/µ 2 )h forced to oscillate by the incident and scattered waves and subjected to the reactive forces from the extractors and from the buoy motion. The potential outside the array must satisfy the radiation condition, hence

Φ 0 = C -Z(z ′ )e -ik ′ x ′ e -iKX e -iΩT , X < 0, (7.1a) Φ 0 = C + Z(z ′ )e ik ′ x ′ e iK(X-L) e -iΩT , X > L, (7.1b)
where the complex coefficients C + and C -are yet unknown. As in the scattering problem, we introduce β ± by defining β ± = β ± e -iΩT and assume the solution inside the array to be of the form

Φ 0 = β + (X)e ik ′ x ′ + β -(X)e -ik ′ x ′ Z(z ′ )e -iΩT . (7.2)
Continuity of pressure and horizontal velocity are required at the edges of the array

Φ 0 | X=0 -= Φ 0 | X=0 + , Φ 0 | X=L -= Φ 0 | X=L + , (7.3a) ∂Φ 0 ∂x ′ X=0 - = ∂Φ 0 ∂x ′ X=0 + , ∂Φ 0 ∂x ′ X=L - = ∂Φ 0 ∂x ′ X=L + . (7.3b)
Using (7.1) and ( 7.2), we get

C -= β + (0) + β -(0), C + = β + (L) + β -(L), (7.4a) -C -= β + (0) -β -(0), C + = β + (L) -β -(L). (7.4b) It follows that β -(L) = 0, β + (0) = 0, (7.5a)
Article submitted to Royal Society which imply C -= β -(0) and C + = β + (L). To complete the analytical solution, let us write the forcing term in (6.9) as (7.6) Using (4.30), we find

α + + α -= C g A + e iKS x + A -e -iKSx .
A + = e -iKSL K S + K 2 (-i(Ω -Ω 0 ) sin(K S L) + Ω 1 cos(K S L))
, (7.7a)

A -= e iKSL K S -K 2 (-i(Ω -Ω 0 ) sin(K S L) + Ω 1 cos(K S L)) . (7.7b)
To find β ± , it is simpler to use linearity and just to solve first for one of the complex exponentials on the right-hand side of (6.9). Let the response to the forcing (e iKSX , e -iKSX ) be denoted by ( β ± 1 , β ± 2 ) respectively, then (7.8b) where K ≡ Ω/C ′ g and K 0 ≡ Ω 0 /C ′ g . By cross-differentiation these two equations can be decoupled to give

-iK β + 1 + ∂ β + 1 ∂X = -iK 0 (1 -G)( β + 1 + β - 1 ) + e iKS X , (7.8a) -iK β - 1 - ∂ β - 1 ∂X = -iK 0 (1 -G)( β + 1 + β - 1 ) + e iKS X ,
d 2 β + 1 dX 2 + K 2 R β + 1 = i (K + K S ) e iKSX , (7.9a) d 2 β - 1 dX 2 + K 2 R β - 1 = i (K -K S ) e iKSX , (7.9b) with K R ≡ K [K -2K 0 (1 -G)] (7.10)
being the natural wave number of the radiated wave. As G is complex for any non zero λ ′ g , K R will be likewise, implying spatial attenuation. Since K R = K S , the spatial period of the forcing term differs from the natural period in (6.9) and no resonance is expected. Variations of the real and imaginary parts are shown in figure 7 as functions of the detuning frequency Ω and the extraction rate λ ′ g . One sees that there is no band gap and that all radiated waves are spatially attenuated. The special case of free buoys without extractors will be treated later.

The general solutions of (7.9) are,

β + 1 = C (1) 1 e iKRX + C (2) 1 e -iKRX + i K + K S K 2 R -K 2 S e iKSX , (7.11a)
β - 1 = C (3) 1 e iKRX + C (4) 1 e -iKRX + i K -K S K 2 R -K 2 S e iKSX . (7.11b)
For brevity we let

B + 1 = K + K S K 2 R -K 2 S , B - 1 = K -K S K 2 R -K 2 S .

Article submitted to Royal Society

We can now solve for

C (i)
1 for all i = 1, 2, 3, 4:

C (1) 1 = -i K R B + 1 (K R + K S ) + B - 1 K G + 1 cos(K R L) -A - 1 K G K R e iKSL +i -B + 1 (K -K G )(K R + K S ) + A - 1 K G K S -K sin(K R L) [2K R (K R cos(K R L) -i(K -K G ) sin(K R L))] -1 , (7.15a)

C

(2)

1 = i K R B + 1 (K S -K R ) + B - 1 K G + 1 cos(K R L) -B - 1 K G K R e iKSL +i -B + 1 (K -K G )(-K R + K S ) + B - 1 K G K S -K sin(K R L) [2K R (K R cos(K R L) -i(K -K G ) sin(K R L))]
-1 , (7.15b)

C (3) 1 = -i B - 1 (K R -K + K G )e iKS L + (K S + K -K G )e -iKRL -e -iKRL 2 (K R cos(K R L) -i(K -K G ) sin(K R L)) , (7.15c)
C (4) 1 = i -B - 1 (K R + K -K G )e iKS L -(K S + K -K G )e iKRL -e iKRL 2 (K R cos(K R L) -i(K -K G ) sin(K R L)) . (7.15d)
Corresponding to the forcing exp(-iK S x), the responses β ± 2 can be treated similarly. Let the solution be of the form

β + 2 = C (1) 2 e iKRX + C (2) 2 e -iKRX + i K -K S K 2 R -K 2 S e iKSX , (7.16a)
β - 2 = C (3) 2 e iKRX + C (4) 2 e -iKRX + i K + K S K 2 R -K 2 S e iKSX . (7.16b)
The coefficients C

2 , j = 1, 2, 3, 4 can be obtained by replacing

K S by -K S in C j 1 , i,e., C (j) 2 (K R , K S ) = C (j) 1 (K R , -K S ), j = 1, 2, 3, 4. (7.17)
It is now easy to see that

β + = iΩ 0 G A + β + 1 + A -β + 2 , (7.18a
)

β -= iΩ 0 G A + β - 1 + A -β - 2 , (7.18b)
are the solutions of the radiation problem, satisfying (6.9) as well as the boundary conditions (7.5). We have also confirmed these formulas by direct numerical solution by the Finite Volume method.

(b) Freely floating buoys

In the limiting case of freely floating buoys, λ ′ g = 0 so that G = 1. The evolution equations for β ± are no longer coupled. After omitting the factor e -iΩT , we get simply

-iK β + + ∂ β + ∂X = -iK 0 ( α + + α -), (7.19) -iK β -- ∂ β - ∂X = -iK 0 ( α + + α -).
(7.20)

Article submitted to Royal Society

These decoupled first order ordinary differential equations have solutions of the form

β + = C 1 e iKX - K 0 A + K -K S e iKS X - K 0 A - K + K S e -iKSX , (7.21)
β -= C 2 e -iKX - K 0 A + K + K S e iKS X - K 0 A - K -K S e -iKSX , (7.22)
with constants C 1 and C 2 to be determined. Note that

K 0 K -K S + K 0 K + K S = 2 K 0 K K 2 -K 2 S = 1 (7.23)
after using the fact that the scattering wave number is given by KX+ΩT ) . This gives the total potential

K S = √ K 2 -2KK 0 , so that Φ 0 = β + + β -e -iΩT = (C 1 -α + )e i(KX-ΩT ) + (C 2 -α -)e -i(
φ 0 + Φ 0 = Z(z ′ ) α + + α -+ β + + β -e -iΩT = Z(z ′ ) C 1 e i(KX-ΩT ) + C 2 e -i(KX+ΩT ) , (7.24a) and thus ζ 0 = η 0 = C 1 e i(KX-ΩT ) + C 2 e -i(KX+ΩT ) , (7.24b)
where η 0 is defined in a similar way as ζ 0 . Matching directly the total potential φ 0 + Φ 0 with the potential outside the array, we find

C 1 = 1, C 2 = 0, (7.25) so η 0 = ζ 0 = e i(KX-ΩT ) . (7.26)
Hence the buoys move with the same amplitude and phase as the surrounding free surface and to the leading order a sparse array of small freely floating buoys does not affect the incoming wave. This result is expected as the consequence of negligible buoy inertia.

The combined effects of scattering, radiation and extraction

The combined effect of scattering, radiation and energy extraction on the free surface and buoy displacements are presented in figures 8 and 9. Two representative values of the detuning are chosen here: one outside the band gap of the scattering problem (Ω = -Ω 0 ) and one inside (Ω = Ω 0 ). One sees in figure 9 that the effect of the length is minor. In contrast, the effect of the extraction force shown in figure 8 has much stronger influence, as shown for the the two limits of fixed and free buoys. The displacements of the free surface and the buoys decrease faster through

Article submitted to Royal Society As shown in §7(b), the free surface amplitude is constant and equal to 1 in the limiting case of free buoys λ ′ g = 0.

the array for frequencies inside the band gap, and as the extraction rate increases.

Understandably the displacement of the buoy is always smaller than that of the water surface. Note that there is no resonance.

In order to characterize the waves outside of the array, let us introduce transmission and reflection coefficients for the complete problem, denoted respectively by T and R -such that the free surface in the open water regions is given by η 0 = e i(KX-ΩT ) + Re -i(KX-ΩT ) , X <0, (8.1a)

η 0 = T e i(KX-ΩT ) , X >L. (8.1b)
In terms of the scattering and radiation problems studied previously, the coefficients are given by

T = α -(0) + β -(0) = C T + C + , (8.2a) R = α + (L) + β + (L) = C R + C -. (8.2b)
Analytical formulas can be obtained from the results in §4 and §7. While scattering is negligible for freely floating buoys when there is no power extracted as shown in §7(b), for moving and energy-extracting buoys |R| is close to 1 when the detuning frequency is within the band gap 0 ≤ Ω ≤ 2Ω 0 , as shown in figures 10 and 11. In particular reflection increases with the energy extraction rate λ ′ g and with the length of the array, and is of course the greatest for fixed buoys, which is equivalent to λ ′ g → ∞. In comparison with fixed buoys, the detuning corresponding to the maximum reflection is shifted slightly from Ω = Ω 0 towards Ω = 0.

Energy extraction

Now the flow and buoy displacements are known. The period-averaged rate of energy extracted by the j th buoy is given, in physical form, by:

E j ≡ 1 2 ω 2 λ g |ζ (j) | 2 . (9.1)
The total energy extracted by N buoys is therefore

E = N j=1 E j = 1 d N j=1 E j d. (9.2)
Article submitted to Royal Society Since N ≫ 1, the above series can be approximated by an integral:

E ≈ 1 2d ω 2 λ g L 0 |ζ| 2 dx = 1 2 ρ gh 3 A 2 π d ′ ω ′ 2 λ ′ g L 0 |ζ 0 | 2 dx ′ . (9.3)
Dividing by the energy flux rate of the incident wave across the width d : 1 2 ρ √ gh 3 A 2 C ′ g d ′ , the efficiency of power absorption is found in terms of the buoy displacement:

E ≡ π d ′2 C ′ g ω ′ 2 λ ′ g L 0 |ζ 0 | 2 dx. (9.4)
As a check, the absorbed energy can also be calculated from the difference of the incoming and outgoing energy flux rates at x ∼ ±∞. From the radiation and scattering coefficients, we can find

E ≡ 1 -α -(0) + β -(0) 2 + α + (L) + β + (L) 2 .
(9.5)

It is shown the appendix that (9.4) and (9.5) are the same.

Article submitted to Royal Society Figure 12 shows how the extracted energy varies with the detuning, the length of the array and the energy extraction rate. Note first that when the detuning frequency is outside the band gap of the pure scattering problem (0 < Ω/Ω 0 < 2), the energy extraction efficiency E tends to the same constant value predicted by the approximate value E 0 given by (3.12) based on Froude-Krylov approximation.

Ω/Ω0 R, T R T (a) λ ′ g = 1/10 (b) λ ′ g = 1/2 -4 -2 0 2 4 6 0 0.2 0.4 0.6 0.8 1 Ω/Ω0 R, T R T -4 -2 0 2 4 6 0 0.2 0.4 0.6 0.8 1 Ω/Ω0 R, T R T (c) λ ′ g = 1 (d) λ ′ g = ∞ (
Article submitted to Royal Society The good agreement shows the robustness of the asymptotic theory here despite its intended realm in the small neighbourhood of Bragg resonance.

It is clear that within the band gap, the extraction efficiency is substantially reduced by Bragg resonance. The reduction is the greatest around Ω = 0 (perfect tuning) and is larger for wider arrays (see figure 12(a)) and greater extraction rate (see figure 12(b)).

Although the maximum efficiency over most of the detuning frequencies is approximately equal to E 0 which increases with the length of the array, it is interesting that shorter arrays can yield more energy for small detuning, as shown in figure 13.

Conclusions

We have developed an asymptotic theory for wave interaction with a periodic array of small buoys in order to gain a better understanding of its potential as a wavepower farm.

We showed that fixed cylindrical buoys of small dimensions produce scattering effect of the same order of magnitude as vertical piles of the same radius but extending across the entire depth. It can easily be shown that the results derived here can be extended to small buoys of arbitrary geometry. In particular, we have shown that hemi-spherical buoys have the same influence as circular cylindrical buoys. We have extended the work of Li and Mei (2007a) and deduced analytically the frequency band gap within which one dimensional wave propagation is inhibited.

By multiple-scale analysis we have further solved the radiation problem of an array of movable buoys partially constrained by energy absorbing devices. The absorber is modeled as a damping force proportional to the velocity of the buoy. We have shown that Bragg resonance reduces the potential for energy extraction, somewhat similar to viscous damping in wall boundary layers studied recently by [START_REF] Tabaei | Viscous effects on Bragg scattering of water waves by an array of piles[END_REF]. While the present theory is designed only for the immediate vicinity of Bragg resonance, it agrees with the approximate theory valid far outside the band gap. Therefore it may be a practical tool for analyzing random incident waves with a broad frequency band. Modifications for shorter waves satisfying the Bragg condition kd = nπ appears straightforward. For practical design of arrays of energy converters with complex phase control of the power-takeoff system, numerical techniques would be necessary. Future extensions may include oblique incidence on a wide array. In this case, the wave physics of multiple scattering is two dimensional. Approximation used here can be extended as in Li and Mei (2007a); [START_REF] Tabaei | Viscous effects on Bragg scattering of water waves by an array of piles[END_REF] for piles.

This research has been supported by a grant from MASDAR Institute of Science and Technology in the program of MIT-Abu Dhabi Alliance. CCM also received partial support from US -Israel Bi-National Science Foundation and from Earth Systems Initiative at MIT.

Appendix A. Energy extraction from a sparse array

To prove that the two expressions of the energy extracted given in §9 are equal, let us recall the governing equations of the scattering and radiation amplitudes:

C ′ g ∂ α + ∂X = -iΩ 0 α + + α -+ iΩ α + , (A 1a) -C ′ g ∂ α - ∂X = -iΩ 0 α + + α -+ iΩ α -, (A 1b) C ′ g ∂ β + ∂X = -iΩ 0 (1 -G)( β + + β -) -G( α + + α -) + iΩ β + , (A 1c) -C ′ g ∂ β - ∂X = -iΩ 0 (1 -G)( β + + β -) -G( α + + α -) + iΩ β -. (A 1d)
The corresponding buoy displacement is given by: Similarly, by taking (A 1a) + (A 1b) + (A 1c) + (A 1d), we find

ζ 0 = G α + + α -+ β + + β -. ( A
C ′ g d ξ dX = -i (2Ω 0 (1 -G) -Ω) ζ 0 .
Denoting their complex conjugates of ζ 0 , ξ by ζ † 0 , ξ † respectively, we get

C ′ g d ζ † 0 ξ dX = -i (2Ω 0 (1 -G) -Ω) | ζ 0 | 2 -iΩ| ξ 2 |, and
C ′ g d ζ 0 ξ † dX = i 2Ω 0 1 -G † -Ω | ζ 0 | 2 + iΩ| ξ 2 |,
which can be summed up to give

C ′ g d ζ † 0 ξ + ζ 0 ξ † dX = 2iΩ 0 G -G † | ζ 0 | 2 .

Article submitted to Royal Society

Recall that

G = 1 1 -iω ′ λ ′ g , so that G -G † = 2iω ′ λ ′ g 1 + ω ′2 λ ′2 g . Since ζ † 0 ξ + ζ 0 ξ † = -2|G| 2 α -+ β -2 -α + + β + 2 ,
we finally get

d dX α -+ β -2 -α + + β + 2 = 2Ω 0 λ ′ g ω ′ | ζ 0 | 2 .
and, after integration,

2Ω 0 λ ′ g ω ′ C ′ g L 0 | ζ 0 | 2 dX = Ω 0 α -+ β -2 -α + + β + 2 L 0 . ( A

Figure 1 .

 1 Figure 1. A sparse array of buoys attached to energy absorbers.

Figure 2 .

 2 Figure 2. Energy extraction efficiency E0 vs. k ′ for a square array with d ′ = 1. (a) Effects of the extraction rates λ ′ g = 0.1, 0.5, 1, 2 for L = 1 (b) Effects of complete optimization for every frequency using (3.13) for L = 0.25, 0.5, 1, 2.

Figure 3 .

 3 Figure 3. A unit cell surrounding a buoy.

Figure 4 .

 4 Figure 4. Dispersion relation in an infinite periodic array. Real values are plotted in solid lines, and imaginary values in dashed line.

Figure 5 .

 5 Figure 5. Coefficients of reflection and transmission for fixed buoys as functions of the detuning frequency, for two array lengths.

Figure 6 .

 6 Figure 6. Wave amplitude inside the buoy array (0 < X < L) around Bragg resonance for d ′ = 1 and Ω = Ω0 (inside the band gap) and Ω = -2Ω0 (outside of the band gap).

Figure 8 .

 8 Figure 8. Free surface elevation η0 and buoy displacement ζ0 in the array. L = 1, d ′ = 1 and λ ′ g = 0.1, 0.5, 1 as marked next to the curves. In figures (a) and (b), the dashed curve represent the limiting case of fixed buoys, i.e. λ ′ g → ∞.As shown in §7(b), the free surface amplitude is constant and equal to 1 in the limiting case of free buoys λ ′ g = 0.

Figure 9 .

 9 Figure 9. Effects of detuning frequencies and array length on free surface elevation η0 and buoy displacement ζ0 for λ ′ g = 0.5 and d ′ = 1. L = 0.25 (diamond), L = 0.5 (square) and L = 1 (triangle).

Figure 10 .

 10 Figure 10. Reflection and transmission coefficients for different L. The energy extraction rate is set at λ ′ g = 1/2 and d ′ = 1.

E

  Figure 12. (a) Efficiency over a range of Ω/Ω0 for d ′ = 1. (a): λ ′ g = 1/2 and L = 0.1, 0.25, 0.5. (b) : L = 1/2 and λ ′ g = 0.1, 0.5, 1.0. The horizontal lines give the results obtained by (3.12).

Figure 13 .

 13 Figure 13. Influence of L and λ ′ g on the the minimum energy extracted due to Bragg resonance. d ′ = 1 and λ ′ g = 0.25, 0.5, 1, 2.

  ξ = G α + + β +α -+ β -.

|

  4)Recall that the boundary conditions (4.29) and (7.5) giveα + (0) = 1, α -(L) = β + (0) = β -ζ 0 | 2 dX = α -(L) + β -

Article submitted to Royal Society

Because of the coupling by (7.8), only two of the four coefficients C (i) 1 , i = 1, 2, 3, 4 are independent. Two relations among them can be found by invoking (7.8) at any X ∈ [0, L], say X = 0,

From the two boundary conditions (7.5) , we get

C

(3) 1 e iKRL + C

(4) 1 e -iKRL = -iB - 1 e iKSL . (7.14d)

Article submitted to Royal Society