Nonlinear Dynamic and Kinematic Model of a Spar-Buoy: Parametric Resonance and Yaw Numerical Instability

: Mathematical models are essential for the design and control of offshore systems, to simulate the ﬂuid–structure interactions and predict the motions and the structural loads. In the development and derivation of the models, simplifying assumptions are normally required, usually implying linear kinematics and hydrodynamics. However, while the assumption of linear, small amplitude motion ﬁts traditional offshore problems, in normal operational conditions (it is desirable to stabilize ships, boats, and offshore platforms), large motion and potential dynamic instability may arise (e.g., harsh sea conditions). Furthermore, such nonlinearities are particularly evident in wave energy converters, as large motions are expected (and desired) to enhance power extraction. The inadequacy of linear models has led to an increasing number of publications and codes implementing nonlinear hydrodynamics. However, nonlinear kinematics has received very little attention, as few models yet consider six degrees of freedom and large rotations. This paper implements a nonlinear hydrodynamic and kinematic model for an archetypal ﬂoating structure, commonplace in offshore applications: an axisymmetric spar-buoy. The inﬂuence of nonlinear dynamics and kinematics causing coupling between modes of motion are demonstrated. The nonlinear dynamics are shown to cause parametric resonance in the roll and pitch degrees of freedom, while the nonlinear kinematics are shown to potentially cause numerical instability in the yaw degree of freedom. A case study example is presented to highlight the nonlinear dynamic and kinematic effects, and the importance of including a nominal restoring term in the yaw DoF presented.


Introduction
The kinematics and dynamics of floating bodies is traditionally related to offshore engineering problems, such as naval applications and the design of large oil and gas platforms [1]. For these applications, the main objective is usually to stabilize the motion of the floating objects, therefore the resulting small amplitude motions are within the limits of where linear theory is sufficiently accurate for modeling the system. However, harsh sea states may induce undesired large motions that do not satisfy linear assumptions, therefore making the model less representative. Note that load and motion prediction in such extreme conditions are imperative for a reliable design that ensures survivability and minimizes the risk of failures. Moreover, even in mild sea states, at certain wave frequencies the buoy may experience dynamic instability [2], threatening the integrity of the structure.
Conversely to conventional offshore applications, wave energy converters (WECs) are designed and controlled with the objective of enhancing the wave-induced motion to maximize power absorption [3,4]. Therefore, linear often models become inapt to accurately predict the behavior of a WEC [5]. The fidelity of mathematical models is crucial for a reliable estimation of the cost of electricity and for the effectiveness of model-based control strategy [6,7], which are essential for achieving economic viability [8,9]. As the wave energy field grows in experience and maturity, the necessity of nonlinear models, for a comprehensive design of most WEC types, becomes increasingly apparent [10][11][12].
For both traditional offshore engineering applications and novel WEC investigations, fully nonlinear hydodynamic models, such as the ones solving Navier-Stokes equations, can achieve a high level of accuracy in a broad range of operational and survival conditions. However, such models are not yet computationally viable for control or optimization applications. Therefore, more computationally efficient partially-nonlinear models, based on potential theory, have been the subject of numerous studies and developments, with research focusing on modeling nonlinear hydrodynamic forces [13][14][15][16]. However, scarcely any effort can be found towards modeling nonlinear kinematics in this area.
The consideration of nonlinear kinematics is usually necessary for systems with large amplitude motion and multiple, coupled degrees of freedom (DoFs). The inclusion of nonlinear kinematics is shown to be important in applications such as biomechanics [17,18], robotics [19,20], transportation [21,22], tracking control [23,24], and design of manipulators [25,26], to name a few. However, for ocean engineering applications, numerical models employed to simulate the dynamic behavior of a floating structure generally assume the motion to be planar, in the direction of wave travel, with up to three DoFs considered (horizontal translation, vertical translation, and rotation in the resulting plane: surge, heave, and pitch, respectively) [27,28]. Moreover, the rotational displacement and velocity are normally assumed to be small. Few nonlinear studies are performed in six DoFs, especially considering roll/pitch parametric resonance or yaw instability [29,30].
This paper presents a 6-DoF nonlinear model, relevant for conventional offshore applications as well as wave energy applications, including both nonlinear kinematics and nonlinear hydrodynamics. Typical modeling approaches are challenged, discussing potential issues arising from employing the usual simplifying assumptions. In particular, potential numerical instability in yaw may arise, since restoring and damping terms in this DoF are typically neglected. Furthermore, the model is able to properly articulate parametric resonance in roll and pitch.
The remainder of the paper is organized as follows. Section 2 provides details of the 6-DoF nonlinear model, derived for a generic axisymmetric spar-buoy. The possible coupling between DoFs and the situations in which they are likely to manifest are discussed in Section 3. The existence of yaw instability and its causes are presented in Section 4. A case study example is then provided in Section 5 to demonstrate the nonlinear dynamic and kinematic effects discussed in the paper. Finally, a number of conclusions are drawn in Section 6. Figure 1 depicts a generic spar-buoy structure, common in oil and gas applications, offshore wind towers, and a popular WEC concept. As one of the design objectives of such structures is to be independent of the incoming wave direction, they are normally axisymmetric. Unlike conventional models, which linearize the hydrodynamics and kinematics around the equilibrium position (shown in Figure 1a), the 6 DoF nonlinear model derived in this section calculates hydrodynamic forces and the kinematics with respect to the actual displaced position and orientation of the buoy (example depiction in Figure 1b), allowing the model to remain valid for large amplitude displacements and rotations. A unified framework for the vectorial parametrization of inertia, Coriolis and centrifugal, and hydrodynamic added mass forces for marine vehicles in 6 DoF was first presented in [31].  Figure 1. Schematic of a generic spar-buoy, depicting the reference frames for the nonlinear kinetic model. The inertial frame, (x, y, z), shown in blue, has its origin at the still water level (SWL) and the body-fixed frame (x,ŷ,ẑ), shown in green, has its origin at the buoy's centre of mass (CoM). The buoy is shown (a) at rest and (b) displaced.

