2. Problem specification and solution

Two-dimensional Cartesian coordinates, $(x,z)$ , are used, $z=0$ corresponding to the rest position of the fluid surface and $z$ directed vertically downwards. In $x \lt 0$ the fluid depth is taken to be a constant, $h$ . At $x=0$ a vertical wall extends upwards from $z=h$ to $z=b$ and from $z=a$ through the surface. The gap in the wall $a \lt z \lt b$ connects the fluid in $x \lt 0$ to an internal fluid surface in $a \lt x \lt b$ via a curved annular duct in $x \gt 0$ . Its inner and outer walls are described in polars coordinates, $(r,\theta )$ with $x = r \sin \theta$ , $z = r \cos \theta$ , by $r=a$ and $r=b$ for $0 \leq \theta \leq \pi /2$ . The duct is divided into $N$ curved channels of uniform width by $N-1$ thin internal curved baffles along $r = r_{\kern-1pt j}$ , $1 \leq j \leq N-1$ $0 \leq \theta \leq \pi /2$ where $r_{\kern-1pt j} = a + j \Delta$ for $0 \leq j \leq N$ . Thus $\varDelta = (b-a)/N$ is the uniform width of each fluid channel and we suppose that $N$ is large enough that $\Delta /a$ is suitably small. This quantity defines the largest aspect ratio of width to length within the array of $N$ fluid channels and its smallness is needed for the method of solution we propose.

Figure 1.

Definition sketch of the OWC device with $N$ curved fluid filled channels with internal free surfaces above which air is pumped from an enclosed chamber through air turbines, indicated by bow ties.

It is supposed that above the internal free surface of each of the $N$ partitioned channels lies an air chamber, which is isolated from the other $N$ channels and connected to the external atmosphere by a duct containing its own self-rectifying Wells air turbine. Thus, when the fluid is set into motion, the rise and fall of the internal fluid surface within each channel drives the air back and forth through the turbine.

The water is assumed to be incompressible, inviscid (although viscous effects are considered later in the paper) and has density $\rho$ ; its motion is assumed to be of small amplitude and irrotational. This allows us to describe the fluid velocity as the gradient of a velocity potential $\varPhi (x,z,t)$ which we express as

(2.1) \begin{align} \varPhi (x,z,t) = {\rm Re} \big\{{-}({\mathrm{i}} g A/ \omega ) \phi (x,z) \mathrm{e}^{-{{i}} \omega t} \big\} , \end{align}

where $A \in \mathbb{C}$ is the incident wave amplitude (encoding a modulus and phase) and $g$ is acceleration due to gravity. We have anticipated the linearity of the governing equations that follow implying that a monochromatic incident wave of angular frequency $\omega$ drives a response of the same frequency thereby allowing us to remove a time-harmonic factor from the motion; $\phi (x,z)$ is consequently a complex-valued velocity potential satisfying

(2.2) \begin{align} {\nabla} ^2 \phi = 0 \end{align}

everywhere in the fluid. The surface elevation (measured upwards from $z=0$ ) is represented by $\zeta (x,t) = {\rm Re} \{ A \eta (x) \mathrm{e}^{-{{i}} \omega t} \}$ , $\eta (x)$ is the complex spatial amplitude of the free-surface elevation, following the standard decomposition of the motion into a time-harmonic factor and a spatial part, such that the linearised kinematic condition can be written as

(2.3) \begin{align} K \eta (x) = -\phi _z(x,0), \qquad \big(K = \omega ^2/g \big) , \end{align}

on all free surfaces: both in $x \lt 0$ , and in each of the channels across $a \lt x \lt b$ . The pressure is defined as $P(x,z,t) = P_a + {\rm Re} \{ A \rho g p(x,z) \mathrm{e}^{-{{i}} \omega t} \}$ where $P_a$ is the atmospheric pressure. Use of the linearised Bernoulli equation

(2.4) \begin{align} P(x,z,t) = P_a - \rho \varPhi _t + \rho g z \end{align}

on the free surface $z = -\zeta (x,t)$ , gives

(2.5) \begin{align} \phi (x,0) - \eta (x) = 0, \qquad x \lt 0 \end{align}

and on each internal free surface $r_{\kern-1pt j-1} \lt x \lt r_{\kern-1pt j}$ ( $j=1,\ldots ,N$ ) we have

(2.6) \begin{align} p_{\kern-1pt j} = \phi (x,0) - \eta (x) , \end{align}

where $p = p_{\kern-1pt j}$ represents the oscillating air pressure in the chamber above the $j$ th channel. Combining (2.3) with (2.4) and (2.5) results in

