Parametric excitation suppression in a floating cylinder via dynamic vibration absorbers: a comparative analysis

Parametric excitation in the pitch/roll degrees of freedom (DoFs) can induce dynamic instability in floating cylinder-type structures such as spar buoys, floating offshore wind or wave energy converters. At certain frequency and amplitude ranges of the input waves, parametric coupling between the heave and pitch/roll DoFs results in undesirable large amplitude rotational motion. One possible remedy to mitigate the existence of parametric resonance is the use of dynamic vibration absorbers. Two prominent types of dynamic vibration absorbers are tuned mass dampers (TMDs) and nonlinear energy sinks (NESs), which have contrasting properties with regard to their amplitude and frequency dependencies when absorbing kinetic energy from oscillating bodies. This paper investigates the suppression of parametric resonance in floating bodies utilizing dynamic vibration absorbers, comparing the performance of TMDs against NESs for a test case considering a floating vertical cylinder. In addition to the type of dynamic vibration absorber utilized, the paper also examines the DoF which it acts on, comparing the benefits between attaching the vibration absorber to the primary (heave) DoF or the secondary (pitch) DoF. The results show that the TMD outperforms the NES and that it is more effective to attach the vibration absorber to the heave DoF when eliminating parametric resonance in the pitch DoF.


Introduction
Dynamic instability in floating bodies can occur due to the parametric excitation of the pitch and roll degreesof-freedom (DoFs). This phenomena was first noted by Froude in the nineteenth century [1], and has since been investigated in applications such as shipping [2][3][4], offshore spar-type structures [5][6][7][8] and wave energy converters [9].

Suppression of parametric resonance in floating bodies
For most offshore applications, parametric resonance is undesirable. The large amplitude motions occurring due to parametric resonance can have an adverse effect on the performance and safety of offshore structures. For example, container ships lost cargo overboard due to large parametric roll motions [10]. As such, meth-ods for suppression and stabilization of parametric pitch/roll in ships have been investigated, see for example [11][12][13][14] and the references therein. Considering spar-type structures, Haslum et al. [15] investigate alternative hull shapes to improve the stability of spar platforms in hostile areas. Rho et al. [16] consider the stability of a spar platform with an approximate 2:1 ratio of pitch to heave natural periods, performing experiments and numerical analysis for regular waves. They consider including damping plates and strakes onto the cylindrical spar and find that the use of damping plates to limit the heave motion has the greatest effect of reducing the parametric pitch motion. Yang and Zhou [17,18] investigate the occurrence of parametric resonance for a spar in irregular waves and means to avoid it. In [17], they consider the effect of hull geometry on the occurrence of parametric resonance in irregular waves. They compare five different spar geometries and conclude those geometries which minimize the heave response also minimize the occurrence of parametric resonance. In [18], some suggestions for designing spar platforms to avoid parametric resonance are presented, such as increasing the damping coefficient, designing the frequency response appropriately and increasing the metacentric height.
A number of studies have investigated suppressing parametric resonance in wave energy converters (WECs). All except for one of these studies, consider parametrically excited pitch/roll motions in spar-type WECs. Nicoll and Wood [19] mitigate the occurrence of parametric resonance in the sway DoF for the Surf-Power WEC, which is tethered to the sea floor and acts like an inverted pendulum. Nicoll and Wood show that the parametric sway motion can easily be detuned by adjusting the mooring line tension, which is controllable through the hydraulic pressure in the power take-off (PTO). Considering a two-body heaving WEC with a central spar, the Wavebob, Villegas and van der Schaaf [20] also utilize the PTO system, applying a notch filter to eliminate the coupling between the two bodies at the resonances responsible for the parametric roll. Also considering a Wavebob-type WEC, Beatty et al. [21] investigate the use of fins to prevent parametric pitch and roll by increasing the drag-induced damping in those DoFs. Follow-up work on the same device, by Ortiz [22], investigates an alternative method which utilizes the mooring dynamics to reduce the occurrence of PR. The use of fins was also investigated by Gomes et al. [23] to suppress parametric pitch and roll in the OWC spar buoy, which was previously observed to be prone to parametric resonance [24,25]. Also considering the OWC spar buoy, Davidson et al. [26,27] investigated the use of a pressure relief valve in the OWC chamber to decouple the motion of the spar and internal water column, analogous to the notch filter used in [20].
One suppression method which has not been thoroughly examined is the use of dynamic vibration absorbers (DVAs), as detailed in the next section.

Dynamic vibration absorbers (DVAs)
DVAs consist of a small mass attached to a primary structure through a visco-elastic connection. Establishing a modal coupling between the primary system and the absorber, part of the vibration energy of the primary system is transferred to the absorber and dissipated, resulting in an overall reduction in vibration amplitude. In terms of frequency response, the use of tuned mass dampers (TMDs) enables the suppression of an essentially undamped high resonant peak, by generating two largely damped smaller peaks, thus avoiding large amplitude oscillation of the primary system [28]. The same principle can be used to avoid selfexcited oscillations, for instance in wing profiles subject to flutter oscillations [29]. This kind of devices were first proposed for reducing rolling motions in ships [30,31] by Watt already in 1883. However, their design was formalized in more rigorous terms only few decades later by Ormondroyd and Den Hartog [28,32]. Nowadays, TMDs are extensively used in engineering practice for a variety of application; between the most famous examples, we cite the one installed in the Taipei 101 skyscraper [33] and the series of TMDs in the Millenium Bridge in London [34]. The main advantages of the TMD are related to its low cost, easy implementation and low maintenance required.
The main limitation of the TMD is related to its narrow operation bandwidth. Indeed, this prevents the TMD to work efficiently in systems undergoing large oscillations at various frequencies or systems with variable natural frequency. This limit motivated the development of alternative DVAs. In particular, the nonlinear energy sink (NES) is analogous to a TMD, but it is connected to the primary system through an essentially nonlinear spring, instead of a linear one. This allows it to resonate virtually at any frequency, enabling broadband energy dissipation [35]. However, an efficient modal coupling between a primary system and an NES is obtained only at a specific amplitude range. In particular, below a specific energy threshold, the NES becomes practically ineffective [36]. Besides many potential applications, its implementation was also proposed for a submerged body oscillating on an elastic support [37], however to the authors' knowledge, has not been investigated for a floating body.

Objective and outline of paper
The main objective of this work is to perform a comparative analysis of the TMD and NES vibration absorbers for the suppression of parametric excitation in a floating cylinder. Considering their different range of operation, the TMD and the NES have somehow complementary properties. The TMD works very efficiently only at its natural frequency, but at virtually any vibration amplitude, while the NES works only at specific amplitude ranges, but over a broad frequency band. A floating body may be subject to excitation from waves that span a large frequency and amplitude range, therefore presenting a challenging application for both the TMD and the NES.
In addition to the type of DVA, the comparative analysis also considers the degree of freedom (DoF) in which the DVA acts on. Parametric resonance results from the transfer of energy from a primary DoF to a secondary DoF. For a floating cylinder, the primary DoF is heave, which can parametrically excite the pitch and roll DoFs. Therefore, the study will compare the performance of the DVAs acting on both the primary heave and the secondary pitch DoFs to analyse which configuration is more effective at suppressing the occurrence of parametric resonance. The comparative analysis is performed via numerical simulations of a floating cylinder which is known to be prone to parametric resonance due to the close proximity of the heave natural frequency to twice the pitch natural frequency.
The remainder of the paper is laid out as follows. Section 2 describes the comparative study, providing details on the floating cylinder, the considered DVA configurations, the modelling and simulation of the system, the input wave conditions considered and the metrics used to assess and compare the performance of the different DVA configurations. The results are then presented in Sect. 3, and conclusions are drawn in Sect. 4.

