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.

Floating cylinders moving in vertical heave motion are a classical model of “point-absorber” wave energy converters1–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 numerical14,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 Manasseh19 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 McIver21 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 model21 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 Ursell22 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 McIver21 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 Meylan25 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.

Consider water of density ρ, constant finite depth h, and infinite horizontal extent. A Cartesian coordinate system $x = ( x , y , z ) ∈ Ω$ 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 Mp, radius ap, and draught dp. The wetted surface of cylinder p is denoted Γp with unit normal $n ̂ 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 $Γ f$ (see Fig. 1).

FIG. 1.

Schematic of a floating truncated cylinder in finite water depth.

FIG. 1.

Schematic of a floating truncated cylinder in finite water depth.

Close modal
The cylinders move in heave only, with the vertical displacement of cylinder p denoted $ζ p ( t )$ and its velocity $ζ ̇ p ( t )$, where the overdot denotes differentiation with respect to time. At t = 0, the cylinders are released from prescribed initial displacements Zp, giving the initial conditions
(1)
Linear wave theory is used to model the resulting motion of the cylinders, in terms of a velocity potential $Φ ( x , t )$.9 Assuming that the water is inviscid, incompressible, and undergoing irrotational motions, the velocity potential satisfies Laplace's equation,
(2a)
and boundary conditions that represent kinematic coupling of water to cylinder motions, free surface motion, no normal flow at the bed, and the far-field condition, respectively,
(2b)
(2c)
(2d)
(2e)
where g is the constant of gravitational acceleration, $∂ / ∂ n ≡ n ̂ p · ∇$ is the normal derivative, linearized Bernoulli pressure is assumed for (2c), and the domain Ω and boundaries $Γ •$ 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 pth cylinder to the surrounding water),
(3)
where $C p = ρ g a p 2 π$, again assuming linearized Bernoulli pressure.
The velocity potential and displacements, $Φ ( x , t ) ∈ ℝ$ and $ζ p ( t ) ∈ ℝ ( p = 1 , … , N )$, are related to their Fourier transforms, $Φ ̂ ( x , ω ) ∈ ℂ$ and $ζ ̂ p ( ω ) ∈ ℂ ( p = 1 , … , N )$, via the expressions
(4a)
and
(4b)
where the transform variable ω represents angular frequency. The frequency-domain potential also satisfies Laplace's equation and the bed condition, respectively,
(5)
It satisfies the frequency-domain free-surface condition,
(6)
and is coupled to the frequency-domain displacements, $ζ ̂ p ( p = 1 , … , N )$, via
(7)
where the derivative with respect to n is in the normal direction ($n ̂ p$). In the far field, the frequency-domain potential satisfies the Sommerfeld radiation condition
(8)
where the wavenumber $k ∈ ℝ +$ satisfies the dispersion relation $k tanh ( k h ) = ω 2 / g$.
The frequency-domain equations of motion are expressed in the matrix–vector form as24
(9)
where the vector
(10a)
contains the frequency-domain cylinder displacements, and
(10b)
contains the initial displacements that provide the frequency-domain forcing on the right-hand side of Eq. (9). The mass and buoyancy matrices are
(11)
The frequency-dependent hydrodynamic matrix is usually expressed as
(12)
where A and B have entries ${ A } p , q = A p , q$ and ${ B } p , q = B p , q$ ($p , q = 1 , 2 , … , N$), which are defined via9
(13)
in which $ϕ q ( x , ω )$ is the frequency-domain potential that satisfies the multiple scattering problem forced by unit displacement in cylinder q. Therefore, the equations of motion (9) can be expressed in the form of a CDHO with frequency-dependent coefficients, i.e.,
(14)
which is why A and B are referred to as the added mass and damping matrices, respectively.

The hydrodynamic coefficients in D are calculated using eigenfunction matching26 and diffraction transfer matrices27 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 Taylor17 in terms of the radiation impedance matrix $− i ω D$.

FIG. 2.

Comparison with digitized results from Siddorn and Eatock Taylor17 of the (a) and (c) imaginary parts and (b) and (d) real parts of the elements of the radiation impedance matrix, (a) and (b) $− i ω D 1 , 1$ and (c) and (d) $− i ω D 1 , 3$, for N = 4 identical cylinders with radius a and draught $d = 2 a$, in a square array with side length $ℓ = 4 a$, in $h = 4 a$ depth water.

