Coexisting attractors in floating body dynamics undergoing parametric resonance

This study pertains to analysing the dynamical behaviour of a floating body undergoing parametric resonances. A simple vertical cylinder, representing a classical spar-buoy, is considered, limiting its motion to heave and pitch degrees of freedom. Its geometry and mass distribution are chosen such that a 2:1 ratio of heave to pitch/roll natural frequency makes the spar-buoy prone to parametric resonance. The system is then studied by the shooting method, combined with a pseudo-arclength continuation, and the harmonic balance procedure. Results show that an extensive bistable region exists, where stable parametric resonance coexists with a regular resonance response. The analysis also unveiled the existence of stable quasiperiodic motions existing in correspondence of both pitch and heave resonance. Results are qualitatively validated using a model based on the explicit nonlinear Froude–Krylov force calculation.


Introduction
Parametric pitch and roll motions can occur in a floating body due to the time variation of the hydrostatic restoring torque coefficients, as wave-structure interaction causes the body submergence level to oscillate. This phenomenon was first documented by Froude in the nineteenth century [9] and has since been observed and studied in a wide range of naval and offshore structures such as ships [28,35,40], spar-buoys [12,22,23,32], floating platforms [37,38,43], and marine renewable energy devices [5,39]. The occurrence of parametric resonance in these structures can pose significant safety risks and performance degradation. For example, container ships lost cargo overboard due to large parametric roll motions [8] and wave energy converters incurring reductions in power output [2,3,17].
As such, methods for suppression and stabilization of parametric pitch/roll in ships have been investigated, see for example [10,11,19,41] and the references therein. A key requirement is the ability to predict the occurrence of parametric resonance, a task that relies heavily on numerical modelling and analysis. While a rich catalog of effective methods exists within the field of nonlinear dynamics for modelling and analysing parametric resonance, relatively few of these methods have been applied to the study of parametric resonance in offshore vessels and structures. In this work, we adopt numerical techniques common in the nonlinear dynamics scientific community but rarely implemented in ocean engineering, aiming at revealing dynamical phenomena hidden in the complex structure of parametric resonances that these kinds of systems experience.

Test case
The floating body used in this study is a simple vertical cylinder, representing a classic spar-buoy, whose particular geometry and mass distribution follow a widely used test case from the literature [4,12,13,20,[24][25][26][27]34,44]. This test case spar-buoy is popular for studies involving parametric resonance, since it has a pitch/roll natural frequency (ω 5 = 0.108 rad/s) which is approximately half the heave natural frequency (ω 3 =0.216 rad/s). The 2:1 ratio of heave to pitch/roll natural frequency makes this spar-buoy prone to parametric excitation in pitch/roll.
Hong et al. [44] performed physical experiments of the spar-buoy at model scale, considering free decay and regular wave inputs. Several studies then utilize the information from Hong et al. in the development of numerical models of the spar-buoy, which are then used to examine the combination resonance response in Jingrui et al. [20], for the study of nonlinear vibration modes and instability in Gavasoni et al. [12], to evaluate the nonlinear dynamical behaviour of the spar in Liu et al. [27], to explore parametric resonance considering wave and vortex exciting loads in Li et al. [24][25][26], to highlight the existence of yaw instability when considering nonlinear kinematics within the numerical model in Giorgi et al. [13], and as a test case for a real-time detection system of parametric resonance in Davidson and Kalmár-Nagy [4].

Outline and objective of the paper
The objective of this study is to analyse the system's response in correspondence with the generation of parametric resonances. In particular, we aim to define the extent of the parametric resonance region by adopting numerical continuation and the harmonic balance (HB) technique. These methods should allow us to define the actual extent of parametric resonant solutions, including bistable regions, which are often overlooked by direct numerical simulations. A partial analysis of the dynamical behaviour of a similar system was already illustrated in [25,26], where bistable regions were first detected exploiting the method of multiple scales [30]. In the present study, different techniques are implemented, complementing the work done in [25,26]. Additionally, stable quasiperiodic motions are disclosed, and the dynamical integrity of coexisting solutions is investigated.
The remaining part of the paper is organized as follows. In Sect. 2, the mechanical model is introduced. In Sect. 3, direct numerical simulations, the shooting method, and the HB procedure are implemented to obtain a general understanding of the system dynamics in correspondence with the onset of parametric resonances, considering constant wave amplitudes. In Sect. 4, the region of existence of parametric resonances is identified in the frequency-amplitude space, its robustness is studied, and coexisting quasiperiodic solutions are also detected. In Sect. 5, the obtained results are qualitatively compared with a high fidelity numerical model of the same mechanical system. Finally, in Sect. 6, concluding remarks are provided.