Reference Frames
Two right-handed reference frames are defined, as schematically shown in Figure 1. The first, (x, y, z), is world-fixed, inertial, with the origin at the still water level (SWL) and on the center of the buoy at rest, the x−axis horizontal in the same positive direction as the wave propagation, the y−axis horizontal and orthogonal to the wave direction, and the z−axis vertical upwards. The inertial frame is used to describe the body displacements (ζ), divided into translations (p) and rotations (Θ): where φ, θ, and ψ are the rotations around the x−, y−, and z−axes respectively. The second right-handed frame of reference is (x,ŷ,ẑ), body-fixed, and therefore non-inertial. The origin is at the CoM of the buoy and (x,ŷ,ẑ) are aligned with (x, y, z) when the body is at its initial position. This is used for writing the dynamic equation of the system, since the inertial matrix remains constant. Therefore, both forces and velocities are represented in the body-fixed frame, along the axis of the buoy. Velocities (ν), divided into linear (v) and translational (ω), are defined as where p, q, and r are the rotational velocities around thex−,ŷ−, andẑ−axes, respectively.

Kinematic Mapping
It is worth remarking that forces and velocities are along time-varying axes, while displacements are along fixed axes. In linear hydrodynamic models there is no difference between such axes, based on the assumption of small displacements. However, in a nonlinear approach, a mapping from body-to world-frame velocities should be applied, at each time step, in order to obtain the correct displacements. One possible mapping is the following, where R Θ is the rotation matrix (applied to linear velocities), depending on the Euler angles Θ, defined according to the 3-2-1 convention as with c and s standing for cos() and sin() trigonometric operators, respectively. The mapping, T Θ , is applied to rotational ones, and is defined as follows, where t stands for the tan() trigonometric operator. Note that the singularity of T Θ in ±π/2 is usually not an issue in wave energy applications, as the amplitude of the pitch angle is, by design, always expected to be smaller than π/2.

Coriolis and Centripetal Forces
Another consequence of using a body-fixed frame are Coriolis and centripetal forces, which are normally neglected under the assumption of small rotational velocities. Let us define, for convenience of notation, the skew-symmetric operator S : R 3 → R 3×3 as S : Using such a notation, it is possible to define Coriolis and centripetal forces as where M is the mass of the body, r g is the vector from the origin of the body-fixed frame (reference point) to the CoM, and I r is the matrix of the moments of inertia with respect to the reference point. If the reference point is coincident with the CoM, then r g is the null vector and I r is a diagonal matrix, with I x , I y , and I z on the diagonal. Consequently, the Coriolis and centripetal force in (7) becomes

