2. The weakly nonlinear wave radiation problem

Wave–structure interactions are often simplified using linear potential flow theory where the input, in this case the forced displacement of the body $\zeta$ , is proportional to the output, the radiation force $F$ and the radiated waves $\eta$ generated. This formulation is ubiquitous and formally derived in most hydrodynamics textbooks such as Newman (Reference Newman2017) and Falnes & Kurniawan (Reference Falnes and Kurniawan2020). For a spherical body undergoing large oscillatory heave motions, representative of an operating point absorber WEC, a limitation is that effects such as the instantaneous change in submerged volume, $V(t) - V_0$ , are linearised such that $V(t) - V_0 \approx - S_0\,\zeta (t)$ , using the mean waterplane area $S_0$ . This prismatic approximation of the change in submerged volume is formally valid only for infinitesimal motion amplitudes, but in practice performs very well over a relatively large range. There is great interest in determining how well linear theory performs, and how large the most important corrections are in the wave frequency range.

Since this study focuses entirely on first-harmonic nonlinearities that are most relevant for power capture, we consider mild to moderate motion amplitudes where the wave–structure interaction remains weakly nonlinear. For such a system, any output such as the radiation force or radiated waves can be represented by a Stokes perturbation expansion

(2.1) \begin{equation} F(t) = F^{(11)}(t) + F^{(20)}(t) + F^{(22)}(t) + F^{(31)}(t) + F^{(33)}(t) + \mathcal{O}(A^4), \end{equation}

where the superscript $(pq)$ indicates a term of $p$ th order and $q$ th harmonic, i.e. a term that scales with the input amplitude $A$ to the $p$ th power, and appears around $q$ times the peak frequency $f_{\!p}$ in the frequency domain. Therefore, $F^{(11)}$ is simply the first-order first harmonic (linear term), and $F^{(31)}$ is the third-order first harmonic, a term that scales cubically with the amplitude but shares a similar frequency range with the linear term.

In this work, we use a broad-banded motion input with phases of all Fourier components deterministically aligned. Such a focused motion signal (analogous to a focused wave; see Tromans et al. Reference Tromans, Anaturk and Hagemeijer1991; Taylor & Williams Reference Taylor and Williams2004) is given by

(2.2) \begin{equation} \zeta _{\psi }(t) = \sum _i \mathcal{A}_i \cos (\omega _i t + \psi ), \end{equation}

where $\mathcal{A}_i$ is the Fourier amplitude corresponding to a component with angular frequency $\omega _i$ , and $\psi$ is the focus phase prescribed to all frequency components. The envelope peak of the displacement input is then $A = \sum _i \mathcal{A}_i$ . If we use two motion inputs $\pi$ phase apart to generate corresponding radiation forces $F_{\psi }$ and $F_{\psi +\pi }$ , and combine them in the form

(2.3) \begin{equation} F_{\textit{odd}}(t) = \frac {1}{2}\left (F_{\psi }(t) - F_{\psi +\pi }(t)\right ) = F^{(11)}(t) + F^{(31)}(t) + F^{(33)}(t) + \mathcal{O}(A^5), \end{equation}

then we can isolate the odd harmonics of the output radiation force (or alternatively the output radiated waves). Since the even harmonics $F^{(20)}$ and $F^{(22)}$ often overlap with the first-harmonic frequency range and contaminate the results, this method allows cleaner isolation of the first harmonic. Applying a simple low-pass filter in (2.3) is then sufficient to isolate the first-harmonic terms $F^{(11)} + F^{(31)}$ from the higher harmonics. Phase-based harmonic separation methods have been applied in many hydrodynamic studies (e.g. Fitzgerald et al. Reference Fitzgerald, Taylor, Eatock Taylor, Grice and Zang2014; Zhao et al. Reference Zhao, Wolgamot, Taylor and Eatock Taylor2017; Chen et al. Reference Chen, Zang, Taylor, Sun, Morgan, Grice, Orszaghova and Tello Ruiz2018; Vyzikas et al. Reference Vyzikas, Stagonas, Buldakov and Greaves2018; Hlophe et al. Reference Hlophe, Adcock, Ding, Zang, Dai, Tang and Taylor2025; Liu et al. Reference Liu, Tang, Mostert and Adcock2025). Further decomposition into $F^{(11)}$ and $F^{(31)}$ is, however, non-trivial for experimental data (Tan et al. Reference Tan, Orszaghova, Wolgamot, Kurniawan and Todalshaug2025), therefore the total first harmonic $F^{(1)} = F^{(11)} + F^{(31)}+\mathcal{O}(A^5)$ is analysed in the remainder of this paper.

3. Experimental method

The experiments were conducted at the Coastal and Offshore Research Laboratory, University of Western Australia. A schematic of the layout and key dimensions can be seen in figure 1(a). The set-up comprised a LinMot actuator (model no. P01-48x240F-SSC) to impose the heaving motion of the sphere, an Interface 3-axis 50 N load cell (model no. 3A120) to measure the force on the sphere, and an array of Edinburgh Designs wave gauges (model no. WG8USB) to measure the radiated waves (gauge locations detailed in figure 1(c)). While the imposed motion of the sphere was recorded by the actuator, we also implemented several external displacement measurement instruments: a Penny & Giles linear displacement transducer (model no. SLS130), a LinMot magnetic strip linear encoder (model no. MS01-1/D-SSI) and a Qualisys motion capture system (model no. 5+) to ensure that the input motions were known accurately – all displacement measurements were found to be consistent to within $0.1\,\%$ . An annotated photo of these instruments is shown in figure 1(b). Most instruments were sampled at 128 Hz or higher, and all of them were synchronised through an electronic trigger. The aluminium sphere model had radius $R = 0.126 \pm 0.001$ m. The water depth was 1.1 m.