Mechanical model
The hydrodynamic model for the test case spar-buoy follows the validated model of the same spar-buoy from [4,12,20,44]. The model captures the occurrence of parametric resonance by the 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. Normalised frequency 5 . Frequency is normalized with respect to ω 5 The inertial forces are given by the product of the inertia and acceleration: where M is the cylinder mass, m 3 is the hydrodynamic added mass in heave, I 5 is the pitch moment of inertia, and m 5 is the hydrodynamic added moment of inertia in pitch. Out-of-diagonal elements of the added mass matrix are neglected. The hydrodynamic restoring force and moment terms have time-varying parameters, which enables the existence of parametric and autoparametric resonance in the numerical model. In [12,20], the hydrodynamic restoring force and moment are given by where ρ is the water density, g is the gravitational constant, A C is the spar-buoy horizontal cross-sectional area, L SC is the distance from the still water level to the spar-buoy centre of mass, L D is the spar-buoy draft, G M is the spar-buoy metacentric height, and η(t) is the wave elevation. Following [12,20], the hydrodynamic damping is represented as where C i are the damping coefficients identified using measured data from the free decay experiments in [44].
For an input wave with amplitude A, frequency ω and phase φ, the wave elevation is The excitation force is given by where H E,i (ω) and φ E,i (ω) are the hydrodynamic excitation force coefficients for the amplitude gain and phase shift, respectively, which are calculated using the open-source linear potential flow boundary element method software Nemoh [2], and plotted in Fig. 1.
The values for the parameters used in the model are listed in Table 1 following the values provided in [12,20,44] for the validated model of the same device.

Preliminary analysis of the system dynamics
The system's dynamics is described by the system of ordinary differential equations provided in Eq. (1). The strongly nonlinear nature of the system makes it challenging to thoroughly study its potentially rich dynamics. Various strategies can be adopted for this purpose. In the following, we illustrate the different outcomes that various numerical techniques provide.

Direct numerical simulations
In ocean engineering, it is common to study the dynamics of a system by performing direct numerical integration of the equations of motion. Typically, to reduce the effect of transient dynamics, the wave forcing signal is introduced by a linear ramp [14]. In this way, the system, starting from the trivial equilibrium position, is gradually excited. A drawback of this approach is that, for each set of parameter values, the system always converges to the same attractor. If other attractors coexist with the identified one, they are overlooked. We apply this method on the system under study, performing simulations for a fixed wave amplitude of 1 m (A/G M = 0.0954), over a range of frequency covering both natural resonance frequencies of the system. We then collect the maximal amplitude of the steady-state oscillations for each frequency, obtaining the frequency response function (FRF) of the system. The result is illustrated in Fig. 2a and b.
From the Figure, we notice the detached branch existing for ω/(2ω 5 ) ∈ (0.951, 1.036). This is a branch of parametric resonances, which has a kind of V-shape, both in heave and in pitch. This parametric resonance, which corresponds to a so-called parametric instability of the regular resonance, is due to the terms η(t)x 5 and 1/2x 3 x 5 in Eq. (4). Both terms contribute to generate an instability in pitch, similarly to a classical Mathieu's equation. η(t)x 5 is properly a Mathieu term generating a parametric instability, while 1/2x 3 x 5 generate a so-called autoparametric instability. The contribution of 1/2x 3 x 5 is much more important than that of η(t)x 5 because, for ω ≈ ω 3 , x 3 is much larger than η(t). The V-shape of the parametric resonance curve is due to the quadratic terms x 2 5 in Eq. (3) and 1/2x 3 x 5 in Eq. (4). From this preliminary analysis, it seems that parametric instabilities are limited to the frequency range ω/(2ω 5 ) ∈ (0.951, 1.036) (for A = 1 m).
For ω ≈ ω 5 , although the system has only one natural frequency, the FRF exhibits two nearby peaks. The reason for this shape is explained in the next Subsection.

