A coupled damped harmonic oscillator model is derived for the long-time motion of an array of floating cylinders moving in heave, forced by a prescribed initial displacement. Eigenmodes of the coupled damped harmonic oscillator model are set to be the complex resonances of the corresponding linear potential theory model of the array. Thus, the solution of the coupled damped harmonic oscillator model has the same form as the singularity expansion method. The derivation of the coupled damped harmonic oscillator model is underpinned by a homotopy method to calculate the complex resonances for arbitrary arrays. A recursive homotopy method is used to study the complex resonances as one of the geometrical parameters is varied. The coupled damped harmonic oscillator model is shown to give accurate approximations of the full linear solutions, except for exceptional cases in which additional complex resonances appear.

## I. INTRODUCTION

Floating cylinders moving in vertical heave motion are a classical model of “point-absorber” wave energy converters^{1–3} and other floating bodies, and are still commonly used in contemporary studies.^{4–6} The cylinders (point absorbers) are designed to resonate in response to surface water waves while being much smaller than typical wavelengths, such that the waves they scatter are dominated by radiated waves. To provide a source of commercially viable offshore renewable energy, they are deployed in arrays with the aim to create wave interactions between the cylinders that give a positive array effect, i.e., to improve the overall energy capture in comparison with the devices in isolation.^{1,7,8}

Full linear (potential flow) theory is the standard benchmark for model predictions of hydrodynamics of floating bodies (such as heaving cylinders)^{4,9,10} and wave interactions between the bodies.^{8,11–13} Full linear interaction theory involves multiple wave scatterings between the bodies so that solution methods are either numerical^{14,15} or semi-analytical.^{7,16–18} Solutions are typically obtained in the frequency domain, where the equations of motion for the bodies can be written in the form of a coupled damped harmonic oscillator (CDHO) model, but with frequency dependent coefficients owing to the coupling with the water (see Sec. II B). A standard CDHO (with constant coefficients) is attractive due to its simplicity and the intuitive nature of its coefficients. This motivated De Chowdhury and Manasseh^{19} to propose a CDHO model for an array of “oscillating water column” wave energy converters, in which the coupling coefficients for wave interactions between the cylinders were found with an approach developed to analyze bubble acoustics.^{20} They used the CDHO model to approximate eigenmodes of the array, although it is unclear how these relate to full linear theory.

Fitzgerald and McIver^{21} developed a DHO approximation for the cognate problem of a “moonpool” within a motionless structure that encloses part of the free surface. The approximation is based on the complex resonances (or scattering frequencies) for the problem, which are the unforced solutions in the frequency domain, analytically extended into complex frequency space. They showed that the DHO model provides accurate approximations for resonant behaviors dominated by a single complex resonance. The DHO model^{21} is closely related to a degenerate (single term) form of the singularity expansion method (SEM), in which a the long-time solution of the full linear problem is expressed as a linear superposition of the complex resonances. The SEM has a long history in water wave problems, dating back to Ursell^{22} and Maskell and Ursell.^{23}

In this article, we outline a novel CDHO model for wave radiation by an array of heaving truncated cylinders based on the SEM. The key advance in comparison to the approach of Fitzgerald and McIver^{21} is to derive coupling coefficients for an array that set the eigenmodes of the CDHO to be the complex resonances. The main computational challenge addressed is calculating the complex resonances, which, in general, relies on searching for the resonances in complex-frequency space. Wolgamot *et al.*^{24} developed an efficient method to calculate the complex resonances when the array satisfies certain symmetries. In contrast, we develop a homotopy method for arbitrary arrays (cylinder dimensions and locations), in which the complex resonances evolve from real-valued resonances of the uncoupled array. The homotopy method extends the method proposed by Bennetts and Meylan^{25} but for the more challenging setting in which the resonant frequencies of the initial guesses (the uncoupled array) overlap, so that the mapping to the complex resonances is not clearly defined *a priori*. We show that the CDHO model accurately approximates the long-time solution to the full-linear problem, but that it breaks down in rare cases where an additional complex resonance appears with a frequency close to one of those of the existing complex resonances.