Hydrodynamic Forces
The wave-structure interaction is modeled using the partially nonlinear hydrodynamic model detailed in [32], which decomposes the total fluid force on the buoy into the Froude-Krylov (FK) force F FK , the diffraction force F d , the radiation force F r , and the viscous force F v . The model is labeled "partially nonlinear" as the diffraction and radiation forces are modeled linearly, whereas the viscous and FK force terms are nonlinear. The viscous force is described by an integral quadratic representation, and the nonlinear FK force is calculated by integrating the undisturbed pressure field from the incident wave over the instantaneous (updated at each time step) wetted surface of the buoy, balancing the gravity force. The FK and diffraction forces provide the external excitation to the system due to the incident waves. Full details of this hydrodynamic model are given in [32]. An open source toolbox for the implementation of the nonlinear FK method is provided in [33].

Mooring Force
To keep the spar-buoy on station, a mooring system is required. While the balance between the FK force and gravity provides a restoring force to the heave, roll, and pitch DoFs, the mooring system is the only mechanism to provide a restoring for to the surge, sway and yaw DoFs. The mooring force, F m , can be modeled using the methods reviewed in [34].

Equation of Motion
Finally, the dynamical equation in six DoFs for the buoy becomes where M is the inertial matrix, which is a (6×6) diagonal matrix with M 1,1 = M 2,2 = M 3,3 = M, M 4,4 = I x , M 5,5 = I y and M 6,6 = I z .

Validation
The validation of this model is presented in [35], performed against experiments of an axisymmetric oscillating water column (OWC) spar-bouy. As discussed in Section 3, the nonlinear dynamics and kinematics included in this model enables coupling and transfer of energy between DoFs. In particular, the OWC spar-buoy exhibits parametric resonance in pitch and roll in the experiments, which the model is shown to accurately reproduce. However, a dynamic instability not present in the experiments (due to the utilized mooring system) but which can manifest in the model, is yaw instability, as discussed in Section 4.

Coupling between DoFs
When employing a linear hydrodynamic model, an incoming unidirectional wave will induce a planar external excitation on an axisymmetric buoy (surge, heave and pitch). However, when considering nonlinear FK forces, coupling between DoFs can manifest under certain conditions, providing excitation to the lateral translational and rotational DoFs (sway and roll) [36]. In particular, when the excitation frequency is about twice the natural frequency in roll, a Mathieu-type of instability induces parametric excitation to the roll DoF [2], which in turn excites the sway DoF, due to a radiation coupling between roll and sway.
In these regions of parametric instability, a nonlinear FK model can provide 5 DoFs of excitation. Note, there is no means of exciting the yaw DoF. However, if nonlinear kinematics effects are introduced, the Coriolis and centripetal forces, as well as the kinematic mapping J Θ , have the mathematical structure to provide a coupling with yaw. The following sections will show that, if these forces and the kinematic mapping are not appropriately taken into account, then the model can exhibit numerical instability.

Kinematic Mapping
The last row of Equation (3) represents the mapping from the body-fixed rotational velocities, ω, to the rate of change of the yaw displacement,ψ: If φ is not exactly zero, Equation (10) shows that a velocity in roll can cause an acceleration in yaw. For the case of 3-DoF excitation (linear FK model or nonlinear FK model away from the parametric instability region), φ is not excited and simply decays from a small initial value φ 0 ; consequently,ψ ≈ 0. On the other hand, when 5-DoF excitation occurs (nonlinear FK model in the parametric instability region), roll is internally excited and eventually the term q sin φ cos θ is non-negligible and neither isψ. Therefore, it results that the yaw DoF is coupled with other rotational DoFs, either weakly (in the 3-DoF case) or strongly (in the 5-DoF case). However, if there is no excitation of the yaw DoF, the yaw displacement is bounded, asψ is of the same order of magnitude of q:ψ = O(q). Nevertheless, under certain conditions, yaw may also be weakly excited by Coriolis and centripetal forces, potentially inducing the model to be unstable and generate an unbounded yaw response.