(2.7) \begin{align} \phi _z(x,0) + K \phi (x,0) = \left \{ \begin{array}{l} 0, \qquad x \lt 0 ,\\[5pt] K p_{\kern-1pt j}, \qquad r_{\kern-1pt j-1} \lt x \lt r_{\kern-1pt j}. \end{array} \right . \end{align}

The oscillating pressure drives a flow through the air turbine connected to the $j$ th channel which is assumed to obey a standard semi-empirical law (see Sarmento & Falcão Reference Sarmento and Falcão1985; Raghunathan Reference Raghunathan1995; Martin-Rivas & Mei Reference Martin-Rivas and Mei2009) relating pressure to the volume flux (per unit length) of fluid, $Q(t)$ , across the mean free surface via

(2.8) \begin{align} Q(t) = \lambda _1 (P- P_a) + \lambda _2 \frac {{\mathrm{d}} P}{{\mathrm{d}} t}. \end{align}

When applied in channel $j$ (2.8) reduces to

(2.9) \begin{align} q_j = {\mathrm{i}} \varLambda _j p_{\kern-1pt j} , \end{align}

where $\varLambda _j$ is a complex turbine coefficient, defined in (2.11), that relates the oscillating chamber pressure to the volume flux through the turbine in channel $j$ .

(2.10) \begin{align} q_j = -\int _{r_{\kern-1pt j-1}}^{r_{\kern-1pt j}} \phi _z(x,0) \, {\mathrm{d}} x \end{align}

represents the volume flux (per unit length) of fluid crossing $z=0$ . In (2.9)

(2.11) \begin{align} \varLambda _j = \rho \omega \lambda _{1,j} - {\mathrm{i}} \rho \omega ^2 \lambda _{2,j}. \end{align}

The law (2.8), applied to channel $j$ , accounts for the compressibility of air (Sarmento & Falcão Reference Sarmento and Falcão1985; Martin-Rivas & Mei Reference Martin-Rivas and Mei2009) with

(2.12) \begin{align} \lambda _{2,j} = \frac {V_{\kern-1pt j}}{\rho _a c_a^2} , \end{align}

where $\rho _a$ is air density, $V_{\kern-1pt j}$ is the volume of air (per unit length) in the $j$ th chamber (with respect to the rest position of the fluid) and $c_a$ is the speed of sound in air. The parameter $\lambda _{1,j}$ can be written in terms of the control and design parameters of a Wells turbine, but we imagine this is a free power control parameter.

No flow conditions apply on all fixed submerged boundaries in the flow and this is expressed as

(2.13) \begin{align} \partial \phi /\partial n = 0, \end{align}

representing the vanishing of the normal derivative of $\phi$ .

Finally, we need to impose a condition in the far field and write

(2.14) \begin{align} \phi (x,z) = \big(\mathrm{e}^{{{i}} k x} + R \mathrm{e}^{-{{i}} k x} \big) \psi _0(z), \qquad \mbox{as $x \to - \infty $}. \end{align}

where $R \in \mathbb{C}$ is the reflection coefficient and

(2.15) \begin{align} \psi _0(z) = \cosh k(h-z)/\cosh kh \end{align}

is the depth dependence associated with propagating waves of wavenumber $k$ , found from separating variables: $k$ is defined to be the positive real root of

(2.16) \begin{align} K = k \tanh kh. \end{align}

The standard method of solution for solving the problem of $N$ independent oscillating pressure distributions is to decompose the problem into $N+1$ potentials: the scattering of incident waves in the absence of pressure distributions plus $N$ radiation potentials defined by a prescribed isolated forcing in each of the $N$ channels in the absence of incoming waves. The individual solutions are then coupled via the turbine law connecting volume fluxes and pressures in each of the $N$ channels. The general theory can be found in Evans (Reference Evans1982). However, each problem is complicated and the complete solution is expressed in terms of $N \times N$ matrices. The advantage of this approach it that it does furnish an explicit condition maximising efficiency, but only by supposing that there is complex-valued control coupling each turbine to every other one, which is not practical.

We take a different approach, the first step of which combines the turbine law (2.9), (2.10) with the surface condition (2.7) to eliminate $p_{\kern-1pt j}$ resulting in

(2.17) \begin{align} \phi _z(x,0) + K \phi (x,0) = {\mathrm{i}} \frac {K}{\varLambda _j} \int _{r_{\kern-1pt j-1}}^{r_{\kern-1pt j}} \phi _z(x,0) \, {\mathrm{d}} x, \qquad r_{\kern-1pt j-1} \lt x \lt r_{\kern-1pt j} \end{align}

( $j=1,\ldots ,N$ ).

2.1. Homogenisation in $x \gt 0$

Next, we aim to simplify the complexity of the duct containing the baffles by approximating the problem using homogenisation methods (see, for example, Mei & Vernescu Reference Mei and Vernescu2010). The basis of this approximation is $\epsilon = \Delta /r_0 \ll 1$ (where $r_0 = a$ ) implying that the width of the channels is small compared with their length. First, we non-dimensionalise the governing Laplace equation (2.2) in the duct by writing $r = r_0 \hat {r}$ , $\phi (r,\theta ) = \hat {\phi }(\hat {r},\theta )$ so that we have

(2.18) \begin{align} \left ( \frac {\partial ^2}{\partial \hat {r}^2} + \frac {1}{\hat {r}} \frac {\partial }{\partial \hat {r}} + \frac {1}{\hat {r}^2} \frac {\partial }{\partial \theta ^2} \right ) \hat {\phi }(\hat {r},\theta ) = 0, \qquad \hat {r}_{j-1} \lt \hat {r} \lt \hat {r}_j,\,0 \lt \theta \lt {\textstyle \frac 12} \pi \end{align}

for $j=1,\ldots ,N$ , where $\hat {r}_j = r_{\kern-1pt j}/r_0$ , noting that the condition (2.13) applied to the channel walls is expressed as

(2.19) \begin{align} \frac {\partial \hat {\phi }}{\partial \hat {r}} = 0, \qquad \hat {r} = \hat {r}_j, \, 0 \lt \theta \lt {\textstyle \frac 12} \pi . \end{align}

Next, we make a multiple-scale approximation in which $\hat {R} = \hat {r}/\epsilon$ where $\epsilon = \hat {r}_{j}-\hat {r}_{j-1} \ll 1$ is the dimensionless distance between adjacent baffles ( $0 \lt \hat {R} \lt 1$ ) is defined as the fast (channel) scale and we write

(2.20) \begin{align} \hat {\phi }(\hat {r},\theta ) = \hat {\phi }_0(\hat {r}, \hat {R}, \theta ) + \epsilon \hat {\phi }_1(\hat {r}, \hat {R}, \theta ) + \epsilon ^2 \hat {\phi }_2(\hat {r}, \hat {R}, \theta ) + \ldots . \end{align}

Consequently, (2.18) is to be written as

(2.21) \begin{align} & \left ( \!\left (\! \frac {\partial }{\partial \hat {r}} + \frac {1}{\epsilon } \frac {\partial }{\partial \hat {R}}\! \right )^2 + \frac {1}{\hat {r}} \left (\! \frac {\partial }{\partial \hat {r}} + \frac {1}{\epsilon } \frac {\partial }{\partial \hat {R}} \! \right ) + \frac {1}{\hat {r}^2} \frac {\partial ^2}{\partial \theta ^2}\! \right ) \left (\! \hat {\phi }_0(\hat {r},\hat {R},\theta ) + \epsilon \hat {\phi }_1(\hat {r},\hat {R},\theta ) + \ldots \! \right ) = 0, \nonumber \\ & \qquad 0 \lt \hat {R} \lt 1,\,0 \lt \theta \lt {\textstyle \frac 12} \pi ,\end{align}

whilst (2.19) becomes

(2.22) \begin{align} \left ( \frac {\partial }{\partial \hat {r}} + \frac {1}{\epsilon } \frac {\partial }{\partial \hat {R}} \right ) \left (\hat {\phi }_0(\hat {r},\hat {R},\theta ) + \epsilon \hat {\phi }_1(\hat {r},\hat {R},\theta ) + \ldots \right ) = 0, \qquad \hat {R} = 0,1,\, 0 \lt \theta \lt {\textstyle \frac 12} \pi . \end{align}

From (2.21), (2.22) the $O(1/\epsilon ^2)$ problem for $\hat {\phi }_0$ is defined by

(2.23) \begin{align} \frac {\partial ^2 \hat {\phi }_0}{\partial \hat {R}^2} = 0, \qquad 0 \lt \hat {R} \lt 1 \end{align}

with $\partial _{\hat {R}} \hat {\phi }_0=0$ on $\hat {R}=0,1$ . The clearly integrates to give $\hat {\phi }_0 \equiv \hat {\phi }_0(\hat {r},\theta )$ (i.e. independent of $\hat {R}$ , the channel scale).

For the problem at $O(1/\epsilon )$ (2.21), (2.22) gives

(2.24) \begin{align} \frac {\partial ^2 \hat {\phi }_1}{\partial \hat {R}^2} = 0, \quad 0 \lt \hat {R} \lt 1, \qquad \mbox{with} \qquad \frac {\partial \hat {\phi }_1}{\partial \hat {R}} = - \frac {\partial \hat {\phi }_0}{\partial \hat {r}}, \quad \mbox{on $\hat {R}=0,1$} \end{align}

and so

(2.25) \begin{align} \hat {\phi }_1(\hat {r},\hat {R},\theta ) = - \hat {R} \frac {\partial \hat {\phi }_0}{\partial \hat {r}} + \chi _1(\hat {r},\theta ) , \end{align}

where $\chi _1$ is an arbitrary function. From (2.21) the governing equation at $O(1)$ turns out to be, after simplifying using (2.25),

(2.26) \begin{align} \frac {\partial ^2 \hat {\phi }_2}{\partial \hat {R}^2} - \frac {\partial ^2 \hat {\phi }_0}{\partial \hat {r}^2} + \frac {1}{\hat {r}^2} \frac {\partial ^2 \hat {\phi }_0}{\partial \theta ^2} = 0, \qquad 0 \lt \hat {R} \lt 1 \end{align}

and the associated boundary conditions are

(2.27) \begin{align} \frac {\partial \hat {\phi }_2}{\partial \hat {R}} = \hat {R} \frac {\partial ^2 \hat {\phi }_0}{\partial \hat {r}^2} - \frac {\partial \chi _1}{\partial \hat {r}}, \qquad \hat {R}=0,1. \end{align}

Integrating (2.26) over $0 \lt \hat {R} \lt 1$ and using (2.27) results in the leading-order governing equation for $\hat {\phi }_0$ , namely

(2.28) \begin{align} \frac {\partial ^2 \hat {\phi }_0}{\partial \theta ^2} = 0. \end{align}

This is what we might have expected as it implies that the leading-order flow is parallel to the channel walls (plug flow); (2.28) integrates to

(2.29) \begin{align} \hat {\phi }_0(\hat {r},\theta ) = B(r_0 \hat {r}) \theta + C(r_0 \hat {r}) , \end{align}

where $B$ , $C$ are arbitrary functions of the macroscale variable, $r$ . Meanwhile on the boundary $z=0$ , equivalent to $\theta = ( 1/2) \pi$ , the transformation applied to (2.17) gives

(2.30) \begin{align} \left ( K r_0 -\frac {1}{\hat {r}} \frac {\partial }{\partial \theta } \right ) \left ( \hat {\phi }_0 \left(\hat {r},{\textstyle \frac 12} \pi \right) + \ldots \right ) = - {\mathrm{i}} \frac {K \Delta }{\varLambda _j} \int _{0}^{1} \frac {1}{\hat {r}} \frac {\partial }{\partial \theta } \left ( \hat {\phi }_0 \left(\hat {r},{\textstyle \frac 12} \pi \right) + \ldots \right ) \, {\mathrm{d}} R \end{align}

and we write $\varLambda _j = K \Delta \tilde {\varLambda }(r_0 \hat {r})$ since it does not depend on $\hat {R}$ . Consequently,

(2.31) \begin{align} \tilde {\varLambda }(r) = \frac {\tilde {\lambda }_1(r)}{\sqrt {K (b-a)}} - {\mathrm{i}} \tilde {\lambda }_2(r) \end{align}

is a dimensionless parameter implying, from (2.11), (2.12), that

(2.32) \begin{align} \tilde {\lambda }_2(r) = \rho g H(r)/\rho _a c_a^2 , \end{align}

where $H(r)$ measures the height the air chamber as a function of position. Additionally,

(2.33) \begin{align} \tilde {\lambda }_1(r) = \frac {\rho \sqrt {g (b-a)}}{\Delta } \lambda _1(r), \end{align}

like $\tilde {\lambda }_2(r)$ , is frequency-independent and related to the geometry and mechanics of the turbine: we regard this as a free parameter, dependent on position $r$ , which controls the damping.

At leading order (2.30) is therefore

(2.34) \begin{align} \left ( 1 - \frac {{\mathrm{i}}}{\tilde {\varLambda }(r_0 \hat {r})} \right ) \frac {\partial \hat {\phi }_0}{\partial \theta } \left(\hat {r}, {\textstyle \frac 12} \pi \right) = K r_0 \hat {r} \hat {\phi }_0 \left(\hat {r}, {\textstyle \frac 12} \pi \right)\!. \end{align}

Now $\phi (r,\theta ) = \hat {\phi }(\hat {r},\theta ) \approx \hat {\phi }_0(\hat {r}, \theta )$ to leading order. So from (2.29) we have

(2.35) \begin{align} \phi (r,\theta ) = B(r) \theta + C(r), \qquad \mbox{for $0 \lt \theta \lt {\textstyle \frac 12} \pi $, $a \lt r \lt b$} , \end{align}

where

(2.36) \begin{align} \left (1 - \frac {{\mathrm{i}}}{\tilde {\varLambda }(r)} \right ) \frac {\partial \phi }{\partial \theta } = Kr \phi , \qquad \mbox{on $\theta = {\textstyle \frac 12} \pi $, $a \lt r \lt b$.} \end{align}

It is useful to define a function $U(r)$ which represents the horizontal fluid velocity along $x=0$ , the vertical opening of the duct, so that

(2.37) \begin{align} \left . \frac {1}{r} \frac {\partial \phi }{\partial \theta } \right |_{\theta = 0} = U(r), \qquad a \lt r \lt b. \end{align}

It follows that $B(r) = r U(r)$ and then using (2.36) we have

(2.38) \begin{align} \left ( 1 - \frac {{\mathrm{i}}}{\tilde {\varLambda }(r)} \right ) U(r) = K \left ( {\textstyle \frac 12} \pi r U(r) + C(r)\right )\! , \end{align}

which defines $C(r)$ and so, finally,

(2.39) \begin{align} \phi (r,\theta ) = r U(r) \left(\theta - {\textstyle \frac 12} \pi \right) + K^{-1} (1 - {\mathrm{i}}/\tilde {\varLambda }(r)) U(r) \end{align}

represents the general solution in $x \gt 0$ in terms of a single unknown function, $U(r)$ .

2.2. Solution in $x \lt 0$ and matching

In $x \lt 0$ we can write the general solution, satisfying (2.2), (2.7) and (2.13) on $z=h$ , using separation of variables (see Linton & McIver Reference Linton and McIver2001), as

(2.40) \begin{align} \phi (x,z) = \big( \mathrm{e}^{{{i}} k x} + R \mathrm{e}^{-{{i}} k x} \big) \psi _0(z) + \sum _{n=1}^\infty \frac {\psi _n(z) \mathrm{e}^{k_n x}}{N_n k_n h} \int _a^b U(z') \psi _n (z') \, {\mathrm{d}} z' , \end{align}

where, for $n \geq 1$ ,

(2.41) \begin{align} \psi _n(z) = \cos k_n (h-z)/\cos k_n h \end{align}

and $k_n$ are the positive real roots of $K = -k_n \tan k_n h$ , ordered in increasing size. The functions $\psi _n(z)$ are orthogonal, satisfying

(2.42) \begin{align} \frac {1}{h} \int _0^h \psi _m(z) \psi _n(z) \, {\mathrm{d}} z = N_n \delta _{mn}, \qquad m,n \geq 0 , \end{align}

which extends the definition of $k_n$ to include $n=0$ where $k_0 = -{\mathrm{i}} k$ and

(2.43) \begin{align} N_n = \frac 12 \left ( 1 + \frac {\sin 2 k_n h}{2 k_n h} \right ) \mbox{sec}^2 k_n h. \end{align}

Note that in the case $n=0$ this becomes

(2.44) \begin{align} N_0 = \frac 12 \left ( 1 + \frac {\sinh 2 k h}{2 k h} \right ) \mbox{sech}^2 k h. \end{align}

In writing (2.40) we have also satisfied the condition

(2.45) \begin{align} \left . \frac {\partial \phi }{\partial x} \right |_{x=0} = \left \{ \begin{array}{l} U(z), \qquad a \lt z \lt b, \\[6pt] 0, \qquad \{ 0 \lt z \lt a\} \cup \{ b \lt z \lt h \} ,\end{array} \right . \end{align}

provided

(2.46) \begin{align} {\mathrm{i}} k h N_0 (1-R) = \int _a^b U(z) \psi _0(z) \, {\mathrm{d}} z. \end{align}

All that remains is to match the pressure across $x=0$ , $a \lt z \lt b$ (noting that $r=z$ on $\theta = 0$ ) which, by (2.4), is proportional to $\phi$ and so from (2.39) and (2.40) we have

(2.47) \begin{align} \frac {\pi }{2} z U(z) - K^{-1} \big( 1- {\mathrm{i}} /\tilde {\varLambda }(z)\big) U(z) + (\mathcal{K} U)(z) = -(1+R) \psi _0(z), \qquad a \lt z \lt b , \end{align}

where

(2.48) \begin{align} (\mathcal{K}U)(z) \equiv \int _a^b U\big(z'\big) \kappa \big(z,z'\big) \, {\mathrm{d}} z' \end{align}

is a real, positive, symmetric integral operator with a kernel defined by

(2.49) \begin{align} \kappa (z,z') = \sum _{n=1}^\infty \frac { \psi _n(z) \psi _n(z')}{N_n k_n h}. \end{align}

Letting $u(z)$ satisfy

(2.50) \begin{align} \left [ \frac {\pi }{2} z - K^{-1}\big( 1- {\mathrm{i}} /\tilde {\varLambda }(z)\big) \right ] u(z) + (\mathcal{K} u)(z) = \psi _0(z), \qquad a \lt z \lt b , \end{align}

implies that $U(z) = -(1+R) u(z)$ satisfies (2.47). Substituting this relation into (2.46) gives

(2.51) \begin{align} {\mathrm{i}} khN_0 (1-R) = - (1+R) \mathcal{A}, \qquad \mbox{where} \quad \mathcal{A} = \int _a^b u(z) \psi _0(z) \, {\mathrm{d}} z. \end{align}

In other words, the reflection coefficient may be determined by

(2.52) \begin{align} R = -\frac {\mathcal{A} + {\mathrm{i}} k h N_0}{\mathcal{A} - {\mathrm{i}} kh N_0}. \end{align}

We remark that, when $\tilde {\varLambda }(z) \to \infty$ , equivalent to switching off the power take-off, the integral equation (2.50) becomes real and so the solution $u(z)$ and, hence $\mathcal{A}$ , is real. Then, $|R| = 1$ , implying no energy is absorbed by the device, as must be the case. The limit $\tilde {\varLambda }(z) \to 0$ is equivalent to replacing the power take-off with a rigid lid with no compression, the integral equation furnishes $U(z) = 0$ and $|R|=1$ . If just $\tilde {\lambda }_1 = 0$ but $\tilde {\lambda }_2 \neq 0$ then the air turbine is shut off but the air compressibility allows a non-zero flow in the duct. In this case the integral equation is real once again and so $|R|=1$ .

3. Efficiency

Under the inviscid assumption, the mean power, $W$ , absorbed at the OWC is calculated in two different ways. In the first method we equate $W$ to the flux of energy crossing a vertical boundary far from the device. That is

(3.1) \begin{align} W = \lim _{X \to -\infty } \frac {\omega }{2\pi } \int _{0}^{2\pi /\omega } \int _0^h \varPhi _x(X,z,t) (P(X,z,t)-P_a) \, {\mathrm{d}} z \, {\mathrm{d}} t. \end{align}

In frequency-domain variables

(3.2) \begin{align} W = \lim _{X \to -\infty } \frac 12 {\rm Re} \left \{ -\frac {{\mathrm{i}} \rho g^2 |A|^2}{\omega } \int _0^h \phi _x(X,z) \overline {\phi (X,z)} \, {\mathrm{d}} z \right \} = W_{\textit{inc}} (1 - |R|^2) \end{align}

after use of (2.14), where the overbar denotes complex conjugate and

(3.3) \begin{align} W_{\textit{inc}} = \frac 12 \rho g |A|^2 c_g \end{align}

is the incident wave energy in terms of the group velocity defined by

(3.4) \begin{align} c_g = \frac {{\rm d} \omega }{{\rm d} k} = \frac {g}{\omega } kh N_0 , \end{align}

where $N_0$ is defined by (2.44). We are interested in the efficiency of power capture $E \in [0,1]$ defined by $E = E_{\textit{ff}}$ where

(3.5) \begin{align} E_{\textit{ff}} = \frac {W}{W_{\textit{inc}}} = 1 - |R|^2 \end{align}

(the subscript $ \textit{ff} $ is used to indicate this is a calculation made in the far field). A second method of calculating mean power comes from measuring the mean rate of working of the flux against the pressure across the turbines. According to this definition in terms of the $N$ discrete channels

(3.6) \begin{align} W = \frac {\omega }{2 \pi } \int _0^{2 \pi /\omega } \sum _{j=1}^N \int _{r_{\kern-1pt j-1}}^{r_{\kern-1pt j}} \zeta _t(x,t) (P(x,0,t)-P_a) \, {\mathrm{d}} x \, {\mathrm{d}} t. \end{align}

Invoking homogenisation replaces discrete variables by continuous variables so that

(3.7) \begin{align} W = \frac 12 \frac {\rho g^2 |A|^2}{\omega } \textrm{Im} \left \{ - \int _a^b \phi _z(x,0) \overline {\phi (x,0)} \, {\mathrm{d}} x \right \}\!. \end{align}

Using (3.4), (3.5) gives the efficiency as $E= E_{\textit{nf}}$ ( $nf$ indicates this is a near-field calculation) where

(3.8) \begin{align} E_{\textit{nf}} = \frac {1}{kh N_0} \textrm{Im} \left \{ \int _a^b \frac {1}{r} \frac {\partial \phi }{\partial \theta } (r, \pi /2) \overline {\phi (r,\pi /2)} \, {\mathrm{d}} r \right \}\!. \end{align}

Inserting the definition (2.39) gives

(3.9) \begin{align} E_{\textit{nf}} = \frac {1}{ kh N_0} \textrm{Im} \left \{ \int _a^b \frac {{\mathrm{i}}}{K \tilde {\varLambda }(r)} |U(r)|^2 \, {\mathrm{d}} r \right \} = \frac {|1+R|^2}{K kh N_0} {\rm Re} \left \{ \int _a^b \frac {|u(z)|^2}{\tilde {\varLambda }(z)} \, {\mathrm{d}} z \right \} \end{align}

after using the relation $U(r) = -(1+R) u(z)$ . This simplifies in the case that $\tilde {\varLambda }$ is constant. We see that $\tilde {\lambda }_1$ controls the power take-off and the efficiency tends to zero when $\tilde {\lambda }_1 \to \infty$ . Note that if $\tilde {\lambda }_2 = 0$ then $E_{\textit{nf}} \to 0$ as $\tilde {\lambda }_1 \to 0$ since $u(z) \to 0$ also in this case. See Evans (Reference Evans1981) for more details on the efficiency of power capture from water waves.

4. Numerical approximation

In order to solve the integral equation (2.50) we use a Galerkin method in which the solution is approximated by the $(2M+1)$ -term expansion

(4.1) \begin{align} u(z) \approx \sum _{m=-M}^M c_m u_m(z) , \end{align}

where $c_m$ are coefficients to be determined and

(4.2) \begin{align} u_m(z) = \mathrm{e}^{{{i}} m \pi (2 z - b - a)/(b-a)} \end{align}

are the complex Fourier series basis for the interval $z \in (a,b)$ satisfying

(4.3) \begin{align} \int _{a}^b u_m(z) \overline {u_q(z)} \, {\mathrm{d}} z = (b-a) \delta _{qm}. \end{align}

This approximation method is applied to a similar geometry in Jan & Porter (Reference Jan and Porter2018). It is justified by the fact that homogenisation of the plate array in the duct means that the governing equation near the corner is no longer Laplace’s equation. Galerkin’s method applied to solutions of integral equations involves substituting (4.3) into (2.50) before multiplying by $\overline {u_q(z)}$ for $q=-M,\ldots ,M$ and integrating over $a \lt z \lt b$ . This results in the linear system of equations

(4.4) \begin{align} \sum _{m=-M}^M c_m \mathcal{K}_{qm} = \overline {F_{q0}}, \qquad q = -M,\ldots ,M \end{align}

defining $c_m$ where

(4.5) \begin{align} \mathcal{K}_{qm} = \frac {\pi }{2} T_{qm} - K^{-1}(b-a)\delta _{qm} + {\mathrm{i}} K^{-1} S_{qm} + \sum _{r=1}^\infty \frac {F_{mr} \overline {F_{qr}}}{N_r k_r h} \end{align}

and we have written

(4.6) \begin{align} F_{mr} = \int _a^b u_m(z) \psi _r(z) \, {\mathrm{d}} z \end{align}

with

(4.7) \begin{align} S_{qm} = \int _{a}^b \frac {u_m(z) \overline {u_q(z)}}{\tilde {\varLambda }(z)} \, {\mathrm{d}} z \qquad \mbox{and} \qquad T_{qm} = \int _{a}^b z u_m(z) \overline {u_q(z)} \, {\mathrm{d}} z. \end{align}

Explicit expressions for $F_{mr}$ and $T_{qm}$ are given in Appendix A. In the simplest case that $\tilde {\varLambda }$ is constant, $S_{qm} = (b-a)\delta _{qm}/\tilde {\varLambda }$ .

We now use (4.1) in (2.51) to get

(4.8) \begin{align} \mathcal{A} \approx \sum _{m=-M}^M c_m F_{m0} , \end{align}

allowing us to approximate $E_{\textit{ff}}$ from (3.5) once $R$ is calculated using (2.52).

The alternative calculation of efficiency uses (4.1) in (3.9) to get

(4.9) \begin{align} E_{\textit{nf}} \approx \frac {4 khN_0}{K|\mathcal{A} - {\mathrm{i}} kh N_0|^2} {\rm Re} \left \{ \sum _{n=-M}^M \sum _{m=-M}^M c_n \overline {c_m} S_{mn} \right \}\!. \end{align}

In the case that $\tilde {\lambda }_{1}$ and $\tilde {\lambda }_{2}$ are both constant we have, on account of (4.3),

(4.10) \begin{align} E_{\textit{nf}} \approx \frac {4 khN_0 (b-a)^{3/2}}{ K^{1/2}|\mathcal{A} - {\mathrm{i}} kh N_0|^2} \boldsymbol{\cdot }\frac { \tilde {\lambda }_1} {\tilde {\lambda }_1^2 + K (b-a) \tilde {\lambda }_2^2} \sum _{n=-M}^M |c_n|^2. \end{align}

In the work of Parry-Barnard (Reference Parry-Barnard2025) (not repeated here as it is too much of a diversion from our main goal) it is shown that approximations, made by implementing Galerkin’s method, to the solution of the integral equations imply the equivalence of $E_{\textit{nf}}$ and $E_{\textit{ff}}$ described previously. Numerical computations bear this out with only small differences due to rounding errors between computed values of $E_{\textit{nf}}$ and $E_{\textit{ff}}$ . This is true for any parameters used including the truncation parameter, $M$ . Thus, the two methods are useful for validating the numerical method, but cannot be used as a means of assessing the accuracy of the numerical results.

5. Viscous losses in the duct

Since the OWC design supposes a cascade of annular baffles separate the curved duct into narrow channels within which the fluid is subjected to resonant excitation, one might therefore be concerned about the role natural dissipative effects have on the device efficiency. Thus, we consider here a simple extension to the previous formulation of the problem which includes damping due to viscous wall effects.

In a full-scale design we reasonably envisage channels within the duct to be large enough that the plug-flow assumption holds across the bulk of each channel ( $\Delta \gg \sqrt {\nu /2 \omega }$ , where $\nu$ is the kinematic viscosity). That is, we assume that $N$ is not so large that the width of the channel is close to the width of the boundary layer; plausible values of $N$ are discussed in § 6.4. We remark that this idealisation of laminar flow with boundary-layer effects will not capture all of the complicated turbulent walls effects relating to oscillating flows in channels, although Batchelor (Reference Batchelor1967, p.353) notes that it is reasonable to assume that vorticity is confined to the walls in oscillating flows. Then the wall stress, per unit length, due to boundary-layer shear associated with the no-slip condition on the walls of the channel is easily calculated (from the solution to the well-known Stokes problem (see, for example, Batchelor Reference Batchelor1967) since wall curvature is negligible under the boundary-layer scaling) to be

(5.1) \begin{align} \left . \sigma _{R\theta } \right |_{R=0,1} = -{\mathrm{i}} \sqrt {\frac {\nu }{\omega }} \mathrm{e}^{-{{i}} \pi /4} U(r) \end{align}

(noting that we have applied the appropriate variable scalings stated in § 2 to the stress and to $U$ ). Therefore, its integrated effect along the length of each channel, from both walls, can be written

(5.2) \begin{align} -{\mathrm{i}} \pi r \sqrt {\frac {\nu }{\omega }} \mathrm{e}^{-{{i}} \pi /4} U(r) = \Delta p_{vis}(r) , \end{align}

which defines the jump in pressure needed to overcome viscous damping ( $\Delta$ is the channel width). This pressure drop needs to be added to the solution in the duct and means that continuity of pressure across $x=0$ is no longer equivalent to matching $\phi$ from both sides resulting in (2.47), but requires matching $\phi$ from $x \lt 0$ with $\phi - p_{vis}$ from $x \gt 0$ . Thus, instead of (2.47) we have

(5.3) \begin{align} \left [ \frac {\pi }{2} z + \frac {\pi z \tilde {\varGamma }\mathrm{e}^{{{i}} \pi /4} }{(K(b-a))^{1/4}} - \frac {\big( 1- {\mathrm{i}} /\tilde {\varLambda }(z)\big)}{K} \right ] U(z) + (\mathcal{K} U)(z) = -(1+R) \psi _0(z), \end{align}

for $a \lt z \lt b$ , where we have written

(5.4) \begin{align} \tilde {\varGamma } = \frac {\nu ^{1/2} N}{(b-a)^{3/4}g^{1/4}}, \end{align}

a dimensionless constant, independent of frequency, which determines the effect of the viscosity on the OWC design. As expected, $\tilde {\varGamma }$ increases linearly with the number of channels, $N$ .

The impact of the modification (5.3) to (2.47) is easy to follow and eventually the first term $(\pi /2) T_{qm}$ in (4.5) needs to be replaced by

(5.5) \begin{align} \frac {\pi }{2} \left ( 1 + \frac {2 \tilde {\varGamma } \mathrm{e}^{{{i}} \pi /4}} {(K (b-a))^{1/4}} \right ) T_{qm} \end{align}

in the numerical method.

The two calculations of efficiency made under inviscid assumptions are now no longer equal and $E_{\textit{nf}}$ is used to calculate the efficiency of power conversion by the device. The larger value of $E_{\textit{ff}}$ includes losses due to viscous damping. The work of Parry-Barnard (Reference Parry-Barnard2025) shows how the two formally equate via the integral equation (5.3) by calculating the rate of viscous dissipation due to the boundary layers in the channels but the lengthy details are omitted from this paper. The two different calculations of $E$ are useful since $E_{\textit{ff}}- E_{\textit{nf}}$ indicates the percentage of incident wave power lost to viscous wall effects.