FIG. 2.

Comparison with digitized results from Siddorn and Eatock Taylor17 of the (a) and (c) imaginary parts and (b) and (d) real parts of the elements of the radiation impedance matrix, (a) and (b) $− i ω D 1 , 1$ and (c) and (d) $− i ω D 1 , 3$, for N = 4 identical cylinders with radius a and draught $d = 2 a$, in a square array with side length $ℓ = 4 a$, in $h = 4 a$ depth water.

Close modal
The time-domain displacement, $ζ ( t )$, is the inverse Fourier transform of the solution in the frequency domain, $ζ ̂ = ζ ̂ ( ω )$, such that24
(15a)
(15b)
where the properties that the time-domain displacements are real-valued and zero for t < 0 have been used.
If the cylinders are uncoupled from the water (and hence each other), the equations of motion (9) degenerate to the system of forced harmonic oscillator equations,
(16)
The uncoupled system has resonances at the (generalized) eigenfrequencies $ω = ν p ∈ ℝ +$ ($p = 1 , … , N$), where
(17)
i.e., non-trivial unforced solutions (for $Z = 0$; also denoted $ζ ̂$) exist at the eigenfrequencies and forced solutions are unbounded. The eigenfrequencies are the zeros of $det { Q unc }$, as is shown in Fig. 3(a). As the cylinders are uncoupled, the corresponding eigenvectors are $ζ ̂ = e p$ ($p = 1 , … , N$), where $e p$ is a standard unit vector, i.e., a length N vector with unitary pth entry.
FIG. 3.

Resonant frequencies for a single cylinder N = 1, with draught $d = 2 a$ in water of depth $h = 4 a$. (a) Absolute determinant values of the Q-matrices versus non-dimensional frequency, $ω ¯ = ω a / g$, for the uncoupled problem ($| det { Q unc } |$; black curve) and the coupled problem ($| det { Q } |$; red), where $a / g = 0.5$ s is used in all calculations. (b) Contour plot of $| det { Q } |$ in complex-frequency space, with zero indicated ($ω 1 ∈ ℂ$; black bullet), and with zero of uncoupled problem shown for reference (ν1; black circle). Values of the non-dimensional complex resonances are given in Table I.

FIG. 3.

Resonant frequencies for a single cylinder N = 1, with draught $d = 2 a$ in water of depth $h = 4 a$. (a) Absolute determinant values of the Q-matrices versus non-dimensional frequency, $ω ¯ = ω a / g$, for the uncoupled problem ($| det { Q unc } |$; black curve) and the coupled problem ($| det { Q } |$; red), where $a / g = 0.5$ s is used in all calculations. (b) Contour plot of $| det { Q } |$ in complex-frequency space, with zero indicated ($ω 1 ∈ ℂ$; black bullet), and with zero of uncoupled problem shown for reference (ν1; black circle). Values of the non-dimensional complex resonances are given in Table I.