The comparative study
The comparative study is performed using a numerical testbed, evaluating the ability of the different DVA configurations to suppress the occurrence of parametric resonance in a floating cylinder. The floating cylinder considered in this study is described in Sect. 2.1. The hydrodynamic model employed to simulate the wave induced motion of the cylinder is then detailed in Sect. 2.2. The attachment of the DVAs to the cylinder and their inclusion in the numerical model is presented in Sect. 2.3, along with a discussion on the selection of the optimal DVA parameters. The methodology for the comparative analysis is then described in Sect. 2.4. Finally, in Sect. 2.5, the simulation details for the numerical testbed are given.

The floating cylinder
The floating cylinder employed in this study is a widely used test case from the literature considering parametric resonance in spar-type structures. The natural period of the cylinder's pitch/roll motion is about twice the heave resonant period, which makes the floating cylinder susceptible to parametric resonance. Hong et al. [38] first introduced this case study, performing modelscale tests on the floating cylinder in a wave tank. Numerous authors then investigated various aspects relating to the dynamic instability and parametric excitation of the cylinder using numerical models [8,[39][40][41][42][43][44][45]. A schematic of the cylinder is shown in Fig. 1, and the relevant parameters are listed in Table 1.

Hydrodynamic modelling
A brief introduction on articulating parametric resonance in hydrodynamic models for floating bodies is first given in Sect. 2.2.1, to provide the reader with a background on the main driver of the nonlinear behaviour herein analysed. The particular model employed in this comparative study is then detailed in Sect. 2.2.2.

Modelling parametric resonance in floating bodies
Parametric resonance is generated by the time-varying changes in the parameters of a system [46]. For floating  bodies, this condition is met due to the wave-induced variation in the pitch and roll restoring stiffness. The pitch and roll restoring stiffness parameters are proportional in magnitude to the distance between the centre of buoyancy (CoB) and the centre of mass (CoM). In the presence of waves, the relative motion between the oscillating water surface and floating body continually alters the position of the CoB, thereby inducing time variation in the pitch and roll stiffness parameters. The mathematical consideration of parametric excitation requires differential equations with time-varying coefficients [47]. Floating structures, like many dynamical systems, can be modelled as a second-order ODE, comprising an inertia term, a resistance/damping term, a restoring term and an excitation term, analogous to a mass-spring-damper system in mechanics or an L-R-C system in electronics [48]. Traditionally, each of these terms consists of a linear time-invariant (LTI) function, allowing computationally efficient frequency domain analysis to be employed when considering the vast range of wave spectra likely to be encountered by the structure in the ocean. The function parameters in each of these terms depend on the geometry of the floating structure's wetted surface relative to its CoM. The LTI approach considers constant parameters derived from the wetted surface at the floating structure's equilibrium position, based on the assumption of small amplitude waves and body motion. However, when considering larger amplitude waves and/or body motion, for example if the wave frequency matches the natural frequency of the floating body and its start to resonate, then the LTI approach breaks down [49].
Modelling parametric resonance in floating structures can therefore be achieved employing a timevarying function for the restoring term, which articulates the variation in the pitch and roll restoring stiffness parameters, due to the evolving wetted surface, as the body submergence oscillates in the presence of waves. For simple geometries, the time-varying function for the restoring terms can be derived analytically [50,51]. Alternatively, system identification techniques can be employed from measured data to derive a suitable function to encapsulate the variation in the restoring parameters [52].
A higher fidelity but more computationally costly approach is the nonlinear Froude-Krylov force (NLFK) method. The Froude-Krylov force is the integral of the static and dynamic pressure fields of the undisturbed wave over the wetted surface of the body, which can be calculated using boundary element methods. Since the Froude-Krylov force considers the undisturbed wave, it excludes the effects of wave diffraction on the body. However, for bodies of small extension compared to the wavelength, diffraction is negligible, and the Froude-Krylov force is a good approximation to the sum of the restoring and excitation forces [48]. The NLFK method, calculates the Froude-Krylov force at each time step based on the instantaneous wetted surface of the body. Therefore, the parameters of the restoring and excitation force terms are continuously recalculated during the simulation and the NLFK method can articulate parametric resonance [53][54][55].
The most high-fidelity, but computationally expensive, method is to employ computational fluid dynamics, or other similar numerical modelling tools, to implement a numerical wave tank (NWT) [56]. The NWT domain is typically discretized spatially into over a million cells/nodes/particles and the Navier-Stokes equations are numerically solved in each of these cells/nodes/particles at each time step to calculate the evolving fluid pressure, velocity and free-surface elevation in the water. The motion of the floating body is calculated by a coupled solver which integrates the fluid pressure over the body surface at each time step.

The hydrodynamic model for the comparative study
The hydrodynamic model follows the validated model from the previous studies of this cylinder in [8,38,39,45], comprising the following coupled equations for heave x 3 and pitch x 5 where F I,i is the inertia, F R,i is the hydrostatic restoring force, F D,i is the hydrodynamic damping force and F E,i is the wave excitation force. In addition to the experimental validation performed in [38], the ability of this hydrodynamic model to correctly predict parametric resonance is also verified in a code-to-code comparison against a NLFK model in [57]. The inertial forces are the product of the acceleration and inertia: where M is the mass of the cylinder, m 3 is the hydrodynamic added mass in heave, I 5 is the moment of inertia in pitch of the cylinder and m 5 is the hydrodynamic added moment of inertia in pitch. The hydrostatic restoring terms include time-varying parameters, allowing the occurrence of parametric resonance to be possible in the hydrodynamic model: where ρ is the water density, g is the gravitational acceleration, A C is the horizontal cross-sectional area of the cylinder, L M S is distance from the still water level to the centre of mass, L D is the draft, G M is the metacentric height and η(t) is the wave elevation. The hydrodynamic damping is identified from the experimental data in [38] and is represented by where C i are the experimentally determined damping coefficients. The excitation force, for an input wave with amplitude A, frequency ω and phase φ, is given by where H E,i (ω) is the excitation force gain and φ E,i (ω) the excitation force phase shift. When considering polychromatic waves, the excitation force consists of a linear combination of N f frequency components.
The excitation force gain and amplitude shift are frequency-dependent. Their values for the test case cylinder are calculated using Nemoh [53], which is an open-source linear potential flow boundary element method software. The calculated excitation force coefficients are plotted in Fig. 2 as a function of frequency. Table 2 lists the parameter values used in the model.