Figure 1.

(a) Schematic of the experimental set-up. (b) Corresponding photo with annotated instrumentation: $\unicode {x2460}$ linear displacement transducer, $\unicode {x2461}$ magnetic strip linear encoder, $\unicode {x2462}$ linear actuator, $\unicode {x2463}$ wave gauge 3, $\unicode {x2464}$ Qualisys marker, $\unicode {x2465}$ load cell, $\unicode {x2466}$ model sphere. (c) The wave gauge locations are relative to the sphere centre, and the coloured numbers indicate the gauges referenced in figure 5. Note that wave gauge 4 is not used in these radiation tests.

Figure 1. Long description

Panel A: A schematic diagram of the experimental setup. The diagram shows a wave flume with a sphere and actuator positioned at the center. Wave gauges are placed at various locations along the length of the flume. A stiff mounting frame supports the sphere and actuator. The still water level and a porous beach are also indicated. The wave gauges are numbered from 1 to 12, with specific positions relative to the sphere center. Panel B: A photo of the experimental setup with annotated instrumentation. The annotations include a linear displacement transducer, magnetic strip linear encoder, linear actuator, wave gauge 3, Qualisys marker, load cell, and model sphere. The photo shows the physical arrangement of these components in the wave flume. Panel C: A table listing the wave gauge locations relative to the sphere center. The table includes wave gauge numbers, their x and y coordinates, and indicates which gauges are referenced in figure 5. Wave gauge 4 is noted as not being used in these radiation tests.

The demand motion signals were compact NewWave-type signals (Tromans et al. Reference Tromans, Anaturk and Hagemeijer1991; Taylor & Williams Reference Taylor and Williams2004) with underlying spectra of JONSWAP type (with peak enhancement factor $\gamma = 2.5$ ) defined across frequencies from 0 to $2.5f_{\!p}$ . Three peak frequencies $f_{\!p} = 0.4, 0.6, 0.8$ Hz, six amplitudes $A = 0.02, 0.03, 0.04, 0.05, 0.06, 0.07$ m (corresponding to the envelope peak of the motion signal), and two focus phases $[\psi ,\psi +\pi ]=[{3\pi }/{2}, {\pi }/{2}]$ (up- and down-crossing focused displacement signals) were tested. However, for brevity, only plots of the experimental data for $f_{\!p} = 0.8$ Hz are presented in the following sections, and the data from the other two peak frequencies are presented in Appendix A.

We note that while it would have been ideal to position the sphere equidistantly from both ends of the flume to maximise the available time prior to reflections, these tests are part of a larger campaign that included tests with incident waves, so a compromise was made to position the sphere further from the beach (and thus closer to the wavemaker, which was not used in these tests). This short distance to the wavemaker inevitably affected the way measurements were analysed, as further discussed in § 7.

Here, we discuss the extraction of the radiation force from the total force measured by the load cell. By nature of the testing set-up, the total force on the load cell $F_{\textit{LC}}$ is dominated by the inertial $F_{\textit{IN}}$ , the hydrostatic $F_{\textit{HS}}$ and weight $F_{G}$ forces on the sphere, whereas the desired hydrodynamic radiation force makes up only $\sim15\,\%$ of the measured force at most, and much less in many cases. The isolated radiation force on the sphere is thus given by

(3.1) \begin{equation} F(t) = F_{\textit{LC}}(t) - F_{\textit{HS}}(t) - F_{\textit{IN}}(t) - F_G. \end{equation}

This leads to two main challenges: first, the estimation of the inertial and hydrostatic forces has to be very accurate to ensure that their contributions are properly removed, and second, the small absolute magnitude of the radiation force means that the signal-to-noise quality is much lower than e.g. a diffraction test where the entirety of the measured dynamic force is used (Tan et al. Reference Tan, Orszaghova, Wolgamot, Kurniawan and Todalshaug2025). Dry tests of every condition were thus conducted to measure the inertial force on the load cell while the set-up was actuated in air. Theoretically, one could also conduct dedicated tests to directly measure the hydrostatic force by slowly changing the draft of the sphere in still water. Unfortunately, the slow drift of the load cell over the long times necessary for such a test was found to substantially affect these measurements, so the hydrostatic force had to instead be indirectly estimated from $F_{\textit{HS}}(t) = \rho gV(t)$ . The volume beneath the mean free surface was calculated using the body displacements and volume profile, with the latter measured by a HandySCAN three-dimensional laser scanner (model no. BLACK Elite) with accuracy $0.025$ mm. The radiated wave measurements, on the other hand, do not require this post-processing, although they are of the order of millimetres and thus are also affected by low signal-to-noise ratio. Repeated tests conducted across three different days were found to have root mean square error 0.1 mm in the wave measurements, and 0.1 N in the force measurements. No wave breaking was observed.

4. Experimental results

