Fundamental formulae for wave-energy conversion

The time-average wave power that is absorbed from an incident wave by means of a wave-energy conversion (WEC) unit, or by an array of WEC units—i.e. oscillating immersed bodies and/or oscillating water columns (OWCs)—may be mathematically expressed in terms of the WEC units' complex oscillation amplitudes, or in terms of the generated outgoing (diffracted plus radiated) waves, or alternatively, in terms of the radiated waves alone. Following recent controversy, the corresponding three optional expressions are derived, compared and discussed in this paper. They all provide the correct time-average absorbed power. However, only the first-mentioned expression is applicable to quantify the instantaneous absorbed wave power and the associated reactive power. In this connection, new formulae are derived that relate the ‘added-mass’ matrix, as well as a couple of additional reactive radiation-parameter matrices, to the difference between kinetic energy and potential energy in the water surrounding the immersed oscillating WEC array. Further, a complex collective oscillation amplitude is introduced, which makes it possible to derive, by a very simple algebraic method, various simple expressions for the maximum time-average wave power that may be absorbed by the WEC array. The real-valued time-average absorbed power is illustrated as an axisymmetric paraboloid defined on the complex collective-amplitude plane. This is a simple illustration of the so-called ‘fundamental theorem for wave power’. Finally, the paper also presents a new derivation that extends a recently published result on the direction-average maximum absorbed wave power to cases where the WEC array's radiation damping matrix may be singular and where the WEC array may contain OWCs in addition to oscillating bodies.


Summary
The time-average wave power that is absorbed from an incident wave by means of a wave-energy conversion (WEC) unit, or by an array of WEC units-i.e. oscillating immersed bodies and/or oscillating water columns (OWCs)-may be mathematically expressed in terms of the WEC units' complex oscillation amplitudes, or in terms of the generated outgoing (diffracted plus radiated) waves, or alternatively, in terms of the radiated waves alone. Following recent controversy, the corresponding three optional expressions are derived, compared and discussed in this paper. They all provide the correct time-average absorbed power. However, only the first-mentioned expression is applicable to quantify the instantaneous absorbed wave power and the associated reactive power. In this connection, new formulae are derived that relate the 'added-mass' matrix, as well as a couple of additional reactive radiation-parameter matrices, to the difference between kinetic energy and potential energy in the water surrounding the immersed oscillating WEC array. Further, a complex collective oscillation amplitude is introduced, which makes it possible to derive, by a very simple algebraic method, various simple expressions for the maximum time-average wave power that may be absorbed by the WEC array. The real-valued time-average absorbed power is illustrated as an axisymmetric paraboloid defined on the complex collective-amplitude plane. This is a simple illustration of the so-called 'fundamental theorem for wave power'. Finally, the paper also presents a new derivation that extends a recently published result on the direction-average maximum absorbed wave power to cases where the WEC array's radiation damping matrix may be singular and where the WEC array may contain OWCs in addition to oscillating bodies.

Introduction
For a general three-dimensional case, the basic linearized theory for conversion of ocean-wave energy by means of one oscillating body was developed in the mid-1970s [1][2][3][4]. The starting point was to consider power input as the product of the net wave force and the body's oscillation velocity. In addition, Newman [1], based on some reciprocity relations, discussed how the absorbed  respectively, where u i is the complex velocity amplitude for oscillation mode i. The complex proportionality coefficients f e,i = f e,i (β, ω), i.e. the excitation-force coefficient, and Z ii = Z ii (ω), i.e. the radiation impedance, as well as the latter's real and imaginary parts, R ii = R ii (ω), i.e. the radiation resistance, and X ii = ωm ii = X ii (ω) = ωm ii (ω), i.e. the radiation reactance, are functions of ω. The coefficient m ii is called 'added mass' although it may be negative in exceptional cases [19]! Moreover, the coefficient f e,i also depends on β, the angle of wave incidence. Observe that, in terms of complex amplitudes, the radiation-force F r,i has two components, an active component and a reactive one, F r,i,act = −R ii u i and F r,i,react = −iX ii u i = −m ii iωu i , (3.2) which are in phase with the velocity u i and the acceleration iωu i , respectively. By inverse Fourier transformation, where products in the frequency domain correspond to convolutions in the time domain, we may find a corresponding decomposition of the general, time-domain, radiation force [20]: F r,i,t (t) = F r,i,t,act (t) + F r,i,t,react (t). (3.3) As shown below, only the active force component contributes to the net time-averaged energy transfer, while the reactive force component serves temporary energy exchange between differently sized stores of kinetic energy and potential energy.
To provide the desired immersed-body motion, the WEC unit needs to be equipped with a machinery for control and power take-off (PTO). This provides an additional force F pto,i,t (t), with corresponding complex amplitude F pto,i for the monochromatic-wave case. Then we may write the equation of motion, in complex-amplitude representation, as (R ii + r loss,i ) + iω m ii + m i − c i ω 2 u i = F e,i + F pto,i = f e,i (β)A + F pto,i , (3.4) where m i is the mass of the immersed body and c i its hydrostatic stiffness coefficient. We have also introduced a coefficient r loss,i to represent linear power loss. Introducing the body's excursion from equilibrium position s i,t (t)-thus u i,t (t) =ṡ i,t (t)-we may, in time-domain representation, write the equation of motion as {−F r,i,t,act (t) + r loss,iṡi,t (t)} + {−F r,i,t,react (t) + m isi,t (t) + c i s i,t (t)} = F e,i,t (t) + F pto,i,t (t). (3.5) In contrast to the frequency-domain model, for the time-domain model we may include possible additional nonlinear forces in the F pto,i,t (t) term of (3.5).
where (d/dt)W i,t,body (t) = m iui,t (t)u i,t (t) + c i s i,t (t)ṡ i,t (t) = (d/dt){m i u 2 i,t (t) + c i s 2 i,t (t)}/2 is the time derivative of the sum of the body's kinetic energy and potential energy, and where (d/dt)W i,t,water (t) is the time derivative of the sum of kinetic energy and potential energy of the water surrounding the body. At instants when u i,t (t) = 0, there is no kinetic energy, and at instants when s i,t (t) = 0, there is no potential energy. Except for conditions of resonance, the r.h.s. of (3.7) does not vanish at all instants. In general, the PTO machinery has to cope with reactive forces and reactive power, because of unequal magnitudes of the kinetic and potential energy stores. The r.h.s.-and hence also the l.h.s.-of (3.7) has, however, a vanishing time average.
For a sinusoidal oscillation with complex velocity amplitude u i = |u i | exp(iϕ u i ), we have {u i,t (t)} 2 = |u i | 2 cos 2 (ωt + ϕ u i ) = |u i | 2 {1 + cos(2ωt + 2ϕ u i )}/2 and {s i,t (t)} 2 = |u i /ω| 2 sin 2 (ωt + ϕ u i ) = |u i /ω| 2 {1 − cos(2ωt + 2ϕ u i )}/2. Using this, we find (3.8) which, together with (3.7), explicitly shows how the reactive power is directly related to the difference between the maximum values of kinetic energy and potential energy. In analogy with (3.8), (d/dt)W i,t,water (t) is related to such an energy difference associated with the water surrounding an array of immersed WEC units. This matter is discussed in §6.2,6.3 and more extensively in appendix B.
In the remaining part of this paper, we consider only a monochromatic wave and the corresponding sinusoidal oscillation of immersed WEC units. Without considering the details of the PTO machinery, we shall rather consider the WEC units' complex oscillation amplitudes-for instance u i -to be independent variables, and a goal of our analysis is to find their optimum values corresponding to the incident wave as given by (2.1). Then the complex-amplitude version of (3.6) is where the asterisk ( * ) denotes complex conjugate. Assuming ideal conditions, we set r loss,i = 0 in the following. Moreover, we shall find it convenient to make the following substitutions: and |U| 2 = UU * = R ii u i u * i 2 = P r , (3.10) where the non-negative quantity P r represents the radiated wave power (caused by any forced oscillation of the immersed body). Although the introduced complex quantity U = √ P r e iδ might, in general, have any arbitrary phase angle δ in the interval −π < δ ≤ π , we shall find it convenient that it is chosen to have the same phase angle as A * E * (β). Then A * E * (β)/U is a real positive quantity, which, notably, is independent of the complex velocity amplitude u i .
We may now simplify (3.9) for the time-averaged absorbed wave power to P a = P e − P r = AE(β) + A * E * (β) − |U| 2 , (3.11) where P e = AE(β) + A * E * (β) (3.12) is the 'excitation power'. An important motivation behind the substitution of f e,i (β)u * i /4 by E(β) and R ii u i u * i /2 by |U| 2 is that the above (3.11) and (3.12) as well as the following (3.13)- (3.17) are applicable also for a WEC array consisting of several WEC units-oscillating bodies and/or OWCs-provided the parameters E(β) and |U| 2 are properly redefined, as explained later in this paper (see (6.20) and (7.1) and (7.2)). For this reason, we propose the terms 'collective excitation-power coefficient' and 'collective oscillation amplitude' for the complex quantities E(β) and U, respectively.
The usefulness of introducing the quantity U is that (3.11) may be rewritten as from which we, simply by inspection, see that the first term equals the maximum possible absorbed power, provided the last term vanishes, that is, if the quantities U and E(β) have optimum values U 0 and E 0 (β) that satisfy the optimum condition (3.14) Hence, at optimum, the three terms on the r.h.s. of (3.11) have the same, real and non-negative, magnitude. The last one of the three terms is the optimum radiated power. It follows that we have several different alternative expressions for the maximum absorbed power, e.g.
We may consider this series of alternative mathematical expressions as reciprocity relations for the maximum absorbed power. For instance, the maximum absorbed power P a,MAX equals the optimum radiated power P r,OPT = |U 0 | 2 . In (3.11), the last term, the radiated-power term P r = |U| 2 , appears to be a power-loss term, but it should, rather, be considered as a necessity, because the radiated wave is needed to extract power from-that is, to interfere destructively with-the incident wave.
As we have chosen A * E * (β)/U to be a real positive quantity, which is independent of u i and, therefore, also of U, we have where we have made use of the optimum condition (3.14). In general, we shall consider U to be an independent complex oscillation-state variable, while the optimum value U 0 (β) is real and positive, because we have chosen U to have the same phase angle as A * E * (β) has. According to (3.16), A * E * (β) = U * 0 U and AE(β) = U 0 U * . If we insert this into (3.11) and also use (3.15), we obtain the simple equation which, for a fixed value P a < P a,MAX , corresponds to the equation of a circle of radius P a,MAX − P a centred at U 0 (β) in the complex U plane. Equation (3.17) may be illustrated as an axisymmetric paraboloid in a diagram where a vertical real P a axis is erected on a horizontal complex U plane, as shown in figure 1.
Assuming that the PTO machinery of the oscillating WEC body contains sufficient control equipment to achieve the desired oscillation, we may consider the complex velocity amplitude u i , as well as U, to be an independent variable. However, the optimum value u i0 , as well as U 0 , depends on the excitation force f e,i (β)A, and is, consequently, dependent on the incident-wave parameters A and β. For an optimum oscillation velocity u i = u i0 (β), say, corresponding to maximum absorbed wave power P a = P a,MAX -cf. e.g. (3.14) and (3.15)-we have optimum collective parameters E(β) = E 0 (β) and U = U 0 . Note that |U 0 | 2 , as well as all the other alternative expressions given in (3.15) for the maximum absorbed power, depends implicitly on β, and, moreover, it is proportional to |A| 2 -remembering that u i0 (β), in contrast to an independent variable u i , is proportional to A.  Figure 1. The wave-power 'island' , illustrating (3.17). Absorbed wave power P a as a function of the complex collective oscillation amplitude U = Re{U} + i Im{U} = |U| e iδ , where the phase δ = arg{U} is chosen such that AE/U * is a real positive quantity, and where |U| is given by (3.10) for the one-mode oscillating-body case, and by (6.20) for the case of a general WEC array. The largest possible absorbed wave power P a,MAX is indicated by a star on the top of the axisymmetric paraboloid, and U 0 is the optimum collective oscillation amplitude. Colour changes indicate levels where P a /P a,MAX equals 0, 1 4 , 1 2 and 3 4 . (a) Side view, (b) top view and (c) inclined view.