Numerical frequency sweep
In order to identify potentially coexisting solutions and to correctly evaluate the extent of the region presenting parametric resonances, we numerically simulate the dynamics of the system performing a so-called frequency sweep. Starting from a certain frequency value, a numerical simulation of the system is computed; the computation is interrupted when the system converges to a periodic solution. Then, the frequency is slightly increased (or decreased), and another numerical simulation is performed, utilizing as initial condition for the simulation the final state of the previous simulation. No ramp is utilized on the forcing signal. If the frequency step is small enough, and if the total simulation time of each simulation is a multiple of the excitation period, the system will interpret the transition from one simulation to the following one as a smooth change. Consequently, the system will follow a branch of solutions as far as they are stable. Performing a sweep-up (increasing frequency) and a sweep-down (decreasing frequency), some of the coexisting solutions can be found. Another advantage of this approach is that each simulation starts from a condition already very close to its steady state (except when jumps from one branch to another one occur); therefore, the simulation time can be significantly reduced. This technique is very well-known and widely utilized in the nonlinear dynamics scientific community [29]; however, it is very rarely implemented in ocean engineering [21,31].
The results of the analysis are shown in Fig. 2c and d. In the Figures, blue arrows indicate the jumps from one branch to the other one occurring during the sweep-up, while the red arrows refer to the jumps of the sweep-down. Jumps occur either because the solutions of the branch become unstable, or because the branch presents a fold. The Figures are almost identical to Fig. 2a and b; however, the branch of parametric resonance solutions is wider. It exists for ω/(2ω 5 ) ∈ (0.924, 1.064). With respect to the frequency band exhibiting parametric resonance previously measured, this one is 62.8 % larger. Besides, the branch of regular solutions exists also in regions where parametric resonances occur, which was not the case in Fig. 2a and b.

