2. Domain decomposition and matching methodology

To develop the axisymmetric hybrid model, the fluid domain is divided into two regions: the core region above the cylinder and an exterior region (figure 1). The flow in the exterior region is modelled using the linear potential flow eigenfunction expansion method, following Reference Siddorn and Eatock TaylorSiddorn and Eatock Taylor (2008), while the flow in the core region is modelled using the shallow water equations. As indicated in figure 2, in the exterior region we decompose the potential into a component due to the heave oscillations of a surface-piercing cylinder and a component due to the flow at the boundary between the exterior and core region, as in Reference YumYum (1985) (and a similar decomposition in the discussion of Reference Newman, Sortland and VinjeNewman et al. (1984)). Since the flow in the core region is assumed to be hydrostatic, we approximate the effect of the flow at the boundary on the exterior region by a radially oscillating ‘piston’ boundary, where the fluid velocity normal to the piston surface at $r=a$ is constant over $-s\leq z\leq 0$.

The potential in the exterior region due to the oscillating surface-piercing cylinder is not described here as it follows directly from Reference Siddorn and Eatock TaylorSiddorn and Eatock Taylor (2008). The potential due to the oscillating piston surface at the boundary between the exterior and core region follows the same method, differing only in the boundary conditions at $r=a$, $-s\leq z\leq 0$ and on the bottom surface of the cylinder.

Figure 2. Solution decomposition of exterior and core region into flow due to a heaving surface-piercing cylinder and radially oscillating piston.

As noted in § 1, the hybrid model developed here is compared with experiments and the full linear potential flow theory (PT) solution of Reference McCauley, Wolgamot, Orszaghova and DraperMcCauley et al. (2018), where the velocity potential and its radial derivative are matched across the boundaries between the exterior and core regions. In order to test the time-domain matching methodology used in the hybrid model, the nonlinear hybrid model and PT solution are also compared with a time-domain linear hybrid model, where the linearised shallow water equations are used in the core region, though the model is otherwise identical to the nonlinear hybrid model.

For all but the PT model, the total surface elevation at the piston boundary, $\eta _{E}$, is composed of $\eta _{SP}$ due to the oscillating surface-piercing cylinder and $\eta _{PE}$ due to the oscillations of the solid piston boundary. For the hybrid models where a time-domain solution is used in the core, $\eta _{E}$ must be computed in the time domain. This is achieved using IRFs based on the pre-computed frequency-domain solutions to the exterior problems.

3. Time-domain modelling

3.1 Shallow water equations

3.1.1 Nonlinear hybrid model

To account for the nonlinearity of the flow in the core region, we use the nonlinear shallow water equations for cylindrical symmetry, given by Reference ToroToro (2001) (p. 246 (13.6)). If we assume the water depth in the core region is the sum of the mean water level $s$ and the surface elevation variation $\eta _{sw}$ about that mean, and the bed is flat, then the equations are

(3.1)\begin{equation} \left[\begin{array}{c}\eta_{sw} \\ u \end{array}\right]_t + \left[\begin{array}{c}(\eta_{sw} + s)u \\ \dfrac{1}{2}u^2 + g_m\eta_{sw} \end{array}\right]_r = {-}\frac{1}{r}\left[\begin{array}{c}(\eta_{sw} + s)u \\ 0 \end{array}\right], \end{equation}

where $u$ is the radial horizontal depth-averaged velocity. Note that this formulation is in the moving reference frame on top of the cylinder, such that the total surface elevation in the core region in the stationary reference frame is $\eta =\eta _{sw}+\zeta$, where $\zeta$ is the time-dependent heave displacement. The moving bed is accounted for by modification of gravity $g$ by the cylinder acceleration $g_m=(g+\partial _t^2\zeta )$. These equations may be solved by a standard splitting scheme (see Reference ToroToro, 2001) to account for the source term. The finite-volume method is used, along with a HLL Riemann solver and the MUSCL-Hancock scheme. The shallow water solver has been verified for small-amplitude cases using the solution of Reference LambLamb (1993) (see p. 275) for canals of varying cross-section. The equations are implemented with a wetting and drying algorithm (see Reference Brufau, Vázquez-Cendón and García-NavarroBrufau, Vázquez-Cendón, & García-Navarro, 2002) allowing dry-bed scenarios to be simulated by the nonlinear hybrid model.