The input displacement and the extracted first-harmonic radiation forces and radiated wave time series are plotted in figure 2. All six amplitudes of the $f_{\!p} = 0.8$ Hz condition are plotted, with darker lines representing smaller amplitude motions, and progressively lighter colours denoting the larger amplitudes. The black lines (figures 2 a,b) show the imposed heave displacement $\zeta (t)$ , the blue lines (figures 2 c,d) show the vector projection of the radiation force onto the displacement $ (F^{(11)}(t) + F^{(31)}(t) )_{\parallel \zeta (t)}$ (‘the added mass portion of the radiation force’), the orange lines (figures 2 e, f) show the corresponding vector rejection $ (F^{(11)}(t) + F^{(31)}(t) )_{\perp \zeta (t)}$ (‘the radiation damping portion’), and the green lines (figures 2 g,h) show the radiated wave elevation from a wave gauge close to the sphere (wave gauge 7 in figure 1). Figures 2(a,c,e,g) are the raw signals, while figures 2(b,d, f,h) contain the same signals normalised by the displacement envelope peak $A$ (the envelope is partially shown by the dash-dotted line in figure 2 b).

Figure 2.

Experimental (a,c,e,g) raw and (b,d, f,h) normalised first-harmonic (a,b) body displacement, (c,d) vector projection and (e, f) rejection of the radiation force onto the heave displacement (analogous to the added mass and radiation damping coefficients), and (g,h) radiated free surface. Six input motion amplitudes are shown, where the darker lines represent the smaller input motion amplitudes. The motion envelope is partially shown by the dash-dotted line in $(b)$ . Peak frequency of the motion input is $f_{\!p} = 0.8$ Hz.

While the normalised motion time series collapse to high precision (by design), the normalised responses (force and free surface) can be seen to exhibit weakly nonlinear behaviour. The inset plots show that there is an ordered separation of the forces similar to that observed in Tan et al. (Reference Tan, Orszaghova, Wolgamot, Kurniawan and Todalshaug2025), where it was confirmed to be the result of third-order first-harmonic terms $F^{(31)}$ . At the largest displacement amplitude tested ( $A/R = 0.56$ ), the normalised forces appear to be $\sim10\,\%$ lower than at the smallest amplitude $A/R = 0.16$ – i.e. the first-harmonic radiation force on a sphere undergoing heaving motions of half of its radius can be 10 % lower than predicted by linear theory. We emphasise again that the exact (to experimental accuracy) hydrostatic contribution relative to still water has been removed here. The level of separation between linear and nonlinear values appears to be comparable between the ‘added mass’ and ‘radiation damping’ portions of the radiation force. The radiated waves also clearly show the nonlinear ordering with motion amplitude.

5. The linear radiation problem in a channel

In subsequent sections, results from linear theory are compared to experimental data. It is therefore useful to briefly discuss features of the linear radiation problem that arise as a consequence of testing between parallel walls. For an oscillating body positioned centrally in a channel of total width $L$ , the far-field radiated waves at a finite frequency consist of a discrete set of propagating wave components travelling along ( $n = 0$ ) and at an angle ( $n \geq 1$ ) to the channel, with their profile across the channel proportional to $\cos {(2n\pi y/L)}$ and visualised in figure 3(b) (note that antisymmetric modes are not excited in linear theory due to the symmetry of the problem). For $n \geq 1$ , each mode can be viewed as a pair of obliquely propagating components with propagation angles $\theta _n$ and $-\theta _n$ .

Figure 3.

(a) Ratio of the linear radiated wave amplitude and displacement amplitude of a sphere in a channel, decomposed into contributions from each channel mode. (b) Their respective profiles visualised across the half-width $(L/2)$ of the channel. The sphere radius is $R = L/12$ .

Figure 3. Long description

Two line graphs depict the hydrodynamic interaction between surface gravity waves and floating structures. Panel A: A line graph shows the ratio of the linear radiated wave amplitude to the displacement of a sphere in a channel, decomposed into contributions from each channel mode. The x-axis is labeled kL, and the y-axis is labeled |h(11)/ξ| (m^-1). The graph includes four lines representing different modes (n = 0, n = 1, n = 2, n = 3), each with distinct colors. The lines show various peaks and troughs, indicating the interaction strength at different kL values. Panel B: Another line graph visualizes the mode profiles across the half-width of the channel. The x-axis is labeled y/L, and the y-axis is labeled Mode profile. The graph includes four lines representing different modes (n = 0, n = 1, n = 2, n = 3), each with distinct colors. The lines show sinusoidal patterns, indicating the distribution of the mode profiles across the channel width.

Solutions to the linear radiation problem for a sphere in a channel were first presented by Kashiwagi (Reference Kashiwagi1989). In this work, we use Hydrostar (2016) to solve the heave radiation problem. To elucidate features of the problem, we compute the far-field linear radiated waves emitted by a heaving half-immersed sphere positioned centrally in a channel. The result is shown as a ratio of the radiated wave amplitude to the body displacement in figure 3(a) (the radiated waves were calculated at $x = 106R$ from the sphere centre). The thin vertical lines indicate the frequencies where an integer number of wavelengths exactly spans the width of the channel – the cut-off frequencies above which additional channel modes can propagate. Each colour shows the contribution from a given mode, found by decomposing the total into modes present at each frequency. This was calculated utilising the complex amplitudes of the total free-surface elevation from a handful of cross-flume locations, and at each frequency, fitting the relevant mode profiles across the flume width at these locations. The modulus of each fitted modal complex amplitude is what is displayed in figure 3(a).