Shooting method
In order to have a better insight into the parametric resonance, we tried to track the branch of periodic solutions by the shooting method, combined with a pseudo-arclength continuation [29,33]. The shooting method consists in identifying initial conditions of the system, such that after one period of excitation (or two periods, for parametric resonances, which present a period double than the excitation) the system reaches again the initial state. This condition can be verified by direct numerical simulations, and, if it is satisfied, a periodic solution is identified.
In a mathematical sense, it corresponds to satisfy the equation where η = x(0). Several approaches can be implemented for finding the correct initial conditions η. Assuming that we have an initial guess η 0 of the correct η value, we aim at identifying δη such that Exploiting a Newton-Raphson approach for solving Eq. (9), η can be found iteratively. The shooting method can also provide unstable periodic solutions, although in this case numerical issues might arise. A by-product of the algorithm is the so-called monodromy matrix, whose eigenvalues can be used to study the stability of the periodic solutions and identify the bifurcations occurring at the stability loss [33]. The shooting method is fully exploited only if it is combined with a continuation technique. In fact, the solution obtained for a given set of system parameter values can be used to identify a proper initial guess for a close-by set of parameter values. This approach enables one to obtain branches of periodic solutions for a varying parameter and significantly speeds up the convergence of the algorithm. For our system, we implemented the pseudo-arclength continuation, which enables to follow branches of periodic solutions also if foldings occur. Details of the algorithm can be found in [33]. Also this technique is commonly implemented in the nonlinear dynamics community, but it is rarely adopted in ocean engineering for solving hydrodynamics problems.
In order to simplify the implementation of the pseudo-arclength continuation, it is convenient to have analytical expressions of the equations of motion. Therefore, we perform a seventh-order polynomial regression of the force coefficients, which are otherwise provided in the form of a Table as functions of ω.
Results of the procedure, for A = 1 m, are illustrated in Fig. 2e and f. Apart from the branches of stable solutions, already identified by the frequency sweep, the analysis also provided some branches of unstable solutions. Regarding regular solutions, the unstable peak of the heave resonance is clearly identified. The left unstable part of the parametric resonance solution's branch, which merges with the branch of regular solutions, is also generated by the procedure. An analysis of the eigenvalues of the monodromy matrix, in correspondence of the loss of stability, illustrates that one eigenvalue leaves the unit circle of the complex plane through -1. This event corresponds to the occurrence of a period doubling bifurcation. In fact, the parametric resonance solutions have a double period with respect to regular resonances. The fact that the periodic solutions of the branch emerging from the bifurcation are unstable, and that they coexist with the stable regular solutions, indicates that the bifurcation is subcritical.
The more complete picture of the branch of parametric resonance solutions provided by the shooting continuation shows that its apparent V-shape is rather given by two peaks with a valley in between. This shape resembles the FRF of systems with internal resonances, such as a single DoF system with an attached dynamic vibration absorber. If the two subsystems (primary single DoF system and dynamic vibration absorber) have the same natural frequency, the frequency response presents a minimum in correspondence of the natural frequency and two peaks at its sides [6]. The local minimum is due to an anti-resonance. In the system under study, the modal interaction is related to an autoparametric resonance in the nonlinear regime. A similar scenario was disclosed by Li et al. [26] for a similar system, by adopting the method of multiple scales.
We notice that also the pitch resonance has a similar trend. In fact, instead of presenting a single peak, it has two peaks divided by a valley. These two close-by peaks are related to a 2:1 internal resonance generated by the quadratic term x 2 5 in the restoring force (Eq. (3)); in fact, quadratic terms generate asymmetry and double harmonics according to the trigonometric relation cos(ω) 2 = 1/2 + 1/2 cos(2ω). This constitutes a sort of forcing term in heave, and, since ω 3 ≈ 2ω 5 , it causes a transfer of energy from pitch to heave. Figure 3 illustrates the time series and the fast Fourier transform (FFT) of the steady-state response for ω = ω 5 = 0.1096 rad/s. The frequency content of the heave response confirms the internal resonance taking place. In fact, the heave response presents a large frequency peak at the double of the excitation frequency, close to ω 3 (Fig. 3b). An analogous dynamical scenario was recently disclosed for microresonators, sharing a similar modal ratio [1].
For comparison, in Fig. 4 the time series response and the FFT of the steady-state solution for ω = ω 3 = 0.2159 rad/s are represented. In this case, the parametric resonance triggers vibrations in pitch with a frequency half of the excitation frequency, close to ω 5 . Looking at the frequency response in Fig. 2e and f, we notice that the shooting method was unable to track the right branch of unstable parametric resonance periodic solutions. While trying to follow that branch, the continuation failed to converge, probably because of its strongly unstable character.

Harmonic balance method
In order to obtain an approximation of the right unstable branch of parametric resonance solutions, we applied the harmonic balance technique. Differently from the methods utilized so far, the HB does not provide an exact solution. It consists of approximating the steady-state solution of the system to a series of harmonic terms. First, we decide to limit the analysis to harmonic terms of order ω and ω/2, which are the minimal terms required for identifying parametric resonances. In practice, we approximate the solution of the system by Higher and lower harmonic terms were neglected. Also constant terms, which exist because the system encompasses even-order nonlinearities, are omitted by the approximation. We substitute Eq. (10) and its time derivatives into Eq. (1); then, we collect equal harmonic terms, reducing the system of ordinary differential equations to an 8- A Newton-Raphson approach is used for this purpose. By coupling the HB method to a pseudo-arclength continuation, we are able to follow branches of periodic solution for a varying parameter, as we did with the shooting method. With respect to the shooting method, the HB is significantly faster, since it does not require time integration, and it is able to track with the same efficiency stable and unstable solutions. The stability analysis of each solution is performed by identifying if the real part of at least one of the Floquet exponents is positive. For this purpose, Hill's method is adopted [18], following the steps outlined in [7,42], which, for the sake of brevity, are not reported here. The results of the computation are illustrated in Fig. 5a and b (black lines). In the same Figures, the red lines reproduce the result obtained by the shooting continuation for comparison. In this case, the procedure tracked the whole branch of parametric resonance solutions, including the right unstable branch.
Referring to the response in correspondence of the pitch resonance (ω ≈ ω 5 ), not surprisingly the HB approximation completely overlooked the peak of the heave response. This occurs because the approximation includes only the harmonic of excitation and its half, while the peak in heave, for ω ≈ ω 5 , has a dominant frequency which is the double of the excitation frequency (Fig. 3b).
Comparing the black and red lines in Fig. 5a and b, we also notice that there is a non-negligible quantitative difference between the two, concerning the parametric resonance branch. Nevertheless, despite the rough approximation, the procedure was able to provide an informative picture about the parametric resonance and its region of existence.
In order to improve the estimation and also attempting at obtaining the correct frequency response for the pitch resonance, we include in the HB additional harmonic and constant terms as well. That is, The system of algebraic equations generated by the HB procedure is now 14-dimensional. Results of the computation are illustrated in Fig. 5c and d (black lines). In this case, the HB procedure was able to correctly track the frequency response for ω ≈ ω 5 . Besides, the difference between the red and black lines is practically negligible along almost all the considered range of excitation frequency. The HB procedure requires to compute the various harmonic contributions of the steady-state solution, which can then be used to obtain additional insight regarding the system dynamics. For the case under study, they are depicted in Fig. 6. In particular, Fig. 6f confirms that the parametric resonance in pitch, for ω ≈ ω 3 , has almost no harmonic component of order ωt, but mainly of order ωt/2 (Fig. 6g). Figure 6d shows that the internal resonance for ω ≈ ω 5 generates harmonics of order 2ωt in heave. Constant terms c 30 and c 50 are present anytime there are large amplitude solutions (Fig. 6a and e).
The accuracy of the estimation and the rapidity of the procedure suggests that the HB method is a potentially valid tool for rapidly performing a parametric analysis of the system dynamics.

