A second-order theory for an array of curved wave energy converters in open sea

We present a second-order nonlinear theory for an array of curved flap-type wave energy converters (WECs) in open sea, excited by oblique incident waves. The flap model, when at rest, has a weak displacement of its wetted surface about the vertical plane. Using a perturbation-harmonic expansion technique, we decompose the set of nonlinear governing equations in a sequence of linearized boundary-value problems that can be investigated through analytical methods. The velocity potential at the leading order is solved in terms of elliptical coordinates and related Mathieu functions. Due to quadratic interactions of the first-order potential, a drift flow, a static angular displacement and a bound wave appear at the second order. At the same order the effect of the flap shape forces the first harmonic solution. We then compare an array of flat devices with an array of curved flaps, assuming several shape functions that could be of practical engineering interest. We show that hydrodynamic interactions between the curved flaps and incident waves can have constructive effects in terms of power absorption at large frequencies. © 2019ElsevierLtd.Allrightsreserved.


Introduction
This paper presents a second-order mathematical model for an array of curved flap-type wave energy converters.Second-order theories on floating and fixed marine structures have an important practical significance in ocean engineering.This is because quadratic products of the first-order solutions generate a drift flow and several harmonics that are related to unexpected forces over the system components (Mei et al., 2005).Analytical models concerning cylindrical structures involved several authors in the past (Molin, 1979;Chau and Eatock Taylor, 1992;Kim and Yue, 1989;Newman, 1996), while the second-order diffraction theory for elliptical cylinders has more recently been developed by Chatjigeorgiou and Mavrakos (2012).Higher-order theories involving floating offshore structures have been developed in several contexts.For example, Sammarco et al. (1997a,b) devised a third-order theory for the Venice gate barrier, while Michele et al. (2018a) developed a nonlinear theory for an array of oscillating gate-type WECs when excited subharmonically in a channel.
Floating bodies such as oscillating gate-type and flap-type devices are capable to absorb wave energy with large efficiency.This is the reason why this kind of devices received much attention by researchers and industry in recent years.Indeed, several theoretical studies and experiments have been carried out.For an extensive review we leave the reader to the works of Babarit et al. (2012), Dias et al. (2017) and Babarit (2018).
The hydrodynamics of a single flap and of arrays of flap-type WECs has been investigated theoretically, numerically and experimentally by several authors in different configurations.In particular, Renzi and Dias (2012) and Sammarco et al. (2013) analysed the behaviour of Oscillating Wave Surge Converters (OWSCs) in a long-infinite channel, while Renzi and Dias (2013), Renzi et al. (2014), Michele et al. (2015), Noad and Porter (2015), Michele et al. (2016b) and Xu et al. (2017) investigated the case of a single flap and of arrays of oscillating flaps in open sea.The effects of a perfectly reflecting straight wall on the OWSC's performance have been analysed by Sarkar et al. (2015) and Michele et al. (2016a), while investigations on viscous dissipations have been carried out by Wei et al. (2015) and Feng et al. (2018).
A large part of the analytical theories concerning the analysis of WECs mentioned above neglect higher order contributions.However, Michele et al. (2018a) showed that nonlinear resonances, including mode-mode competitions, can potentially increase power extraction.Moreover, recent analysis on curved gates showed that curved structures can both increase energy extraction efficiency (Michele et al., 2018b(Michele et al., , 2019;;Hodge et al., 2017) and affect significantly the exciting torque (Nokob and Yeung, 2015).
Motivated by these recent results, in this paper we investigate the combined effects of flap surface curvature and nonlinear contributions on the energy extraction by an array of WECs.We use a similar approach to the analysis of the Bragg scattering of incident waves (Mei et al., 2005(Mei et al., , 1988;;Alam et al., 2010), by introducing a small deviation of the flap wetted surface about the vertical plane.We show that flap curvature can change significantly the behaviour of the system.This is more evident at large frequencies when short waves interact with the flap shape function and when oblique waves excite odd and even out-of-phase modes.A nonlinear set of governing equations is derived and decomposed in a sequence of linearized boundary value problems through the use of a standard perturbation-harmonic expansion technique (Jordan and Smith, 2011).
At the first order, the oscillating flap-type device is excited by oblique incident waves.The solution in terms of velocity potential can be found via integral equations (Noad and Porter, 2015), hypersingular integral equations based on the Green function method (Renzi and Dias, 2013;Michele et al., 2015Michele et al., , 2016a) ) and analytical approaches in terms of elliptical coordinates and Mathieu functions (McLachan, 1951;Morse and Feshbach, 1981;Michele et al., 2016b).In this paper we use the analytical approach.This method allows us to reduce significantly the computational cost of numerical evaluations and to obtain explicit expressions for the added inertia, radiation damping and exciting torque.Furthermore, we derive the Haskind-Hanaoka relation in elliptical coordinates to check the numerical accuracy of the radiation and diffraction velocity potentials.
At the second order we show that the drift flow is forced by quadratic products of the first harmonic solution.The corresponding flap angular displacement and bound wave can be found without solving the velocity potential, because they only depend on the leading order terms.We show that in the case of a single flap, the drift flow makes the WEC rotate landward, i.e. towards the direction of propagation of the incoming waves.This angular displacement represents the mean value of the WEC oscillations and can be of engineering interest for design purposes.In the case of an array of several flaps, the static displacement shape can include flaps rotating either landward or seaward.When the incident waves frequency is small, the radiation potential is related to in-phase motions and the diffraction potential sets the rotation of each flap.In this case, the angular rotation axis tends to be orthogonal to the incident waves direction.At large frequencies, the effect of even and odd natural modes manifests and the hydrodynamics affecting the static angular displacement increases in complexity.
Finally, we investigate the first-harmonic problem at the second order forced by the array shape function and its derivatives.Solution is then found again via Mathieu functions.We remark that this problem does not occur for flat devices, because in that case the boundary value problem would be damped, unforced and identical to the first-order radiation problem.We show that the array shape yields either constructive or disruptive effects in terms of power extraction efficiency, depending on the concavity of the flap surface.This aspect becomes more relevant in the case of multiple flaps, oblique waves and large frequencies.

Mathematical model
Referring to Fig. 1, consider an in-line array of Q neighbouring flap-type OWSCs hinged on a rigid bottom.Let us define a Cartesian reference system with the x ′ and y ′ axes lying on the undisturbed free surface level z ′ = 0 and the z ′ axis pointing upward.The physical variables are indicated by primes.The horizontal rigid bottom has coordinate z ′ = −h ′ and the OWSCs oscillate about a horizontal axis located at z ′ = −h ′ , x ′ = 0. Let a ′ be the width of each flap, then w ′ = a ′ Q represents the entire width of the array.Let t ′ be the time, G q , q = 1, . . ., Q , denote the qth flap and Θ ′ q ( t ′ ) be the angular displacement of G q positive anticlockwise.Each flap G q spans a y ′ -width given by (1) as the angular displacement piecewise function of the array, defined in . The array has negligible thickness if compared to w ′ , thus the hypothesis for thin obstacles can be applied to evaluate the wave field (Linton and McIver, 2001;Renzi and Dias, 2012).The fluid is assumed inviscid and incompressible and the flow irrotational, hence the velocity potential satisfies the Laplace equation in the entire fluid domain Ω ′ ( x ′ , y ′ , z ′ ) . Let us introduce for convenience a local reference system which is integral with the oscillating flap G q .In this non-inertial reference frame the wetted flap surfaces are (2) where δ ′ denotes the deviation of the flap surface from the rotating plane X ′ = 0, while X ′ and Z ′ are respectively: (3) Let A ′ be the amplitude of the incident waves making an angle Ψ with the x ′ -axis, ω ′ the incident wave frequency, λ ′ the wavelength and g ′ the acceleration due to gravity.Then introduce the following non-dimensional quantities (Michele et al., 2018a(Michele et al., , 2019;;Sammarco et al., 1997a,b): where ζ ′ represents the free surface elevation, δ ′ g is the length scale for δ ′ and G the non-dimensional frequency.Let us now introduce the following two length ratios: (5) The Laplace and Bernoulli equations in non-dimensional variables are then where the subscripts denote differentiation with respect to the relevant variable, p ′ is the fluid pressure and ρ ′ the fluid density.The dynamic and mixed boundary conditions on the free surface are given by while the no-flux condition at the horizontal bottom requires The qth flap surface in the reference frame (X, y, Z ) can be written as hence the corresponding kinematic boundary condition is given by Finally, the equation of motion of the curved flap with a linear generator at the hinge reads IϵΘ q,tt + GC sin ϵΘ q + ν pto ϵΘ q,t = where the superscript + (−) indicates that the corresponding function is evaluated on the right (left) side, is the non-dimensional inertia of the flap about the hinge, C = C ′ /ρ ′ λ ′4 g ′ is the non-dimensional buoyancy restoring torque and ν pto = ν ′ pto /ρ ′ λ ′5 ω ′ is the non-dimensional power take off (PTO) coefficient.The free-surface boundary conditions ( 8)-( 9) can be Taylor expanded about z = 0 up to O ( ϵ 2 ) : Similarly, Taylor expansion of the flap surface expression (11) and of the kinematic boundary condition ( 12) about x = 0 yields for y ∈ [−w/2, w/2]: The equation of motion ( 13) can be written as: where ∆(•) = (•) + − (•) − denotes the difference of (•) on two sides of the flap.Let us assume the following perturbation expansion for the unknowns up the second order O (ϵ): Substitution of the above expansions into the set of governing equations and boundary conditions yields for n = 1, 2: Laplace equation: Free surface dynamic condition: where Free surface mixed condition: where No-flux boundary condition at the bottom: Kinematic condition on the array surface: where Equation of motion of the qth flap: where The combination of harmonic forcing terms at higher orders implies that the response of the array includes several harmonics (Jordan and Smith, 2011).In the next section we apply harmonic-expansion techniques to the unknowns in order to remove the time dependence and simplify the mathematics.

First and second-order solutions
Because higher order solutions imply higher harmonics, we return to physical variables, omit the primes for convenience and assume the following harmonic expansions for the real unknowns: where i denotes the imaginary unit while the symbol * indicates the complex conjugate of the first terms on the righthand side.Substitution of the latter expansions into the governing equation ( 23) and boundary conditions ( 24)-( 30) gives a sequence of linear boundary value problems of order n and harmonic m: The latter must be solved together with the equation of motion for each flap G q (31), which upon substitution of (34) becomes In the latter, the forcing terms F nm , B nm , G nm and D nm 's are defined for each order n and harmonic m in the following subsections.
Since the angular response of the qth flap is the corresponding average extracted power becomes The contribution of the second-harmonic term θ q,22 is of order O ( ϵ 2 ) , and so its effect on the generated power ( 42) is negligible at the second order O (ϵ).Therefore, the second-harmonic problem is not considered here because of its small effects.On the contrary, the first-harmonic terms θ q,11 and θ q,21 are dominant and can be optimized by requesting that both have the same phase in the complex plane.In this case, the double products between real and imaginary parts of the flap amplitude give an extra term of order O (2ϵ), which is not negligible at the second order.
Once the generated power is obtained, the efficiency of the system is quantified by considering the capture factor C F , defined as the ratio between P( 42) and the incident wave energy flux per array width w: where the term ) . (44) In the latter expression, E represents the averaged total energy per unit surface, while C g denotes the group velocity (Mei et al., 2005).
We are now in a position to solve the boundary value problems at the leading order O (1) and second order O (ϵ).

Leading order problem -first harmonic response to incident waves
The incident wave field is assumed to be at the leading order O (1), hence it forces the first harmonic problem.Let the total solution be where φ I is the velocity potential of the monochromatic incident waves of amplitude A and angle of incidence φ S is the scattered velocity potential, φ R q denotes the radiation velocity potential due to the qth moving flap while all the others are at rest, while the term k 0 is the real root of the dispersion relation (see Mei et al., 2005) Both φ R q and φ S satisfy Eqs. ( 35)-( 38) with null forcing terms F 11 = B 11 = 0, while the kinematic boundary conditions on the array (39), read elsewhere on the array surface Finally, it is required that φ R q and φ S are outgoing for This boundary-value problem has been extensively studied in the past.Semi-analytical solutions based on integral and hypersingular integral equations are already available in the literature (Renzi and Dias, 2013;Michele et al., 2015).Recently, Michele et al. (2016b) showed that the problem at the leading order admits separable solutions expressed in terms of elliptic coordinates (χ, ϑ) defined as and Mathieu functions.Michele et al.'s theory has been successfully validated against existing analytical, numerical and experimental data.
Mathieu functions have the advantage of decreasing significantly the numerical tasks and reducing the algebraic manipulations with respect to methods based on integral and hypersingular integral equations.For this reason, their applications can be found in several surface wave diffraction phenomena of the first (Chen and Mei, 1973;Williams and Darwiche, 1988;Chatjigeorgiou and Mavrakos, 2010;Chatjigeorgiou, 2011) and second order (Chatjigeorgiou and Mavrakos, 2012) and other fields involving elliptic structures (Stamnes and Spjelkavik, 1995;Ruby, 1995;Gutiérrez-Vega, 2000;Gutiérrez-Vega et al., 2003).Hereafter we summarize the method of solution adopted in Michele et al. (2016b).
Since separation of the z-dependence use of both separation (53) and transformation from rectangular coordinates (x, y) to elliptical coordinates (χ, ϑ) modifies the Laplace equation ( 35) to the following two-dimensional Helmholtz equation in elliptic coordinates: where τ n = w 2 k 2 n /16 is a real non-dimensional parameter that can be greater (propagating waves) or smaller (evanescent waves) than zero.The general solution of (54) bounded for χ → ∞ corresponds to Michele et al. (2016b), McLachan (1951) and Morse and Feshbach (1981) ϕ R qn (χ , ϑ) ϕ S (χ , ϑ) where ce (ϑ) and se (ϑ) are called respectively cosine-type and sine-type Mathieu functions of the first kind (or angular Mathieu functions), He (1) (χ) and Ho (1) (χ) are the even and odd Hankel-Mathieu functions of the first kind (or radial Mathieu functions), while A, B, C and D denote the unknown coefficients to be determined by applying the boundary conditions on the array.Mathieu functions are equivalent to Bessel functions in cylindrical coordinates and share similar interesting properties.We leave the reader to the book of McLachan (1951) for an extensive analysis.Unfortunately, angular and radial Mathieu Functions are not present in standard software and there are few textbooks regarding their numerical evaluation.For a detailed numerical scheme we refer to the works of Michele et al. (2016b) and Gutiérrez-Vega (2000).Substitution of the expansion ( 53) into ( 50)-( 51) yields the boundary conditions expressed in terms of elliptical coordinates: elsewhere on the array surface ∂ϕ S ∂χ = Awgk 0 cos Ψ e −ik 0 (w/2) sin Ψ cos ϑ sin ϑ where while S q defines the boundary of the moving flap G q expressed in the plane (χ, ϑ): In the latter expression Since the Mathieu Functions satisfy the orthogonality condition (McLachan, 1951) integration of ( 56)-( 57) yields: in which the subscript (•) χ denotes the derivative of (•) with respect to the radial coordinate χ, B (m) j is the jth coefficient for the angular Mathieu function se m (see Michele et al., 2016b, for details) and J j is the Bessel Function of the first kind and order j.To determine the flap response, we now use the equation of motion (40): where { θ 11 } is a column vector of length Q that includes all the angular displacements of the flaps θ q,11 , q = 1, . . ., Q , I is the square identity matrix of size Q × Q , M and N are respectively the added inertia matrix and the radiation damping matrix of size Q × Q as well: in which represent, respectively, the added inertia and the radiation damping of the flap G q due to the unit rotation of the flap G p .The forcing term { F } is a column vector of size Q including the values of the exciting torque acting on each flap G q : We now focus the attention on the accuracy of the numerical computations for the hydrodynamic parameters F and ν by applying the well-known Haskind-Hanaoka identity.Let us consider the simple case of Q = 1.The asymptotic behaviour of the radiation potential in the far field for unit rotational velocity of the flap is given by where r = we χ /4 is the radius expressed in terms of the radial elliptic coordinate χ, while: represents the angular variation of the radially spreading wave.The term s 2m+1 in (73) has the following expression (McLachan, 1951;Michele et al., 2016b): ) . ( The three-dimensional Haskind-Hanaoka relation therefore reads ) . (75) Unlike (Michele et al., 2016b), we disregard the assumption of sin ϑ ∼ se m (ϑ), hence the general identity (75) can be also applied for large values of τ 0 .Eqs. ( 71) and ( 75) will be used in Section 4 to check the accuracy of the numerical evaluations.