After each successive cut-off frequency, another propagating mode is permitted across the channel, and the total radiated wave measurement at a given location in the flume is therefore subject to the interference of these modes with each other. This causes highly irregular features in the frequency-domain transfer functions at any point in space when $kL \gt 2\pi$ since contributions from multiple modes are totalled (even though the individual modes follow smooth curves). Some of these irregular features are seen in the linear results in figure 5(d) below.

The far-field radiated waves are related to the damping, so the decomposed free-surface amplitudes can also be used in the following equation, adapted from Srokosz (Reference Srokosz1980) and Linton & Evans (Reference Linton and Evans1993):

(5.1) \begin{equation} \begin{gathered} b(\omega ) = \frac {2\rho g L c_g(\omega )}{\omega ^2 } \sum _{n}\frac {\cos (\theta _n(\omega ))}{\varepsilon _n}\left |\frac {\hat {\eta }^{(11)}_n(\omega )}{\hat {\zeta }(\omega )}\right |^2 \!, \\ \varepsilon _n = \begin{cases} 1, & n = 0, \\ 2, & n \geq 1, \end{cases}\quad \theta _n(\omega ) = \arcsin \left (\frac {2n\pi }{k(\omega )\,L}\right )\!, \end{gathered} \end{equation}

where $b$ is the radiation damping, $\omega$ is the angular frequency, $\rho$ is the water density, $g$ is the gravitational acceleration, $c_g$ is the group velocity, $n$ is the symmetric mode number, $k$ is the wavenumber, $\hat {\eta }^{(11)}_n$ is the complex amplitude of the linear radiated free surface in the $n$ th channel mode propagating at angles $\pm \theta _n$ , and $\hat {\zeta }$ is the complex amplitude of the body displacement. Thus the radiation damping can be computed in modal components and compared to the total (directly from Hydrostar) and checked for consistency, thus validating the modal decomposition.

6. The volume flow radiation model

Work on fluids is generally governed by volume flow and pressure. In an attempt to approximately account for the nonlinearities observed in the experimental data, we consider the volume flow in and out of the fluid due to the relative motion between body and fluid surface rather than the body heave motion itself. It is our hope that such a change of variables, partly inspired by small-body theory, may be used to extend the applicability of the linear potential flow solution for the problem at hand.

Accordingly, a partially nonlinear heave radiation model based on volume flow is now derived to model the observed nonlinear changes. We start from the time-domain formulation for the linear radiation potential (e.g. McIver & McIver Reference McIver and McIver2011):

(6.1) \begin{equation} \phi ^{(11)}(\boldsymbol{x},t) = \varOmega (\boldsymbol{x})\, \ddot {\zeta }(t) + \int _0^t K(\boldsymbol{x},t-\tau )\, \dot {\zeta }(\tau ) \,\mathrm{d}\tau , \end{equation}

where $\varOmega$ is the infinite-frequency radiation potential, $K$ is the heave radiation impulse response function for the potential, and the accents $\dot {\Box }$ and $\ddot {\Box }$ represent the conventional time derivatives. Heaving bodies are associated with oscillatory volume flow and radiated waves. In the classical literature for a heaving sphere (in open water, Havelock Reference Havelock1955; Hulme Reference Hulme1982) this link is clear, in that the radiated field far from the body is directly proportional to the strength of the wave source term used in constructing the solution. We propose to scale the radiated potential using the ratio of the instantaneous nonlinear volume flux associated with the sphere’s immersed volume to the version from linear theory. Taking (6.1) and multiplying and dividing by the equilibrium waterplane area $S_0$ gives

(6.2) \begin{equation} \phi ^{(11)}(\boldsymbol{x},t) = \frac {1}{S_0} \varOmega (\boldsymbol{x})\, [S_0\, \ddot {\zeta }(t)] + \frac {1}{S_0} \int _0^t K(\boldsymbol{x},t-\tau )\, [S_0\, \dot {\zeta }(\tau )] \,\mathrm{d}\tau , \end{equation}

and ‘correcting’ the products $S_0\, \ddot {\zeta }(t)$ and $S_0\, \dot {\zeta }(\tau )$ to the time derivatives of the exact nonlinear submerged volume (as discussed in § 2), we obtain the partially nonlinear model

(6.3) \begin{equation} \phi (\boldsymbol{x},t) = -\frac {1}{S_0} \bigg ( \varOmega (\boldsymbol{x})\, \ddot {V}(t) + \int _0^t K(\boldsymbol{x},t-\tau )\, \dot {V}(\tau ) \,\mathrm{d}\tau \bigg ). \end{equation}

The free-surface elevation can be computed using the linear free-surface boundary condition

(6.4) \begin{equation} {\eta (t)} = -\frac {1}{g}\frac {\partial \phi }{\partial t}\bigg |_{z=0}, \end{equation}