3.1.2 Linear hybrid model

The linear hybrid model uses the same matching at the boundary as the nonlinear hybrid model, but solves the linearised shallow water equations. Assuming that $u$, $\zeta$ and $\eta _{sw}$ are small quantities, whose squares and products can be ignored in comparison with the linear terms, the equations simplify to

(3.2)\begin{equation} \left[\begin{array}{c} \eta_{sw} \\ u \end{array}\right]_t + \left[\begin{array}{c} su \\ g\eta_{sw} \end{array}\right]_r = {-}\frac{1}{r}\left[\begin{array}{c} su \\ 0 \end{array}\right]. \end{equation}

These are the linearised shallow water equations for axisymmetric flow which are equivalent to the equations used by Lamb for modelling tides in channels of varying width (see Reference LambLamb, 1993, p. 275 (1)). These equations are used in the linear hybrid model.

3.2 Impulse response functions

The surface elevation in the hybrid models is matched at the boundary between the exterior and core regions by use of IRFs for the potential flow solution. These allow the frequency-domain potential flow solution to be converted to a time-domain surface elevation at the boundary which is matched to the time-domain shallow water equations solver. We consider only forced heave motion, so the surface-piercing cylinder time-domain surface elevation may be pre-computed, while the piston convolution must be performed at every time step.

3.2.1 Heaving surface-piercing cylinder

The complex surface elevation coefficient $\hat {\eta }_{SP}(\omega )$ is Hermitian, implying the real part is an even function and the imaginary part is an odd function, such that the IRF can be written as

(3.3)\begin{equation} h_{SP}(t) = \frac{2}{\rm \pi}\int_{0}^{\infty}\boldsymbol{Re} \{\hat{\eta}_{SP}(\omega)\}\cos (\omega t)\, \mathrm{d}\omega. \end{equation}

(Note that $h_{SP}(t)$ has no relation to the water depth $h$.) Thus only the real part of the transfer function is required in computation of the IRF. The real part approaches zero for high frequency faster than the imaginary part, reducing numerical error due to truncating the above integral at a cut-off frequency. The time-domain surface elevation at the cylinder edge can then be calculated by the convolution of the IRF with the time series of body heave velocity:

(3.4)\begin{equation} \eta_{SP}(t) = \int_{-\infty}^{t}u_h(\tau)h_{SP}(t-\tau)\, \mathrm{d}\tau. \end{equation}

3.2.2 Piston

The IRF for the surface elevation at the boundary due to the radially oscillating piston requires special attention. Due to the presence of the oscillating piston at the free surface, the amplitude of the free-surface motion at that boundary does not decay quickly with frequency as in the case of the heaving surface-piercing cylinder. The transfer function may be calculated in the same way as for the heaving surface-piercing cylinder (the solution for this follows directly from the methods outlined in Reference Siddorn and Eatock TaylorSiddorn and Eatock Taylor (2008) and is omitted for brevity) and is plotted in figure 3. It can be seen that even at $ka=20$ the transfer function is relatively large in amplitude. We therefore derive a simpler approximation for the transfer function for more efficient computation. Assuming that at very high frequency the fluid below the cylinder has little effect on the propagating wave generated, we extend the vertical sides of the cylinder to the seabed (such that the boundary condition on the cylinder sides becomes $\partial \phi /\partial r=0|_{r=a, -h\leq z\leq -s}$ such that the fluid beneath the cylinder is ignored, where $\phi$ denotes the velocity potential). Then the complex surface elevation at $r=a$ considering only the propagating mode is

(3.5)\begin{equation} \hat{\eta}_{PE}^0 = \frac{-2\,\mathrm{i}\omega}{g} \left(k_0h + \frac{\sinh(2k_0h)}{2}\right)^{{-}1}\mathrm{cosh}(k_0h) [\mathrm{sinh}(k_0h)-\mathrm{sinh}(k_0(h-s))] \frac{H_0^{(2)}(k_0a)}{H_0'^{(2)}(k_0a)}, \end{equation}

where $H_0^{(2)}$ is a Hankel function of the second kind in standard notation and the prime denotes the radial derivative, as in Reference HavelockHavelock's (1929) classical cylindrical wavemaker. Using Bessel function approximations for large argument, it can be shown that ${H_0^{(2)}(x)}/{H_0'^{(2)}(x)}\cong \mathrm {i}$ for $x\rightarrow \infty$. Approximating $\mathrm {sinh}(x)$ and $\mathrm {cosh}(x)$ by $\mathrm {e}^x/2$ and wavenumber $k_0\approx \omega ^2/g$, we find that