Second order problem O (ϵ) -zeroth harmonic drift
The quadratic products of the first-order solution give the following forcing terms of the drift flow, respectively, on the free surface and on the flap: ] . ( The forcing terms B 20 and D 20 allow us to evaluate directly the real expression for the bound wave independent on time (or wave set up/set down): and the expression for the static angular displacement: This latter term indicates that the flap G q oscillates about θ q,20 and not about θ q = 0 as in the linearized theory.Moreover θ q,20 is a real number, indeed all the integrands above are referred to absolute values of complex functions.

Second order problem O (ϵ) -first harmonic and flap shape effects
This problem represents the effects of the flap shape on the total wave field.Indeed, the corresponding forcing terms are: The boundary conditions on z = 0 and z = −h remain like those of the first-harmonic O (1) problem.This suggests that the solution can be expanded in terms of hyperbolic cosine eigenfunctions: To simplify the algebra, we decompose ϕ 21,n as follows: where Γ n satisfies the following boundary condition on the flap elsewhere on the array surface while the boundary condition for Λ n reads The boundary value problem for Γ n is identical to that for ϕ R qn , thus, we do not need to solve it again.The velocity potential Λ n represents the forcing term for the flap oscillations θ 21 .Its general solution can be found by applying elliptical coordinates: where each constant is now given by Note that, unlike the leading order problem analysed before, the constants A Λ nm and B Λ nm that multiply the even Mathieu functions in (87) differ from zero.This is because the integrals in ( 88)-( 91) do not have orthogonal properties in the interval ϑ ∈ [0, 2π ].
Once Λ n is given, we must apply the equation of motion to the flap G q to find the response θ q,21 .The equation of motion in matrix form is in which the matrices I, M and N are equal to those introduced in (68), while each element F q,21 , q = 1, . . ., Q , inside the column vector The latter is the complex exciting torque for the first harmonic motion at the second order, and represents the action due to the flap profile.In the absence of flap curvature, i.e. δ = δ y = δ z = 0, the exciting torque becomes F q,21 = 0 and consequently θ q,21 = 0, q = 1, . . ., Q .Therefore, the flap curvature has a direct effect on the first-harmonic solution at the second order.On the other hand, the terms on the left side of (92) do not depend on the flap curvature.Hence the angular displacement θ q,21 is maximum if the exciting torque function F q,21 is maximized by choosing appropriate shape functions.

Results and discussion for the case Q = 1
In this section we discuss the simple case of a single flap in open sea.The solution for the leading order radiation potential (55) simplifies as follows while the solution for the leading order scattering potential is given by se m (ϑ) Ho ( 1) m (χ) We remark that the first-order terms are not affected by the small curvature of the flap.
Let us assume water depth h = 10 and flap width w = a = 20 m.For these flap dimensions, buoyancy and inertia are considered respectively C = 3.53 × 10 7 kg m 2 s −2 and I = 2.6 × 10 6 kg m 2 .The ratio between the incident wave amplitude A and the sea depth h must be of order O (ϵ), hence we assume A = 1 m.The corresponding behaviour of added inertia µ, radiation damping ν and magnitude of the exciting torque F for Ψ = 0 rad and Ψ = π/6 rad versus the frequency ω ∈ [0.25, 4] rad s −1 is shown in Fig. 2. Expressions ( 71) and ( 75) can now be used to check the numerical evaluations for both radiation and diffraction potentials.As in Renzi and Dias (2013), let us define the relative error ε as where 'l.h.s.' represents the exciting torque given by (71), while 'r.h.s.' denotes the exciting torque given by the Haskind-Hanaoka relation (75).Taking max m = 5 in the expression for the velocity potentials (55), we obtain a maximum relative error ε = O ( 10 −15 ) . This is a single-degree-of-freedom system, thus it has one eigenfrequency at the leading order (Sammarco et al., 2013).Indeed, the eigenfrequency satisfies the following implicit equation in ω: with µ being dependent on ω as well.Solution of the latter expression yields ω = 0.57 rad s −1 .Maximization of the power output for the resonance condition requires that the ν pto coefficient is equal to the radiation damping for ω = ω (Mei et al., 2005).Numerically, we have ν pto = 9 × 10 6 kg m 2 s −1 (see also Fig. 2(a)).
Once the parameters are fixed, we can evaluate the static displacement due to the second order drift, which again is not affected by the small curvature of the flap.Fig. 3 shows the behaviour of θ 20 versus the frequency ω for Ψ = 0 rad and Ψ = π /6 rad.Both curves correspond to positive values, i.e. the flap equilibrium position is rotated anticlockwise in the same direction of the incoming waves.The curvature of the flap affects the second-order solution of the first harmonic, which in turn has an impact on the extracted power (42).To investigate such effects, we compare the flat device with null curvature (δ = δ y = δ z = 0) with six different configurations in terms of power extraction.We remark that δ g /w ∼ O (ϵ), thus we assume the maximum value of the array surface deviation from X ′ = 0 be equal to w/10 m.Let us introduce the following flap configurations of practical engineering interest described by the following expressions: in which configurations δ 1,2 depend on the y-coordinate, δ 3,4 depend on the vertical coordinate z while δ 5,6 depend both on y and z.Figs. 4(a) and 4(b) show the behaviour of the capture factor C F of each configuration versus ω, respectively, for normal incident waves with Ψ = 0 rad and oblique waves with Ψ = π/6 rad.In all cases, the efficiency decreases with the incidence angle Ψ .The effect of the curvature disappears for oblique waves, thus the efficiency is practically the same as that of a flat body.Note also that the maximum value of C F correspond to the theoretical maximum value predicted by Michele et al. (2016b) when sin ϑ ∼ se m (ϑ): in which k 0 ≃ 0.061 m −1 is the wavenumber related to the eigenfrequency ω = ω.
At low frequencies the effect of the flap surface shape is small.Its constructive or destructive effects arise mainly at large frequencies.This is because long waves do not interact with a small δ and the array appears to them as a flat oscillating body.Configurations 4 and 6 are more efficient than the flat device for large frequencies ω > 1.5 rad s −1 .These configurations correspond to concave surfaces in (x, y, z), hence a phenomenon similar to the focusing of radiated and diffracted waves occurs (Stamnes, 1986).The least beneficial effects occur in Configuration 3 and 5.We now turn to investigate an array of Q = 5 flaps.m 2 .Non-trivial solutions of the following unforced and undamped nonlinear system det

Results and discussion
give one in-phase mode N (ω 1 ), 2 even out-of-phase modes (N 1 , N 3 ) and 2 odd out-of-phase modes (N 2 , N 4 ) with related eigenfrequencies (see Sammarco et al., 2013;Michele et al., 2015Michele et al., , 2016a,b, for details),b, for details).As in the previous section, let us consider normal incident waves (Ψ = 0 rad) and oblique waves with Ψ = π/6 rad.In the first case, only the even modes can be excited, whereas in the second case both even and odd modes can be excited because symmetry of the exciting torque is lost.The numerical values of the eigenfrequencies and periods for each mode are listed in Table 1.Note that the in-phase solution has the same eigenfrequency of the single flap analysed in Section 4, thus we expect that the array behaves as a single flap when N (ω 1 ) is resonated by incident waves.Now, let us maximize power extraction for ω = 0.57 rad s −1 by imposing that the PTO coefficient is ν pto = 1.8 × 10 6 kg m 2 s −1 and evaluate the value of the static angular displacement θ q,20 for each flap G q .Figs. 5(a) and 5(b) show the behaviour of θ q,20 , q = 1, . . ., Q , versus the frequency of the incident waves for Ψ = 0 rad and Ψ = π/6 rad, respectively.Unlike the case of a single flap, here θ q,20 can be negative also for normal incidence.At small frequencies the exciting torque plays a predominant role in the flap displacement.This is the reason why the flaps located on the left side are rotated towards x > 0 when Ψ = π/6 rad (see again Fig. 5(b)).When the frequency ω increases, the effects of the natural out-of-phase modes start to arise and the dynamics increases in complexity.
We now consider the same six configurations δ i , i = 1, . . ., 6, analysed in the previous section and investigate the variation of the capture factor for different incident wave frequencies.The behaviour of C F versus ω for each configuration is shown in Fig. 6(a) for normal incidence (Ψ = 0 rad) and Fig. 6(b) for oblique incoming waves (Ψ = π/6 rad).
Constructive and disruptive interactions appear mainly at larger frequencies characterized by short waves.As in the case analysed before, configurations δ 4 and δ 6 are the most efficient when ω > 1.5 rad s −1 ., shows that significant energy extraction occurs at large frequencies when the array is excited by inclined waves with respect to the case of a single flap.Furthermore, differences between flat and curved arrays are appreciable even in the case of oblique incidence (see again Fig. 6(b)).On the other hand, at small frequencies the overall array motion returns in-phase, the wavelength starts to increase and the array behaves as a single large flat body.Indeed, the maximum values of the capture factor almost coincide with the maxima of the single flap showed in Figs.4(a) and 4(b).