in which case the first term of (6.3) does not enter. This volume flow model is consistent with the local nonlinearity at the body being of greater importance than the nonlinearity associated with the free-surface boundary condition, similar to body-exact models noted earlier. In contrast to those models, however, this method does not enforce satisfaction of the body boundary condition. Due to the complexity of the far-field solution in a channel (discussed in § 5), construction of the kernel $K$ is unappealing in this case (as also noted by Cohen & Troesch (Reference Cohen and Troesch1988) for the case of a sphere in a channel). Using the Fourier transform of the nonlinear volume flux time series, an equivalent calculation can be carried out in the frequency domain using the free-surface transfer functions shown in figure 3. The volume flow model is physically intuitive, with smaller radiated waves generated when the volume flux is smaller than its prismatic approximation, and vice versa.

For a heaving sphere with an initial draft of one radius (i.e. nominally half-immersed), the submerged volume (beneath the still free surface) is

(6.5) \begin{equation} V(t) = \frac {1}{3}\pi \Big (2R^3 - 3R^2\,\zeta (t) + \zeta ^3(t)\Big ), \end{equation}

where $R$ is the radius of the sphere. We note that this expression holds for any heave displacement time series $\zeta (t)$ , including the specific ‘focused’ signals $\zeta _\psi$ used in this work. With such an analytical expression in hand, there is essentially no extra computational cost in using the volume flow model compared to the linear. Importantly, since (6.5) only contains first- and third-order terms, the proposed nonlinear model (with contributions proportional to $\dot {V}$ ) yields only these terms. No even-ordered terms are produced due to the even geometry about the waterline. In fact, one can further refine the model to only simulate the first-harmonic terms (i.e. removing the third-harmonic $\eta ^{(33)}(t)$ terms) by making the substitution $\zeta ^3(t) \rightarrow ({3}/{4})\,\zeta (t)\, [\zeta ^2(t) + H(\zeta (t))^2 ]$ in (6.5), where $H$ is the Hilbert transform operator (as per Walker et al. Reference Walker, Taylor and Eatock Taylor2004). This allows the model to be directly compared with the experimentally extracted first-harmonic data. Another useful observation is the opposing signs of the first- and third-order terms, which would reduce the amplitude scaling of the first harmonic to be slightly less than linear, consistent with the results shown. Due to manufacturing limitations, our spherical model is very slightly elongated vertically ( $R_z/R_x = 1.008$ ), but a spheroidal volume also shares the same equations (down to some constants), thus all of the traits discussed above are valid.

Figure 4.

Experimental (a–c) ‘added mass’ and (d–f) ‘radiation damping’ coefficients for six input motion amplitudes, where the darker lines represent the smaller input motion amplitudes. Thin solid lines are the experimental data calculated from (7.1). Dotted lines are results from linear potential flow theory calculated from Hydrostar (2016). Bold solid lines are also from linear potential flow theory, truncated in the same way as the experimental data. Thin black vertical lines are cut-off frequencies. The smaller plots are used to better show regions of interest. Peak frequency of the motion input is $f_{\!p} = 0.8$ Hz.

7. Comparison and discussion

Experimental ‘added mass’ and ‘radiation damping’ coefficients are now derived in the conventional manner using the total first-harmonic force

(7.1) \begin{equation} \tilde {m}(\omega ) = \frac {1}{\omega ^2}\, \textrm{Re} \left(\frac {\hat {F}^{(1)}(\omega )}{\hat {\zeta }(\omega )} \right)\!, \quad \tilde {b}(\omega ) = -\frac {1}{\omega }\, \textrm {Im} \left(\frac {\hat {F}^{(1)}(\omega )}{\hat {\zeta }(\omega )}\right)\!, \end{equation}

where the $\tilde {\Box }$ accent is used to differentiate them from the theoretical added mass $m$ and radiation damping $b$ calculated from linear theory. This distinction is important, because the experimental ratios presented cannot be used in the same way as the linear coefficients: they are functions of the spectral shape and amplitude of the input signal. The experimental values are plotted in figure 4 with the same colouring conventions as in figure 2. Additionally, we also plot the linear result in dotted and bold solid lines, where the dotted lines are the direct frequency-domain calculations from Hydrostar (2016) (thus assuming a steady-state system), and the bold solid lines are processed in the same way as the experimental data, i.e. a time series truncated to the same duration (prior to experimental reflections), and then transformed to the frequency domain via a Fourier transform. The thin vertical black lines are the cut-off frequencies described in § 5, where the effects of channel modes are discussed. The smaller plots provide magnified views for regions of interest.

Naturally, the progressive reductions in the amplitude-normalised responses observed in the time series in figure 2 are also apparent in the frequency domain, where the ‘transfer functions’ appear to gradually reduce with increasing input motion amplitudes. Additionally, the smallest amplitude lines appear to be consistent with the theoretical linear results calculated from Hydrostar, as expected. The reductions in both added mass and damping are in contrast to the body-exact computations of Lin & Yue (Reference Lin and Yue1991) where the ‘added mass’ was found to decrease considerably, while the ‘damping’ remained close to constant. However, those results were only for a single monochromatic motion input (at a frequency of approximately $1.75f_{\!p}$ here).

Figure 5 shows the ratio of the first-harmonic amplitudes of the radiated free surface and body displacement. In addition to the experimental and linear potential flow theory curves, results from the volume flow model are displayed. The model appears to effectively capture the nonlinearity in the radiated waves, both in the near field (figures 5(a–c), green colour) and in the far field (figures 5(d–f), purple colour), corresponding to wave gauges 7 and 8, respectively, in figure 1. These observations largely hold true for the other peak frequencies shown in Appendix A.