Dynamic vibration absorbers
where M a is the mass of the absorber and M c the mass of the rest of the cylinder. Similarly, the addition of the DVA is assumed to not alter the moment of inertia of the cylinder, I 5 .

Heave
When acting in the heave DoF, the DVA is attached to the cylinder through a spring and a resistive damper and oscillates along the central axis of the cylinder, as shown in Fig. 3a. The force applied by the DVA on the cylinder is given by: where x a3 is the translational displacement of the DVA mass (m a3 ) along the cylinder axis from its rest position at the CoM, b a3 is the damping coefficient, k a3 is the linear spring coefficient and k nl3 is the nonlinear spring coefficient. For the TMD, k nl3 = 0, and for the NES, k la = 0. Following the derivation in "Appendix A", we obtain the following system of equation governing the dynamics of the system: F 3 and F 5 account for all hydrodynamic forces and inertial effects related to the water, as detailed in Sect. 2.2.

Pitch
For the pitch DoF, a rotational DVA is attached to the cylinder at its CoM, as shown in Fig. 3b. Examples of rotational DVAs can be found in applications such as rail [58] and vibration isolation [59]. The attachment of a rotational DVA acting in pitch has the following The DVA is coupled to the pitch DoF through the spring-damper torque: where x a5 is the rotation of the TMD with respect to the inertial reference frame. The equation of motion for the TMD is therefore given by: where I a5 is the DVA mass moment of inertia in the pitch direction. We thus obtain the following system of equation governing the dynamics of the system:

Parameter optimization
The performance of DVAs is dependent on the mass, spring and damping properties, thus these parameters must be optimized when assessing a particular DVA configuration. To aid in the optimization and provide a better understanding of the physical meaning of the results obtained, the following non-dimensional parameters are introduced, as listed in Table 3. The optimization procedure involves first setting the DVA's mass/moment of inertia as a fraction of the cylinder's mass/moment of inertia, denoted by the mass ratio parameter, ε. For a realistic and viable DVA, the size and mass should ideally be as small as possible to reduce its cost and obtrusiveness on the rest of the system. Therefore, the mass ratio value is an important assessment criteria in the DVA comparison, with smaller ε values being more favourable. Once the mass ratio is set then the natural frequency ratio, γ a , for the TMD or the pseudo-natural frequency ratio, κ a , for the NES, and the damping ratio, ζ a , for both, must be optimized, according to a specific objective, by tuning the damping and spring parameters.
Considering a single DoF primary structure, several studies provide the optimal damping and spring parameters, depending on different objectives that optimize for different metrics, such as: (a) minimizing the maximum displacement of the primary mass [32,60], (b) minimizing the total displacement of the primary mass over all frequencies [61,62], (c) minimizing the transient vibration of the system [63,64], (d) minimizing the relative displacement (i.e. the DVA stroke) [65], (e) minimizing the total kinetic energy of the primary mass Damping ratio ζ a ba3 2ma3ω3 b a5 2I a5 ω 5 over all frequencies [61], or (f) maximizing the power dissipated by the absorber [66]. The different optimization methods are summarized in [66]. Optimal performance of nonlinear energy sinks in multiple-degree-offreedom systems was studied by, for example, Tripathi et al. [67].
However, for the comparative study in the present paper, the system under consideration has two coupled DoFs, and does not have existing analytical solutions for the optimal parameters. Therefore, herein the optimal parameters will be found through a brute-force approach by sweeping through the parameter value combinations. This approach requires two key components: • A well-defined metric that encapsulates the objective of the optimization. This is discussed in Sect. 2.4. • A computationally fast method to simulate the system and calculate the metric for all of the parameter value combinations. This is discussed in Sect. 2.5.

The comparative analysis
The performance of the DVAs is assessed on their ability to suppress the occurrence of parametric pitch motion in the cylinder for a wide range of wave conditions. The input wave conditions are described in Sect. 2.4.1, and the performance metrics on which the DVA parameters are optimized and the different DVA configuration are compared are discussed in Sect. 2.4.2.

Input wave conditions
Parametric resonance is both frequency-and amplitudedependent, thus the performance of the DVAs are assessed in a range of input wave frequencies and amplitudes. In addition, the comparative analysis considers the performance of the DVAs in both singlefrequency regular waves and multi-frequency wave spectra. The regular waves encompass 8000 different wave conditions, parameterized by their frequency and amplitude, spanning an array of 80 frequencies and 100 amplitudes. Parametric resonance is expected at twice the pitch natural frequency; therefore, the wave frequencies range from 0.8 to 1.2 times this parametric resonant frequency (2ω 5 ). The wave amplitudes span 0.001 to 0.5 times the metacentric height of the floating cylinder (G M).
The irregular waves are derived from the JONSWAP spectra, comprising N f = 300 harmonic components (spanning 0.01 to 3 rad/s) with randomly assigned phases. The irregular waves also encompass 8000 different wave spectra, parameterized by the peak frequency and significant wave amplitude (half the significant wave height). The 80 peak frequencies range from 0.2 to 1.7 times the parametric resonant frequency. The 100 significant wave amplitudes span 0.001 to 5 times the metacentric height of the floating cylinder.

Performance criteria
The performance of the DVAs are assessed on their ability to suppress the occurrence of parametric resonance. Each particular DVA configuration, with a certain mass, spring and damping parameter combination, is simulated over the array of 8000 different wave conditions. The calculated cylinder motion for each wave condition is then analysed to detect the presence or absence of parametric resonance. From these results, the following two performance criteria are considered: • u * : the parametric resonance region ratio.
This criteria refers to the region in the frequencyamplitude space for which parametric resonance occurs compared to the total frequency-amplitude space spanned by the 8000 wave conditions. It is essentially the ratio of simulations for which parametric resonance was present to the total number of simulations. For example, if parametric resonance is present in 4000 of the 8000 simulations, then u * = 0.5. • A * : the minimum normalized wave amplitude which generates parametric resonance. This criteria refers to the wave amplitude (normalized by the cylinder's metacentric height) with the lowest amplitude across all frequencies for which parametric resonance is present.
These two criteria are used in the optimization of the DVA parameters and as a metric in the overall assessment of the DVA performance. The brute-force parameter optimization procedure, described in Sect. 2.3.3, involves sweeping through the parameter values and identifying where u * is minimum and A * is maximum. Some other criteria, which are also important in the overall assessment of the DVA performance but not explicitly used in the optimization of the DVA parameters, are the DVA mass and the maximum DVA displacement. The smaller the size of the DVA, in terms of weight and required displacement, the more attractive it is as a viable solution. Therefore, the DVA mass and maximum DVA displacement should ideally be as small as possible.

Simulation details
To identify the optimal performance of the four different DVA configurations, many combinations of DVA stiffness, damping and inertia values must be explored, each requiring 8000 simulations. Therefore, it is imperative that the simulations be computationally efficient. For this reason, to perform the optimization of the DVA parameters, we adopted the harmonic balance method (HBM) [68], as described in Sect. 2.5.1. Additionally, although significantly slower to compute than the HBM, direct numerical simulations are also necessary for the optimization and analysis for the following two reasons: 1) to verify the results produced by the HBM, and 2) because for certain conditions the HBM could not converge on a result, thus direct numerical simulations are needed to compute the results for these cases. The details of the direct numerical simulations are presented in Sect. 2.5.2.