Close modal
Coupling to the water domain allows energy to be carried away from the cylinders in the form of radiated waves (known as radiation damping), so that the homogeneous version of Eq. (9) ($Z = 0$), i.e.,
(18)
(in general) does not possess real-valued eigenfrequencies [Fig. 3(a)]. Taking the analytic extension of $Q ( ω )$ to complex-valued frequencies, the eigenfrequencies, i.e., the zeros of $det { Q }$, are found to be shifted into the lower-half of the complex-frequency plane [Fig. 3(b)]. In the context of water wave theory, the eigenfrequencies are referred as complex resonant (or scattering) frequencies. The small negative imaginary components of the complex resonant frequencies are associated to decay with increasing time, so that causality is obeyed. The corresponding eigenvectors, $ζ ̂ = ζ ̂ p$, are non-trivial as the cylinders are coupled.
The complex resonances are used to capture the long-time motion of the cylinders. Following Refs. 23, 24, 29, and 30 and others, the integral on the right-hand side of Eq. (15a) is approximated by closing the integration contour in the fourth quadrant of the complex frequency plane. Applying Cauchy's residue theorem, and assuming that integrals along complex contours tend rapidly to zero in time, gives the singularity expansion method (SEM),24
(19a)
(19b)
The prime superscript denotes a complex derivative, which is calculated by applying a central differencing scheme on the real and the imaginary parts of D over real part of a complex frequency, assuming the matrix D is analytic at ωp. The form of the residuals used in (19b) is based on the Laurent series of $ζ ̂$ around the complex zeros of $det { Q }$, assuming simple poles, i.e., (19b) approximates the frequency-domain solution close to $ω = ω p$, as derived by Ref. 30. Note that, in general, $ζ SEM ( 0 ) ≡ Re { Z } ≠ Z$, where
(20)
due to the components of the full linear solution neglected in the SEM.
In principle, finding the eigensolutions to the coupled problem involves a costly search in the lower half of the complex frequency plane for the zeros of $det { Q }$, followed by retrieval of the associated eigenvectors. Methods are not typically discussed in the existing literature. Here, a homotopy method is proposed to connect the pairs of eigensolutions (complex resonant frequencies $ω p ∈ ℂ$ and eigenvectors $u p ∈ ℂ N × 1$) with the real-valued resonances of the corresponding uncoupled problem ($ν p ∈ ℝ$ and $e p ∈ ℝ N × 1$). Thus, Eq. (18) is modified to
(21)
where $ℏ ∈ [ 0 , 1 ]$ is the homotopy parameter, such that $ℏ = 0$ is the uncoupled problem and $ℏ = 1$ is the fully coupled problem. The aim is to evolve the uncoupled eigensolutions into the fully coupled eigensolutions, i.e., to find complex resonant frequencies, denoted $ω = w p ( ℏ )$, and corresponding eigenvectors $ζ ̂ = u p ( ℏ )$ ($p = 1 , … , N$) of Eq. (21), at each homotopy step, such that
(22)
The homotopy is implemented in finite steps of the parameter $ℏ$ by uniformly sampling the unit interval at the points $ℏ = ℏ j$ ($j = 1 , 2 , … , J$), where
(23)
The initial step (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 Jth step, which is the fully coupled problem ($ℏ = 1$). An iterative method is applied to make the mapping at each step, such that
(24)
(25)
where the superscript q denotes the iterations. The iterations are based on the assumption that20,25
(26)
Substituting Eq. (26) into Eq. (21) and performing a Taylor series expansion around $w p ( q )$ up to first order give
(27)
Equation (27) can be rewritten as
(28)
so that
(29)
Hence,
(30)
where $p p ( q ) = − inv { Q ℏ j ( w p ( q ) ) } Q ′ ℏ j ( w p ( q ) ) u p ( q )$.
When there are repeated resonant frequencies of the uncoupled problem, i.e., when some of the cylinders are identical, the homotopy between the uncoupled and coupled problems defined above is not necessarily one-to-one. In this case, the initial eigenvectors are set to be orthogonal to any already obtained eigenvectors for the coupled problem. Thus, suppose $ν 1 = ν 2 = … = ν q ≡ ν$ and that $w 1 ( 1 ) = ω 1$ and $u 1 ( 1 ) = ζ 1$ are obtained from the homotopy by using $w 1 ( 0 ) = ν$ and $u 1 ( 0 ) = e 1$. Then, set $w 2 ( 0 ) = ν$ and
(31)
where $c p ( 2 )$ are chosen to give orthogonality of $u 2 ( 0 )$ and $ζ ̂ 1$, such that
(32)
to obtain $w 2 ( 1 ) = ω 2$ and $u 2 ( 1 ) = ζ ̂ 2$. Similarly, set
(33)
where $c p ( 3 )$ are chosen to give orthogonality, such that
(34)
and so on. For example, for two identical cylinders ($ν 1 = ν 2 ≡ ν$), say the homotopy for p = 1 is such that
(35a)
(35b)
then the $ℏ = 0$ case for p = 2 is set as
(36)

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 a$ in water of depth $h = 4 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 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.

FIG. 4.

Examples of homotopy steps (symbols) in complex-frequency space, superimposed on corresponding contour plots of $| det { Q } |$, for (a) N = 1 and (b) N = 2 identical cylinders, with $d = 2 a$ and $h = 4 a$. The contours denote relatively large values (yellows) to relatively small values (blues), similar to Fig. 3.

FIG. 4.

Examples of homotopy steps (symbols) in complex-frequency space, superimposed on corresponding contour plots of $| det { Q } |$, for (a) N = 1 and (b) N = 2 identical cylinders, with $d = 2 a$ and $h = 4 a$. The contours denote relatively large values (yellows) to relatively small values (blues), similar to Fig. 3.

Close modal
TABLE I.

Non-dimensional complex resonant frequencies shown in Fig. 4(a) [row one; N = 1; also in Fig. 3(b)] and Fig. 4(b) (row two; N = 2).

 $ν ¯ 1 = 0.7071$ $ω ¯ 1 = 0.6225 − 0.0111 i$ $ν ¯ 1 , 2 = 0.7071$ $ω ¯ 1 = 0.6197 − 0.0067 i$ $ω ¯ 2 = 0.6260 − 0.0191 i$
 $ν ¯ 1 = 0.7071$ $ω ¯ 1 = 0.6225 − 0.0111 i$ $ν ¯ 1 , 2 = 0.7071$ $ω ¯ 1 = 0.6197 − 0.0067 i$ $ω ¯ 2 = 0.6260 − 0.0191 i$
In certain situations where the array contains symmetry, it is possible to diagonalize Q using a frequency-independent transformation matrix $T = [ τ 1 , τ 2 , … , τ N ]$, such that24
(37)
where $Λ ( ω ) = diag { λ 1 ( ω ) , λ 1 ( ω ) , … , λ N ( ω ) }$ is a diagonal matrix. Therefore, it is possible to find the eigenfrequencies as the zeros $ω = ω p$ of the scalar equations
(38)
The columns of the transformation matrix are the corresponding eigenvectors, i.e., $τ p = ζ ̂ p$ ($p = 1 , … , N$).
The simplest example is an array with two identical cylinders, for which the transformation matrix is
(39)
The first column of T, $τ 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, $τ 2 = [ 1 ; − 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, $ℓ$. A recursive homotopy method is used, such that the complex resonances at the smallest separation, $ℓ = 4 a$, are calculated using the direct homotopy method outlined in Sec. III B. These are used as initial guesses for the second smallest separation, $ℓ = 4 a + Δ ℓ$, where $Δ ℓ = 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 $λ 1 ( ω ) = λ 2 ( ω ) = 0$.

FIG. 5.

Validation of homotopy methods (the direct homotopy method described in Sec. III B, and the recursive homotopy method described in Sec. III C) to calculate (a) real and (b) imaginary parts of complex resonant frequencies, for N = 2 identical cylinders, with $d = 2 a , h = 4 a$ and center-to-center spacing $ℓ$.

FIG. 5.

Validation of homotopy methods (the direct homotopy method described in Sec. III B, and the recursive homotopy method described in Sec. III C) to calculate (a) real and (b) imaginary parts of complex resonant frequencies, for N = 2 identical cylinders, with $d = 2 a , h = 4 a$ and center-to-center spacing $ℓ$.

Close modal
Consider an array of N = 4 identical cylinders with radius a, arranged in a square array with side lengths $ℓ = 4 a$. The complex resonances are calculated by the direct homotopy method and verified using the symmetry approach, using the transformation matrix24
(40)
The eigenvectors $τ p$ have corresponding complex resonant frequencies ωp ($p = 1 , 2 , 3 , 4$), where $ω 3 = ω 4$.

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 $ε ℓ [ − 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 ($ε = 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 $ω 1$ and $ω 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 $ω 3$ and $ω 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 ($ω 2$), which contrasts with the findings of Ref. 24.

FIG. 6.

Deviations in the locations of the complex resonant frequencies (a) $ω 1$, (b) $ω 2$, (c) $ω 3$, and (d) $ω 4$, in the square array due to perturbation strengths ε.

FIG. 6.

Deviations in the locations of the complex resonant frequencies (a) $ω 1$, (b) $ω 2$, (c) $ω 3$, and (d) $ω 4$, in the square array due to perturbation strengths ε.

Close modal
FIG. 7.

Variation versus radius ratio of (a) and (b) the CDHO coupling coefficients and (c) and (d) complex resonant frequencies, for two cylinders with radii a1 and a2, spacing $ℓ = 4 a 2$, and where $d = 2 a 2$ and $h = 4 a 2$, where single cylinder and identical cylinder references are given in (c) and (d).

FIG. 7.

Variation versus radius ratio of (a) and (b) the CDHO coupling coefficients and (c) and (d) complex resonant frequencies, for two cylinders with radii a1 and a2, spacing $ℓ = 4 a 2$, and where $d = 2 a 2$ and $h = 4 a 2$, where single cylinder and identical cylinder references are given in (c) and (d).

Close modal
For an isolated cylinder, assume a DHO model of the form
(41)
where $ζ ( t )$ now denotes an approximation of the full linear solution. Motivated by the form of the SEM, let
(42)
where ω1 is the complex resonant frequency for the specified single cylinder and $ζ ̃$ 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
(43)
In the case of arrays (N > 1), the DHO (41) is applied to each body, along with coupling between the bodies, to give the CDHO model
(44)
where
(45)
and the real-valued coupling matrices d, $e ∈ ℝ N × 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
(46)
where the eigenmodes are set as the complex resonances by defining the matrices X, $Y ∈ ℂ N × N$ as
(47a)
(47b)
The vector $ζ ̃ ∈ ℂ N × 1$ depends on the initial conditions (see Sec. IV D).
Substituting Eq. (46) into (44) gives
(48)
where
(49)
Thus, Eq. (48) holds for arbitrary $ζ ̃$ if
(50)
(50a)
(50b)
and the coupling matrices, d and e, are found from (50a), as
(50c)
(50d)

The homotopy method may yield $Im { 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.

Consider the problem of two cylinders (N = 2) with the same draught but with different radii (ap 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 $λ j = ω − ω 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 / a 2$, where a1 is varied and a2 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 / a 2 < 0.75$), waves radiated by the smaller cylinder do not excite motion of the larger cylinder, such that d21 and $e 21 ≈ 0$, although the larger cylinder affects the smaller cylinder by reflecting waves, which is expressed through relatively large absolute values of d11 and e11. 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 $≈ e 1$ (not shown). Waves radiated by the larger cylinder excite motion of the smaller cylinder (d12 and e12 are non-zero) but the resulting waves created by the smaller cylinder do not influence the larger cylinder (d22 and $e 22 ≈ 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 / a 2 ∈ ( 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 / 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 $≈ 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.

For a given initial displacement vector for the full linear solution, Z [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
(51a)
(51b)
so that, similar to the SEM, the CDHO satisfies different initial conditions to the full linear solution.

The CDHO solutions for two identical cylinders with radius a and spacing $ℓ = 4 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 ¯ ≡ g / a 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 $ℓ = 4 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 ¯ < 15$ for the cylinders initially at rest.

FIG. 8.

Displacements ζp of cylinders (a) p = 1 and (b) p = 2, normalized by Z1, versus non-dimensional time, $t ¯ = g / a t$ ($g / a = 2$ Hz in computations), with draught $d = 2 a$, water depth $h = 4 a$ and spacing $ℓ = 4 a$, and with initial conditions $Z 1 = 1$ and $Z 2 = 0$, leading to $ζ 1 ( 0 ) / Z 1 = 1.030$ and $ζ 2 ( 0 ) / Z 1 = 0.001$ for the CDHO.

FIG. 8.

Displacements ζp of cylinders (a) p = 1 and (b) p = 2, normalized by Z1, versus non-dimensional time, $t ¯ = g / a t$ ($g / a = 2$ Hz in computations), with draught $d = 2 a$, water depth $h = 4 a$ and spacing $ℓ = 4 a$, and with initial conditions $Z 1 = 1$ and $Z 2 = 0$, leading to $ζ 1 ( 0 ) / Z 1 = 1.030$ and $ζ 2 ( 0 ) / Z 1 = 0.001$ for the CDHO.

Close modal
FIG. 9.

As in Fig. 8 but for cylinders (a) p = 1, (b) $p = 2 ( = 4 )$, and (c) p = 3, in an N = 4 square array, with initial displacements $Z 1 = 1$ and $Z 2 = Z 3 = Z 4 = 0$, where cylinders p = 1 and 3 are in opposite corners and $ζ 2 = ζ 4$ due to symmetry. The initial displacements of the full linear solution lead to $ζ 1 ( 0 ) / Z 1 = 0.975 , ζ 2 , 4 ( 0 ) / Z 1 = − 0.056$, and $ζ 3 ( 0 ) / Z 1 = − 0.073$ for the CDHO.

FIG. 9.

As in Fig. 8 but for cylinders (a) p = 1, (b) $p = 2 ( = 4 )$, and (c) p = 3, in an N = 4 square array, with initial displacements $Z 1 = 1$ and $Z 2 = Z 3 = Z 4 = 0$, where cylinders p = 1 and 3 are in opposite corners and $ζ 2 = ζ 4$ due to symmetry. The initial displacements of the full linear solution lead to $ζ 1 ( 0 ) / Z 1 = 0.975 , ζ 2 , 4 ( 0 ) / Z 1 = − 0.056$, and $ζ 3 ( 0 ) / Z 1 = − 0.073$ for the CDHO.

Close modal

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 $ℓ = 33.08 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 $ζ ̂ 1$, i.e., the same as ω1. It first appears for $ℓ ≈ 33.08 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.

FIG. 10.

(a) As in Fig. 9(a) but for spacing $ℓ = 33.08$ a. (b) The corresponding modulus of $det { Q }$, where the symbols (black bullet and red cross) denote complex resonant frequencies. (c) As in (b) but for spacing $ℓ = 32$ a.

FIG. 10.

(a) As in Fig. 9(a) but for spacing $ℓ = 33.08$ a. (b) The corresponding modulus of $det { Q }$, where the symbols (black bullet and red cross) denote complex resonant frequencies. (c) As in (b) but for spacing $ℓ = 32$ a.

Close modal

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.

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.

The authors have no conflicts to disclose.

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).

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

1.
K.
Budal
, “
Theory for absorption of wave power by a system of interacting bodies
,”
J. Ship Res.
21
,
248
253
(
1977
).
2.
J.
Falnes
, “
Radiation impedance matrix and optimum power absorption for interacting oscillators in surface waves
,”
Appl. Ocean Res.
2
,
75
80
(
1980
).
3.
D. V.
Evans
, “
Some analytic results for two and three dimensional wave-energy absorbers
,” in
Power from Sea Waves
, edited by
B. M.
Count
(
Academic
,
1980
), pp.
213
249
.
4.
H.
Chen
,
Q.
Xu
,
X.
Zheng
,
L. G.
Bennetts
,
B.
Xie
,
Z.
Lin
,
Z.
Lin
, and
Y.
Li
, “
Viscous effects on the added mass and damping forces during free heave decay of a floating cylinder with a hemispherical bottom
,”
Eur. J. Mech. B
98
,
8
20
(
2023
).
5.
N.
Sergiienko
,
B. S.
Cazzolato
,
B.
Ding
,
P.
Hardy
, and
M.
Arjomandi
, “
Performance comparison of the floating and fully submerged quasi-point absorber wave energy converters
,”
Renewable Energy
108
,
425
437
(
2017
).
6.
P.
McIver
and
R.
Porter
, “
The motion of a freely floating cylinder in the presence of a wall and the approximation of resonances
,”
J. Fluid Mech.
795
,
581
610
(
2016
).
7.
S. A.
Mavrakos
, “
Hydrodynamic coefficients for groups of interacting vertical axisymmetric bodies
,”
Ocean Eng.
18
,
485
515
(
1991
).
8.
M.
Göteman
, “
Wave energy parks with point-absorbers of different dimensions
,”
J. Fluid Struct.
74
,
142
157
(
2017
).
9.
C. M.
Linton
and
P.
McIver
,
Handbook of Mathematical Techniques for Wave/Structure Interactions
(
CRC Press
,
2001
).
10.
J.
Sun
,
S.-L. J.
Hu
, and
H.
Li
, “
Computing free motion of floating structures by pole-residue operations
,”
Appl. Ocean Res.
109
,
102558
(
2021
).
11.
C.
Stavropoulou
,
J.
Engström
, and
M.
Göteman
, “
Two-body, time domain model for a heaving point absorber
,” in
Advances in the Analysis and Design of Marine Structures
(
CRC Press
,
2023
) pp.
213
220
.
12.
Q.
Zhong
and
R. W.
Yeung
, “
Wave-body interactions among energy absorbers in a wave farm
,”
Appl. Energy
233-234
,
1051
1064
(
2019
).
13.
J. C.
McNatt
, “
A novel method for deriving the diffraction transfer matrix and its application to multi-body interactions in water waves
,”
Ocean Eng.
94
,
173
185
(
2015
).
14.
F.
Kara
, “
Time domain prediction of power absorption from ocean waves with wave energy converter arrays
,”
Renewable Energy
92
,
30
46
(
2016
).
15.
F.
Fàbregas Flavià
,
C.
McNatt
,
F.
Rongère
,
A.
Babarit
, and
A.
Clément
, “
A numerical tool for the frequency domain simulation of large arrays of identical floating bodies in waves
,”
Ocean Eng.
148
,
299
311
(
2018
).
16.
H.
Kagemoto
and
D. K. P.
Yue
, “
Interactions among multiple three-dimensional bodies in water waves: an exact algebraic method
,”
J. Fluid Mech.
166
,
189
(
1986
).
17.
P.
Siddorn
and
R.
Eatock Taylor
, “
Diffraction and independent radiation by an array of floating cylinders
,”
Ocean Eng.
35
,
1289
1303
(
2008
).
18.
A. N.
Williams
and
A. G.
Abul-Azm
, “
Hydrodynamic interactions in floating cylinder arrays-II. Wave radiation
,”
Ocean Eng.
16
,
217
263
(
1989
).
19.
S.
De Chowdhury
and
R.
Manasseh
, “
Behaviour of eigenmodes of an array of oscillating water column devices
,”
Wave Motion
74
,
56
72
(
2017
).
20.
A.
Ooi
,
A.
Nikolovska
, and
R.
Manasseh
, “
Analysis of time delay effects on a linear bubble chain system.
J. Acoust. Soc. Am.
124
,
815
826
(
2008
).
21.
C. J.
Fitzgerald
and
P.
McIver
, “
Approximation of near-resonant wave motion using a damped harmonic oscillator model
,”
Appl. Ocean Res.
31
,
171
178
(
2009
).
22.
F.
Ursell
, “
The decay of the free motion of a floating body
,”
J. Fluid Mech.
19
,
305
319
(
1964
).
23.
S. J.
Maskell
and
F.
Ursell
, “
The transient motion of a floating body
,”
J. Fluid Mech.
44
,
303
313
(
1964
).
24.
H. A.
Wolgamot
,
M. H.
Meylan
, and
C.
Reid
, “
Multiply heaving bodies in the time-domain: Symmetry and complex resonances
,”
J. Fluid Struct.
69
,
232
251
(
2017
).
25.
L. G.
Bennetts
and
M. H.
Meylan
, “
Complex resonant ice shelf vibrations
,”
SIAM J. Appl. Math.
81
,
1483
1502
(
2021
).
26.
R. W.
Yeung
, “
Added mass and damping of a floating vertical cylinder in finite-depth waters
,”
Appl. Ocean Res.
3
,
119
133
(
1981
).
27.
S. A.
Mavrakos
and
P.
Koumoutsakos
, “
Hydrodynamic interaction among vertical axisymmetric bodies restrained in waves
,”
Appl. Ocean Res.
9
,
128
140
(
1987
).
28.
M.
Peter
and
M.
Meylan
, “
Infinite-depth interaction theory for arbitrary floating bodies applied to wave forcing of ice floes
,”
J. Fluid Mech.
500
,
145
167
(
2004
).
29.
C.
Hazard
and
F.
Loret
, “
The singularity expansion method applied to the transient motions of a floating elastic plate
,”
ESAIM
41
,
925
943
(
2007
).
30.
M. H.
Meylan
and
M.
Tomic
, “
Complex resonances and the approximation of wave forcing for floating elastic bodies
,”
Appl. Ocean Res.
36
,
51
59
(
2012
).
31.
L. G.
Bennetts
,
M. A.
Peter
, and
F.
Montiel
, “
Localisation of Rayleigh–Bloch waves and damping of resonant loads on arrays of vertical cylinders
,”
J. Fluid Mech.
813
,
508
527
(
2017
).