## II. PRELIMINARIES

### A. Geometry and governing equations

Consider water of density *ρ*, constant finite depth *h*, and infinite horizontal extent. A Cartesian coordinate system $ x = ( x , y , z ) \u2208 \Omega $ defines locations in the water, where (*x*, *y*) is the horizontal coordinate and *z* is the upward-pointing vertical coordinate, with origin set to coincide with the undisturbed water surface. An array of *N* circular cylinders, with arbitrary (non-overlapping) configuration, float on the surface of the water, where cylinder *p* has mass *M _{p}*, radius

*a*, and draught

_{p}*d*. The wetted surface of cylinder

_{p}*p*is denoted Γ

_{p}with unit normal $ n \u0302 p$ directed into the water from the cylinder. The free surface of the water, i.e., the surface not occupied by the array, is denoted $ \Gamma f$ (see Fig. 1).

*p*denoted $ \zeta p ( t )$ and its velocity $ \zeta \u0307 p ( t )$, where the overdot denotes differentiation with respect to time. At

*t*= 0, the cylinders are released from prescribed initial displacements

*Z*, giving the initial conditions

_{p}^{9}Assuming that the water is inviscid, incompressible, and undergoing irrotational motions, the velocity potential satisfies Laplace's equation,

*g*is the constant of gravitational acceleration, $ \u2202 \u2009 / \u2009 \u2202 n \u2261 n \u0302 p \xb7 \u2207$ is the normal derivative, linearized Bernoulli pressure is assumed for (2c), and the domain Ω and boundaries $ \Gamma \u2022$ are unchanged from their equilibrium states, according to linear theory. Cylinder displacements satisfy the equations of motion, i.e., the dynamical cylinder-coupling conditions (coupling of motion of the

*p*th cylinder to the surrounding water),

### B. Spectral solution

*ω*represents angular frequency. The frequency-domain potential also satisfies Laplace's equation and the bed condition, respectively,

*n*is in the normal direction ( $ n \u0302 p$). In the far field, the frequency-domain potential satisfies the Sommerfeld radiation condition

^{24}

**A**and

**B**have entries $ { A} p , q = A p , q$ and $ { B} p , q = B p , q$ ( $ p , q = 1 , 2 , \u2026 , N$), which are defined via

^{9}

*q*. Therefore, the equations of motion (9) can be expressed in the form of a CDHO with frequency-dependent coefficients, i.e.,

**A**and

**B**are referred to as the added mass and damping matrices, respectively.

The hydrodynamic coefficients in **D** are calculated using eigenfunction matching^{26} and diffraction transfer matrices^{27} for the individual cylinders, and Graf interaction theory for wave interactions between the cylinders.^{16,28} Figure 2 validates the computation of the hydrodynamic coefficients with those of Siddorn and Eatock Taylor^{17} in terms of the radiation impedance matrix $ \u2212 i \u2009 \omega \u2009 D$.

^{24}

*t*< 0 have been used.

## III. COMPLEX RESONANCES AND HOMOTOPY METHOD

### A. Complex resonances

*N*vector with unitary

*p*th entry.

^{24}

**D**over real part of a complex frequency, assuming the matrix

**D**is analytic at

*ω*. The form of the residuals used in (19b) is based on the Laurent series of $ \zeta \u0302$ around the complex zeros of $ det { Q}$, assuming simple poles, i.e., (19b) approximates the frequency-domain solution close to $ \omega = \omega p$, as derived by Ref. 30. Note that, in general, $ \zeta SEM ( 0 ) \u2261 Re \u2009 { Z} \u2260 Z$, where

_{p}### B. Direct homotopy method