Before discussing these results in more detail, we note some discrepancies between the experimental and numerical radiated free-surface amplitudes in the form of oscillations seen in the low-frequency tails of the experimental data in figure 5. These are due to reflections of the radiated wave from either the wavemaker or the beach, and therefore do not appear in the numerical data, which are simulated assuming a channel of infinite length. The time series used to generate these plots are thus truncated early to minimise the impact of these reflections, but ultimately a compromise is required since we also try to retain as much of the true behaviour as possible. This is further exacerbated in the lower peak frequencies (shown in Appendix A) where the waves travel even faster and more reflections inevitably contaminate the results, incurring larger spurious oscillations (especially in the far-field wave gauge 8).

Figure 5.

Ratio of the first-harmonic amplitudes of the radiated free surface and body displacement for wave gauges (a–c) 7 and (d–f) 8 for six input motion amplitudes, where the darker lines represent the smaller input motion amplitudes. Thin solid lines are the experimental data. Thin dashed lines are derived from the nonlinear numerical volume flow model derived in § 6. Dotted lines are results from linear potential flow theory calculated from Hydrostar (2016). Bold solid lines are also from linear potential flow theory, truncated the same way as in the experimental data. Thin black vertical lines are cut-off frequencies. The smaller plots are used to better show regions of interest. Peak frequency of the motion input is $f_{\!p} = 0.8$ Hz.

One of the notable features of the observed nonlinear behaviour is that the amplitude separations are more pronounced at the high-frequency tails, which was also observed and discussed in Tan et al. (Reference Tan, Orszaghova, Wolgamot, Kurniawan and Todalshaug2025). The nonlinear terms present, in addition to their different amplitude scaling properties, also possess broader bandwidths than the linear term, a property known as ‘spectral smearing’ (Fitzgerald et al. Reference Fitzgerald, Taylor, Eatock Taylor, Grice and Zang2014), caused by the sum–frequency and difference–frequency interactions. For example, given a linear term with frequency content that spans $0.6f_{\!p}$ to $2.5f_{\!p}$ , the third-order response from the $(++-)$ frequency interaction will span from 0 to $4.4f_{\!p}$ . As a result of this difference in bandwidths, the first-harmonic $F^{(11)} + F^{(31)}$ is usually dominated by the linear term $F^{(11)}$ around the peak frequency where most of its energy is concentrated. However, the linear term $F^{(11)}$ subsides more rapidly than the third-order term $F^{(31)}$ towards the tails of the spectrum, therefore the first-harmonic lines deviate more from the linear prediction as the nonlinear contribution grows larger. This effect is, by nature, limited to broad-banded inputs, thus the influence of these terms is much more difficult to observe in monochromatic inputs (such as those used by Vantorre Reference Vantorre1986b ). The free surface (experimental and numerical) and the ‘radiation damping’ results exhibit this spectral smearing property, while the experimental ‘added mass’ appears to be the sole exception. Since radiation damping and the far-field radiated free surface are directly linked, it is encouraging to see analogous nonlinear behaviour in these two quantities.

Added mass, on the other hand, is dominated by local effects and is more sensitive to experimental uncertainties. One of these uncertainties is the removal of the hydrostatic force discussed in § 3. While every care was taken to estimate the volume of the model as accurately as possible, small local defects due to the manufacturing method made this very challenging, even with the three-dimensional scan. It is apparent that the experimental ‘added mass’ estimates exhibit different trends in amplitude change as a function of frequency compared to the other results. Unfortunately, the Kramers–Kronig relations assume a linear system, so cannot be used on the nonlinear quantities studied here.

Since the observed nonlinear changes in the radiated wave data are largely consistent with the volume flow model (which only contains a cubic nonlinear term), this implies that these effects are caused by the third-order first-harmonic $F^{(31)}$ and $\eta ^{(31)}$ terms, as initially proposed. While there is a method to estimate the order of the nonlinear terms from the amplitude scaling of the experimental data, as shown in Tan et al. (Reference Tan, Orszaghova, Wolgamot, Kurniawan and Todalshaug2025), we do not use it here since the smaller amplitude cases are somewhat noisy.

While we have considered that viscous effects potentially contribute to the observed nonlinearities, due to the small size of the spherical model and low Keulegan–Carpenter number ( $\textit{KC} \lt 2.1$ ) it is likely that the boundary layer is laminar and does not separate. This is known to linearly affect the added mass and radiation damping (Stokes Reference Stokes1851; Molin & Etienne Reference Molin and Etienne2000), and thus would manifest in the frequency-domain plots (figures 4 and 5) as a single line offset from the linear potential flow solution, as opposed to the systematic separation of lines observed. We also estimate the magnitude of the laminar drag force later, in § 8, and find that this viscous force is approximately 0.5 % of the measured force at most.

The lower peak frequencies in Appendix A show similar nonlinear behaviours, where the ordered reductions in the ‘transfer functions’ can be observed clearly for the $0.6$ Hz results. While the $0.4$ Hz results (with correspondingly smaller hydrodynamic forces) are more corrupted by noise, they also show a general decreasing trend with increasing motion amplitude. Interestingly, the volume flow model appears effective in modelling the radiation behaviour at all peak frequencies, suggesting that the volume corrections are the dominant source of nonlinearity across the tested frequency range. Therefore, while there are different features in the transfer functions across the peak frequencies (largely influenced by the channel modes discussed in § 5), the nonlinearities appear to be largely independent of the peak frequency, depending only on the instantaneous submerged volume.