Harmonic balance method
The HBM consists of approximating the steady-state solution of the system by a series of harmonic terms. The approximated solution is then substituted into the system of differential equations under study. Subsequently, all harmonic terms are collected, reducing the nonlinear system of differential equations to a nonlinear system of algebraic equations, which can be solved with root finding procedures, such as the Newton-Raphson algorithm. Harmonic terms not included in the trial solution are neglected. In order to solve the system for a given frequency (or amplitude) range, the solution obtained for a specific frequency can be used as trial solution for the following one, significantly reducing the computational time, which, overall, is several order of magnitude faster compared to regular numerical simulations.
In practice, we approximate the solution of the system by x a = c a1 cos(ωt) + s a1 sin(ωt) where x a indicates either x a3 or x a5 for the cases of the absorber acting in heave or in pitch, respectively. This formulation implies that we consider only the frequency of excitation ω and its half, ω/2, which is the frequency of the parametric excitation if triggered. Then, we plug this approximate solution into the equations of motion of the system [Eqs. (13) or (16)] for the absorber set either in heave or in pitch, respectively). For the case of the absorber in heave, before inserting the approximate solution, trigonometric terms of x 5 are expanded in Taylor series up to the fifth order.
Collecting terms with the same trigonometric order, we obtain a 12-dimensional nonlinear system of algebraic equations whose unknowns are c 31 , s 31 , c 32 , s 32 , c 51 , s 51 , c 52 , s 52 , c a1 , s a1 , c a2 and s a2 , in the form where k is a vector including all system parameters.
The algebraic equations f i are extremely lengthy and can be practically handled only with computer algebra. The solution of this system with the Newton-Raphson method provides an approximation of the steady state solution of the original system. In order to predict when parametric resonances are triggered, the stability of the obtained solution is investigated. The stability analysis is performed by determining if the real part of at least one of the Floquet exponents is positive. For this purpose, Hill's method is adopted [69], following the steps outlined in [70,71]. In practice, if the analysis shows that for certain parameter values the response is unstable, we assume that a parametric resonance is triggered. On the contrary, parametric resonance responses, are not directly investigated with the HBM in this work. This approach is rather standard for the analysis of this kind of dynamical behaviour.
We remark some limitations of the implemented HBM for the system under study. First of all, neglecting higher harmonics and constant terms (harmonics of order zero) might have an effect on the system dynamics, especially at high amplitude [57]. However, for the analysed cases, the considered harmonic components seemed sufficient. Besides, as mentioned above, the procedure does not directly study parametric resonances, information about the existence of parametric resonances is extrapolated by the instability of the regular resonance solution; clearly, a similar approach is not necessarily correct. Furthermore, the procedure studies only periodic solutions, while nonlinear systems, and especially NESs, typically exhibit quasiperiodic motions, relaxation oscillations [72] or even chaotic motions [73], completely overlooked by the HBM adopted. Additionally, the procedure is unable to identify coexisting branches of periodic solutions [57]; it would be possible to implement control schemes to partially analyse if more branches coexist; however, considering the already high computational cost of the analysis, we decided not to implement such techniques. Finally, the HBM is not effective for studying multifrequency sea spectra and cannot be used for the irregular wave simulations. Therefore, the usage of direct numerical simulations will also play a key role in the analysis where necessary. A code-to-code verification between the HBM and direct numerical simulation is performed in "Appendix B".

Direct numerical simulations
The direct numerical simulations involve integrating the governing system of equations using the ode45 solver in MATLAB. Each simulation is 10,000 s in length; however, if the pitch angle exceeds 90 degrees, the simulation is aborted. At such extreme angles, the validity of the model is diminished and predictions of the cylinder continuing to oscillate past its horizontal orientation are unrealistic. Thus, for these cases, the simulation is terminated and flagged as a capsizing event due to parametric resonance.
The same method to identify parametric resonance in regular waves as used in [45] is employed here. For normal operating conditions, in single-frequency regular waves, the frequency of the pitch motion and the input waves is the same. However, when parametric resonance occurs, the frequency of the pitch motion becomes half of the input wave frequency. Therefore, the occurrence of parametric resonance in regular waves can be easily identified during post processing by comparing the frequency of the pitch motion to the input wave frequency. If the peak amplitude of the power spectral density (PSD) occurs at a frequency less than 75% of the input wave frequency, then it is clear that parametric resonance has occurred for that input wave condition.
For irregular waves a different method must be adopted to identify parametric resonance, as discussed in [45]. For an irregular wave spectrum, there is energy in the wave frequencies below the peak of the spectrum; thus, if there is pitch motion at half the peak frequency of the input waves, it cannot be instantly identified whether this pitch motion response is due to parametric excitation or direct excitation from the spectral wave components at that lower frequency. Similarly, there is energy at wave frequencies above the peak of the spectrum, so once again it cannot be instantly identified whether large pitch motion at low frequencies is due to direct excitation from the wave components at that frequency, or from parametric excitation caused by the wave components at twice that frequency. Therefore, a different metric is applied for the identification of parametric resonance in irregular waves, based on the energy increase in the pitch motion observed during the simulation. The energy in the pitch motion (cal- culated from the root-mean-square value of x 5 ), measured at the start of the simulation is compared to that at the end of the simulation. Based on the findings in [45], if the energy increase from the start to the end of the simulation is more than double, then the simulation is added to the parametric resonance region. Following [45], the energy is measured over a duration equal to 5T p , where T p is the peak period of the wave spectrum. The energy at the start of the simulation is measured from 5T p ≤ t ≤ 10T p , so that the motion of the buoy is not measured from a cold start and the energy at the end of the simulation is measured from 10,000 [s] − 5T p ≤ t ≤ 10,000 [s].

Results
The results of the comparative study are presented in the section, first considering the regular waves and then the irregular wave spectra. In Sect. 3.1, to create a baseline against which the performance of the different DVA configurations can be compared, the occurrence of parametric resonance for the floating cylinder without any DVA attached is analysed. In Sect. 3.2, the DVA parameters are optimized, considering a fixed mass ratio of ε = 0.05. The mass ratio is initially fixed, considering that typically the larger the mass ratio, the better a DVA performs, as also discussed later in Sect. 3.3. The performance of the different DVA configurations are then compared in Sect. 3.2.5. Next, in Sect. 3.3, the effect that varying the mass ratio has on the performance of the DVAs is analysed. Lastly, in Sect. 3.4, the performance of the DVAs in irregular wave spectra are analysed and compared.