*j*= 0) has the known uncoupled solution (Sec. III A). The initial solution is mapped to the solution at the first step (

*j*= 1), and so on, until the

*J*th step, which is the fully coupled problem ( $ \u210f = 1$). An iterative method is applied to make the mapping at each step, such that

*q*denotes the iterations. The iterations are based on the assumption that

^{20,25}

*p*= 1 is such that

*p*= 2 is set as

Figure 4(a) shows the homotopy steps to recover the complex resonant frequency (first row of Table I) for a single cylinder (draught $ d = 2 \u2009 a$ in water of depth $ h = 4 \u2009 a$) superimposed on a contour plot of $ | det { Q} |$. Four steps (*J* = 4) are shown, although the complex resonant frequency (and associated eigenvector) can be obtained with only two steps (*J* = 2). Similarly, Fig. 4(b) shows the homotopy steps to recover the two complex resonant frequencies for a pair of identical cylinders with spacing $ 4 \u2009 a$ (second row of Table I). As the cylinders are identical, the resonant frequencies for the uncoupled system are the same, and the choice of the initial eigenvectors is important to evolve to the correct complex resonant frequencies of the coupled system.

$ \nu \xaf 1 = 0.7071$ | $ \omega \xaf 1 = 0.6225 \u2212 0.0111 i$ | |

$ \nu \xaf 1 , 2 = 0.7071$ | $ \omega \xaf 1 = 0.6197 \u2212 0.0067 i$ | $ \omega \xaf 2 = 0.6260 \u2212 0.0191 i$ |

$ \nu \xaf 1 = 0.7071$ | $ \omega \xaf 1 = 0.6225 \u2212 0.0111 i$ | |

$ \nu \xaf 1 , 2 = 0.7071$ | $ \omega \xaf 1 = 0.6197 \u2212 0.0067 i$ | $ \omega \xaf 2 = 0.6260 \u2212 0.0191 i$ |

### C. Recursive homotopy method and validation using symmetry analysis

**Q**using a frequency-independent transformation matrix $ T = [ \tau 1 , \tau 2 , \u2026 , \tau N ]$, such that

^{24}

**T**, $ \tau 1 = [ 1 ; 1 ]$, is the symmetric eigenfunction, i.e., where the two cylinders move with the same amplitude and in-phase with one another. The second column, $ \tau 2 = [ 1 ; \u2212 1 ]$, is the anti-symmetric solution, i.e., where the cylinders have the same amplitude but are out of phase.

Figure 5 shows the complex resonant frequencies for an array with two identical cylinders versus the center-to-center cylinder separation, $\u2113$. A recursive homotopy method is used, such that the complex resonances at the smallest separation, $ \u2113 = 4 \u2009 a$, are calculated using the direct homotopy method outlined in Sec. III B. These are used as initial guesses for the second smallest separation, $ \u2113 = 4 \u2009 a + \Delta \u2113$, where $ \Delta \u2113 = 0.01$, and so on. The recursive homotopy method, which is used to find the complex resonances as functions of a parameter (here the cylinder separation), is more efficient than using the direct homotopy (Sec. III B) at each point on the abscissa. Values calculated using the method in Sec. III B directly are shown to validate the calculations, along with values found from the symmetry method, i.e., by solving $ \lambda 1 ( \omega ) = \lambda 2 ( \omega ) = 0$.

### D. Example: Random positional perturbations

*N*= 4 identical cylinders with radius

*a*, arranged in a square array with side lengths $ \u2113 = 4 \u2009 a$. The complex resonances are calculated by the direct homotopy method and verified using the symmetry approach, using the transformation matrix

^{24}

*ω*( $ p = 1 , 2 , 3 , 4$), where $ \omega 3 = \omega 4$.