Coriolis and Centripetal Forces
As shown in Equation (8), the surge component of F Cor is (introducing an informal index notation): Notice that in the simple 3-DoF case, the product rv ≈ 0, so that F Cor (1) ≈ −Mqw. Therefore, the mean of F Cor (1) depends on the phase difference between pitch and heave, which are both externally excited. In particular, a zero mean is obtained if the phase difference is 90 • , while strongly negative or positive means are obtained for phase differences of 0 or 180 • , respectively. In a linear hydrodynamic model, the surge exciting force has zero mean, so that the resulting surge displacement is bounded to have the same sign of the mean of F Cor (1), and magnitude depending on the mooring restoring force [37], as no hydrostatic force is present in surge. On the other hand, if a nonlinear hydrodynamic model is used, second-order drift effects shift the mean of the surge exciting force to positive values, so that the resulting mean displacement is a combination of both the wave and the Coriolis and centripetal forces.
The yaw component of the F Cor around the CoM, as shown in (8), is the following, Therefore, F Cor (6) is exactly zero if and only if I x = I y .

Numerical Yaw Instability
Generally, both intuition and experience teach that no significant yaw response is expected from an axisymmetric system. Physically, the only restoring force in yaw for an axisymmetric spar is provided by moorings. In many cases, the restoring is provided by the torsional stiffness of the mooring line, which is normally negligible. Thus no yaw restoring term is usually implemented in these types of numerical models [34]. In addition, when modeling yaw, no dissipative mechanism are usually implemented, because hydrodynamic radiation damping is ideally zero and viscous losses are reasonably negligible, due to the smooth axisymmetric geometry. While neglecting restoring and dissipative terms in the yaw DoF for conventional models is commonplace and without consequence, this section details that when such models are extended to include nonlinear kinematics unexpected yaw responses may arise, potentially generating conditions for numerical instability.

Initial Conditions
To study the possible existence and effects of instability, it is common practice to consider very small nonzero initial conditions in nonlinear hydrodynamic models, in order to provide some initial energy to all DoFs [38,39]. In the absence of external-internal excitation or strong coupling, the small initial conditions rapidly decay. Nonzero initial conditions are also consistent with the common application where a mathematical model is coupled with a physical system, taking measured displacement and/or velocity signals as inputs (either in an experimental test rig [40] or in real sea deployment [41]), due to the noise inherent in any measured signal. Therefore, let us take the initial conditions, (ζ 0 ), as nonzero and a small fraction of their expected steady state response, so that ζ 0 can be considered "almost" zero. In the following discussion, the initial roll displacement, (φ 0 ), is slightly greater than zero (say, 0.5 • ), so that φ,φ, and p are nonzero.
Furthermore, let us assume that the initial yaw displacement ψ 0 is zero. Although this is an unnecessary assumption, it will highlight that a response in yaw (with no external nor internal hydrodynamic excitation) can appear solely due to the nonlinear kinematics.

Transversal Moments of Inertia
Theoretically, for an axisymmetric system, the two transversal moments of inertia, I x and I y , should be identical. However, physically and numerically, it is often possible that the two transversal moments of inertia are not identical. In a real physical spar-buoy, while the hull geometry may practically be axisymmetric, the location of system components within the spar, and thus the distribution of mass, may not be perfectly axisymmetric. In a numerical model, the geometrical discretization of the buoy will be reproduced with finite accuracy, thus I x and I y may be not exactly the same as each other. Therefore, often the "ideal case" of I x = I y will not exist. Instead an "almost ideal case" (I x ≈ I y ) is more likely, whereby even the smallest difference between the two transversal moments of inertia will enable an internal excitation in yaw due to the Coriolis force, as discussed in Section 3.2.