We now briefly mention some wider implications of these results. Since the volume flow model is an extremely simple model to implement, only requiring the linear hydrodynamic coefficients and basic geometric parameters (the equilibrium waterplane area and instantaneous submerged volume), it can be evaluated almost instantly and could ultimately be used in real-time control of WECs. In fact, if we inspect the linear equation (6.1) and the volume flow nonlinear equation (6.3), then there is practically no extra computational cost if the submerged volume can be directly obtained from the draft.

We note that spherical or near-spherical geometries, as studied here, excite these third-order wave-structure interactions particularly strongly given the cubic relationship between their draft and submerged volume. As such geometries are common in point absorber designs, the effects are worthy of study.

8. Near- and far-field radiated power

To further test the consistency of the measurements, we also compare the input power in the near field (transferred from the moving body to the fluid), and the radiated power in the far field (transported by the radiated waves). Since the nonlinearities appear most strongly in the high-frequency tail, we focus on the frequency band from $1.4f_{\!p}$ to $1.7f_{\!p}$ between the first and second cut-off frequencies (avoiding frequencies very close to the cut-offs due to the extremely slow propagation of the nearly trapped waves).

Two factors prevent simply using the measured radiated free-surface measurements directly. First, the inline propagating mode $(n = 0)$ and the obliquely propagating channel mode $(n = 1)$ must be explicitly separated since they have different relationships between modal amplitudes and power, which need to be accounted for when calculating the transported power. Second, dissipation of energy during propagation from the body to the far-field gauge should be accounted for. To separate the total measured wave elevation into modes, we therefore fit to the following basic equation, which exploits the inline propagation speed of each mode to differentiate their contributions:

(8.1) \begin{equation} \hat {\eta }^{(1)}(\omega , x) = \sum _{n=0,1} \hat {\eta }^{(1)}_n(\omega , x_0) \exp {[\left ( -\mathrm{i}\,k(\omega ) \cos {(\theta _n(\omega ))} - \alpha _n(\omega ) \right ) (x - x_0)]}, \end{equation}

where $x_0$ is a reference point, and $\alpha _n$ is the (complex) dissipation coefficient of the $n$ th mode accounting for laminar boundary layer friction on the flume walls and free surface (discussed later and derived in Appendix B). This fitting procedure therefore takes the measured wave elevation from several inline wave gauges and, at each frequency, finds the best-fit complex amplitudes $\hat {\eta }^{(1)}_0$ and $\hat {\eta }^{(1)}_1$ . This approach (including only linear dispersion with laminar dissipation) is appropriate if $k\delta \gg (ka)^2$ , where $\delta = (\nu /\omega )^{{1}/{2}}$ , and $a$ is a representative wave amplitude (Mei & Liu Reference Mei and Liu1973). Here, the largest total wave elevation measured anywhere in the flume (for the largest motion amplitude) was $9$ mm (figure 2) – yielding $(k_{\!p} \eta _{\textit{max}}^{(1)})^2 = 6 \times 10^{-4}$ , while $k_{\!p}\delta = 1 \times 10^{-3}$ . Given that this is the largest free-surface elevation measured (close to the body) for the largest body motion, and that the groups disperse rapidly with distance, this is a reasonable assumption. We note in passing that if transverse measurements of the far-field free-surface were available at a single distance, then one could fit the $n = 1$ mode profile since the $n = 0$ mode is flat across the channel (as done in § 5 for the theoretical linear results). However, the wave gauge layout in our tests did not contain transversely positioned gauges in the far field (see figure 1).

Figure 6 shows the measured free-surface time series from the four far-field wave gauges used for the fitting, which clearly demonstrates the different propagation speeds of the two modes, and their progressive separation from wave gauges 8–11. The shaded regions show the range of propagation speeds accounting for frequency dispersion and the dispersion due to the different propagation angles in the $n = 1$ group. Wave gauge 12 is not used for this analysis since it is very close to the beach, and the uncompromised time window is too short. After truncating the reflections and using a conventional least squares fitting method, the separated wave elevations at each gauge were obtained, and plotted for wave gauge 8 in figure 7 for all six input motion amplitudes, alongside a linear solution that was processed identically (i.e. accounting for the laminar dissipation over the distance from the sphere to wave gauge 8, and using the same truncations).

Figure 6.

Extended experimental first-harmonic wave elevation time series for $A = 0.07$ m and $f_{\!p} = 0.8$ Hz across wave gauges (WG) 8–11 (see figure 1), band-pass filtered from $f_l/f_{\!p} = 1.4$ to $f_u/f_{\!p}=1.7$ . Vertical axis offsets for each gauge are proportional to their respective $x$ locations in the flume. The shaded regions show the dispersion of the radiated wave group accounting for both frequency and propagation angle of all components. For completeness, waves radiated from the sphere in the $-x$ direction, and subsequently reflected from the wavemaker (WM), are also shown.

Figure 6. Long description