_{p}Random perturbations are introduced to the cylinder locations in the horizontal plane, where the perturbations in the *x*- and *y*-directions for each cylinder are chosen randomly from a uniform distribution over the interval $ \epsilon \u2009 \u2113 \u2009 [ \u2212 0.5 , 0.5 ]$ for a prescribed perturbation strength *ε*, similar to Ref. 31. Therefore, the array loses symmetry and the symmetry approach cannot be used to find the complex resonances. The direct homotopy method (Sec. III B) can be used to find the complex resonances of the perturbed array, but it is more efficient to use the complex resonances of the unperturbed array as the initial guess for a single homotopy step.

Figure 6 shows the complex resonant frequencies for four different perturbation strengths ( $ \epsilon = 0.01 , 0.03 , 0.07$, and 0.1), with twenty-five realizations of the perturbed array for each perturbation strength. As expected, larger perturbation strengths lead to broader ranges of the resonant frequencies. For $ \omega 1$ and $ \omega 2$, the deviations in the real and imaginary components of the complex resonances appear to be linearly correlated, which can also be seen for certain resonances in the results of Wolgamot *et al.*^{24} for a perturbed array of nine cylinders (their Fig. 7). The complex resonance are more scattered for $ \omega 3$ and $ \omega 4$. The random perturbations cause some complex resonances to move closer to the real axis, even for the complex resonance closest to the real axis ( $ \omega 2$), which contrasts with the findings of Ref. 24.

## IV. COUPLED DAMPED HARMONIC OSCILLATOR METHOD

### A. Damped harmonic oscillator model for the isolated cylinder

*ω*

_{1}is the complex resonant frequency for the specified single cylinder and $ \zeta \u0303$ depends on the (as yet unspecified) initial conditions. Substituting Eq. (42) into (41), and enforcing the roots of the characteristic equation to be the conjugate pairs of the complex resonant frequency

*ω*

_{1}, gives the coefficients of the ODE as

### B. Coupled damped harmonic oscillator model for arrays

*N*> 1), the DHO (41) is applied to each body, along with coupling between the bodies, to give the CDHO model

**d**, $ e \u2208 \mathbb{R} N \xd7 N$ are to be determined. In order to derive the coupling matrices, let the solution to Eq. (44) be a superposition of eigenmodes, such that

**X**, $ Y \u2208 \u2102 N \xd7 N$ as

**d**and

**e**, are found from (50a), as

The homotopy method may yield $ Im \u2009 { X} = 0$ when the array is symmetric, e.g., two identical cylinders or four identical cylinders in a square array, so that Eq. (50) are undefined. These cases can be rectified by scaling the columns of **X** by an arbitrary (non-real) complex number. In practice, the need to scale the eigenvectors is avoided by constructing the coefficients of the DHO [Eq. (44)] based on the complex eigenfrequencies and then invoking in-built function polyeig in MATLAB^{®} to find the suitable complex eigenvectors.

### C. Example: Non-symmetric array

Consider the problem of two cylinders (*N* = 2) with the same draught but with different radii (*a _{p}* for

*p*= 1, 2), such that the array is non-symmetric, i.e., it cannot be transformed using the method outlined in Sec. III C. In fact, results not shown demonstrate that a diagonalization of

**Q**for this problem, in the form (37) with $ \lambda j = \omega \u2212 \omega j$ (

*j*= 1, 2), leads to a frequency dependent matrix

**T**, indicating that a suitable transformation does not exist. Figure 7 shows the coupling coefficients and complex-resonant frequencies for a non-symmetric, two cylinders array versus the radius ratio $ a 1 \u2009 / \u2009 a 2$, where

*a*

_{1}is varied and

*a*

_{2}is held fixed. The results are found using the direct homotopy method in Sec. III B for $ a 1 = a 2$, and using these complex resonances as a starting point for a recursive homotopy method to find complex resonances over a surrounding interval of radius ratios. The complex resonant frequencies for the non-symmetric array are shown alongside the complex resonant frequencies for the cylinders in isolation (i.e.,

*N*= 1).