Excitations
Let us consider the two hydrodynamic excitation conditions:  Table 1 summarizes all possible conditions that can arise in the absence of restoring and damping terms in yaw. In the ideal case (I x = I y ), there is no forcing term in yaw, so that the yaw response will follow roll and pitch angles, according to Equation (10). In particular, in the 3-DoF excitation condition, yaw will follow the decay of roll; in the 5-DoF excitation condition, the oscillatory part of yaw will follow the pitch sustained response, modulated by the frequency of the roll response. A slowly increasing mean of yaw is also present, due to the absence of a restoring force.
However, in the almost-ideal case (I x ≈ I y ), Equation (12) shows that there is a forcing term of the yaw DoF, much smaller in a 3-DoF scenario than a 5-DoF scenario, but never exactly zero. Consequently, due to the lack of viscous and restoring terms, the yaw DoF is not restrained and becomes unstable, so that its response diverges at a rate proportional to the difference between I x and I y . Therefore, when implementing Coriolis and centripetal forces in a 6-DoF model, it is important to include a yaw restoring term, which prevents the numerical instability from appearing. Table 1. Characteristic of the yaw oscillatory response, ψ. Considering no restoring nor damping terms in yaw and a small perturbation of the initial condition in roll and pitch. Two transversal moment of inertia cases are considered: ideal (I x = I y ) and almost-ideal (I x ≈ I y ), in combination with two hydrodynamic excitation conditions: 3-DoF (linear FK model or nonlinear FK model away from the parametric resonance region) and 5-DoF (nonlinear FK model close to the parametric resonance region).

Case Study Example
In this section example results are shown, highlighting the coupling between DoFs and dynamic instabilities, discussed in Sections 3 and 4, which can result when considering both nonlinear hydrodynamics and kinematics in the 6-DoF model of a floating spar-type structure. Section 5.1 provides a description of the test case and then the results are shown in Sections 5.2 and 5.3.

Test Case Description
The spar-buoy considered in this test case, has geometrical dimensions and mass distribution following the cylindrical spar studied in [42][43][44], which is a renown example of a parametrically unstable buoy, due to the approximate 2:1 ratio of the natural periods in pitch and roll, T n,4 and T n,5 , compared to the natural period in heave, T n,3 . Hong et al [42] performed physical experiments on the cylindrical spar in regular waves, using four wire-spring mooring lines to keep the buoy on station, and measured the responses in 5-DoF (surge, sway, heave, roll, and pitch). The studies in [43,44] perform numerical simulations in 2 DoF (heave and pitch), to investigate the parametric excitation in pitch, and do not include the mooring system, as surge and sway are not considered in the 2 DoF model. The present test case considers all 6 DoFs and includes a simplified mooring system to keep the buoy on station and to provide a mechanism for a restoring force in yaw. The relevant parameters for the system are listed in Table 2.
A simplified representation of the mooring force (F m ) is adopted in this study, to avoid additional complexity that would make the discussion more challenging and the cause-effect relationships less transparent. In particular, it is convenient that the mooring system does not provide coupling between DoFs, so that any coupling is only provided by kinematics and hydrodynamics. Thus, herein F m is simply considered as linearly proportional to displacements, with a coefficient K m,1 in surge/sway, K m,3 in heave, and K m,6 in yaw. The values for these mooring coefficients are listed in Table 2 and were selected to tune the natural periods for these DoFs to values comparable to those measured for the physical experiments in [42] and the numerical simulations in [32] (which investigates the effect of the mooring line parameters in inducing parametric resonance in a spar-buoy). The natural periods for the 6 DoFs are listed in Table 3, normalized against the pitch and roll natural periods.
The hydrodynamic forces are modeled following [32]. This model includes a viscous drag force, which provides a mechanism for damping to be applied to the yaw DoF. Modelling the viscous drag force utilizes an integral formulation (detailed in [32]), which requires the selection of a drag coefficient in each DoF. The drag coefficients employed here are the same as those used in [32] and are listed in Table 3, where C d,c is the drag coefficient around a circular cylinder and is applied to surge, sway, pitch, and roll, C d,3 is the drag coefficient in heave and C d,6 is the drag coefficient in yaw.
Following the discussion in Section 4.2, a 0.1% perturbation in one of the two transverse moments of inertia is considered (I y = 0.999I x ).