A line graph showing wave elevation time series across different wave gauges, with shaded regions indicating the dispersion of radiated wave groups. The graph includes multiple lines representing different wave gauges (WG 8, WG 9, WG 10, WG 11) and their respective offsets. The horizontal axis represents time normalized by the peak period (t/Tp) ranging from 0 to 100. The vertical axis represents wave elevation (eta) in meters, with offsets for each gauge proportional to their locations in the flume. The shaded regions highlight the dispersion of the radiated wave group, accounting for frequency and propagation angle. Additionally, waves radiated from the sphere in a specific direction and reflected from the wavemaker are also shown. Each line exhibits oscillatory behavior with varying amplitudes and phases, indicating the wave elevation patterns at different gauges. The trends show how the wave elevations change over time and how they are influenced by the radiated wave groups and reflections.

Figure 7.

Filtered (1.4 $f_{\!p}$ to 1.7 $f_{\!p}$ ) and amplitude-normalised first-harmonic wave elevation time series for $f_{\!p} = 0.8$ Hz from wave gauge 8 decomposed into the $n = 0$ and $n = 1$ modes through fitting the experimental data to (8.1). Six amplitudes are shown, where the darker lines represent the smaller input motion amplitudes. The thin lines are the experimental data, and the bold line is the corresponding linear solution for the radiated field from the sphere, with linear dissipation applied. Inset plots show magnified views of regions of interest.

The fitted $n = 0$ and $n = 1$ modes appear to be consistent with the linear solution, suggesting that the separation method is effective. The fitting was substantially improved by including the dissipation coefficient from side walls only, which reduced the root mean square error between the fitted and measured total amplitude by approximately $40\,\%$ . (Friction from the bottom of the flume is negligible as the waves in this frequency band are on deep water.) However, the observed dissipation rate was still noticeably higher, so we also included dissipation at the free surface, which is affected by surface contamination and ‘ageing’ effects commonly experienced in laboratory tank tests (Van Dorn Reference Van Dorn1966; Mei Reference Mei1989). This free-surface dissipation coefficient was scaled by a single factor $\beta$ , common to all frequencies, between 0 (no contamination) and 1 (inextensible film present) to best fit the measurements, as noted in Appendix B. The final value used was $\beta =0.31$ , which reduced the error by a further $55\,\%$ . With the amplitudes for each mode now separated, we can perform the power calculation using the far-field radiated waves, and compare it to the results in the near field.

The nonlinear input mean power in the near field $\overline {P}_{\!F}$ follows from the product of the experimental force and velocity summed over the frequencies present in our band of interest:

(8.2) \begin{equation} \overline {P}_{\!F} = \frac {1}{2} \sum _i \textrm{Re}\left (\mathrm{i} \omega _i\, \hat {\zeta }(\omega _i)\,\hat {F}^{*(1)}(\omega _i)\right )\!, \end{equation}

where $\Box ^*$ denotes the complex conjugate. Consistent with (5.1), for the radiated waves we obtain

(8.3) \begin{equation} \overline {P}_\eta = \rho g L \sum _i c_g(\omega _i) \sum _{n} \frac {\cos {(\theta _n(\omega _i))}}{\varepsilon _n}\left |\hat {\eta }_n^{(1)}(\omega _i,0)\right |^2 \!, \end{equation}

where $\overline {P}_\eta$ represents the mean power calculated from the (far-field) wave elevation measurements. Note that in the above, the mode amplitudes $\hat {\eta }_n^{(1)}$ are taken at the sphere location (i.e. at $x=0$ , using (8.1)) such that the dissipation does not contaminate the far-field power estimate. We also note that in both of the above expressions for mean power estimates, only products of the same frequency terms are needed, as any other terms would be oscillatory in nature and therefore not contribute to the time-averaged value.

Figure 8.

Experimental (circle and square markers), volume flow numerical (dashed black line) and linear (solid black line) mean power from the near-field (orange) and far-field (purple) calculations over the frequency band from 1.4 $f_{\!p}$ to 1.7 $f_{\!p}$ . Linear viscous effects have been removed from the experimental data as described in the text. The vertical axis is normalised by the power from linear theory.

We also attempted to remove the laminar drag force on the body from the experimental result to achieve a more accurate comparison to the numerical potential flow results. If we assume that the boundary layer thickness is much smaller than the body radius, and thus the body surface can be locally treated as a flat plate, then we can use the classical solution from Stokes’ oscillating plate to estimate the shear stress. This is calculated for 65 discretised area elements over the hemisphere, using the local tangential velocities calculated from Hydrostar, and then integrated to obtain the total laminar drag. Consistent with Stokes (Reference Stokes1851) and Molin & Etienne (Reference Molin and Etienne2000), this force was linear and at $\pi /4$ phase to the displacement. The magnitude was, however, extremely small (less than 0.5 % of the measured force), and thus hardly made a difference.

Using (8.2) and (8.3) on the linear and experimental data, and (8.3) on the volume flow numerical data, then gives figure 8, which shows the mean power in the near and far fields across the range of tested displacement amplitudes. The vertical axis is normalised by the linear power. The nonlinear changes in power are not the same in the near- and far-field results, with the power in the near field registering smaller changes with amplitude than in the far field, and a small offset between the smallest amplitudes. The offset is suggestive of additional linear damping, and the slightly different shapes of the curves are possibly indicative of quadratic drag – effects that could increase the near-field power but decrease the far-field power. Nevertheless, a systematic, monotonic decrease in radiated power is observed in all curves, with the dominant trend appearing to be proportional to $(A/R)^2$ , as would be expected as the leading-order effect of a third-order contribution. The volume flow model gives a reasonable match to the far-field data.