Regular waves: baseline case with no vibration absorber
The baseline case, which considers the floating cylinder without any DVAs, is simulated using the HBM, and the results are shown in Fig. 4 which illustrates the frequency/amplitude range for which parametric resonances are triggered. The colours indicates the heave/pitch amplitude of the regular wave solutions. The thick black line marks the stability boundary of the regular solution, meaning that inside the boundary the system undergoes parametric resonances. Regular resonance solutions were generated through the HBM, while parametric resonances through direct numerical simulations. We remind that the HBM enables to detect parametric resonances directly from instabilities of regular resonances, without requiring to actually obtain the parametric resonance response. The figure shows that parametric resonances are generated already for a very low wave amplitude ( A * 0 = 0.011) when the wave frequency matches 2ω 5 . The frequency region for which parametric resonance occurs widens as the wave amplitude increases and overall 36 % of the area in the considered frequency/amplitude range is subject to parametric resonance (u * 0 = 0.36).

Fig. 5
Referring to the system with a TMD acting in heave, for each couple of absorber's natural frequency ratio γ a and damping ratio ζ a (ε = 0.05), a illustrates the results for the u * metric, while b shows the results for the A * metric. The red and blue dots mark the optimal absorber configuration according to the given criterion The optimization shall initially consider the mass ratio fixed at ε = 0.05. For each DVA configuration, the optimal values for the natural frequency ratios, γ a or κ a , and the damping ratio, ζ a , are determined from the brute-force approach of sweeping through the parameter values.

TMD acting in heave
The first configuration considered is the TMD acting in the heave direction. The result of the optimization is summarized in Fig. 5a, b, which represent the parametric resonance region ratio u * and the minimal wave amplitude for generating parametric resonance A * , respectively. The optimal conditions obtained with the first criterion are γ a = 0.9125 (k a = 4.44 × 10 5 N/m) and ζ a = 0.074 (b a = 376,100 Ns/m), while for the second criterion, they are γ a = 0.9125 (k a = 4.44×10 5 N/m) and ζ a = 0.0786 (b a = 388,000 Ns/m). Optimal conditions are marked by red and blue dots, respectively, in the figure. The optimal values are very similar in the two cases. The natural frequency ratio is the same, while damping ratio differs by only 6 %. In terms of performance, the TMD operating in heave seems extremely efficient. Utilizing a mass of only 5 %, it is able to reduce the parametric resonance area ratio u * from 0.36 to 0.015, while A * is increased from 0.011 to 0.44. Figure 6 shows the response of the system with the attached TMD optimally tuned according to Fig 5a. In the three subplots, the response of the system in heave, pitch and the relative displacement of the absorber are represented. In Fig. 6a, the two resonant peaks generated by the TMD are clearly recognizable. The unsta-ble area, which generates the parametric resonance, is limited to two small regions existing for high wave amplitudes. Interestingly, the right unstable area does not appear where the heave response is maximal, but on the left of the right peak. This is probably related to the frequency content, which should be as close as possible to 2ω 5 . The TMD manages to suppress parametric resonance instabilities by significantly reducing the heave oscillation amplitude for ω ≈ ω 3 . If the heave response is reduced, the term ρg A C L D x 3 x 5 /2, which is responsible for the generation of the parametric resonance, is also reduced, suppressing this kind of behaviour. The reduced amplitude of the heave response can be noticed comparing Fig. 6a with 4a. A comparison of Fig. 6b  and 4b shows that, apart from changing the stability properties of the response, the amplitude of oscillation in pitch (without considering the parametric resonance) is practically unvaried. This is expected, considering that the absorber operates in heave and it does not directly affect oscillations in pitch direction (except for nonlinear terms, which depends on x 5 and are practically negligible if x 5 is small). Figure 6c highlights an important technological aspect of the system. Namely, in the considered frequency range, the TMD has a relatively large vertical displacement, which poses a challenging aspect in the practical design of this system. In fact, although the considered mass is relatively small, the space required for its vertical displacement is not negligible. In Fig. 6, the parametric resonance solutions are obtained through direct numerical simulations, differently from the regular resonance solutions which are obtained through the HBM. We notice that large pitch values are exhibited only in the left parametric resonance region, which suggest that the TMD acting in heave might be even more efficient than what is predicted by the HBM.

TMD acting in pitch
The second configuration investigated is the TMD acting in pitch. The result of the optimization is shown in Fig. 7. For the minimal parametric resonance region criterion, the optimal value is u * = 0.28, provided by the parameter values: γ a = 0.4427 (k a = 8.418×10 8 Nm) and ζ a = 0.0487 (b a = 1.938 × 10 9 Nms). Instead, by maximizing the wave amplitude required for generating parametric resonance, the optimal parameter values are γ a = 0.5464 (k a = 1.2821 × 10 9 Nm) and ζ a = 0.024 (b a = 9.57 Nms), which provide A * = 0.086. The absorber in pitch seems to have performance significantly worse than the absorber in heave. The minimal achievable u * with the TMD in pitch is 0.28 and the maximal A * is 0.086, compared against u * = 0.015 and A * = 0.44 provided by the TMD in heave. Also, the optimal values for the TMD's natural frequency ratio γ a and damping ratio ζ a are significantly different for the two optimization criteria, which was not the case for the TMD in heave. The curvy trend of the optimal region in Fig. 7b cannot be easily interpreted. It suggests that if the natural frequency of the TMD is exactly as the pitch frequency (γ a = 0.5), absorber damping should be slightly higher than if γ a is slightly lager or smaller than 0.5. Figure 8 illustrates the response of the system with the TMD optimally tuned according to Fig. 7a. Comparing it with Fig. 4, we notice that the improvement of stability is limited to a relative thin strip for ω slightly smaller than 2ω 5 . Such performance is probably insufficient to justify the implementation of the TMD in pitch in any engineering system.