Parametric Resonance and Coupling between DoFs
The coupling between DoFs, and the ability of the model to capture the existence of parametric resonance, is highlighted in Figure 2. In this figure, a contour of the displacement amplitude in response to an input unidirectional regular wave, for each of the 6 DoFs, is plotted as a function of wave period, T w , and wave height, H w . As discussed in Section 3, there is no external excitation in the sway, roll nor yaw DoFs from a unidirectional wave input, due to the symmetry of the system, thus any motion in these DoFs is a consequence of the internal excitation arising from the coupling provided by the nonlinear dynamics and kinematics. In particular, parametric resonance, a Mathieu-type dynamic instability, is expected around an external wave period of 1 2 T n,5 . Finally, note that all parameters and results here presented are normalized, enabling application to structures of varying size, such as the large spars in [43,44], and to smaller WEC-like structures, as in [45]. The relevant common feature is to realize a 2:1 ratio between pitch and heave natural periods.  Figure 2 shows that at the roll and pitch natural period (denoted by the dash-dotted red line), the pitch DoF experiences a resonant amplitude increase, whereas the roll DoF exhibits zero motion. However, at half the roll and pitch natural period (denoted by the dashed red line), both the pitch and roll DoFs exhibit a very large amplitude increase, due to parametric resonance, as expected. In addition, at the same period, motion is detected in the sway DoF and there is a clear reduction of heave response, due to an internal exchange of energy between DoFs. Finally, under the 5-DoF excitation condition occurring with the existence of parametric resonance around half the roll and pitch natural period, there is a small response in the yaw DoF, made possible by the nonlinear kinematics and the perturbation of the transverse moment of inertia. However, due to the restoring term in yaw, numerical instability is avoided and the yaw response is contained below 1 degree (explored further in the next subsection).
To further highlight the effect of nonlinearities for different frequencies, Figure 3 shows the response amplitude operator (RAO) for different DoFs, with the color code proportional to the normalized wave height. As in a fully-linear system the response is proportional to the input, the RAO is insensitive to changes in wave height and all curves would overlap. Therefore, variations in RAOs for different H w are due to nonlinear effects. Figure 3 clearly shows that such deviations from linear response are around the two parametric resonance periods, i.e., at T w = 1 2 T n,5 and T w = T n,5 . Consistently with Figure 2, there is a sharp drop of heave response at T w = 1 2 T n,5 . The region of parametric instability, particularly evident in the roll and pitch RAOs, becomes wider as H w increases; however, the severity of parametric response grows less than linearly with H w , since the RAO decreases for larger H w . Conversely, yaw response, although small, is the only DoF showing a superlinear response, i.e., faster than the growth of H w . Furthermore, it is worth describing the trend of the internally-excited-only DoFs (sway, roll, and yaw). In linear conditions, the RAO would be null due to the absence of external excitation. However, due to the initial energy provided in all DoFs, nonzero RAOs are found across the whole frequency range, also in linear conditions (for small waves and/or away from parametric resonance). Moreover, the initial condition (the same regardless of the wave height) divided by smaller waves appears as slightly higher RAO. Focusing on the parametric instability region, the sustained nonlinear response is only due to dynamic instability related to parametric resonance, and the final amplitude depends on both the wave height and the dissipation of the system. As viscous drag grows quadratically with the wave height, the amplification of parametric response is weaker than the increase in wave height, therefore partially explaining why RAO(4) is smaller for larger H w . Furthermore, Figure 3 shows that, due to nonlinear kinematic coupling, as the wave height increases there is additional energy transferred from RAO(4) and RAO (5) to RAO (6). For further detailed discussion of the nonlinear dynamics of parametric resonance, please refer to the work in [35,46].

Yaw Instability
As discussed in Section 4, a consequence of including nonlinear kinematics in the 6-DoF model is that unexpected yaw responses may arise, potentially generating conditions for numerical instability. In particular, while the common approach in conventional models of neglecting restoring and damping forces in yaw is generally valid, due to their insignificant contribution, for the present model with nonlinear hydrodynamics and kinematics, this can lead to numerical instability. The effects of neglecting restoring and damping forces in yaw is demonstrated in Section 5.3.1. Next, in Section 5.3.2, the influence of including the restoring and damping forces in yaw on the numerical stability of the model is shown.