When cylinder *p* = 1 is much smaller than *p* = 2 (approximately $ a 1 \u2009 / \u2009 a 2 < 0.75$), waves radiated by the smaller cylinder do not excite motion of the larger cylinder, such that *d*_{21} and $ e 21 \u2248 0$, although the larger cylinder affects the smaller cylinder by reflecting waves, which is expressed through relatively large absolute values of *d*_{11} and *e*_{11}. It follows that one complex resonant frequency (*ω*_{1}) is similar to the complex resonant frequency of the smaller cylinder (*p* = 1) in isolation, and the corresponding eigenvector $ \u2248 e 1$ (not shown). Waves radiated by the larger cylinder excite motion of the smaller cylinder (*d*_{12} and *e*_{12} are non-zero) but the resulting waves created by the smaller cylinder do not influence the larger cylinder (*d*_{22} and $ e 22 \u2248 0$). Thus, the other complex resonant frequency (*ω*_{2}) is similar to the complex resonant frequency of the larger cylinder (*p* = 2) in isolation, but the corresponding eigenvector ≉ $ e 2$.

As cylinder *p* = 1 increases in radius towards the radius of cylinder *p* = 2, and the complex resonant frequencies become closer, the coupling coefficients are all non-zero and no particular coupling coefficients are dominant. In this intermediate regime, approximately $ a 1 \u2009 / \u2009 a 2 \u2208 ( 0.85 , 1.25 )$, the real parts of the complex resonant frequencies for the non-symmetric array become near-locked, i.e., they closely follow one another, whereas the imaginary parts veer away from one another before coming back together. The corresponding eigenvectors vary rapidly with respect to the amplitude ratio, passing through the symmetric and anti-symmetric solutions for equal radii [cf. Eq. (40)].

For radius ratios above the intermediate regime, i.e., when cylinder *p* = 1 is much larger than the other, the complex resonances are again close to those of the isolated cylinders. Complex resonant frequency *ω*_{1}, i.e., the one associated to motion of the smaller cylinder (*p* = 1) in isolation for small $ a 1 \u2009 / \u2009 a 2$, is also close to the complex resonant frequency of the smaller cylinder (*p* = 2) for large values of the radius ration with corresponding eigenvector $ \u2248 e 2$. Similarly, complex resonant frequency *ω*_{2} is close to the complex resonant frequency for the larger cylinder for both small and large radius ratios. Therefore, as the radius ratio moves from small to large, the complex resonances for the array switch between cylinders. If the direct homotopy method is applied at the different radius ratios, the initial guesses corresponding to the uncoupled problems for the individual cylinders switch between the continuously evolving solutions obtained from the recursive homotopy method in the intermediate regime.

### D. Initial conditions

**[Eq. (10b)], initial conditions for the CDHO are set such that the solution of Eq. (44) is accurate for large times, i.e., the CDHO is consistent with the SEM. Thus, using the SEM (19), the initial conditions are**

*Z*### E. Assessment against full linear solution

The CDHO solutions for two identical cylinders with radius *a* and spacing $ \u2113 = 4 \u2009 a$ (whose complex resonances were studied in Sec. III C), with initial displacement vector $Z$ such that $ Z = [ 1 , 0 ] T$, are shown against the full linear solution (Fig. 8), where the full linear solution is given in Eq. (15a) and derived from potential-flow theory including multiple wave scattering, as described in Sec. II B. The CDHO solutions are identical to the SEM (not shown). As expected, the CDHO does not satisfy the same initial displacements as the full linear solution. However, even at early times, the CDHO is almost exact with only very small differences from the full linear solution for $ t \xaf \u2261 \u2009 g \u2009 / \u2009 a \u2009 t < 5$. As expected, the CDHO tends towards the full linear solution as time increases. Moreover, the CDHO is ∼300 times faster to calculate than the full linear solution. (The full linear solution took 1123.74 s and the CDHO took 3.8, 3.2 s for the complex resonances and 0.6 s for the time integrator, on an ASUS-Zephyrus-G502D laptop with 24 GB RAM and a 2.3 GHz AMD Ryzen processor with linux ubuntu 18.04 lts.) Similarly, the CDHO for an array of four identical cylinders in a square with side lengths $ \u2113 = 4 \u2009 a$ (whose complex resonances were studied in Sec. III D) and, with $Z$ such that $ Z = [ 1 , 0 , 0 , 0 ] T$, is identical to the SEM (not shown) and almost identical to the full linear solution (Fig. 9), although the differences from the full linear solution are greater and extend up to approximately $ t \xaf < 15$ for the cylinders initially at rest.