Conclusion
We investigated the first and second order dynamics of an array of flap-type curved WECs in open sea of constant depth.The theoretical solution has been pursued in terms of angular Mathieu functions and Hankel-Mathieu functions of the first kind.This method is fast converging and becomes more useful than existing semi-analytical models when extensions to higher orders are carried out.Moreover, hydrodynamic parameters such us added inertia, radiation damping and exciting torque can be evaluated explicitly.
Perturbation-harmonic expansion up to the second order allowed us to find the forcing terms that depend on the array shape function.First, we assumed monochromatic incident waves of constant amplitude and arbitrary direction to be present at the leading order.Then, we optimized power extraction by fixing the PTO coefficient to be equal to the radiation damping when synchronous resonance occurs.Numerical accuracy has been checked by deriving the Haskind-Hanaoka identity in elliptical coordinates.
Next, we obtained the static angular displacement of each flap and the bound wave at the second order.We remark that the corresponding expressions depend on the leading order terms only and that can be found without passing through the drift flow solution.We compared the single flap and the array of flaps and found that the static angular displacement has a strong dependence on the wave direction and the occurrence of out-of-phase motion.
Finally we analysed the effects of the flap shape by solving the first-harmonic problem at the second order.We found that the corresponding velocity potential cannot exist for flat devices because forcing terms depend on interactions between the flap profile and the wave field at the leading order.We investigated array configurations that can be of interest for engineering and design purposes and compared them with a standard flat device of the same size.
Our results show that curved arrays can have beneficial effects at large frequencies when short waves interact with the array curvature.We found that in the case of normally incident waves the differences between the single flap and the array are minimal in terms of power extraction.On the contrary, for oblique incident waves the array is more efficient than the single flap at large frequencies.Indeed, differences in the capture factor behaviour of a flat array and a system made of curved flaps are appreciable.
We also found that the vertical shape of each element of the array is more important than the horizontal curvature in maximizing power extraction.In particular, vertical concave configurations analysed in this work showed the maximum efficiency, while configurations that varies with the horizontal coordinate did not give significant beneficial effects.
The results obtained in this work shows the relevance of higher-order contributions, combined with curved flap surfaces, for energy extraction purposes.The zeroth-harmonic effects at the second order clearly do not influence energy extraction but can be of engineering interest in order to optimize the system layout and to predict unwanted behaviours such as the excessive inclination of the flaps.Moreover, the presented analytical model is numerically fast and can be implemented in algorithms able to seek the optimal shape function that maximizes energy extraction efficiency for a particular configuration, depending on a large number of parameters.For this reason, a parametric analysis is necessary and will be addressed in future works.This will allow us to obtain a set of optimized configurations that can be of interest to the wave energy research community and industry.
We assumed that the flap surface function is small when compared to the spatial coordinate scales.Thus the ratio between the wavelength and the length scale for the flap shape must be small.When this condition stops to be valid, the array system at the first order cannot be approximated as a flat vertical plane and the theory fails.Solution to this problem can be potentially found by following semi-analytical approaches based on Green's function methods and integral equations, however the analysis at higher orders may result algebraically complex and cumbersome.This aspect is rather important and should be taken into account to evaluate the nonlinear response of floating flaps in open sea when the curved surface profile is not small anymore.
Finally, waves in real seas require representations that include broad banded spectra.Since the problem analysed here is nonlinear and the incident frequency in irregular waves is not fixed, the harmonic expansion applied in this work cannot be extended to the case of random seas.Indeed, complex interactions arise and a new approach based on the Fourier formulation is needed (Mei et al., 2005).

Fig. 1 .
Fig. 1.Plan geometry of the array and side view of the oscillating flap.

Fig. 2 .
Fig. 2. (a) Added inertia µ, radiation damping ν and (b) magnitude of the exciting torque |F | versus incident wave frequency ω.The exciting torque diminishes when the angle of incidence Ψ increases.

Fig. 3 .
Fig. 3. Static rotation of a single flap θ 20 due to quadratic products of the velocity potential φ 11 at the second order.For both configurations, the static rotation is positive and anticlockwise.

Fig. 4 .Fig. 5 .
Fig. 4. Behaviour of the capture factor C F of each configuration versus incident wave frequency ω for two different values of the angle of incidence Ψ .(a) Ψ = 0 rad.(b) Ψ = π/6 rad.The power take-off coefficient is optimized for the eigenfrequency ω = 0.57 rad s −1 .

Table 1
Eigenfrequencies and period for each natural mode.