Region of existence of parametric resonance solutions
With the objective of identifying the parametric resonances' region of existence, we repeat the direct numerical simulations performed for obtaining Fig. 2c and d over a wide range of wave amplitude. The results of the analysis are illustrated in Fig. 7. Figure 7a and b depicts the stable branches of the heave and pitch responses, respectively. The amplitude of oscillation increases almost linearly with the wave amplitude for both heave and pitch. However, for wave amplitude larger than 2 m ( A/G M = 0.191), the maximal heave and pitch amplitudes remain constant. This kind of saturation of the amplitude is due to capsizing of the buoy, occurring approximately for x 3 /G M > 2.5. In this event, which corresponds to a physical capsizing of the mechanical system, the model fails, and simulations are interrupted. Clearly, no parametric resonance can be identified after capsize occurs.  (ω, A) space. According to the numerical simulations performed, the system converges towards regular solutions in the clear area. Conversely, in the black area the system experiences a parametric resonance. In the grey area the system can either display a regular solution, or can undergo a parametric resonance, depending on the initial conditions. We notice that the grey area in Fig. 7c has an extent equal to 68.8 % of the black area, which makes it particularly relevant from an engineering perspective. The red lines in Fig. 7c indicate the boundary of the three regions obtained by the HB method. The two methods provide results in almost perfect agreement regarding the boundary of the black region, while the boundaries of the grey region present a slight mismatch, especially at high amplitude. The mismatch is related to the increased relevance of harmonics neglected in the HB approach.
For A ≈ 2 m (A/G M = 0.191), the right boundary of the grey area has a non-smooth deviation. This trend can be explained by analysing the frequency response function around those values. Figure 8 illustrates that, for an amplitude between 1.8 and 1.9 m, the branch of parametric solutions experiences a fold, which results in a non-smoothness in the continuation performed for increasing ω values (curves obtained by the HB method). In fact, once the fold is reached, the system necessarily leaves the branch it was following. Probably, the transient dynamics caused by the system leaving the branch makes it reach conditions causing capsizing. Numerical results in Fig. 7 show that the change of trend occurs for wave amplitude of 2 m, while Fig. 8 suggests that the fold appears for A between 1.8 and 1.9 m. The difference between the two is probably related to the approximation of the HB solution. In the case of more coexisting stable solutions, the system converges towards one of the attractors, depending on its initial conditions. Since the system is 5-dimensional (including the phase of the excitation), it is computationally very demanding to identify the basin of attraction of the system, even for a single parameter configuration. We, therefore, limited the analysis to sections of it, for a fixed value of the wave amplitude (A = 1 m) and for various frequencies of excitation. The result of the analysis is depicted in Fig. 9.
The Figure illustrates sections of the system's basin of attractions for x 3 (0) =ẋ 3 (0) = 0. The force is provided with no initial ramp. The brown area indicates initial conditions leading to regular solutions, while the white area marks the basin of attraction of the parametric resonance solution. The system tends to converge towards regular solutions for small initial pitch values (initial heave values are always zero for the computation). Conversely, to trigger parametric resonances, one of the variables must be significantly different from zero. This phenomenon is consistent with the fact that direct numerical simulations, performed with the initial ramp on the wave amplitude, completely overlooked parametric resonances in the bistable ranges ( Fig. 2a and b).
For A = 1 m, the system is bistable for ω/(2ω 5 ) ∈ (0.924, 0.951) and ω/(2ω 5 ) ∈ (1.037, 1.064). Figure 9a and b refers to the first interval; in the plots ω/(2ω 5 ) = 0.93 and 0.94, respectively. While for ω/(2ω 5 ) = 0.93 (Fig. 9a) the regular solution seems relatively robust, increasing the frequency ratio to 0.94 (Fig. 9b), the white region gets much closer to the origin of the axis. These are two steps of the transition that goes from only regular solutions for ω/(2ω 5 ) < 0.924, where the white area would disappear, to only parametric resonances for ω/(2ω 5 ) > 0.951, where there would be no brown area at all. The black area, visible near the corners of all the plots in Fig. 9, indicates initial conditions which caused the capsizing of the system. Referring to the second interval with bistable solutions, Figs. 9c (ω/(2ω 5 ) = 1.05) and 9d (ω/(2ω 5 ) = 1.06) show how the basin of attraction of the parametric resonance recedes as ω increases.  Fig. 10 correspond to the so-called global integrity measure (GIM) [36], which is the hyper-volume of the basin of attraction of a given solution. The trend of the curves in the Figure shows that the mutual erosion of the basin of attraction is gradual, and the variation of the global integrity measure is practically linear in the region of bistability. Figure 11 depicts the frequency response of the system for various values of the wave amplitude, from 0.5 to 4 m (obtained by the HB method). From the Figure, we can recognize that the left part of the branch of parametric resonance solutions has always much smaller amplitude than the right part of the branch; both regarding heave and pitch. For excessively high oscillation amplitudes, even small perturbations can make the system undergo capsizing. Therefore, even if solutions are theoretically stable, their practical relevance is questionable. For large wave amplitudes, we also notice that some parts of the branches related to parametric resonances become unstable even in the absence of folds. These instabilities are related to Neimark-Sacker bifurcations, which generate quasiperiodic motions; quasiperiodic solutions cannot be captured by the HB procedure adopted. Most of these instabilities occur for very large values of the response, where capsizing often occurs; therefore, it is challenging to identify them even by numerical simulations. For the same reason, they are probably practically irrelevant, and they exist in a phase space region where the validity of the model is arguable.