(3.6)\begin{equation} \hat{\eta}_{PE}^0\cong \frac{2}{\omega}(1-\mathrm{e}^{-\omega^2s/g}), \end{equation}

for large $\omega$. This approximate transfer function solution for the propagating mode is plotted as a dashed line in figure 3, and shows good agreement with the real part of the exact solution for high frequency. The propagating mode approximation is used for computation of the piston IRF, i.e. $\boldsymbol {Re}\{\hat {\eta }_{PE}\}\approx (2/\omega )(1-\mathrm {e}^{-\omega ^2s/g})$. A similar approximation may also be derived for the evanescent mode contribution to the transfer function (shown in figure 3). Calculation of the IRF for the piston solution by (3.3) also requires special care. We split the integral into three parts:

(3.7)\begin{align} h_{PE}(t) & = \frac{2}{\rm \pi}\left(\int_{0}^{\omega_L}\boldsymbol{Re}\{\hat{\eta}_{PE}(\omega)\}\cos (\omega t)\, \mathrm{d}\omega + \int_{0}^{\infty}\frac{2}{\omega}(1-\mathrm{e}^{-\omega^2s/g})\cos(\omega t)\, \mathrm{d}\omega\right.\nonumber\\ & \quad - \left.\int_{0}^{\omega_L}\frac{2}{\omega}(1-\mathrm{e}^{-\omega^2s/g})\cos(\omega t)\, \mathrm{d}\omega\right). \end{align}

Figure 3. Frequency transfer function $\hat {\eta }_{PE}(\omega )$ for surface elevation at $r=a$, piston model. Here $s/a=0.16$, $T/a = 0.4$, $h/a=2.4$.

The first integral, which uses the exact transfer function, is evaluated numerically over a low-frequency range $\omega \leq \omega _L$, where $\omega _L$ is a suitable cut-off frequency. Similarly, the third integral is evaluated numerically. The second integral may be computed exactly, and the closed form solution is given by

(3.8)\begin{align} \int_{0}^{\infty}\frac{1}{\omega}(1-{\rm e}^{-\omega^2s/g})\cos(\omega t)\, \mathrm{d}\omega & = [gt^2 \,_2F_2(1,1;3/2,2;-gt^2/4s)\nonumber\\ & \quad -2s(\gamma-\ln(s/g) + \ln(t^2))]/4s, \end{align}

given $t\in \mathbb {R}$, $\boldsymbol {Re}\{s/g\}>0$. Function $\,_2F_2(1,1;3/2,2,x)$ is a hypergeometric function in standard notation and $\gamma$ is the Euler–Mascheroni constant. Note that (3.8) results in $h_{PE}(t)\rightarrow \infty$ as $t\rightarrow 0$. This singularity has been observed in similar problems (see e.g. Reference RobertsRoberts, 1987 or Reference CointeCointe, 1988) and is an artefact of the ‘wavemaker-like’ piston problem used for matching. For best convergence properties, this singularity should be removed. However, for the purpose of initial comparisons to determine the feasibility of the hybrid model approach, it has been found that sufficiently accurate results can be achieved by calculating the convolution integral at a small time increment away from the singularity. In future work further effort should be devoted to removing the singularity. An example of the IRF is shown in figure 4, where it can be seen to quickly decay to zero. Comparison between numerically evaluated convolution of the IRF with harmonic piston motion, and the corresponding frequency-domain solution, confirm the validity of this approach, provided a sufficient number of data points is used near $t=0$. In practice the convolution integrals for the heaving surface-piercing cylinder and the radially oscillating piston are evaluated using the Prony method (see Reference MarpleMarple, 1987). The number of exponentials used in the Prony approximation is based on checking the convolution result with a harmonic input velocity, so that the result may be compared with the exact frequency-domain solution. The IRF for the piston surface elevation at the edge requires a large number of exponentials due to the steep gradient near $t=0$.