Neglecting Restoring and Damping Forces in Yaw
Considering no restoring nor damping terms in yaw, and a small perturbation of the initial condition in roll and pitch, four different scenarios can be considered which dictate the resulting stability of the yaw response, as tabulated in Table 1. Each of these scenarios are exemplified in Figures 4-7, showing the time traces for the 6-DoF response to representative incoming regular waves. The wave height (H w ) is chosen to be about equal to the metacentric height (0.95GM in particular) and two different wave periods are considered. To demonstrate the two scenarios with the 3-DoF excitation, the wave period is set to 3 4 T n,5 , exactly in between T n,5 and the parametric resonance period ( 1 2 T n,5 ). For the other two scenarios, with the 5-DoF excitation, the wave period is set to 1 2 T n,5 , so that parametric instability is triggered. The time traces are computed for a time-window 200 times long the wave period, with an initial sigmoid ramp five wave periods long, using a Runge-Kutta time integration scheme of the second order, with 75 time steps per wave period. Initial positions are defined according to Section 4.1. Figure 4 shows the case of the ideal transversal inertia conditions (I x = I y ) with 3-DoF hydrodynamic excitation. A transient response is seen for each of the 6-DoFs, due to the nonzero initial conditions imposed in the model. The motion of the 3 DoFs with external excitation (surge, heave, and pitch), is seen to reach a steady state response, whereas the motion of the other 3 DoFs, which includes yaw, decays with time. The nonzero initial conditions and the following transient response sees the buoy rotate less than 0.1 degree in yaw, and then the energy in this DoF begins to decay, consistent with the response characterized in Table 1.    The other case with the ideal transversal inertia conditions (I x = I y ), but considering 5-DoF hydrodynamic excitation is shown in Figure 5. Here, the hallmark exponential increase in oscillation amplitude, associated with parametric resonance, can be seen in the pitch and roll DoFs. When this occurs the heave motion amplitude is seen to decrease, highlighting the internal transfer of energy from this DoF to the others. Consistent with Table 1, the yaw response is seen to grow linearly with time for this scenario.
Considering the scenarios where there is a tiny difference between transversal inertia conditions (I x ≈ I y ), Figures 6 and 7 show the results for which Section 4 predicts the yaw response to be unstable. The case with 3-DoF hydrodynamic excitation is displayed in Figure 6, for which the yaw response is seen to be growing exponentially, though only rotating a total of 20 degrees after 200 wave periods. The instability is much more pronounced for the case with 5-DoF hydrodynamic excitation in Figure 7, where the buoy has approximately completed two and a half full rotations after 200 wave periods and the responses for the other DoFs have grown unstable (i.e., a pitch displacement of 90 degrees corresponds to a capsized buoy, with the spar oriented horizontally).

Including Restoring and Damping Forces in Yaw
As shown in Section 5.2, including restoring and damping forces in yaw prevents the occurrence of yaw instability. In this section, the influence of these two forces in stabilizing the yaw response, will be investigated separately. The three scenarios from Table 1, for which the yaw response did not decay, that were demonstrated in Section 5.3.1 and shown in Figures 5-7, are again simulated here but with the yaw restoring and/or damping forces included in the numerical model. Figure 8 considers the cases in which the restoring force is not included in the model, for which the yaw response with and without the damping force included in the model is plotted. For Figure 8a,b, corresponding to the cases in Figures 5 and 6, respectively, almost no difference in the yaw response can be observed by introducing the damping force. For Figure 8c, corresponding to the case in Figure 7, the inclusion of damping is seen to slow the large increase in yaw response compared to the undamped case, but does not eliminate the instability. Figure 9 shows the cases in which the restoring force is included. Here it can be seen that the influence of the restoring force eliminates the instability in yaw, with the maximum yaw response for all three cases only reaching around 1 degree. Again in this figure, the effect of including damping in yaw is shown to have negligible effects on the response.

Conclusions
This paper proposes a model in six degrees of freedom for axisymmetric buoys, including nonlinear kinematics, Coriolis and centripetal forces, and nonlinear Froude-Krylov forces. The model is shown to enable energy transfer between degrees of freedom not present in more simplistic models for such systems, due to the couplings introduced by the nonlinear dynamics and kinematics. The nonlinear FK enables the model to predict the existence of parametric resonance, a Mathieu type of instability, which is observed in physical experiments. The inclusion of nonlinear kinematics is shown to make the occurrence of numerical instability a possibility, due to a coupling between the pitch and roll DoFs with the yaw DoF, which can be avoided by the inclusion of a restoring term in yaw.