Quasiperiodic motions
One of these unstable branches exists for ω/(2ω 5 ) ≈ 0.95 and A = 4 m where the oscillation amplitudes are not very large. Therefore, with the help of numerical simulations, we are able to track these quasiperiodic motions. A close-up of the FRF around the second resonant peak is illustrated in Fig. 12a and b. The red curve in the Figures, obtained by direct numerical simulations, marks the maximal amplitude of the quasiperiodic solutions. The rest of the curve was obtained by the HB method. Figure 12c and d represents the phase portrait of the solution existing for ω/(2ω 5 ) = 0.955. The depicted curve seems to move on a torus, which confirms its quasiperiodic nature. The red lines in Fig. 12c and d are stroboscopic Poincaré sections of the solution, obtained by picking one point for each excitation period. Their closed shape further proves that the motion is quasiperiodic. In pitch (Fig. 12d), the red curve is broken into two parts because of the doubled period of oscillation. Finally, Fig. 12e and f illustrates the FFT of the response, highlighting the multiple frequency components of the motion.
The yellow curve in Fig. 11 shows that also in correspondence of the pitch resonance the frequency response presents an instability. This can be better recognized in the enlargement of the first resonant peak provided in Fig. 13a and b.
Also this instability is generated by a couple of Neimark-Sacker bifurcations, which, therefore, generate quasiperiodic motions. The same analysis performed regarding the branch of quasiperiodic motions at heave resonance is repeated for this case, around pitch resonance. The red curves in Fig. 13a and b, obtained by direct numerical simulations, mark the maximal amplitude of the quasiperiodic solutions in that region. These curves complete the FRF obtained by HB and represented in black. Figure 13c and d represents the phase portrait of the solution existing for ω/(2ω 5 ) = 0.4926. In heave (Fig. 13c), the torus makes a double loop in one period because it is dominated by a frequency double than the excitation frequency. Poincaré sections (red

Validation of the observed phenomena against high fidelity numerical simulations
Although the model implemented was already validated in previous studies [44], in this Section we verify if the most important findings of this research are confirmed by more sophisticated numerical models describing the same physical system. For this purpose, we implement a more representative model of the system, based on the explicit nonlinear Froude-Krylov (NLFK) force calculation: both the static and dynamic pressure fields of the undisturbed incident wave are integrated over the instantaneous wetted surface of the floater. This approach, available for both axisymmetric and prismatic floaters [16], has been validated against experimental tests in [14], and has already been effectively used to investigate parametric resonance for a floating wave energy converter in [15]. Moreover, in [13], this model has been applied to the same cylindrical floater considered in this paper, highlighting parametric resonance regardless of the simple archetypal geometry. In this paper, in order to ensure the most meaningful comparison, a 2-DoF NLFK model (heave and pitch) is considered, with parameters and structure equal to the model in Sect. 2, apart from the nonlinear description of Froude-Krylov forces. We excite the system by a linear chirp waveform of constant amplitude (A = 1 m). First, we consider increasing, then decreasing, frequencies. In this way, we can provide a situation similar to the frequency sweep illustrated in Fig. 2c and d, except that the excitation frequency is continuously increased, without waiting for the system to reach its steady state. The frequency ratio ranges from 0.8 to 1.2, which covers the whole bistable region. The response of the system to the excitation is depicted in Fig. 14 for the two cases. Figure 14a and b refers to the increasing frequency. We notice that pitch is suddenly excited for ω/(2ω 5 ) ≈ 0.98, marking the parametric resonances' triggering. The system persists in its parametric resonance state until ω/(2ω 5 ) = 1.142.
Looking now at the decreasing frequency case (Fig. 14c and d), we notice that parametric resonances are triggered for a frequency ratio of 1.035. The system converges again towards regular solutions for ω/(2ω 5 ) ≈ 0.945. Since the system takes some time to settle on the new solution at each transition, it is not easy to mark the transition frequency precisely. Nevertheless, the analysis confirms that a frequency range presenting both stable solutions exists. Additionally, the envelopes of the time series confirm the double peak shape of the frequency response in parametric resonance, both in heave and in pitch, as exhibited by the simplified model.

Conclusions
The dynamical behaviour of a pitch-heave spar-buoy was investigated. The floating body, having a heavepitch natural frequency ratio of 2:1, undergoes parametric resonances when the wave frequency is close to the heave natural frequency. The obtained results showed that a bistable region exists where either regular or parametric resonances can occur, and its extent is significant. This result is in contrast with numerical observations relative to similar systems [15]. Conversely, it confirms recent analytical studies about analogous mechanical models [25]. Interestingly, the overall dynamical scenario obtained closely resembles the internal resonance of a nonlinear asymmetric microbeam resonator [1]. Furthermore, the performed analysis unveiled the existence of stable quasiperiodic motions. However, their practical relevance for floating bodies is probably marginal.
In a more general framework, this study illustrated how relatively simple numerical or semi-analytical techniques can be implemented to analyse these kinds of systems. Although this might sound obvious for most scientists dealing with nonlinear dynamics, they are rarely adopted in ocean engineering. The set of analytical and numerical methods utilized in this paper partially complements those implemented in [25].
The results were validated qualitatively using a model based on the explicit nonlinear Froude-Krylov force calculation. Future developments of this work should comprise a thorough analysis of all the dynamic phenomena unveiled, utilizing this more accurate model. Besides, the effect of including other degrees of freedom, such as roll and sway, should also be investigated.