Figure 4. The IRF $h_{PE}(t)$ for surface elevation at $r=a$ for piston solution. Here $s/a=0.16$, $T/a = 0.4$, $h/a=2.4$.

4. Wave basin tests

Forced oscillation tests of a submerged buoyant cylinder were conducted in an experimental campaign at the COAST laboratory, University of Plymouth, UK, in 2018. The cylinder, of radius $a=0.625$ m, was positioned in the centre of the 35 m by 15.5 m wave basin. In the forced heave tests, motions of the cylinder were imposed using a taut three-cable arrangement with each cable controlled by a separate winch and passing through a pulley on the basin floor before connecting to the underside of the cylinder. Sinusoidal cable motions were imposed (displacement control), but due to the inclined cable geometry the resulting heave motion profiles differed slightly from sinusoidal. Forces in the cables were measured with load cells connected between the cable and the cylinder. Free-surface elevations were also measured at $r=0$. The motions of the cylinder were captured in six degrees of freedom by an optical tracking system (Qualisys). Further experimental details are given in Reference Orszaghova, Wolgamot, Draper, Taylor and RafieeOrszaghova et al. (2020).

A summary of the experimental cases is given in table 1, though results are not shown for all cases for brevity. The model was tested at two submergences, $s/a = 0.08$ and $0.16$, with the majority of tests at $s/a=0.16$ (hence results for $s/a = 0.08$ are shown only in the supplementary material available at https://doi.org/10.1017/flo.2023.32). A range of amplitudes were tested for each submergence and frequency. The model geometry and water depth were constant for all tests, with $T/a = 0.4$ and $h/a=2.4$. Figure 5 shows the frequency–submergence space of interest for the tests. The colour indicates the (normalised) heave radiation damping coefficient $b_h$ (predicted by the linear PT flow model). The first heave resonance, which corresponds to a peak in the radiation damping, can be clearly seen. The heave frequencies in the experiments (shown by the crosses) were chosen to span this heave resonance. The vertical lines indicate wave periods (at full scale) of 8 s and 16 s, the expected incident wave period range if the CETO device were to be deployed off the southwest coast of Australia. Figure 6 shows two snapshots of the free surface in experiments for $ka=0.35$, $A/s=0.33$ and $ka=0.60$, $A/s=0.31$, both at $s/a=0.16$. A large water jet at the centre of the cylinder is evident at $ka=0.60$.

Table 1. Experiment parameters.

Figure 5. Frequency-dependent heave radiation damping coefficient $b_h$ (from linear PT) for the experimental cylindrical geometry at different submergence depths. The experiment test frequencies are shown by the crosses.

Figure 6. Experiment snapshots: (a) $ka=0.35$, $A/s=0.33$; (b) $ka=0.60$, $A/s=0.31$. Both cases for $s/a=0.16$ and taken just after maximum positive heave displacement. The wave gauge at $r=0$ is fixed to the Qualisys stalk seen in the centre of the cylinder.

5. Comparison with experiments

Force and surface elevation data from experiments are compared to the nonlinear hybrid model, the hybrid model with linearised shallow water equations and a linear PT model for a submerged cylinder as described in Reference McCauley, Wolgamot, Orszaghova and DraperMcCauley et al. (2018) (see § 2 for explanation of each model). The linear potential flow model time-domain results are generated using IRFs for the surface elevation and force given input heave body motion from the experiments.

The heave displacement measured in the experiments by Qualisys is used to force the hybrid model, after a low-pass filter with cut-off at $6.5\omega$ is used to remove noise from the velocity and acceleration time series (which are calculated from the displacement). The IRFs for the heaving surface-piercing cylinder and the piston are pre-calculated. The surface elevation at the matching boundary due to the heaving surface-piercing cylinder may be calculated for all times as the body velocity is known. However, the convolution integral for the surface elevation at the cylinder edge due to the piston must be evaluated at every time step, as it is dependent on the fluid velocity at the boundary from the shallow water equations solver.

The heave force may be calculated once the piston velocity and surface elevation in the core are known. Impulse response functions for the force on the lower cylinder surface due to motion of the heaving surface-piercing cylinder and the piston are calculated similarly to the surface elevation IRFs shown in § 3.2. The force on the lower cylinder surface also includes a correction for the change in hydrostatic pressure due to the cylinder motion (since the force on the top face is evaluated based on the exact position of the cylinder in time). The force acting on the top cylinder surface is given by integrating the surface elevation from the shallow water solver over the area (i.e. the hydrostatic force).

Figure 7 shows the surface elevation measured above the centre of the cylinder (top row of each panel), cylinder heave displacement measured in the experiments and used to force the hybrid and linear models (middle row of each panel) and the heave hydrodynamic force (bottom row of each panel). Force is defined as positive upwards. Results from three models are shown and the position of the top surface of the cylinder is shown in the top row by the solid black line (except in cases where the free-surface motions are small and the top surface is out of the range plotted). The vertical lines in this figure indicate the predicted arrival time of reflections from the basin sidewalls (calculated using the group velocity of the fundamental frequency); comparison to the left of these lines is therefore most relevant. It can be seen that at high frequency and near the heave resonance (around $ka=0.6$), the heave motion deviates significantly from the desired sinusoidal signal. As the heave motion was controlled via a taut three-cable arrangement, which could not apply an upward force (the only upward force on the cylinder was from net buoyancy), at high frequencies the required upward velocity was not always reached and the tethers became slack in some cases. Even though the non-harmonic cylinder motion is undesirable, the experimental data nevertheless provide valuable validation since the calculated free-surface and heave force responses are based on the measured motion signals input into the numerical models. We also note that noticeable motion instabilities eventually developed in yaw, pitch and roll at some forcing frequencies (see Reference Orszaghova, Wolgamot, Draper, Taylor and RafieeOrszaghova et al., 2020). However, these built up over longer time scales and therefore did not affect the initial oscillation periods before reflections from the basin walls arrived (maximum rotation in any mode at the time of reflections from the tank walls is 0.64$^{\circ }$).

Figure 7. Surface elevation at the cylinder centre (top), heave displacement (middle) and heave hydrodynamic force (bottom): (a) $ka=0.35$, $A/s=0.33$, (b) $ka=0.42$, $A/s=0.34$, (c) $ka=0.60$, $A/s=0.31$, (d) $ka=0.77$, $A/s=0.29$. All plots: $s/a=0.16$, $h/a=2.4$, $T/a=0.4$.

Figure 7 shows larger-amplitude cases ($A/s = 0.29\unicode{x2013}0.34$) for increasing frequency ($ka=0.35\unicode{x2013}0.77$) at submergence $s/a=0.16$ (a similar plot for $s/a=0.08$ is given in the supplementary material). The nonlinear hybrid model performs significantly better than the linear PT model in predicting the force in all cases. The linear hybrid model crashes in all cases except at $ka=0.35$, due to the water level (in one or more cells) on the top of the cylinder dropping to zero. This cannot be simulated by the linear shallow water equations model, but can be handled by the wetting and drying scheme implemented in the nonlinear hybrid model. At frequencies above $ka=0.42$ the linear PT model is seen to predict the surface elevation above the cylinder (at $r=0$) to oscillate below the top face of the cylinder, a clear indicator that the model is outside its range of applicability. In all cases the hybrid model provides an improvement over the linear PT model. An estimate of the error in the nonlinear hybrid model and linear PT model total heave force when compared with experimental measurements (for $s/a=0.16$, $50< U<877$) is given in figure 8. The error is calculated as the root-mean-square difference between the experimental force and the linear PT or hybrid model force in the first five oscillations. The value is normalised by the maximum of the experimental force in the same range. It is clear that the nonlinear hybrid model more accurately predicts the heave force for all cases, and that with increasing amplitude the hybrid model performs increasingly better than the linear PT model. Though the nonlinear hybrid model performs significantly better than the linear model, some discrepancy with the experiment remains. Some of this discrepancy may be expected to be associated with viscous losses due to flow around the lower corner of the cylinder and the ‘jet-like’ separation of flow coming off the top of the cylinder. The former problem is common in offshore engineering (e.g. Reference HuseHuse, 1990) while the latter is more unusual. The results indicate that the shallow water flow is indeed the dominant nonlinearity, with acceptable results achieved without attempting to quantify viscous losses. Model–experiment comparisons for the shallower submergence $s/a = 0.08$ are shown in the supplementary material and demonstrate even better performance in this shallower regime.

Figure 8. Root-mean-square (RMS) error in the linear PT and hybrid model force when compared with experiments, $s/a=0.16$.

It may be seen that the linear hybrid and linear PT models give very similar results in all cases. The deviation between these two models is small compared with the difference between either model and the nonlinear hybrid, and can be attributed to the matching process between the shallow water equations and potential flow components of the hybrid model. Thus, the error introduced by the matching method is not significant compared with the shallow-flow effects in the range tested. The discrepancy between the linear hybrid model and linear PT model is also small compared with the difference between the linear hybrid model and experiments.

Large amplification of the free-surface motion is seen at the centre of the cylinder, observed in experiments as a vertical water ‘jet’ which occurs when a radially converging bore reaches the centre of the core region (see figure 6b). This feature is most pronounced at high frequency (see figure 7c,d), where the wave gauge measurements of the surface elevation are over twice the submergence. The nonlinear hybrid model cannot capture this striking, but very localised, feature of the flow as the shallow water equations do not account for vertical accelerations. However, the hybrid model does capture the arrival of a shock or bore front at a similar time to the jet. There is little evidence from the experiments that the jet significantly affects the total force on the top surface of the cylinder, and at the time at which it occurs the heave force from the hybrid model and the experiments are generally in good agreement.

6. Analysis with harmonic input forcing

Having validated the nonlinear model with experimental data, we now use it (within the range tested) to systematically explore the trends in the heave force for harmonic motion with increasing amplitude. As shown in § 1, calculation of the added mass and damping coefficients has been the focus of much previous work on submerged oscillating structures. It is possible to extract nonlinear heave ‘added mass’ and ‘damping’ coefficients at the fundamental frequency from the hybrid model by isolating the first harmonic of the heave force. In these simulations the cylinder is forced to oscillate harmonically. Once the resulting heave force has reached a steady state, the signal is trimmed (to contain a minimum of eight oscillations) and subsequently band-pass filtered (using $0.75 \omega$ and $1.25 \omega$ cut-offs) to isolate the first harmonic. The signal is converted back to the time domain and a sine function fitted to determine the linear amplitude and phase of the force. The resulting force coefficients are plotted in figure 9 for submergence $s/a=0.16$ and $0.3< ka<0.8$ against the linear potential flow solution (solid black line). The nonlinear hybrid model results are plotted as solid coloured lines and those of the linear hybrid model as a dashed line, with the colour of the line indicating the amplitude of the heave motion. The linear hybrid model results are independent of amplitude (except at large amplitude where the shallow water equations ‘go dry’ and the model breaks down); thus the result is only plotted for the smallest amplitude.

Figure 9. Nonlinear heave added mass (a) and damping (b) at the fundamental frequency extracted from hybrid and linear PT models. Here $s/a=0.16$, $h/a=2.4$, $T/a=0.4$.

At low and high frequency both hybrid model solutions tend towards the linear PT solution (but do not reach it due to approximations introduced in the matching; see supplementary material). For small amplitude $A/s=0.01$ the nonlinear hybrid model agrees well with the linear hybrid model, with small discrepancy in the added mass near the resonance. With increasing amplitude of motion the added mass is reduced for $ka<0.60$ and increased for $ka>0.65$ (though it remains negative in this region). The damping is reduced for $ka>0.55$ and increased for $ka<0.5$, such that we see a shift in the damping peak to lower frequency. The damping peak has shifted from $ka=0.59$ for $A/s=0.01$ to $ka=0.54$ for $A/s=0.4$. For the smaller submergence case where $s/a=0.08$ (not shown for brevity) the effects of increasing amplitude have greater influence on both hydrodynamic coefficients, with a reduction of over 50 % in the damping peak for $A/s=0.2$ as compared with the linear model. Although these plots ignore the higher harmonics in the force, it is clear that the linear model is insufficient for predicting the added mass and damping coefficients where the amplitude of the heave motion is greater than $A/s=0.1$. While obvious, it should be noted that these curves are valuable for illustrating nonlinear changes in behaviour, but cannot be interpreted and used in the way in which linear added mass and damping coefficients are used – the nonlinear time domain model should be used instead.

While examining the nonlinear coefficients at the fundamental frequency helps to relate the hybrid model results to conventional understanding, the full force time histories provide more information; results for four frequencies are plotted in figure 10 for $A/s=0.01\unicode{x2013}0.5$ and $s/a=0.16$. Changes in the volume of water in the core region (away from the mean) are also plotted for each frequency and amplitude below the corresponding force plot. The coloured vertical lines in each plot correspond to the times of maximum and minimum force for the corresponding amplitude. Free-surface profiles for two frequencies ($ka=0.35$ and 0.55) and two amplitudes ($A/s=0.2$ and 0.5) are shown in figure 11. The force on the bottom cylinder surface is generally small relative to the force on the top face, except at low frequency where they are 180$^\circ$ out of phase, making the second harmonic force component more prominent. Hence, the force on the top face dominates for most frequencies. Three nonlinear effects are particularly noticeable in figure 10:

1. (i) At each frequency, a delay in the phase of maximum force as the amplitude of oscillation increases, a large effect most important at low frequency.

2. (ii) At each frequency, a smaller delay in the phase of minimum force as the amplitude of oscillation increases, most apparent around the linear resonance frequency.

3. (iii) Significant reductions in the amplitude of the maximum and minimum forces as the amplitude of oscillation increases, for frequencies around (linear) resonance.

Figure 10. Total heave force (top row) and change in volume in core region (bottom row): (a) $ka=0.35$, (b) $ka=0.45$, (c) $ka=0.55$, (d) $ka=0.65$. Heave displacement is plotted on the right-hand axis. The vertical lines indicate the times of maximum and minimum force for each amplitude. All plots: $s/a=0.16$, $h/a=2.4$, $T/a=0.4$.

Figure 11. Free-surface snapshots showing the flow on (top row) and the flow off (bottom row) phases for two amplitudes at (a) $ka=0.35$ and (b) $ka=0.55$. The arrow represents the direction and amplitude of the velocity at the boundary and the cylinder position for each amplitude is shown by the dashed lines. All plots: $s/a=0.16$, $h/a=2.4$, $T/a=0.4$.

To explain these observations, we divide a cylinder oscillation into two phases: the ‘flow on’ phase (where the fluid velocity at the piston boundary is onto the cylinder) and the ‘flow off’ phase (where the fluid velocity is off the cylinder). The flow on phase starts with a shallow flow propagating onto the cylinder and ends when the reflection from the cylinder centre propagates off the edge, at the same time as (or just after) the minimum in force occurs (corresponding to maximum fluid on the cylinder). For large-amplitude oscillations, the shallow flows on the cylinder dominate during the flow off phase. The strong flow off acts to delay the reversal of flow at the boundary and results in a phase delay in the maximum force with increasing amplitude. This partly explains observation (i), which is also influenced by the second harmonic at low frequency. Observation (ii) suggests that despite starting later, during the flow on, the bores at larger amplitude reflect from the centre and back to the edge more quickly, a nonlinear effect. At higher frequencies ($ka>0.35$) the cylinder is already moving upwards when flow off starts. The surface elevation on the cylinder near the centre drops to zero (or close to) as it starts to move downwards (see figure 11b), resulting in a minimum in the fluid volume on top of the cylinder and a maximum in the force. For larger amplitudes and at higher frequency the cylinder is almost completely dry on top, such that the non-dimensional peak in force must be reduced with amplitude, as the volume of fluid on the cylinder does not scale. Similarly, at higher frequency the cylinder is closer to the mean free surface (relative to the amplitude) when the maximum volume occurs, resulting in a reduction in the volume of fluid on top of the cylinder (see figure 11b). Hence the force on the top cylinder surface cannot scale linearly with amplitude. Thus the maximum and minimum force reductions in observation (iii) are explained.

The phase shift in the peak force at low frequency results in the increase in nonlinear damping with amplitude observed earlier, since the force is closer to the phase of the velocity. However, at higher frequencies ($ka>0.45$) the decrease in relative magnitude of the force with increasing amplitude overcomes the phase shift and the damping decreases with increasing amplitude.