The previous examples are indicative of the general accuracy of the CDHO. However, when the spacing of the square array becomes large, e.g., when $ \u2113 = 33.08 \u2009 a$, the CDHO loses accuracy [Fig. 10(a)]. The loss of accuracy is caused by another complex resonant frequency appearing close to *ω*_{1} [Fig. 10(b)], which is not one of the complex resonant frequencies included in the CDHO, i.e., *ω _{p}* (

*p*= 2, 3, 4). The additional complex resonant frequency is associated with the eigenvector $ \zeta \u0302 1$, i.e., the same as

*ω*

_{1}. It first appears for $ \u2113 \u2248 33.08 \u2009 a$ but is initially not close to the other complex resonant frequencies [Fig. 10(c)], so that it does not noticeably affect the accuracy of the CDHO. The motion associated with the additional complex resonance can be incorporated in the SEM by summing over all five complex resonances in (19) to generate an accurate approximation [Fig. 10(a)]. The CDHO framework does not permit the additional complex resonance to be incorporated, and, thus, its accuracy is limited in these rare cases.

## V. CONCLUSIONS

An efficient CDHO model has been derived for the long-time motions of heaving cylinder arrays forced by prescribed initial displacements. The model extends cognate previous work (i) to arrays using coupling coefficients between the cylinders and (ii) by linking the model to full linear potential flow theory. The coefficients of the CDHO were set such that the eigenmodes are the complex resonances of the array, and the coupling coefficients were shown to imply the behavior of the array. The CDHO initial conditions are set to coincide with the SEM initial displacements, so that the CDHO solution is tends to that of full linear theory for large times, although it was shown that the CDHO often gives reasonable approximations for early times.

A direct homotopy method was developed to calculate the complex resonances for arbitrary array, in which the starting points are the resonances of the uncoupled arrays that exist at real-valued frequencies. The direct homotopy method was adapted from a previous study, such that it could be applied when the resonant frequencies of the uncoupled problem overlap, i.e., for identical cylinders.

Further, based on the direct homotopy approach, a recursive homotopy method was developed for efficient calculation of the complex resonances as a parameter (e.g., cylinder separation) is varied. The method provides a benchmark for future developments of automated calculation of complex resonances.

An example involving an array with relatively large cylinder separations was also given to show that the CDHO model loses accuracy if an additional complex resonance appears that is excited by the initial displacements. However, the CDHO model typically gives accurate approximations. It is available as the basis for designing point absorber arrays with optimal energy capture properties.

## ACKNOWLEDGMENTS

This work was supported by the Australian Research Council (Nos. LP180101109 and FT190100404). L.G.B. thanks the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the Theory and Application of Multiple Wave Scattering programme when work on this paper was undertaken and supported by EPSRC under Grant No. EP/R014604/1. The first author acknowledges Australian Renewable Energy Agency (ARENA) who supported a post-doctoral fellowship with the Emerging Renewable Program Grant No. A00575 during 2015–2018.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Swapnadip De Chowdhury:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Luke G. Bennetts:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Richard Manasseh:** Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

## REFERENCES

*Power from Sea Waves*

*Handbook of Mathematical Techniques for Wave/Structure Interactions*

*Advances in the Analysis and Design of Marine Structures*