NES acting in pitch
The next configuration considered is the NES acting in the pitch DoF. For this case (and for the NES acting in heave), the HBM was ineffective. Often, during the computation, the continuation procedure failed and was unable to converge to any solution. A more in-depth analysis revealed that the failure of the continuation was due to the presence of folds in the frequency response function. The adopted continuation technique consists of discretizing the continuation parameter (usually the wave frequency), finding the solution for the first point by utilizing as initial guess the solution of the linearized problem, then, for each subsequent point, by utilizing the solution of the previous point as initial guess, finding the solution for the new point. This technique is extremely simple and has the advantage of providing the solution on a given grid of values in the desired variables space, for instance the frequency, amplitude space (ω, A). However, the presence of folds makes the branch of solutions turn back, which is something that the algorithm cannot follow. Furthermore, the solution of the system after the fold is not necessarily close the previous point, therefore the initial guess utilized might be unable to converge towards any solution. Therefore, the optimization procedure for the NES cases (pitch and Fig. 7 Referring to the system with a TMD acting in pitch, for each couple of absorber's natural frequency ratio γ a and damping ratio ζ a (ε = 0.05), a illustrates the results for the u * metric, while b shows the results for the A * metric. The red and blue dots mark the optimal absorber configuration according to the given criterion The results of the optimization procedure revealed that the NES acting in pitch is unable to reduce the occurrence of parametric resonance. Sweeping a wide range of κ a and ζ a values, spanning many orders of magnitude, the values of A * and u * remained unchanged from the baseline values.

NES acting in heave
The results of the optimization for the NES acting in heave is shown in Fig. 9. A very large range of spring and damper coefficient values were swept (note the log-scale on both axes in the contour plots). The results in Fig. 9a have a minimum value of u * = 0.30, revealing that the NES is unable to reduce the parametric resonance region much below the baseline value of u * = 0.36. In fact for many spring and damper values, the NES in heave actually increases the parametric resonance region. Figure 9b shows that the for 0.1 ≤ κ a ≤ 0.2, the NES is able to increase the normalized minimal wave amplitude generating instability to above A * = 0.05, which is an improvement on the baseline value of A * = 0.011 but well below the values achieved by the TMD acting in heave.
Further analysing the results raises an interesting question towards the metric used to evaluate the performance of the DVA: whether the DVA should be judged on its ability to suppress parametric resonance when it occurs or on its ability to mitigate the onset of parametric resonance before it occurs. This question arises after noting that in some instances, parametric resonance would appear at the start of the simulation and then be suppressed by the NES, such that the second half of the time simulation shows no unstable motions. An example is shown in Fig. 10a. The PSD for the entire time series is plotted in Fig. 10b, showing a much larger response around half the input wave frequency than at the input wave frequency (ω = 1.86ω 5 ). Therefore, the identification algorithm classes this case as a positive occurrence of parametric resonance, and the NES is deemed ineffective. On the other hand, if only the second half of the time series is considered then the corresponding PSD is plotted in the Fig. 10c, which shows a single peak that corresponds with the input wave frequency. The identification algorithm would class this case as being stable and that the NES was effective in suppressing the occurrence of parametric resonance. For many applications, it may not be acceptable to incur any large amplitude, parametrically induced oscillations before they are, i.e. parametric resonance should be mitigated before it occurs. However for the sake of completeness and out of interest a comparison of these parametric resonance suppression results against the original parametric resonance mitigation results are presented.
Reanalysing the results, considering only the second half of the time series is shown in Fig. 11. Here, it can be seen that the values of the u * and A * metrics are greatly improved in the region around ζ a ≈ 10 0 for very small spring coefficients, κ a ≤ 10 −4 . Compared to Fig. 9, the minimum u * value is improved from to 0.3 to 0.13 and the maximum A * value increased from 0.085 to 0.21. At such low spring coefficient values, the spring force is basically negligible and the NES is behaving similar to a Lancaster damper (which is a mass damper system without any spring).
To examine the cause of the improved performance for the Lancaster damper type NES configurations with very low spring coefficients, the simulation results for the parameter values producing the optimal A * result (marked by Point 1 on Fig. 11) are investigated in detail. Figure 12 plots the response of the system over the considered frequency-amplitude range and compares the difference between the identified regions of parametric resonance considering the entire time series (inside the red dashed line) and considering the second half of the time series (inside the black dashed line). The results show that when considering the second half of the time series only, parametric resonance is identified almost only when the system goes catastrophically unstable and the buoy capsizes (denoted by the white region in the graphs, where the simulation is aborted after the pitch angle exceeds 90 degrees). The region between the red and black dashed curves corresponds to cases Fig. 11 Considering the second half of the pitch displacement time series only, results of the NES acting in heave, for each couple of absorber's pseudo-natural frequency ratio κ a and damping ratio ζ a (ε = 0.05), a illustrates the results for the u * metric, while b shows the results for the A * metric where parametric resonance occurs at the start of the simulation but then the NES is able to suppress it.
The results in Fig. 12b show that Lancaster damper type NES requires very large displacements in order to suppress the occurrence of parametric resonance, with stroke distances 5 -25 times G M. Another interesting result can be observed when comparing the parametric resonance region for the Lancaster damper type NES against the baseline parametric resonance region for the case with no DVA, denoted by solid grey line in Fig. 12. It can be seen that although the NES significantly decreases the region of parametric resonance for frequencies higher than ω/2ω = 1, it actually slightly increases the parametric resonance region for lower frequencies. The reason for the very large displacements and for the shift of the parametric resonance region to lower frequencies is due to centripetal acceleration on the DVA mass, from the pitch motion, forcing the DVA mass to drift away from its equilibrium position at the buoy's CoM.
Centripetal acceleration on the DVA mass is represented by the term −m a3ẋ 2 5 x a3 in Eq. (13), which therefore acts as negative restoring force proportional to the pitch velocity squared. For the cases when the NES spring coefficient is very small (seen on Fig. 11, for κ a < 10 −4 ) the centripetal acceleration term dominates the spring restoring force, k nl x 3 a3 , and the DVA mass begins to drift from x a3 = 0. The effect of the drifting DVA position on suppressing parametric resonance can be gleaned from inspection of the inertia term for the pitch motion in Equation (13); I 5 + m a3 x 2 a3 ẍ 5 , which shows that the pitch inertia increases with the square of the DVA displacement. An increase in pitch inertia results in a decrease in the pitch natural frequency, which explains why the parametric resonance region is slightly shifted to lower frequencies in Fig. 12. This also explains why the large amplitude DVA displacements found between the red and black dashed lines in Fig. 12b corresponds to cases where parametric resonance occurs at the start of the simulation but then gets suppressed by the second half of the time series. The large amplitude parametric pitch motions result in a large centripetal force which causes the DVA mass to drift far from x a3 = 0 and therefore detunes the parametric resonant frequency from the frequency of the input waves. These effects are highlighted in Figs. 13 and 14, which plot the simulation results from Points 1 and 2 on Fig. 12, respectively. Figure 13 corresponds to Point 1, which is in between the red and black dashed curves in Fig. 12, where parametric resonance is identified when considering the whole pitch displacement time series but not when considering only the second half of the time series. The pitch displacement time series is shown in Fig. 13b, where at the start of the simulation an exponential increase in the pitch displacement oscillation amplitude occurs. The DVA displacement is plotted in Fig. 13c, showing a large displacement from equilibrium occurring when the pitch oscillations increase in amplitude. When the DVA displacement increases, the large parametric resonance pitch motions are seen to be suppressed and the frequency of the pitch oscillations is seen to decrease, highlighting the detuning effect the drifting DVA displacement has on the pitch natural frequency. This can also be seen in the PSD plot in Fig. 13d, corresponding to the Fourier transform of the pitch displacement for the entire time series, which shows a very large broadband region that starts around half the input wave frequency and then decreases to a main peak at a much lower frequency. By comparison, the PSD corresponding to only the second half of the pitch displacement time series, plotted in Fig. 13d, shows a single small amplitude, narrow banded, peak at the input wave frequency, indicating that the occurrence of parametric resonance has been suppressed by the DVA. The simulation results for Point 2, which is out of the parametric resonance region, are shown in Fig. 14, where the PSD for both the entire and the second half of the pitch displacement time series contains only a single peak at the input wave frequency. The effect of the centripetal acceleration causing the DVA mass to drift from its equilibrium can be seen in Fig. 14c.

Regular waves: DVA performance comparison
The performance of the different DVA configurations are directly compared in Fig. 15, which depicts the parametric resonance boundary computed for each of the DVA configurations, as well as for the baseline case with no DVA. Note that for the NES acting in heave, the two different optimal cases are shown: (1) considering the pitch displacement from the full simulation ( Fig. 9) and (2) only considering the second half of the simulation (Fig. 12) are shown. The figure summarizes the results presented in the previous subsections, namely: • the absorber in heave outperforms the absorber in pitch, and • the TMD outperforms the NES.
Accordingly, the TMD acting in heave is the most effective, almost completely eliminating the parametric resonance region in the area of interest. Conversely, the NES acting in pitch is the least effective, providing no distinguishable reduction in the parametric resonance region. Therefore, coupled with the fact that the NES can only be assessed utilizing the computationally expensive direct numerical simulations, the NES acting in pitch shall be omitted from the further investigations in the following sections.

Regular waves: variation of the mass ratio
Here, the effect of varying the DVA mass is examined. As mentioned in Sect. 2.4.2, the DVA mass is one of the most important parameters since the smaller the DVA mass the more feasible the device. The optimization procedure for the spring and damper parameters is performed utilizing various values of the mass ratio ε, and the results are shown in Fig. 16. The TMD acting in heave (blue line) is seen to completely eliminate the parametric resonance region once the mass ratio exceeds 0.06. By comparison, the NES acting in heave does not completely eliminate the parametric resonance region until the mass ratio approaches 0.2. (Note that the resolution of the considered mass ratios for the NES cases is more sparse than for the TMD cases, since the NES cases require more computational effort using the direct numerical simulations.) For the NES acting in heave, the two different optimal cases that consider the pitch displacement from either the full simulation Fig. 15 Comparison of the stable region provided by the optimal TMD acting in heave (blue lines), the optimal TMD acting in pitch (yellow line), the optimal NES acting in heave when considering the second half of the simulation only (red line), the optimal NES acting in heave when considering the full simulation (green line) and the optimal NES acting in pitch (hidden beneath the black line), with respect to the case of no absorber (black line); ε = 0.  Fig. 12 where the buoy capsizes and the simulation is aborted, which means that for mass ratio values less than 0.15 the NES acting in heave cannot stop the buoy from capsizing but so long as the simulation does not abort it can suppress the parametric resonance by the second half of the simulation (albeit with very large stroke displacements).
The TMD acting in pitch is unable to completely eliminate the parametric resonance region for the considered mass ratio values up to 0.2. Here, an important point should be highlighted: practical considerations of the mass ratio reveals that the TMD in pitch is unfeasible. The inertia in pitch is I 5 + m 5 = 1.85 × 10 12 kgm. Considering the design of the DVA in pitch shown in the schematic in Fig. 3, the upper limit for the DVA inertia is physically constrained by the radius of the cylindrical buoy, which is the maximum distance the DVA mass can be located from the centre of rotation. Therefore, the maximum inertia is I a = M a5 (D/2) 2 , which means that in the extreme upper limit case where the DVA mass equals in the entire mass of the buoy, Fig. 17 a Maximum pitch displacement for the range of considered irregular wave spectra. The parametric resonance region is denoted inside the white dashed line. b The same information as in a but presented using contour lines with a 0.01 rad step size. The parametric resonance region is denoted inside the black dashed line M a = M, then I a = 7.4 × 10 10 , which equates to a mass ratio value of ε = 0.04. So, the analysis performed in Sect. 3.2, which considered a mass ratio of ε = 0.05, is infeasible, requiring a DVA mass greater than the entire buoy mass. Therefore, it can be stated that the TMD acting in pitch is significantly inferior to the DVAs acting in heave. As such, the TMD acting in pitch will be omitted from further investigations. Further explanation of the superior performance of the DVA acting in heave compared to pitch is given in the comparative analysis in Sect. 3.5.

Irregular wave spectra
While the results from the regular waves showed the superior performance of the TMD over the NES, it is possible that the NES may possess some advantages over the TMD for multi-frequency wave spectra due to the intrinsic broadband nature of the NES. In this section the performance for the TMD and the NES acting in heave are compared considering external excitation from irregular wave spectra. A baseline case with no DVA is given in Sect. 3.4.1, and then the performance of the TMD and NES are compared in Sect. 3.4.2. Figure 17a plots the maximum pitch displacement for each sea state and identifies the parametric resonance region within the white dashed line. Figure 17b shows the same information, but uses contour lines (each contour line represents 0.01 rad) in order to better highlight the results. In particular, the synchronous location of the identified parametric resonance region and the area where the density of the contour lines increases, provides confidence in the procedure employed to identify parametric resonance in irregular wave spectra.

Baseline case with no vibration absorber
Regarding the baseline values for the performance criteria, the parametric resonance region ratio 1 u * 0 = 0.417 and the minimal normalized wave amplitude to generate parametric resonance A * 0 = 0.07.

Comparison of TMD and NES acting in heave
The performance of the DVAs will be optimized and assessed considering ε = 0.05, which is the same mass ratio value used for the regular waves in Sect. 3.2.
The results for the TMD acting in heave are shown in Fig. 18, revealing the TMD is able to completely eliminate the occurrence of parametric resonance in the area of interest for a relatively wide range of spring and damper values. In comparison, Fig. 19 shows that the NES acting in heave is unable to completely eliminate the occurrence of parametric resonance, thus the TMD has displayed the superior ability to suppress the occurrence of parametric resonance in both regular and irregular waves. Further analysis and comparison of the results is given in Fig. 20, which shows contour plots of the maximum pitch displacement over the range of considered sea states for the optimal DVA parameter combinations. Figure 20a plots the TMD results, using the parameters from point 1 in Fig. 18, where the maximum pitch displacement is almost entirely below 0.03 rads. The maximum pitch displacement values over the entire frequency-amplitude space in Fig. 20a are comparable to the values for the baseline case that fall outside of the parametric resonance region in Fig. 17.
For the NES, two optimal values are considered and shown in Fig. 20b, c, corresponding to points 1 and 2 Fig. 18 Results of the TMD acting in heave, for each couple of absorber's natural frequency ratio γ a and damping ratio ζ a (ε = 0.05), a illustrates the ratio of unstable area in the frequency/amplitude space, while b indicates the minimal normalize amplitude at which instability occurs Fig. 19 Results of the NES acting in heave, for each combination of pseudo-natural frequency ratio κ a and damping ratio ζ a (ε = 0.05), a illustrates the ratio of unstable area in the frequency/amplitude space, while b indicates the minimal normalize amplitude at which instability occurs in Fig. 19, respectively. Similar to the case for the regular waves, the NES has an optimal performance region around κ a = 0.1 and then another region for very low spring coefficient values, κ a ≤ 0.001, which behaves similar to a Lancaster damper. Although the parametric resonance region ratio value is similar for the two NES configurations, Fig. 20b shows that the Lancaster damper type configuration, increases the region, relative to the baseline case, where the pitch motion is catastrophically unstable and the buoy capsizes.
Additional comparison between the DVA performances in irregular waves is given in Fig. 21, which shows the maximum DVA displacements for the corresponding simulations in Fig. 20. Here, it can be seen the Lancaster damper type configuration of the NES with a very low spring value, Fig. 21b, requires extremely large, infeasible stroke displacements, similar to the case for the regular waves, whereas the other NES configuration in Fig. 21c requires much smaller stroke displacements. Comparing the maximum DVA displacements for the TMD in Fig. 21a with the results for the regular waves in Fig. 6c shows slightly smaller DVA displacements in the irregular waves; however, the TMD stroke displacement required for the TMD in heave are still quite large.

Comparative analysis
The results from the comparative analysis clearly show that DVAs acting in heave outperform the DVAs acting in pitch. This agrees with the findings in [16,17] that reducing the heave oscillations is the most beneficial means to mitigate parametric resonance in pitch/roll. Parametric resonances are a sort of selfexcited oscillations, generated by a harmonically varying stiffness term, which for the system under study is ρg A C L D x 3 x 5 /2. The DVA acting in heave directly reduces the heave oscillation amplitude x 3 acting against forced vibrations, therefore reducing the amplitude of the varying stiffness term. On the contrary, the TMD acting in pitch has practically no effect on heave oscillations, but it increases the damping of pitch oscillations, thus preventing, to some extent, the onset of parametric resonances.
Another significant benefit of applying to the DVA to the heave rather than pitch DoF is that DVA acting in heave would also be able to suppress parametric resonance in the roll DoF. Although this study only considered a 2 DoF system, heave and pitch, parametric resonance is activated in roll by the same means as in pitch, therefore the DVA acting in heave can suppress both Fig. 20 Comparison of the maximum pitch displacement over the a TMD acting in heave using the parameters from point 1 in Fig. 18, b NES acting in heave using the parameters from point 1 in Fig. 19 and c NES acting in heave using the parameters from point 2 in Fig. 19. The contour lines represent a 0.01 rad increase in max pitch displacement Fig. 21 Comparison of the maximum DVA stroke displacement. a The TMD acting in heave using the parameters from point 1 in Fig. 18, b the NES acting in heave using the parameters from point 1 in Fig. 19 and c the NES acting in heave using the parameters from point 2 in Fig. 19 parametric resonance in pitch and roll. By comparison, the DVA in pitch would have no effect on the roll DoF; therefore, a second DVA would need to be implemented in roll. Thus, the single DVA acting in heave is much more efficient and practical than requiring two separate DVAs acting in the rotational DoFs. Similarly, although out of the scope of the present paper, the implementation of DVAs acting in both the heave and pitch DoFs at the same time can be considered, in order to investigate if utilizing multiple DVAs can enhance the suppression of parametric resonance enough to justify the added complexity to the system.
The results from the comparative analysis also demonstrate that the TMDs outperformed the NESs. This is probably due to the nature of parametric res-onances, which are related to a relatively narrow frequency band, despite waves having a broad frequency content. Therefore, the TMD can work efficiently for the whole frequency range of interest. Indeed, previous studies revealed that the TMD frequency band of operation can be enlarged if its damping and natural frequency are properly tuned [74]. This observation practically eliminates the advantage of the NES, of being able to resonate in a broad frequency range. Additionally, the NES requires a minimum energy threshold in the primary system to be activated, probably limiting its effectiveness. Conversely, the TMD, being a linear device, is not affected by the energy level, at least within the limits of validity of the mechanical model adopted. TMD displacement poses a practical limit on its feasibility. In the irregular wave spectra, the results in Fig. 21a showed the maximum TMD displacement values approximately increase linearly with the significant wave amplitude, requiring an 8-9 times greater stroke displacement than the significant wave height. For the regular wave case, in Fig. 6c, the TMD displacement is over 10 times greater than the wave amplitude. In a realistic system, it would be necessary to constrain the physical stroke length of the DVA system, resulting in end-stop collisions whose effects should be added to the model and assessed in future work, following methods such as those used for the vibro-impact WEC [75].
An interesting result emerging from the wide sweep of the DVA parameter values was the identification of a Lancaster damper-type system, which is a mass and damper without a spring, being a theoretically good means of suppressing parametric resonance after it occurs (but not at mitigating it before it starts). The Lancaster damper provides a sort of self-detuning mechanism, where the pitch rotations cause a centrifugal force on the DVA mass, which accelerates the DVA mass away from its resting position. The displaced DVA mass alters the mass distributions within the cylinder, thereby changing the moment of inertia which shifts the pitch natural frequency and detunes the parametric resonance. Although this is a theoretically interesting application, the results showed that such a system requires unrealistically large stroke displacements thus is not practically feasible.

Conclusion
The ability of different DVA configurations to suppress the occurrence of heave induced parametric pitch motion in a floating cylinder was assessed via simulation using a 2 DoF coupled heave-pitch numerical model. Two prominent DVA methods, the TMD and the NES, were considered, in combination with the DoF in which they acted, either on the heave motion or the pitch motion. The results demonstrated that suppressing parametric resonance in the pitch DoF is far more efficiently achieved by applying a DVA to the heave DoF than applying a DVA to the pitch DoF. In addition, the TMD demonstrated superior performance to the NES in both DoFs and in single frequency regular waves as well as multi-frequency irregular waves. Therefore, the TMD acting in the heave DoF was most effective at suppressing the occurrence of parametric resonance, however the maximum displacement/stroke length required by the TMD acting in heave became quite large as the wave amplitude increased, therefore further work is required to understand its performance in stroke-limited situations involving end-stop collisions.
We then define the potential energy function of the system, including the gravitational potential energy of the absorber mass and the elastic potential energy. The gravitational potential energy of the buoy is neglected, since it is compensated by hydrodynamic forces. We have +m a3 g (x 3 + (x a3 + l 0 ) cos(x 5 )) .
We then choose k 2 = 3k 3 l 0 , such that the absorber's spring restoring force has no quadratic component.
If the absorber under consideration is a TMD, then we should have k 2 = k 3 = 0 and m a3 g = k 1 l 0 , such that k nl3 = 0. Conversely, in the case of an NES, we should impose k 1 = 3k 3 l 2 0 , such that k a3 = 0. We notice that realizing a spring which satisfies the requirements of an NES is very challenging. However, the practical realization of such a system is out of the scope of this research.

B Verification of harmonic balance method against direct numerical simulation
The results for the code-to-code verification between the HBM and direct numerical simulations is shown in Fig. 22, which considers the baseline case with no attached DVA presented in Sect. 3.1. The results look qualitatively the same and quantitatively the values for the two performance criteria metrics are also the same. The minimum normalized wave amplitude which generates parametric resonance for the direct numerical simulations is: A * = 0.011, and the parametric resonance region ratio is: u * = 0.36, which agree with the values obtained by the HBM in Sect. 3.1.