Destructive far-field wave interference
power by analysing wave interference in the far-field region, that is, many wavelengths away from the immersed WEC body, or more generally, WEC array. We then have to assume that the wave propagation takes place in water that we may consider to be an ideal loss-free fluid. In the following, we shall apply such a GPV, and then compare the results with the above LPV results. If we assume that, within a localized area near our chosen origin (x, y, z) = (0, 0, 0), the body-or, more generally, a WEC array-is installed, but not oscillating, then an incident plane, monochromatic wave, with wave elevation η 0 = η 0 (r, θ ) = A e −i{kr cos(β−θ)} , as given by (2.1), produces a diffracted wave, for which the wave elevation has a complex amplitude η d = η d (r, θ), say. (This may also include diffraction effects of possible reefs and rocks.) Next, imagine that the immersed body, or the array, is performing forced oscillations with no incident wave, that is, with A = 0. Then a wave will be radiated from the body, or the array. Let η r = η r (r, θ ) denote the complex elevation amplitude of this radiated wave.
When there is an incident wave, and the immersed WEC array is oscillating, then the complex amplitude of the wave elevation is η = η 0 + η g = η 0 + η d + η r . We have here introduced η g = η d + η r as the complex elevation amplitude of the total 'outgoing' wave. Depending on the geometrical details of the array, the outgoing waves may have a complicated mathematical structure near the array, in the so-called near-field region. We shall, however, only need far-field mathematical details, which, in general, are asymptotically valid several wavelengths away from the WEC array. The far-field diffracted/radiated/outgoing wave elevations may thus be expressed in the form η d/r/g = − iω g C d/r/g (θ )(kr) −1/2 e −ikr + · · · , as kr → ∞, (4.1) where the complex functions C d (θ ), C r (θ ) and C g (θ) are the far-field coefficients for the diffracted wave, the radiated wave and the outgoing wave, respectively. It is convenient to express these coefficients in terms of the so-called Kochin functions Note that the diffracted wave is linearly related to A, the complex amplitude of the incident-wave elevation amplitude at the origin, r = 0, while the radiated wave is linearly related to all WEC units' oscillation amplitudes. WEC arrays will be discussed in §6. At present, we shall consider the simpler case of only one immersed, single-mode oscillating, body. Introducing complex Kochin function coefficients of proportionality by corresponding lower case symbols, we may write the Kochin functions as for the diffracted wave and the radiated wave, respectively. The total generated wave's Kochin function is Note that for an optimum oscillation vector u i = u i0 (β), there corresponds optimum Kochin functions H r (θ ) = H r0 (θ ) and H g (θ ) = H g0 (θ ), which depend, implicitly, also on β. In particular, u i0 (β) and, thus, H r0 (θ ) are linearly related to the excitation force F e,i (β) = f e,i (β)A. However, the coefficient h i (θ) does not depend on β, in contrast to the coefficient h d (θ ), which depends implicitly on β, since the diffracted wave is a response to the incident wave.
In correspondence with our derivation of (3.9) for the wave power absorbed by an immersed oscillating body as the product of the net wave force and the body's oscillation velocity, Newman [1, §10] expressed the power P a absorbed by an oscillating immersed body as an integral over the body's wetted surface, where the integrand is the hydrodynamic pressure multiplied by the normal component of the fluid velocity. Then, applying Green's theorem, he expressed P a as an integral over an, envisaged, cylindrical control surface in the far-field region, a surface that encloses the immersed body and all water between the body's wetted surface and the control surface. In this way, Newman [1, eqns 58 and 59] moved from the LPV to the GPV, and expressed the absorbed wave power in terms of Kochin functions.
Accordingly, following Newman, we may write the time-averaged absorbed wave power as P a = P i − P g = I(β) + I * (β) − |G| 2 (4.5) where P i = I(β) + I * (β) is the 'input power' and is the (non-negative) total 'outgoing power'. Here v p = ω/k and v g = dω/dk are the phase velocity and the group velocity, respectively. Moreover, and we may write the input power as An approach corresponding to (4.5)-(4.8) has been applied by Farley [7,15] and Rainey [17]. Their approach shows the physical details of wave-interference energy removal in the far-field region. By wave interference in the far-field region, wave-energy removal takes place where the outgoing wave η g travels in the same direction as the incident wave η 0 , that is, for direction θ coinciding with the incident-wave direction β. As there is no energy exchange between two plane waves propagating in different directions, there is no contribution to far-field wave-energy removal by the outgoing wave in directions where θ = β (or, more precisely, outside a small θ interval around θ = β, an interval that tends to zero as kr → Noting that (4.5) has a similar mathematical structure as (3.11), it might appear that (3.14)- (3.17) are valid also if we replace the LPV parameters P e , P r , E(β) and U by the GPV parameters P i , P g , I(β)/A and G, respectively. However, note that the optimum LPV parameters are more directly related to the optimum WEC body oscillations than the GPV ones are. If we compare LPV equations (3.10)-(3.12) with GPV equations (4.5)-(4.8), we may note that P e is proportional to A and linearly related also to the WEC body's oscillation amplitude, while P r is quadratically related to this amplitude, but independent of A. By contrast, P i , as well as P g , is related in a more complicated way to A and to the WEC body's oscillation amplitude. Equations (3.14)-(3.17) therefore do not apply for the GPV parameters.
We may mitigate this drawback by rearranging the GPV equations (4.5)-(4.8) as follows. Firstly, we observe that if the single-mode oscillating WEC body does not oscillate, i.e. u i = 0, then no wave energy is being absorbed, i.e. P a = 0. Moreover, the radiated wave's Kochin function vanishes, i.e. H r (θ) = 0. Then the GPV equations (4.4)-(4.8) agree with the following reciprocity relation for the diffracted wave's Kochin function [1, eqn 33]: Secondly, from the same GPV equations (4.4)-(4.8), we then find, for the oscillating-body case (i.e. u i = 0), that the power P a , which is removed by the far-field wave interference, is as given by the LPV(!) equations (3.10)-(3.12), but now with collective parameters |U| 2 and E(β) expressed in terms of far-field quantities, namely, and Note that |U| 2 , in contrast to |G| 2 , is independent of the wave amplitude A, and quadratic in the WEC body's oscillation amplitude. Further, E(β)A, in contrast to I(β), is linearly related to the oscillation amplitude, and proportional to the incident wave amplitude A. Equations (3.10) and (4.10) present two different expressions for the radiated power |U| 2 . Physically, this means that the power which is radiated from the WEC body's wave-interacting surface into a lossless fluid equals the power that is associated with the radiated wave in the far-field region of the fluid. From a physical point of view, what may be observed in the far-field region is a superposition of the plane incident wave and the outgoing wave. As observed in the far-field region, one may not know whether the outgoing wave originates from one single-mode oscillating body or from an array consisting of many WEC units. For this reason, all equations in the present section are valid for this latter WEC system, provided the two-factor product h i u i that appears in (4.3) is generalized to a sum of such products, one product for each of the WEC array's oscillating modes. Details are given in §6.

Relationships between radiated and diffracted waves
Among many water-wave reciprocity relations, there are two relations, which relate diffraction and radiation parameters, and about which Newman [1, eqns 45 and 48] remarked that the corresponding two involved physical problems that 'are not physically related to each other in any obvious manner'. Newman used these two relations to convert the formula for absorbed wave power from the version of (4.5) to the version of (3.11), but with the collective parameters |U| 2 and E(β) expressed solely in terms of radiation Kochin function coefficients, that is, without the diffraction Kochin function-which we still need to eliminate from (4.11).
The second one is a relation between H d (θ ) and h i (θ), a relation which Newman [1, eqn 61] used to simplify (4.11) to  In this way, Newman succeeded to eliminate, mathematically, the diffracted wave's Kochin function H d = h d A, that appears for example in (4.5)-(4.8) and in (4.11). Referring to (3.10), we may, however, arrive at the same result (5.2) without referring to the second one of the above-mentioned two, not very obvious, reciprocity relations. The result simply follows by applying the Haskind relation (5.1) to the excitation-force coefficient f e,i (β) in (3.10).
Although reciprocity relations between diffraction and radiation parameters connect different physical problems 'which are not physically related to each other in any obvious manner', as admitted by Newman [1, §7], the Haskind relation (5.1) may be supported by the following physical argument. Imagine a non-symmetric WEC, e.g. the well-known nodding-duck device [23], which is installed with an optimum orientation to absorb waves arriving from west, thus incident waves propagating eastwards. If, in a case with no incident wave, the device is performing forced oscillations, the device will primarily radiate waves propagating westwards. Thus, the addition of an angle π in the argument on the r.h.s. of (5.1) seems reasonable. Moreover, it is reasonable that the excitation-force coefficient f e,i (β) of the incident-wave problem (diffraction problem) is proportional to the radiation-ability coefficient h i (β + π ) of the forced-oscillation problem (radiation problem). Admittedly, however, the second one of the two above-mentioned reciprocity relations, which directly connects (4.11) and (5.2), is less obvious from a physical point of view, namely the reciprocity relation presented by Newman [1, eqn 48]: Earlier, the LPV quantities |U| 2 and E(β) appearing in (3.11) were given by the two equations (3.10). However, when we now have, alternatively, expressed E(β) by equation (5.2) and |U| 2 by equation (4.10), which are far-field, or global, equations, this corresponds to a mixed, or hybrid, global-local point of view (GLPV), because we have now expressed the LPV parameters |U| 2 and E(β) in terms of radiation Kochin functions, which are far-field parameters. Although it is not easy to give the GLPV version a direct physical interpretation, it has the advantage that it may be a basis for several reciprocity relations [1] and, moreover, also for certain mathematical derivations below, as exemplified later on in this paper; see §6.1,6.2.

Generalization to wave-energy converter arrays
We consider a case of wave-energy absorption by an array of immersed oscillating rigid bodies and of OWCs, as indicated in figure 2. Let us assume that the number of wave-interacting oscillators is N = N u + N p , where N p is the number of OWCs and N u is the number of used body modes, whose number may be up to six times the number of bodies. The oscillation state and the excitation due to an incident plane wave may be described by N-dimensional column vectors v and x, respectively, where Here we have introduced two N u -dimensional column vectors u = [u 1 u 2 u 3 · · · u N u ] T and F e = [F e,1 F e,2 F e,3 · · · F e,N u ] T , where u i and F i are the complex amplitudes of the oscillation velocity and of the excitation force for rigid-body oscillation mode i. Correspondingly, we have introduced two N p -dimensional column vectors p = [p 1 p 2 p 3 · · · p N p ] T and Q = [Q e,1 Q e,2 Q e,3 · · · Q e,N p ] T , where p i and Q i are the complex amplitudes of the oscillating dynamic air pressure and of the excitation volume flow for OWC i. We may think of v and x as vectors in an N-dimensional complex space. The superscript 'T' denotes the transpose of a matrix, and the complex conjugate transpose of a matrix is correspondingly denoted by the dagger symbol ( †). Applying linear theory, we have also introduced the following N u -and N p -dimensional vectors: the complex vectorial proportionality excitation vector coefficients f e = f e (β) and q e = q e (β), respectively. The complex excitation vector x = x(β), acting on the WEC array, depends on the angle β of wave incidence, and it is proportional to the complex elevation amplitude A of the undisturbed incident wave at the origin (x, y) = (0, 0). We shall, however, here consider the complex oscillation vector v to be an independent variable, assuming that we have an ideal machinery for PTO and motion control.

Global point of view
For this, rather general, WEC array, we may extend the second equation of (4.3) to write the radiated wave's Kochin function as We may then generalize H r correspondingly in, for example (4.4), (4.6)-(4.8), (4.10) and (4.11). Moreover, the Haskind relation (5.1) is generalized to and (5.2) to Note that this latter equation may be considered as a generalized Haskind relation for the collective excitation-power coefficient E(β). Moreover, on the basis of the general equations (4.11) and (6.5), we easily see that (5.4)-the least obvious one of the reciprocity relations presented by Newman [1, eqn 48]is still valid with the general radiation Kochin function H r , as given by (6.2). Equations (4.11) and (6.5) provide two different mathematical relations between the LPV collective excitation-power coefficient E(β) and the N GPV Kochin function coefficients h i for the radiated wave. Equation (4.10), with (6.2), provides a mathematical relation between these h i coefficients and the LPV collective amplitude |U|, and also the corresponding complex amplitude U if we remember that we have chosen the phase of U to equal the phase of A * E * (β), in accordance with (3.16). Let us next consider the optimum case for maximum absorbed power. Algebraic procedures for determining the optimum value v 0 = [v 10 v 20 v 30 · · · v N0 ] T of the complex oscillation-state vector v are treated in more detail in appendix A. Correspondingly, according to (6.2), there exists an optimum Kochin function for the radiated wave. Note that, even if we, in general, consider v to be an independent variable, the optimum value v 0 = v 0 (β), as well as the β-dependent optimum Kochin function H r0 (θ), is linearly related to the incident wave amplitude. From the optimum condition (3.14), we have AE 0 (β) = |U 0 | 2 = A * E * 0 (β), which, in combination with (4.10) and (4.11), gives the condition which the optimum radiated wave's Kochin function H r0 (θ) needs to satisfy. Combining this condition with the reciprocity relation (5.4)-see also (6.4)-yields That this is real and positive corresponds to the radiated wave having optimum phase. If we choose A to be real and positive, then alsoH r0 (β + π ) has to be real and positive. In (3.15), we presented several different expressions for the maximum power P a that is possible to be absorbed by the WEC array. We shall find it convenient to add also the following expressions: Applying the last one of the fractions shown in the LPV equations (6.9), and then inserting from the GPV equations (4.10) and (6.5), we get where J w is the wave-power level, as given by (2.2), and d a ≡ P a /J w is the 'absorption width'. Moreover, we have introduced the-at optimum-gain function which is an extension of a formula presented, independently, by Newman [25] and by Evans [8] for the single-body, one-mode case. It is remarkable that we here have been able to express the maximum absorbed power in terms of optimum far-field Kochin functions for the radiated wave only. It should be emphasized that this gain function G g0 (β) applies only to the optimum case for maximum absorbed wave power. The Haskind relation (5.2) and the collective Haskind relation (6.5), as well as the optimum gain function as given in (6.11), indicate that it is important for a WEC system to possess the ability to radiate a wave propagating in a direction opposite to the direction of the incident wave. This ability is represented, quantitatively, by the N coefficients h i (β + π ) in (6.5). If N = 1, we see from (4.3) and (5.3) that, with i = 1, we have |H r (β + π )| 2 = |h 1 (β + π )v * 1 | 2 = |h 1 (β + π )v 1 | 2 = |H r (β + π )| 2 . However, with (4.3) and (5.3) generalized to (6.2) and (6.4), we should note that, in general, |H r (β + π )| 2 = |H r (β + π )| 2 for N ≥ 2. For instance, referring to (6.2)-(6.4) for N = 2, we have, for any θ, which, in general, deviates from zero for arbitrary as well as for optimum values of the complex velocity amplitudes v 1 and v 2 . For the circularly oscillating Evans Cylinder [2], we may, as shown below in §7.1see (7.9)-replace |H r0 (β + π )| 2 by |H r0 (β)| 2 in the numerator of (6.11). A corresponding replacement may be made if diffraction is negligible or, otherwise, in cases where the integral on the r.h.s. of (5.4) vanishes when H r (θ ) = H r0 (θ ).
Concerning the GPV discussion, presented in § §4 and 5, for the case of one single-mode body WEC unit, we have, so far, here in §6, extended results to the case of an array of WEC units. Our next task will be to generalize some of the LPV matter discussed in §3.

Local point of view
For a single-mode oscillating body, the complex amplitude of two wave-force contributions, the excitation force F e,i and the radiation force F r,i , are given by (3.1). For our WEC array, the excitation vector x, introduced by the second equation of (6.1), is an extension of F e,i , while [10] is an extension of F r,i , where Z and Y are the N u × N u radiation-impedance matrix for the oscillating bodies and the N p × N p radiation-admittance matrix for the OWCs, respectively. These matrices are symmetric, that is, between the oscillating bodies and the OWCs, which compose the WEC array. It is convenient to split these complex matrices into real and imaginary parts: where the radiation resistance matrix R, the radiation reactance matrix X, the radiation conductance matrix G, the radiation susceptance matrix B, as well as the matrices C and J, are real. All these matrices are frequency dependent. Further comments concerning these matrices are given in appendix B; see (B 24)-(B 28) and related text.
Let us now, for the WEC array, extend (3.1) and (3.2), where we defined the radiation force F r,i and split it into active and reactive components. The extension reads where We may note that D = D † , (6.16) which means that the radiation-damping matrix D is hermitian. Also the matrix (D reactive /i) is hermitian. From this it follows that, for any N-dimensional complex column vector v, the scalar matrix products v † Dv and v † D reactive v are real and purely imaginary, respectively. (We may observe that the matrices D and (D reactive /i) are real and symmetric if the WEC array contains no OWCs or no oscillating bodies, that is, in cases where D complete = Z = R + iX or D complete = Y = G + iB, respectively.) If we premultiply (6.14) by −v † /2, we get the 'complex radiated power' where the last term, the reactive-power term, v † D reactive v/2 is purely imaginary, while the first term, the radiated-power term, v † Dv/2 ≡ P r is real and non-negative-see (6.21). Moreover, if we premultiply by v † /2 the excitation vector x, defined by equation (6.1), we get the 'complex excitation power' We may note that the imaginary part Im{P e } represents reactive power (see (B 49)). For any oscillation vector v = [u − p] T , the time-average wave power absorbed by the array is P a = P e − P r , where the 'excitation power' P e and the radiated power P r are given by [10] We may express this in the form of (3.11) provided we define the collective excitation-power coefficient E(β) and the collective oscillation amplitude U by which is an extension of (3.10). We still choose the phase angle of U such as to make A * E * (β)/U a real and positive quantity. For a case with no incident wave, x = 0 (which means that P e = 0), energy conservation requires that the absorbed wave power P a = P e − P r = −P r = −v † Dv/2 cannot be positive. Thus, for all possible finite oscillation-state vectors v, we have v † Dv ≥ 0. (6.21) Thus, in general, the radiation damping matrix D is positive semidefinite. It is singular in cases when its determinant vanishes, It is well known [10] that the maximum wave power that can be absorbed by the array is where v 0 = v 0 (β) is an optimum value of the oscillation-state vector v that has to satisfy the optimum condition For a fixed value of the absorbed wave power P a , where P a < P a,MAX , equation (6.24) represents an 'ellipsoid' in the complex N-dimensional v space, C N -but reduced to an r D -dimensional v space, C r D , in cases where the radiation damping matrix D is singular and of rank r D < N-see (A 15). The centre of the 'ellipsoid' is at the point v = v 0 . The elliptical semi-axes are 2(P a,MAX − P a )/λ i for i = 1, 2, 3, . . . , r D , where λ i are the positive definite (non-zero) eigenvalues of the matrix D-cf. (A 2). The 'ellipsoid' that corresponds to P a = 0 runs through for example points v = 0 and v = 2v 0 . The degenerate 'ellipsoid' that corresponds to P a = P a,MAX is just one point, which represents the (unconstrained) optimum situation. Choosing smaller P a results in increased size of the 'ellipsoid'. If N = 1, then the 'ellipsoid' simplifies to a circle in the complex v 1 plane. Then, as v 10 = x 1 /2D 11 , we may, from the general equation (6.24), derive Note that a similar simple circle equation may be derived for cases where the radiation damping matrix D is singular and of rank r D = 1, although N ≥ 2; as exemplified by (7.25), for an axisymmetric system [26, eqn 37]. Considering how the absorbed power P a varies with v, the relationship (6.24) may be thought of as a 'paraboloid' in the complex N-dimensional v space, C N . The top point of this 'paraboloid' corresponds to the optimum, (v 0 , P a,MAX ). Here, N should be replaced by r D in cases where the radiation matrix D is singular.
The simple equation (3.17), which for a fixed P a represents a circle in the complex U plane, can be shown to be equivalent to (6.24) above, which represents an 'ellipsoid' in the complex v space, by making use of (3.16), (6.20) and (6.23). Starting from (3.17), we have (6.25) noting that the collective excitation-power coefficient E(β) is a scalar, and thus equals its own transpose, and recalling the hermitian property (6.16) of the radiation-damping matrix D.
The proof (6.25) also serves to demonstrate that, with the generalizations (6.20), equations (3.11)-(3.17) are valid not only for a single, one-mode oscillating body, but even for an array consisting of several WEC units-oscillating bodies, as well as OWCs. In particular, equation (3.17), as illustrated in figure 1, is applicable even to the general case of wave energy absorption by an array of oscillating bodies as well as OWCs. The involved physical quantities refer to the array objects' wave-interacting surfaces. Thus, if inverse Fourier transformation is applied to, for example (3.11) and (6.20), they may be applied to analyse the WEC array's wave-power absorption also in the case of non-sinusoidal time variation.
As long as we have not taken any equation of motion into account, we may here consider the components v i of the vector v to be independent variables-assuming that the WEC array contains sufficient control equipment to achieve the desired oscillations. As in the case of a single WEC unit oscillating in one mode (cf. §3), however, all optimum values v i0 , and thus the optimum column vector v 0 , depend on the excitation vector x(β), and are, consequently, dependent on the incident-wave parameters A and β.

Reactive radiation parameters
Let us now return to consider the reactive power corresponding to the radiated wave in a forced oscillation case, that is, without any incident wave. The imaginary part of the complex r.h.s. term in where we, in the last step, made use of (B 25), which applies for a case with no incident wave. Here, T − V is the time-average difference between kinetic energy and potential energy of the water surrounding the WEC array, or more precisely, the water in the near-field region. (Note that, in the far-field region, such an energy difference averages to zero.) Details are discussed in appendix B. For the case of a single-mode one-body WEC unit, (3.8) corresponds to a time-domain analogue of the frequency-domain equation (6.26), but with energy difference in the mechanical oscillating system itself, rather than in the surrounding water. For this one-mode case, (6.26) simplifies to where m ii = m ii (ω) = X ii (ω)/ω is the so-called 'added mass', a term which may appear confusing in particular cases when it shows up to be negative, that is, in potential-energy-dominating cases where T − V < 0 [19]. From (3.8), we may conclude that the time-average difference between kinetic energy and potential energy of the oscillating body itself is (m i − c i /ω 2 )|u i | 2 /4, which is positive only if ω > √ c i /m i . Resonance may occur for angular frequencies ω = ω 0 , which satisfy the equation ω 2 0 {m i + m ii (ω 0 )} = c i . At resonance, the WEC unit's PTO machinery need not exchange reactive power with the oscillating system.
The values of m ii or, more generally, of D reactive at infinite frequency are important for time-domain models. It is well known that the elements of matrix m = X/ω, in general, tend to finite, non-zero, constants at infinite frequency. The infinite-frequency behaviours of the other two matrices which make up the matrix D reactive are less well known, although Evans & Porter [27] observed that the radiation susceptance matrix B is zero at infinite frequency and Kurniawan et al. [28] reported that the real part, C, of the radiation coupling matrix has elements which tend to finite, non-zero, constants at infinite frequency. A physical explanation for these behaviours is given in the following.
Consider first an array of oscillating bodies with no OWCs, where one of the bodies is forced to oscillate harmonically with a unit velocity corresponding to mode i, in the absence of incident waves, while the other bodies are held fixed. As the oscillation frequency is increased to infinity, the acceleration also increases to infinity. The force required to move the body will necessarily also be infinite. There is therefore sufficient force to accelerate the fluid, which on the wetted body surface needs to move with the same velocity as the body. As the potential energy is zero in this limiting case of infinite frequency (as there are no radiated waves at infinite frequency), while the kinetic energy is positive (as the velocity of the fluid is finite), m ii (∞) is necessarily positive, according to (6.27). From (B 1), it also follows that m(∞) is positive definite and that the off-diagonal elements of m(∞), i.e. m ij (∞), are generally non-zero.
Next, consider an array of OWCs with no oscillating bodies, where an oscillating finite pressure is applied on the internal free surface of OWC i, in the absence of incident waves, while the other OWCs are open to the atmosphere. As the oscillation frequency is increased to infinity, the force on the free surface remains finite since the pressure is finite. There is therefore insufficient force to accelerate the fluid, and hence the kinetic energy of the fluid is zero. Since the potential energy is also zero at infinite frequency, B ii (∞) must be zero according to (6.26). It follows that all B ij (∞) are also zero.
The fact that C(∞) has non-zero elements may be explained by recalling that, in an array of oscillating bodies and OWCs, the radiation coupling coefficient H ij relates the velocity of rigid-body oscillation mode j to the resulting volume flow across the internal free surface of OWC i, when it is open to the atmosphere. Since the fluid is assumed to be incompressible, we cannot avoid creating a volume flow by moving the body, even at infinite frequency.

Two-mode wave-energy converter examples
In agreement with (6.20) we may, for the case of N = 2 oscillation modes, write the collective excitationpower coefficient as where the complex wave excitation variables x i and the complex oscillation amplitudes v i are defined by (6.1). Corresponding to the two r.h.s. terms in (7.1), we may, with reference to (6.4) and (6.5), note that also the adjoint radiation Kochin functionH r (θ) is, for these examples, composed of two terms. Moreover, the complex collective oscillation amplitude U is determined, firstly, by the phase requirement that A * E * (β)/U is a real and positive quantity, and secondly, by a modulus (magnitude) requirement that According to (6.23), the column vector v 0 = v 0 (β), of the optimum complex oscillation amplitudes, has to satisfy the algebraic equation Dv 0 = x(β)/2. This optimum vector v 0 determines the optimum collective parameters E 0 (β) and U 0 . Referring to (3.15), the corresponding maximum absorbed power may be expressed as, for example P a,MAX = P e,OPT /2 = AE 0 (β) = P r,OPT = |U 0 | 2 = v † 0 Dv 0 /2. Considering v = [v 1 v 2 ] T as an independent variable in a two-dimensional complex space C 2 , the relationship Moreover, the radiation damping matrix is real and symmetric, D = R = R T . In the first example, with one symmetric body in heave and surge, the matrix is diagonal, that is R 21 = R 12 = 0. Then there is no hydrodynamical coupling between the two modes. In the second example, R 21 = R 12 = 0, and, moreover, R 22 = R 11 , as we have, for convenience, chosen two equal bodies oscillating in the heave mode, only.

One symmetric body in heave and surge
We shall consider an example with only one immersed body, which has two vertical symmetry planes, one perpendicular to the x-axis and one to the y-axis. This body is assumed to oscillate in just N = 2 modes, surge (i = 1) and heave (i = 2), with complex velocity amplitudes v 1 = u 1 and v 2 = u 2 , and excitation forces x 1 (β) = f e,1 (β)A and x 2 (β) = f e,2 (β)A, respectively. The radiated-wave Kochin coefficient is antisymmetric, h 1 (β + π ) = −h 1 (β), for surge, and symmetric, h 2 (β + π ) = h 2 (β), for heave. Because of the body symmetry, the radiation damping matrix is diagonal, i.e. D = R = diag(R 11 , R 22 ); thus there is no hydrodynamical coupling between the two modes, i.e. R 12 = R 21 = 0. One or both of the diagonal matrix elements R 11 and R 22 for the body may become zero for certain frequencies but are otherwise positive. Let us, however, restrict the following discussion to sufficiently low frequencies to ensure that R 11 and R 22 are never zero, but only positive. Then the radiation damping matrix D = R is non-singular in the frequency interval of interest.
According to (7.1) and (7.2) we now have, for this symmetric-body example, It is interesting to note that the last equation here contains two terms, in contrast to the four terms in (7.2). Thus, for this example, (7.4) appears simply as a two-term extension of (3.10). Consequently, because there is no hydrodynamic coupling between the surge and heave modes, the maximum absorbed power may, in agreement with the alternative equations (3.15), be written simply as where and u i0 (β) is the optimum value of u i . Here, in the last step, we made use of the Haskind relation (5.1) or (6.3). Using (6.5), we note that P a,MAX = AE 0 (β), which is one of the alternative expressions in (3.15). For this symmetric body, where the radiation-resistance matrix is diagonal, that is R = diag(R 11 , R 22 ), the last line in (7.3) vanishes, and thus 2(P a,MAX − P a ) = (u 0 − u) † R(u 0 − u) = R 11 |u 10 − u 1 | 2 + R 22 |u 20 − u 2 | 2 . (7.7) For a fixed value of the absorbed wave power P a , where P a < P a,MAX , this equation represents an 'ellipsoid surface' in the complex two-dimensional u space, C 2 . The centre of the 'ellipsoid' is at the point u = u 0 . The elliptical semi-axes are 2(P a,MAX − P a )/R ii for i = 1, 2. Considering how the absorbed power P a varies with u, the relationship (7.7) may be thought of as a 'paraboloid surface' in the complex twodimensional u space, C 2 . The top point of this 'paraboloid' corresponds to the optimum, (u 0 , P a,MAX ). were previously presented by Evans [29]. A particular case of a symmetric body, for which (7.5)-(7.7) apply, is an axisymmetric body, which was analysed by Newman [1, §10], who found P a,MAX1 = 2J w cos 2 β k and P a,MAX2 = J w k (7.8) for this case.
Another particular case of a surging and heaving symmetric body is the famous Evans Cylinder [2]. It is a two-dimensional WEC device, a horizontal circular cylinder, which is submerged below the free water surface. Let the cylinder axis be aligned in the y-direction, and let the incident wave propagate in the positive x-direction, that is β = 0. For this submerged cylinder, there is no reflected wave, that is, no wave diffraction in the up-wave direction. This means that the Kochin function for diffraction, as introduced by (4.1)-(4.3), vanishes in the up-wave direction. Hence, it follows from the principle of conservation of energy, that the transmitted wave has the same amplitude |A| as the incident wave. Consequently, a non-zero diffracted-wave Kochin function coefficient in the down-wave direction cannot contribute to the amplitude, but only the phase of the transmitted wave.
Another feature of the Evans Cylinder is that the Kochin function coefficients h i , for radiation, as introduced by (6.2), have the property that h 1 (β) = ih 2 (β). Thus if we choose u 2 = iu 1 , which corresponds to a circularly polarized oscillation in the clockwise direction if the x-axis is pointing to the right, then H r (0) = H r (β) = h 1 (β)u 1 + h 2 (β)u 2 = (−ii + 1)h 2 (β)u 2 = 2h 2 (β)u 2 , while H r (π ) = H r (β + π ) = h 1 (β + π )u 1 + h 2 (β + π )u 2 = −h 1 (β)u 1 + h 2 (β)u 2 = (ii + 1)h 2 (β)u 2 = 0. With this circularly polarized oscillation, the radiated waves due to surge and heave are equally large, and they cancel each other in the up-wave direction, but add together constructively in the down-wave direction. Thus, there is neither wave diffraction nor wave radiation in the up-wave direction. All incident wave energy will be absorbed by the Evans Cylinder, provided the circularly polarized oscillation has an optimum amplitude and an optimum phase, in such a way that the down-wave radiated wave exactly cancels the abovementioned transmitted wave. Half of the incident wave energy is absorbed by each of the two modes, surge and heave.

Two equal heaving bodies
Let us consider a system of two equal, semisubmerged, axisymmetric bodies with their vertical symmetry axes located at horizontal positions (x, y) = (∓d/2, 0). We shall assume that they are oscillating in the heave mode only. With this assumption, the excitation-force vector is of the form x(β) = F e (β) = [F e,1 (β) F e,2 (β)] T . Further, the radiation damping matrix may be written as Note that the diagonal entry R d is positive, while the off-diagonal entry R c , which represents hydrodynamical coupling between the two bodies, may be positive or negative, depending on the distance d between the two bodies. As explained in appendix A.2, we may assume that the matrix R is non-singular, and hence R 2 c < R 2 d . We observe that (7.11) has a main algebraic structure similar to that of (7.5) and (7.6), which concern example §7.1, where the resistance-damping matrix is diagonal, R = diag(R 11 , R 22 ). Before we, in appendix A, derived (7.11) and (7.12) we carried out a similarity transformation in order to diagonalize our given radiation damping matrix (7.10); see the similarity-transforming equations (A 5)-(A 11). From a wave-body-interaction point of view, it is interesting to note that the first r.h.s. term in (7.11) and the first equation of (7.12) correspond to a sub-optimum situation when the two, equal, heaving bodies cooperate as a source-mode (monopole) radiator, that is, when the constraint u 2 = u 1 is applied. Then the two bodies are constrained to heave with equal amplitudes and equal phases. By contrast, the last r.h.s. term in (7.11) and the last equation of (7.12) correspond to a sub-optimum situation when the two bodies are constrained to cooperate as a dipole-mode radiator, that is, when the constraint u 2 = −u 1 is applied. In general, (7.11) and (7.12) may be considered to quantify the optimum situation for this combined monopole-dipole wave-absorbing system.

According to (A 10), (A 14), (A 19) and (A 21)-(A 23), the maximum wave power absorbed by the two optimally heaving bodies is
If the maximum radius of each body is sufficiently small, say less than 1 30 of a wavelength, it may be considered as a point absorber, for which the heave excitation force F e is dominated by the Froude-Krylov force, and the diffraction force may be neglected. If, moreover, the centre-to-centre distance d between the two bodies is large in comparison with the maximum body radius, then where F 0 = σρgπ [a(0)] 2 A. Here a(0) is each body's water-plane radius, and A is the complex amplitude of the incident-wave elevation at the chosen origin (x, y) = (0, 0). Further, σ ≤ 1 is a factor that corrects for the diminishing of hydrodynamic pressure with distance below the water surface. (In many cases of practical interest, this correction factor may be approximated to σ ≈ 1.) For this point-absorber case, the entries in the radiation-resistance matrix R in (7.10) are approximately given by [9, eqns 43-44] (see also (A 8)) (7.14) where J 0 denotes the Bessel function of the first kind and zero order. We observe that the matrix R is non-singular, and moreover, are positive, since −1 < J 0 (kd) < 1 for kd > 0. Using formulae (7.11)-(7.15), we find and P a,MAX (β) ≈ |F 0 | 2 4R 0 . (7.17) Note that, in general, this maximum absorbed power is not equally divided between these two bodies [9, eqn 51]. We may note from the point-absorber approximation (7.13) that, since F 0 /A is real, f e,i (β) = F e,i (β)/A = F * e,i (β + π )/A * = f * e,i (β + π ) for i = 1, 2. Correspondingly, we then find from (6.2)-(6.4) thatH r (β) = h 1 (β)u * 1 + h 2 (β)u * 2 = h * 1 (β + π )u * 1 + h * 2 (β + π )u * 2 = H * r (β + π ) and, similarly,H r (β + π ) = H * r (β). Thus, for the two considered heaving bodies, we have here explicitly demonstrated that the term containing the integral on the r.h.s. of (5.4) is, as expected, negligible in the point-absorber limit, because diffraction effects are then negligible.

One single-mode body and one oscillating water column
We consider one single floating body that contains one OWC, and we make the simplifying assumption that only one rigid-body oscillating mode is involved. It could be, for instance, a BBDB device structure [30], in a case where the OWC-containing body is restricted to oscillate in the pitch mode only. Otherwise, we shall also discuss an axisymmetric system where the rigid-body structure is restricted to oscillate in the heave mode only.
With this example, the two N-dimensional column vectors v and x, as well as the N × N radiation damping matrix D, introduced by (6.1), as well as (6.15), reduce to the following two two-dimensional vectors: v = u −p and x = F e −Q e , (7.18) as well as the 2 × 2 matrix (7.19) respectively. In order to determine the maximum absorbed power and the corresponding optimum oscillation, it is convenient to apply similarity transformation as shown in some detail in appendix A. For the present example, the eigenvalues λ 1 and λ 2 of the radiation damping matrix (7.19) are solutions of the seconddegree algebraic equation |D − λI| = λ 2 − (R + G)λ + RG − J 2 = 0. Thus, λ 1 and λ 2 are given by The corresponding two eigenvectors, which satisfy (A 2) and (A 4), are In terms of similarity transformed excitation amplitudes x i (β) and corresponding optimum oscillation amplitudes v i0 (β), both of which are given below, the maximum absorbed power may, according to (A 10) and (A 11), be written as corresponding to the optimum condition We note that the main algebraic structure is similar here and in (7.5)-(7.6) and (7.11)-(7.12). According to (A 4)-(A 6), the similarity transformed complex amplitudes are given by is the similarity transforming matrix-see (A 5). By means of the similarity transformation, (7.3) may be simplified to see also (A 15). We note that (7.24) has an algebraic structure, as well as a geometrical interpretation, similar to that of (7.7). For the particular case of a heaving axisymmetric body that contains an axisymmetric OWC, we have J 2 = RG, and thus, from (7.20), we see that λ 1 = R + G and λ 2 = 0, which means that matrix D, in this case, is singular and of rank r D = 1 [10, eqn 73]. In this case, there is only one term in the sum on the r.h.s. of (7.24), which simplifies to which represents a circle in the complex v 1 plane. The centre of the circle is at v 1 = v 10 (β) = x 1 (β)/(2λ 1 ) = x 1 (β)/(2R + 2G), and the radius is 2(P a,MAX − P)/(R + G). While figure 3 may serve to illustrate (7.24), figure 1 is more relevant as an illustration of (7.25). Because of the singularity of the radiation damping matrix, the similarity-transformed variable v 2 is irrelevant, and may have any arbitrary value, without influencing the absorbed power. The physical reason for the singularity is that both modes, the heavingbody mode and the OWC mode can radiate only isotropic outgoing waves. To realize maximum absorbed wave power, the optimum isotropically radiated wave may be realized by any optimum combined wave radiation from the axisymmetric OWC and the heaving axisymmetric body. The transformed oscillation v 2 corresponds to a situation where the heave mode and the OWC mode cancel each other's radiated waves in the far-field region.

Discussion
In this section, we first compare two versions of the so-called 'fundamental theorem of wave power' [17]. We shall discuss, secondly, the direction-averaged maximum absorbed wave power for an array of WEC units, and also, thirdly, the physical interpretation of the absorbed-wave-power surfaces. Finally, we shall comment on a disputed formula applied to the optimum performance of the Evans Cylinder.

The 'fundamental theorem of wave power'
In this paper, by considering the physical process of wave-power absorption at the wetted surface of an oscillating immersed body, and, more generally, at a WEC array's wave-interacting surfaces, we derived, in § §3 and 6.2, respectively, an LPV version of the 'fundamental theorem of wave power', equation (3.11) : P a = AE(β) + A * E * (β) − |U| 2 . Moreover, we presented, in §4, a GPV version, equation (4.5): P a = I(β) + I * (β) − |G| 2 , where I(β) + I * (β) is the wave-power input through an envisaged surface enclosing all WEC units, and |G| 2 is the outgoing wave power through the same envisaged surface, which, for mathematical convenience, is chosen in the far-field region of the generated waves. In § §5 and 6.1 we introduced a mixed, or hybrid, GLPV version, where the LPV parameters |U| 2 and E(β), by means of (4.10), (5.2) and (6.5), are expressed in terms of global far-field quantities. For these versions of the 'fundamental theorem of wave power', the r.h.s. has three terms, the sum of two complex conjugate terms minus a real, non-negative, term. The third term of the LPV versionincluding the GLPV version-contrary to the GPV version, is quadratically dependent on the oscillation amplitudes, but independent of the incident wave amplitude, while the first and second terms are linear in both kinds of amplitudes. For the LPV version, the first two terms, the excitation power, P e = AE(β) + A * E * (β), represent the gross power input from the incident wave, while the third term, P r = |U| 2 , is the necessary, outward-propagating, radiated power.
With the LPV/GLPV and GPV versions, the third terms |U| 2 and |G| 2 , which represent energy associated with the radiated waves and the outgoing waves, as given by (4.10) and (4.6), respectively, should be considered as a necessity rather than a power loss. In order to absorb wave energy, it is necessary, firstly, to have wave-diffracting WEC units immersed in the sea, and, secondly, to let the WEC units oscillate and thus produce radiated waves, which interfere destructively with the incident wave. The WEC units need to oscillate in order to receive wave energy.
Before comparing the LPV/GLPV and GPV versions applied to a point absorber, let us consider a two-dimensional 100% absorbing WEC unit, such as an optimally run Evans Cylinder [2] or a hinged oscillating flap in the down-wave end of a wave channel, a flap that we may consider as an ideal nodding-duck device [23]. For these examples, the LPV equation (3.15) shows that the optimum values of the excitation power P e and the radiated power P r correspond to 200% and 100%, respectively. For the GPV version, which does not discriminate between radiated waves and diffracted waves, the optimum values of the input power, P i = I(β) + I * (β) and the outgoing power P g = |G| 2 correspond to 100% and 0%, respectively. (Note that for a real nodding-duck WEC that absorbs less than 100%, the optimum outgoing power is not zero, and the optimum input power is larger than the maximum absorbed power.) As we shall see below, the two versions show a less drastic difference when applied to a point absorber.
In a case where the WEC array is not oscillating, there is no absorbed wave power, P a | v=0 = 0, and also no radiated power, P r | v=0 = 0. Then it follows from (3.11) that P e | v=0 = 0, and, moreover, from (4.5) that P i | v=0 = P g | v=0 ≡ P d , where P d is the outgoing power associated with the diffracted wave alone. From (4.8), we may note that P i | v=0 = P d = ρv p v g Re{H d (β)A * } ≥ 0.
In cases of rather weak diffraction, as with a wave-power-absorbing very small point absorber, P d may be negligibly small. We may note that, if H d (θ) is small for all θ (including θ = β), then the r.h.s. of (4.9) is small of second order. Thus, in cases of very weak diffraction, H d (β) is, approximately, purely imaginary, if we choose A to be real. This matter has been discussed in more detail by Farley [15]. By comparing (4.7) with (4.11) and (4.6) with (4.10), we observe that, for cases where the diffracted wave is negligible compared to the radiated wave, I(β) ≈ E(β)A and |G| 2 ≈ |U| 2 . Thus, for such cases, there is no great difference between corresponding terms of the LPV/GLPV and GPV versions of the 'fundamental theorem of wave power'.
In §3, oscillations, wave forces, power and energy were quantitatively discussed in the time domain, but elsewhere, in this paper, only in the frequency domain. In the case of non-sinusoidal waves, it may be desirable to carry out analyses in the time domain. In this situation, a time-domain type of the 'fundamental theorem of wave power' may be desirable. This type should correspond to an inverse Fourier transform of the LPV version derived in §6.2-or §3 for the one-mode case. It should neither be the GPV version nor the GLPV version, which are derived and discussed in § §4, 5 and 6.1. These versions cannot represent the instantaneous power absorbed by the WEC, but only the time-average power. With a time-domain analysis, also the reactive-power terms-see §3 and §6.2-need to be taken into account.

Direction-averaged maximum absorbed wave power
For the case of only one immersed WEC unit oscillating in a single mode i = 1, we have, in agreement with (4.3) and (5.3), that H r0 (θ ) = h 1 (θ )v 10 (β), and thatH r0 (β + π ) = h 1 (β + π )v * 10 (β). Then the optimum Since the Kochin function coefficient h 1 (θ ) is a function of geometry and mode of motion, this means that the optimum gain function G g0 for this case depends on geometry and mode of motion only, and not on the WEC velocity. However, to maximize power absorption, the WEC unit needs to move with optimum amplitude and phase. For an isotropically radiating system, such as a heaving axisymmetric body, the optimum gain function is G g0 = 1 and independent of the wave-incident angle β. Then, (6.10) gives the maximum absorption width d a,MAX ≡ P a,MAX /J w = G g0 /k = 1/k, a well-known result since the mid-1970s. From (8.1), we find that the direction-averaged optimum gain function is as averaged over all directions of wave incidence, a result reported by Fitzgerald & Thomas [31]. However, in some singular cases, we may find that G g0,average = 0. For instance, any axisymmetric body oscillating only in the yaw mode can, in an ideal fluid, neither radiate nor absorb wave energy, for any frequency. Then H r (θ ) ≡ 0. For a floating semi-submersible platform, as well as for a floating bottle-shaped axisymmetric body that has a relatively small water-plane area, the heave excitation force vanishes at a certain frequency [32, p. 77]. Hence, according to the Haskind relation (5.1), h 1 (β + π ) = 0, and then G g0 (β) = 0 at this particular frequency. For a general WEC array oscillating in N modes, with N ≥ 2, it is not convenient to apply (6.11) to determine G g0,average . By means of another mathematical procedure, involving similarity transformation, as applied in appendix A, it is found that G g0,average , in general, equals the rank r D of the radiation damping matrix D-see (A 17). In cases where this matrix is non-singular, the rank of the matrix equals its dimensionality N. Thus, in general, the direction-averaged value of the optimum gain function is equal to an integer in the interval 0 ≤ G g0,average ≤ N. For instance, an immersed body may oscillate in N = 6 different modes. However, if the body has a vertical axis of symmetry, then G g0,average = 3 in general, or less in exceptional cases [1, §10]. These results extend the findings of Wolgamot et al. [18], who found the result G g0 (β) = N for cases where the, general, hermitian radiation-damping matrix D specializes to a, non-singular, real radiation-resistance matrix.
We considered in §7.2 an array consisting of two heaving point absorbers, and we derived formula (7.17) for the maximum absorbed power. To find the direction-averaged maximum absorbed power, we need to integrate cos{kd cos β} from β = 0 to β = 2π . Since this integral equals J 0 (kd) (e.g. [33, formula 9.1.18, p. 360]) we see, from (7.17) and the first equation of (7.14), that the direction-average of the absorbed power is |F 0 | 2 /(4R 0 ) = 2J w /k. This result was, according to (A 1) or (A 17), to be expected for this non-singular system of two heaving bodies.

Absorbed-wave-power surfaces
The absorbed wave power relative to its maximum may be expressed as two equally general functions of the WEC oscillation amplitudes relative to the optimum amplitudes. The first of these expressions is given in (3.17), which may be illustrated as an axisymmetric paraboloid on the complex collectiveamplitude U plane (figure 1). The second expression is given in (6.24), which may be thought of as a 'paraboloid' in the complex N-dimensional v space, C N (see figure 3 for an example with N = 2).
In popular terms, it might be useful to think of this 'paraboloid', for P a > 0, as a single 'mountain island' in a C N 'world', where there is otherwise, for P a < 0, only an infinite 'ocean' (cf. figures 1 and 3). The top of the absorbed-power 'mountain' corresponds to optimum, and the 'shore' of the 'island' to P a = 0. This 'mountain' top can be reached only if no technical or practical constraint prevents the required complex amplitudes v i0 from being realized for all i = 1, 2, 3, . . . , N.
For practical reasons, it may not be possible to realize the optimum condition v = v 0 . Note that all components of the excitation vector x(β), and hence also of the optimum oscillation amplitude vector v 0 (β), are proportional to the complex amplitude A of the wave elevation of the undisturbed incident wave, and that P a,MAX is proportional to A * A = |A| 2 . As oscillation amplitudes cannot exceed their design-specified limits, it will not be possible to realize the described optimum situation, if the amplitude of the incident wave exceeds a certain critical value. With such constraints, or for other technical reasons preventing realization of the optimum condition (6.23), the practical, constrained-case, maximum absorbed power P a,max will be less than the ideal P a,MAX . (In such a case, it will not be practically possible to 'climb' to the 'top' of the above envisaged 'mountain island'.)

A disputed 1979 formula
For the optimum case, the maximum absorption width may be expressed as according to (6.9)-(6.11). The single-mode version of this formula was derived by Newman [25] and, independently, by Evans [8]. For the case of more than one oscillation mode, we have found it necessary to introduce the adjoint Kochin functionH r in the numerator (see (5.3) and (6.4)).
It should be emphasized that formula (8.3) applies only to the optimum case, as it is based on the fact that the optimum radiated power is equal to the maximum absorbed wave power. Neither any wave force nor the incident wave amplitude is explicitly present in (8.3). However, as each component v i0 of the optimum oscillation vector is proportional to the incident wave amplitude, the numerator, as well as the denominator, of fraction (8.3) is proportional to the square of the incident wave amplitude.
The controversy [15,16]-concerning the numerator in formula (8.3)-is mainly related to the Evans Cylinder, which we, in the last three paragraphs of §7.1, discussed in some detail. Let this cylinder be aligned perpendicular to the incident wave direction. At optimum oscillation, this submerged horizontal cylinder absorbs all incident wave energy. Then the optimum radiated wave has to propagate only down-wave in order to cancel the transmitted wave, as there is no reflected wave to cancel up-wave. In agreement with this physical observation, it is reassuring to observe that we, according to (7.9), which is valid for the Evans Cylinder, may replace |H r0 (β + π )| by |H r0 (β)| in the numerator of (8.3). Such a replacement seems to resolve the controversy, because H r0 (β) = h 1 (β)u 10 + h 2 (β)u 20 is the Kochin function for the down-wave radiated wave when optimum wave-power absorption is actually taking place, that is, when the rotating Evans Cylinder's surge and heave modes' complex velocity amplitudes have their optimum values, u 10 and u 20 , respectively. In comparison,H r0 (β + π ) = h 1 (β + π )u * 10 + h 2 (β + π )u * 20 is the adjoint Kochin function, which corresponds to a wave that is radiated in the opposite direction if the rotating Evans Cylinder is oscillating with opposite sense of rotation. Observe that to replace u i0 by u * i0 corresponds to time reversal, since (e iωt ) * = e iω(−t) . It seems that we have to include the relation (8.3) among the reciprocity relations between physical quantities about which Newman [1, §7] expressed that they 'are not physically related to each other in any obvious manner'.

Conclusion
After the petroleum crisis in 1973, the basic theory for primary wave-energy conversion was developed during the mid-and late 1970s and the early 1980s. Different versions-the LPV version (cf. § §3 and 6.2), the GLPV version (cf. § §5 and 6.1) and the GPV version (cf. §4)-of the so-called 'fundamental theorem of wave power' have given rise to some controversy even during recent years. Comparative discussion of these different versions has been presented in §8.1. The GLPV version, in particular, is mathematically convenient when proving some useful reciprocity relations, as applied, for instance, by Newman [1]. It is, however, difficult to give a physical interpretation of some of these relations and of the GLPV version. This may be the cause of recent controversy concerning the GLPV version. All of these versions provide, however, the correct value of the time-average absorbed wave power. Hopefully, the discussion in §8.4 helps to do away with some of this controversy.
The LPV, the GPV and the GLPV versions express, respectively, the absorbed wave power P a in terms of the WEC units' complex oscillation amplitudes, in terms of the outgoing (diffracted plus radiated) wave, and in terms of the radiated wave alone. For mathematical convenience, the outgoing and radiated waves in the far-field region are considered, explicitly.
For a general WEC array consisting of oscillating immersed bodies and OWCs, we have found it convenient to introduce complex collective parameters, the collective oscillation amplitude U and the collective excitation-power coefficient E(β) (see (3.10) and (6.20)). Then it is, even for a WEC array, a rather simple algebraic exercise to derive expressions for the maximum absorbed wave power P a,MAX and the corresponding optimum values U 0 and E 0 (β) (see (3.11)-(3.15)). Moreover, we may illustrate the dependence of the absorbed wave power P a versus U as an axisymmetric For an N-mode WEC array, we may, in greater detail than figure 1, consider the real-valued P a as represented by a paraboloid in an N-dimensional complex v space C N . Cross sections of such a paraboloid are, as an example, illustrated in figure 3. Mathematically, the mentioned paraboloids are represented by rather simple mathematical expressions (3.17) and (6.24), which may be considered as alternative variants of the LPV version of the 'fundamental theorem of wave power'.
In contrast to the GPV version and the GLPV version, only the LPV version is applicable for the purpose of quantifying the instantaneous absorbed wave power. Then it is necessary to take also the reactive power into account. When deriving the LPV version (cf. § §3 and 6.2), we also discussed the reactive power that is associated with wave-power absorption. In appendix B, we have derived expressions that relate the 'added-mass' matrix, as well as a couple of additional reactive radiationparameter matrices, to the difference between kinetic energy and potential energy in the water surrounding the immersed oscillating WEC array. To the best of the authors' knowledge, some of these derived relations are new results (e.g. (B 25), (B 27) and (B 28)). In appendix B.5, we have also derived new relations concerning reactive power associated with the incident wave.
In appendix A, we applied similarity transformation of the radiation damping matrix to derive a formula for the direction-average maximum absorbed wave power P a,MAX,average . Correspondingly, as discussed in §8.2, we found that the direction-average value G g0,average of the optimum gain function G g0 (β)-defined by (6.11)-equals an integer in the interval 0 ≤ G g0,average ≤ N, where N is the WEC array's number of modes of oscillation (number of degrees of motion). Only when the radiation damping matrix is non-singular, we have G g0,average = N, as derived by Wolgamot et al. [18, eqn 21] for an N-mode WEC array consisting of oscillating bodies only. Thus, our result is an extension of theirs, to WEC arrays that may contain OWCs and also may have a singular radiation damping matrix. In general, G g0,average equals the rank of this N × N matrix.
Acknowledgements. This work has been initiated following extensive discussions with Prof.  Competing interests. We declare we have no competing interests.

Appendix A. Direction-average maximum absorbable wave power
Consider wave-energy absorption by an array of immersed oscillating rigid bodies and OWCs. The number of wave-interacting oscillators is N = N u + N p , where N p is the number of OWCs and N u is the number of used body modes, which may be up to six times the number of bodies.
Observing, e.g. from (6.10) and (6.22), that the maximum absorbed power P a,MAX = P a,MAX (β), our present aim is to prove the following conjecture [34, p. 17]: where r D is the rank of the radiation damping matrix D, as given in (6.15). For cases when matrix D is non-singular, its rank equals its dimensionality, r D = N. For an array without OWCs and a non-singular radiation resistance N u × N u matrix R, which is real, symmetric and positive definite and has a rank of r R = N u , Wolgamot et al. [18, eqn 21] derived a result in agreement with our conjecture (A 1). However, in many cases of interest, the radiation matrix is singular, that is, its determinant vanishes, |D| = 0. In such cases, it is possible to find a linear combination of two different oscillation modes, in such a way that the resulting far-field radiated wave vanishes. This is possible, for instance, with certain linear combinations of surge and pitch modes for an axisymmetric body [35, p. 175]. A two-dimensional example, a symmetric body, which has N u = 3 modes of motion (surge, heave and pitch), was discussed by Newman [1], who found that two modes (heave and surge or heave and pitch) are sufficient to absorb 100% of the incident wave energy; thus a third mode is not needed. Newman also studied optimum oscillation for maximum wave power absorbed by one body which has a vertical symmetry axis and oscillates in the three modes surge, heave and pitch-an example where either surge or pitch was found not necessary. Including the three remaining modes (sway, roll and yaw) for this axisymmetric body, we find that its radiation resistance matrix R has a dimensionality N u = 6, while its rank is only r R = 3 (cf. [

A.1. Similarity transformation
To derive a result in agreement with our conjecture (A 1), Wolgamot et al. [18] applied Cholesky decomposition of the real symmetric radiation resistance matrix R, provided this matrix is positive definite and thus non-singular. It would be straightforward to apply the same method for our complex hermitian radiation damping matrix D in cases where this matrix is non-singular and thus positive definite. However, we wish to include the general case when matrix D may be singular and thus need to have at least one of its eigenvalues equal to zero. In this case, Cholesky decomposition is less convenient.
Instead we shall apply similarity transformation (e.g. [36, ch. V]) to prove our conjecture (A 1). Then our first step is to determine the eigenvalues λ i and the corresponding, normalized and mutually orthogonal, eigenvectors e i for i = 1, 2, 3, . . . , N − 1, N. They have to satisfy the following system of linear homogeneous equations: Because matrix D is hermitian, all its eigenvalues are real, and since it is positive semidefinite, they are all non-negative. They are the N solutions of the N-degree equation |D − λI| = 0, where I = diag(1, 1, 1, . . . , 1) is the N × N identity matrix, for which all diagonal entries equal 1, while all other entries equal 0. Further, |D − λI| denotes the determinant of matrix (D − λI). We may wish to arrange the N eigenvalues in descending order, that is, Note that if e i is a possible solution when λ = λ i , then also C i e i is a solution, where C i is an arbitrary complex scalar. It is convenient to normalize the eigenvectors e i in such a way that the condition and conversely, v = Sv and x = Sx . Then, the collective excitation-power coefficient E(β) and the collective oscillation amplitude U defined in (6.20) may be expressed as We have here made use of (A 2), (A 4) and (A 5) to obtain the last equality. Note that matrix D has the same eigenvectors and eigenvalues, and hence also the same rank, as matrix D. While D is complex and hermitian in the general case, the diagonal matrix D is real. Moreover, the complex collective parameters E(β) and U remain unchanged through the similarity transformation. By means of a well-known reciprocity relation [10], we may express the radiation damping matrix D in terms of the wave excitation vector x as follows: where k is the angular repetency (wave number) and J w is the wave power level (incident wave power transport per unit width of the wave front), defined in (2.2). Using (A 6)-(A 8), we obtain Moreover, the optimum condition (6.23) may be reformulated as ) and the maximum absorbed power, as given by (6.22), may be rewritten as In the general case, the N × N matrix D may be singular, and then of rank r D < N. For convenience, let its N eigenvalues be arranged to satisfy the following condition: where the integer r D satisfies 1 ≤ r D ≤ N. The matrix is non-singular if r D = N. From (A 7) and (A 9), we may now conclude that the only non-zero matrix-element entries of our transformed radiation damping matrix D are the following r D diagonal entries: Note from (A 10) that at optimum, if λ i = 0, the oscillation component v i0 (β) has to be in phase with the excitation component x i (β). In general, because matrix D may have non-zero off-diagonal entries, unless N = 1, such a simple optimum-phase relationship does not apply between v j0 (β) and x j (β), components of vectors v 0 (β) and x(β), where j = 1, 2, . . . , N. The smallest positive eigenvalue is λ r D . Observe that if it is very small, then there is a risk that the required optimum amplitude |v 0r D (β)| may be impractically large! Note that for r D < i ≤ N, we have no particular optimum requirement for v i0 (β). Then the optimum value of v 0 (β) is ambiguous, although the maximum absorbed power is unambiguous, as is made clear by the following equation (A 14).
Using (A 10) and (A 12), the maximum absorbed power (A 11) can now be written as This result represents a generalization of a previous expression [37, eqn 135] for the maximum wave power absorbed by an array consisting of oscillating rigid bodies to a more general situation when the array contains also OWCs. Also, by using (A 6), (A 7) and (A 12), equation (6.24) may be expressed as which, for any fixed positive l.h.s. value, represents an 'ellipsoid' in the complex r D -dimensional v space, C r D .
In agreement with (A 14), the maximum absorbed wave power averaged over all angles of wave incidence may be written as Using (A 13) to eliminate the integral, we finally obtain in agreement with our conjecture (A 1), which we have now proved to be correct!

A.2. Two-body example
As a simple example of similarity transformation, consider a system of two equal axisymmetric bodies with their vertical symmetry axes located at (x, y) = (∓d/2, 0). We shall assume that they are oscillating in the heave mode only. With this assumption, the excitation-force vector is of the form x = F e = [F e,1 F e,2 ] T . Further, the radiation resistance matrix may be written as the 2 × 2 matrix Note that the diagonal entry R d is positive, while the off-diagonal entry R c may be positive or negative, depending on the distance d between the two bodies. Since waves radiated from the two distinct bodies cannot cancel each other in all directions in the far-field region-which may be proved, by summing two asymptotic expressions (e.g. [35, formulae 4.266 and 4.270])-the matrix R is non-singular, and hence R 2 c < R 2 d . (We assume that the body does not have such a peculiar shape that we may find R d to vanish for some particular frequency.) The eigenvalues λ = λ 1,2 satisfy the algebraic equation It is easy to show that, for this case, the corresponding eigenvectors e 1  Appendix B. Difference between kinetic energy and potential energy in the near-field region It is well known that, for a propagating plane wave, the time-average energy stored per unit horizontal area E = E k + E p is divided into equal amounts of kinetic energy and potential energy. Thus, E k = E p = E/2, that is, E k − E p = 0. As will be evident through a derivation below, the same is true for a radiated, outwards propagating, circular wave in the far-field region. However, generally, in the near-field region of a group of oscillating immersed bodies, E k − E p = 0. For a case with no incident wave, this difference is related to the added-mass matrix m for immersed radiating bodies performing forced oscillations, through the relationship [19,38] T − V = 1 4 u † mu = 1 4ω u † Xu.

( B 1 )
Here u = [u 1 u 2 u 3 · · · u N u ] T is the vector of complex velocity amplitudes of the N u used body modes. The l.h.s. of (B 1) denotes the total near-field value of the difference between kinetic and potential energy, namely, where S 0 is the air-water interface. We have, in (B 1), also introduced the radiation-reactance matrix X = ωm. (This is a less confusing terminology than 'added mass' for cases where V > T.) Our aim is to generalize (B 1) to a situation where the array contains N p OWC oscillating modes in addition to the rigid bodies, as indicated in figure 2. Let the column vector p = [p 1 p 2 p 3 · · · p N p ] T denote the oscillating state of the N p OWCs, where p k is the complex amplitude of the dynamic air pressure of OWC k. For such a more general case, the integral on the r.h.s. of (B 2) needs to be replaced by integrals derived below. We shall find, as a result, that a generalized version of (B 1) has to contain additional terms on the r.h.s., including terms associated with the incident wave.

B.1. Preparatory basic equations
In accordance with linear theory for an ideal fluid, the complex amplitudes p = −iωρφ of the hydrodynamic pressure and v = ∇φ of the fluid (water) velocity are derivable from the velocity potential φ, which, everywhere in the ideal fluid, has to obey the Laplace equation ∇ 2 φ = 0. Moreover, φ has to satisfy certain boundary conditions as shown below. The complex amplitude of the water-surface elevation η on S 0 , as well as η k on S p,k = S k (figure 2) is given by respectively. Here the water-air interfaces, outside the OWC structures and inside the kth OWC structure, are denoted by S 0 and S p,k , respectively. (It is assumed that the ambient air pressure is static.) Moreover, ρ is the mass density of water, g is the acceleration of gravity and ω is the angular frequency of the considered waves and oscillations. The potential energy, associated with the wave, results from water being lifted against gravity, and inside the kth OWC also against the dynamic air pressure p k . On S 0 and S p,k , the time-average of the potential energy, per unit horizontal surface, is respectively.
For convenience, we decompose the velocity potential φ into three terms, φ 0 , φ d and φ r , corresponding to the incident, the diffracted and the radiated waves: where φ s = φ 0 + φ d is the scattered potential. Further, we have introduced column vectors ϕ u = [ϕ u,1 ϕ u,2 · · · ϕ u,N u ] T and ϕ p = [ϕ p,1 ϕ p,2 · · · ϕ p,N p ] T , which are composed of space-dependent, complex, proportionality coefficients ϕ u,i and ϕ p,k , that quantify how the radiated potential φ r is linearly related to u i and p k , components of the oscillation-state column vectors u and p, respectively. It is required that every term in (B 5) satisfies the Laplace equation. Furthermore, it is required that every term in (B 5), except φ 0 -and, consequently, also except φ s and φ-satisfies the radiation condition of outgoing waves at infinity, that is where r is the radius of the envisaged control cylinder S ∞ (figure 2). It may be remarked that, for large values of kr, all these complex amplitudes, φ d , φ r , ϕ u and ϕ p , diminish asymptotically as 1/ √ kr when kr → ∞, and their directional variation may, quantitatively, be described by certain complex factors s d (θ), s r (θ ), s u (θ ) and s p (θ ), respectively.
Furthermore, it is required that the components of the velocity potential in (B 5) satisfy certain boundary conditions. The following homogeneous boundary conditions (B 7) apply on the ambient-airto-water interface S 0 and on the totality of immersed stationary rigid-body surfaces S b , including the sea bed (figure 2): Next, for completeness, we shall state the boundary conditions on the wave-radiating wetted surfaces S i = S u,i of the immersed oscillating bodies and on the wave-radiating internal air-water interfaces S k = S p,k of the OWCs (figure 2). It is convenient to define the two unions S p = N p k=1 S p,k and S u = N u i=1 S u,i (B 8) as the totality of the two types of wave-radiating surfaces. On all OWC wave-radiating internal air-water interfaces S p,k , and on all oscillating rigid-body wetted surfaces S u,i , the following boundary conditions apply: to the difference between kinetic and potential energy stored in the water surrounding the WEC array. This contribution is adding to the energy difference which is quantified by (B 25), and which is applicable to a situation when the WEC array is performing forced oscillation in the absence of an incident wave.