Modelling of Parametric Resonance for Heaving Buoys with Position-Varying Waterplane Area

: Exploiting parametric resonance may enable increased performance for wave energy converters (WECs). By designing the geometry of a heaving WEC, it is possible to introduce a heave-to-heave Mathieu instability that can trigger parametric resonance. To evaluate the potential of such a WEC, a mathematical model is introduced in this paper for a heaving buoy with a non-constant waterplane area in monochromatic waves. The efﬁcacy of the model in capturing parametric resonance is veriﬁed by a comparison against the results from a nonlinear Froude–Krylov force model, which numerically calculates the forces on the buoy based on the evolving wetted surface area. The introduced model is more than 1000 times faster than the nonlinear Froude–Krylov force model and also provides the signiﬁcant beneﬁt of enabling analytical investigation techniques to be utilised.


Introduction
Renewable energy sources are an attractive form of energy supply for the growing human population, since these resources are inexhaustible and clean. Ocean wave energy presents an enormous renewable energy resource, spurring significant research towards the development of wave energy converters (WECs) in recent decades [1,2]. To maximise energy absorption, WECs are typically designed to resonate with the dominant frequency of the input waves. Thus, the concept of resonance is very well known and applied within the field of wave energy conversion. In contrast, the phenomenon of parametric resonance has received far less attention, and its potential to increase WEC performance is mostly unexplored.
Parametric resonance is a dynamic instability caused by the time-varying changes in the parameters of system [3]. The occurrence of parametric resonance results (at least initially) in an exponential growth of oscillations of a system with time. In comparison, normal resonance only causes the oscillation amplitude to increase linearly with time. Therefore, exploiting the dynamic instability from parametric resonance may prove beneficial for WEC performance by enabling large oscillations for increased power output.
A major barrier in studying parametric resonance for WECs is the complexity of the mathematical models required to capture this nonlinear phenomenon, compared to the computationally efficient linear/frequency domain hydrodynamic models traditionally utilised in WEC research and analysis [4,5]. However, with the rise in availability of highperformance computing, the use of nonlinear hydrodynamic models is becoming more prevalent in the WEC design process [6], enabling parametric resonance to be captured by the simulations. In other fields and applications, parametric resonance is well understood, and a wealth of analysis tools and techniques have been developed by the nonlinear dynamics community [7][8][9][10][11][12]. These tools enable the conditions for the existence of parametric resonance to be efficiently examined across large parameter spaces; however, they require an analytical model of the system under study.
The work in this paper first proposes a method for exploiting parametric resonance in heaving WECs via a non-constant waterplane area, which, to the authors' knowledge, has not been previously reported in the literature. As a first step in examining the potential of this method, an analytical model for the system is developed, which enables the tools from nonlinear dynamics to be employed in the design and assessment of the proposed concept.

Parametric Resonance in Heaving Wave Energy Converters
Heaving WECs are one of the most common and widely proposed WEC concepts, designed to react to wave-induced vertical motions of a buoy against: the sea floor, an inertial damping plate, a second less dynamic body or an internal water column. Like other floating bodies, heaving WECs are susceptible to the existence of large unstable parametric pitch and roll motions due to a Mathieu-type instability. The heave motion of the buoy causes a parametric excitation in the pitch and/or roll degrees of freedom (DoFs) due to the time variation of the restoring torque parameters. The restoring torque coefficient depends on the position of the centre of buoyancy with respect to the centre of mass, and since the centre of buoyancy moves as the buoy heaves into and out of the water, the heave motion causes a time variation in the pitch and roll restoring torque parameters.
This type of parametric resonance in the pitch and/or roll DoFs of floating bodies has been studied since the 19th century, following Froude's observations that large roll motions occur when the roll natural period of a ship is twice the heave or pitch natural period [13]. The parametric coupling between heave, pitch and roll DoFs has been investigated in floating structures, such as ships [14][15][16][17][18] and spars [19][20][21][22][23]. The large amplitude oscillations resulting from parametric resonance can jeopardise the safe operation of floating structures and vessels, such as cargo being lost overboard due to large parametric roll motions on container ships. The main focus of research on parametric resonance within naval architecture and offshore engineering has therefore been towards control methods for the mitigation and stabilisation of parametric pitch/roll [24][25][26][27][28].
Heave-to-pitch and heave-to-roll Mathieu-type instabilities have been observed in heaving WECs [29][30][31][32][33][34][35]. This parametric transfer of energy from the primary heave mode of motion to the pitch and roll DoFS leads to reduced power capture from the WEC. Therefore, several studies have investigated control methods to mitigate the onset of parametric resonance, which can be classified as passive or active. The passive control methods consider the inclusion of strakes/fins [36], utilising the mooring system [37] or modifying the initial design of the geometry and mass distribution [38]. The active control methods utilise systems, such as the power-take off [39], variable inertia mechanisms [40] or a pressure relief valve, in an oscillating water column [41,42]. A real-time detection system for the onset of parametric resonance in WECs has also been developed to give early warning to such active control systems [43].
The present paper proposes a new method for exploiting parametric resonance in heaving WECs utilising a heave-to-heave Mathieu instability. To the authors' knowledge, this type of instability has not previously been examined for the case of wave energy conversion. The parametric excitation in the heave DoF is facilitated through a time variation in the heave restoring and excitation force parameters, which is physically realised through a geometry with a non-constant cross-sectional area. The occurrence of a heaveto-heave Mathieu instability for this type of geometry has been investigated for Arctic spars [23], which have a sloped geometry at the free surface interface to break sea ice. The application of this type of spar geometry to exploit parametric resonance in wave energy conversion has not been explored. Several heaving WECs already utilise sloped geometries at the free surface, such as the 'Free Floating Clam' [75], the 'WedgeTop WEC' [76] and the 'MoonWEC' [77], for example. The advantage of such geometries, as stated by Francis Farley in [75], is that they "act like a vertically driven wedge which is strongly coupled to the waves". However, none of these studies have considered the occurrence of heave-to-heave parametric resonance, which is largely due to the conventional linear hydrodynamic modelling methods utilised in their design and analysis.

Modelling and Analysis of Parametric Resonance
The conventional linear modelling methods for WEC simulation and analysis do not predict the occurrence of parametric resonance. Therefore, more complex, high-fidelity modelling techniques are required. Computational fluid dynamics (CFD) enable a fully nonlinear hydrodynamic modelling and provide the highest fidelity tool available for WEC simulation. While the computational requirements of CFD traditionally rendered its use for WEC analysis infeasible, the increasing availability of high-performance computing resources over time has led to a steady growth in the use of CFD for WEC simulation. The review of CFD simulations of WECs in [78] collects over 200 publications from 2004 to 2017, none of which report or consider parametric resonance. This might be due to the number of wave periods necessary for the onset of parametric resonance to occur, requiring longer simulation durations than typically utilised in CFD. However, in 2018, Eskilsson et al. [79] observed the occurrence of parametric resonance in their CFD simulations of a moored heaving point absorber, leading them to produce a follow-up study applying CFD to analyse the parametric excitation of moored, full-scale WECs in six DoF [80]. The simulations required 1600 CPU hours per wave period, highlighting the large computational expense of this modelling method.
Due to the high computational demand of CFD, more computationally efficient partially nonlinear models are typically preferred, which are a compromise between the CFD approach of solving the Reynolds Averaged Navier-Stokes equations and the conventional modelling approach of linearising the solution to the velocity potential at the equilibrium position of the WEC. The nonlinear Froude-Krylov force method (NLFK) extends the conventional modelling approach by applying the solution to the velocity potential to the instantaneous position of the WEC at each time step, which corresponds to integrating the pressure from the undisturbed incident wave field over the evolving instantaneous wetted body surface of the WEC throughout the simulation. An NLFK model, therefore, incorporates both nonlinear excitation and restoring forces, whereby the hydrodynamic parameters vary at each time step, being recalculated for the exact wetted body surface at that instant. This modelling approach has been utilised to successfully capture the parametric pitch/roll motion in heaving point absorber WECs, such as the Wavebob [32,37,40,81,82] , the CorPower WEC [83] and the Oscillating Water Column Spar Buoy [84][85][86][87]. A comparison of a NLFK against a CFD simulation for capturing parametric resonance in a heaving point absorber is performed in [88], reporting that the NLFK model is more than 5000 times faster than the CFD simulation with comparable accuracy.
While CFD and the NLFK methods have been demonstrated to capture parametric resonance in WECs, the analytical consideration of parametric excitation requires differential equations with time-varying coefficients [89]. The classic mathematical example of parametric resonance is the Mathieu (or Hill) Equation [90], which is a second-order ODE with no external forcing and a parameter, a(t), that varies harmonically with time: Since the conventional linear hydrodynamic models used in WEC simulation and analysis are time-invariant differential equations, several studies have extended these conventional linear hydrodynamic models with a time-varying restoring term in the pitch/roll DoFs. This approach yields a type of damped Mathieu equation with external forcing due to the inclusion of the hydrodynamic damping and excitation force terms from the conventional WEC model.
Gomes et al. [36,38] derive a Mathieu equation to predict the occurrence of parametric resonance observed in the Oscillating Water Column Spar Buoy WEC. The restoring torque is derived based on the instantaneous position of the metacentric height, which is calculated using an approximation from shipping [91]. Gomes et al. do not perform simulations with the model but utilise it to identify the input wave conditions for which parametric resonance is likely, using the Mathieu diagram. Considering a simple cylinder, Zou et al. [45,47] and Abdelkhalik et al. [44,46] use geometry and trigonometry to describe the co-ordinates of the centre of buoyancy and the centre of gravity as a function of the heave and pitch displacement, from which the pitch restoring moment can be calculated. They apply a Taylor expansion to linearise the resulting expression and include it in a conventional hydrodynamic model to provide a position-dependent/time-varying restoring torque that is capable of capturing parametric resonance. Similarly, Kurniawan et al. [34] develop a time domain model with a pitch-restoring torque coefficient that is dependent on the instantaneous position of the centre of buoyancy, motivated by the observation of parametric resonance in their physical experiments that was not captured by their original frequency domain models. Davidson et al. [92] utilise system identification techniques [93] to model the restoring force/torque as a function of heave and pitch displacement. The hydrostatic restoring force/torque is measured in a series of CFD numerical wave tank experiments, where the WEC is slowly submerged whilst orientated at different pitch angles. A finite set of basis functions is then fitted to the data as a function of the heave and pitch displacement and used to extend a conventional linear hydrodynamic model to include a time-varying restoring torque, which is demonstrated in a simulation to capture parametric resonance. Davidson and Kalmár-Nagy [43] utilise a model derived for a spar buoy [94][95][96] that describes the pitch-restoring torque as a function of the heave displacement and free surface elevation based on the instantaneous position of the centre of buoyancy.
While all of these studies consider the heave-to-pitch/roll Mathieu instability, none of these methods are able to cater for the case of self-excited parametric resonance in the heave direction from a heave-to-heave Mathieu instability due to a time-varying waterplane area. The objective of the present paper is to fill this gap in the literature by developing an analytical model for heaving WECs with a time-varying waterplane area.

Outline and Objectives of Paper
In this paper, a mathematical model for a non-constant waterplane area heaving WECs that is able to capture parametric resonance is introduced. The model describes the motion of a heaving buoy in monochromatic waves using a nonlinear parametrically excited differential equation. The developed model not only provides a much faster simulation than numerical methods, such as NLFK force and CFD simulations, but also enables analytical investigations using methods that are extremely efficient in demonstrating, predicting, and describing phenomena caused by nonlinear effects, such as perturbation or harmonic balance methods.
The remainder of the paper is organised as follows. In Section 2, the analytical model of parametric resonance for a heaving WEC with a non-constant cross-section is presented. A test case to evaluate the performance of the analytical model is introduced in Section 3. The results for the test case are summarised in Section 4, and some conclusions are drawn in Section 5.

Mathematical Model
In this section, the mathematical model of a heaving WEC, considering a singlefrequency (monochromatic) input wave, is discussed. A schematic view of an axisymmetric WEC with a center of mass G is shown in Figure 1a. The motion of the WEC is restricted to the vertical z (heave) direction only. The origin of the world coordinate system z is at the still water level (SWL)-see Figure 1a. The wave elevation η(t) is also measured with respect to the SWL and is considered along the vertical line through the position of the buoy's CoM, as illustrated in Figure 1a. For a harmonic wave, the elevation is given by where H and ω are the wave amplitude and frequency, respectively. The free-body diagram in Figure 1b shows the forces acting on the buoy, namely the weight F g and the hydrodynamic force F h , resulting from the integral of the fluid pressure over the surface of the buoy. Therefore, the equation of motion of the buoy can be written as In conventional hydrodynamic models for WECs, F h is described using linear theory, as detailed in Section 2.1. The present paper extends this conventional linear hydrodynamic model to include parametric resonance, as detailed in Section 2.2.

Conventional Linear Hydrodynamic Model
Linear hydrodynamic modelling [4,5,97] assumes that the total force of the fluid can be decomposed into three separate forces: hydrostatic restoring force F s , the wave radiation force F r and the wave excitation force F e : Each of these forces are discussed in Sections 2. The hydrostatic restoring force equals the difference between the buoyancy force from the fluid and the weight of the WEC: where ρ is the density of the sea water, g is the gravitational constant and V(z) is the submerged volume of the WEC, which increases/decreases as the WEC heaves into/out of the water. At equilibrium, z = 0, the buoyancy force balances the weight of the WEC, ρgV(0) = F g . Therefore, the hydrostatic restoring force equals the weight of the fluid, which is displaced by the change in the submerged WEC volume: where ∆V(z) = V(z) − V(0). For a linear hydrodynamic model, the change in the submerged WEC volume is assumed to be linearly proportional to the waterplane area A s of the WEC's cross-section at equilibrium: This assumption is strictly true for geometries with constant cross-sections, such as vertical cylinders, and is otherwise approximately true for very small heave displacements from the equilibrium. However, for geometries with non-constant cross-sectional areas, this linearising assumption breaks down with the increase in heave amplitude.

The Wave Radiation Force
The wave radiation is the reaction force due to the radiated wave created by a body moving in the fluid, comprising an inertia term due to the fluid mass surrounding the body that must be accelerated and a resistance term due to the energy carried away by the wave: where m a (ω) and c a (ω) are the frequency-dependent added mass and damping coefficients, respectively. The added mass and damping parameters depend on the WEC geometry and the properties of the surrounding fluid.

The Wave Excitation Force
The wave excitation force is the force exerted on a body that is held fixed in the presence of an incident wave. For linear hydrodynamic models, the excitation force is assumed to be linearly proportional to the wave amplitude and is also phase shifted with respect to the position of the wave: where f e (ω) is the excitation force coefficient, and θ e (ω) is the phase angle between the wave and the excitation force. Both of these parameters are frequency-dependent and also depend on the geometry of the WEC and the properties of the surrounding fluid.

The Linear Hydrodynamic Model
Combining the expressions for the three hydrodynamic forces in in Sections 2.

Analytical Model for Parametric Resonance
Here the linear hydrodynamic model in Section 2.1.4 is modified to more accurately describe the case of a heaving WEC with a non-constant waterplane area. In Section 2.2.1, a nonlinear hydrostatic restoring force is introduced to better represent the actual, rather than linearised, change in the submerged WEC volume as a function of the heave displacement. Section 2.2.2 derives a parametric forcing term by introducing position dependence to the excitation force coefficients. The nonlinear hydrostatic restoring force and parametric excitation terms are incorporated in the conventional hydrodynamic model in Section 2.2.3, which is rearranged to cast it in a form similar to the Mathieu equation.

Nonlinear Hydrostatic Restoring Force
The hydrostatic restoring force is proportional to the change in the displaced water volume due to the heave motion of the WEC, ∆V(z). Rather than assuming ∆V(z) is linearly proportional to the heave motion, as in Equation (7), a nonlinear function can be utilised to more accurately represent ∆V(z) in order to increase the fidelity of the hydrostatic restoring force term.
Nonlinear hydrostatic restoring force terms have previously been implemented in WEC models, whereby the restoring force is expressed as a polynomial function of the displacement [98][99][100]. In these models, the polynomial expressions for the restoring force are identified using data of the hydrostatic force versus the WEC displacement, as measured from static tests in the laboratory [98] or CFD-based numerical wave tanks [99,100]. In the present paper, a polynomial expression is instead identified for ∆V(z), rather than the hydrostatic force directly, using data generated from an analytical function describing the WEC geometry.
Considering an arbitrary shaped buoy, with a cross-sectional area A w (ζ) that can vary along its vertical axis, ζ. The change in the displaced volume of water due to a heave displacement of z can be written as: For axisymmetric buoys, Equation (11) can be simplified to: where the axisymmetric geometry of the buoy is defined by the revolution of a generic function f (ζ), ζ ∈ [−σ 0 , σ 1 ], see Figure 2, where the origin of the body fixed coordinate system ζ is denoted by G, and it is located at the still water level when the buoy is at rest. The integral expression for ∆V(z) in either Equation (11) or (12) can be approximated by using an nth degree polynomial expression, i.e., where the coefficients c 1 , . . . , c n can be determined by the least squares method. The hydrostatic restoring force can then be expressed as

Parametric Excitation
Here the wave excitation force is recast into a parametric excitation term, which to the authors' knowledge has not been proposed before. This is achieved by introducing time-varying parameters to the excitation force model (Equation (9)). While the conventional linear hydrodynamic model determines the excitation force coefficient, f e (ω), and phase angle, θ e (ω), based on the wetted surface area of the submerged buoy geometry at equilibrium, the proposed method here determines these two parameters based on the instantaneous heave position. Therefore, rather than being constant, the excitation force coefficient, f e (z, ω), and phase angle, θ e (z, ω), now depend on the position of the buoy and convert the wave excitation force into a parametric forcing term: Similar to the derivation of the nonlinear hydrostatic restoring force in Section 2.2.1, the excitation force coefficient and phase angle can be expressed as a nonlinear function of the heave displacement by fitting polynomial functions to data where the coefficients a 1 , . . . , a n and b 1 , . . . , b n can be determined (e.g., by the least squares method) from a dataset comprising the values of f e (ω) and θ e (ω) for several different submergence levels of the buoy (rather than just at the equilibrium position as is the case for the linear hydrodynamic model). The general form of the wave excitation force then becomes:

Mathieu Type Equation
Incorporating Equations (14) and (17) into Equation (10), the analytical model for a heaving buoy with a non-constant cross-sectional area is This can be recast into a Mathieu/Hill-type form

Test Case
To evaluate the proposed analytical model, a test case is presented comparing the performance of the model against two other models. First, to investigate the influence of the novel parametric excitation term (introduced in Section 2.2.2), the results of the analytical model are compared against the results from a model that only extends the conventional linear model with a nonlinear hydrostatic restoring force term. Next, the performance of the analytical model is compared and verified against a high-fidelity benchmark produced herein by the NLFK model implemented using an open source toolbox [101]. The WEC geometry is introduced in Section 3.1, and the derivation of the analytical model for this particular geometry is presented in Section 3.2. The nonlinear hydrostatic and NLFK force models and their implementation are detailed in Sections 3.3 and 3.4, respectively.

Test Case Geometry
The test case considers a heaving point absorber WEC with a non-constant waterplane area, as shown in Figure 3, whose geometry comprises a truncated cone-shaped buoy with a cylindrical extension. The parameters of the buoy geometry are listed in Table 1. The shape of this buoy is defined by where γ = tan α. To be consistent with Figure 2, we note that for the test case, σ 0 = h 0 /2 + h 1 and σ 1 = h 0 /2.

Parameter Value Units
Opening angle of the conical part α 11.31

The Analytical Model
The calculation of the parameters in the analytical model (Equation (19)) for the test case WEC geometry is detailed here. Section 3.2.1 describes the calculation of the hydrostatic restoring force parameters, Section 3.2.2 the parametric excitation and Section 3.2.3 the radiation force.

Hydrostatic Restoring Force
To calculate the nonlinear hydrostatic force, the displaced water volume has to be determined. Using Equation (11), the volume of the displaced water due to variations in heave position z can be calculated for the test case geometry. Substituting the displaced water volume into Equation (6), the nonlinear hydrostatic force of the test case geometry can be written as The water density and the gravitational acceleration were set to ρ = 1025 [kg/m 3 ] and g = 9.806 [m/s 2 ], respectively. The resulting nonlinear hydrostatic force F s (z) for the test case geometry is illustrated in Figure 4. The piecewise function of the hydrostatic force has a strong nonlinear characteristic in the −h 0 /2 < z ≤ h 0 /2 domain.

Parametric Excitation
The wave excitation force amplitude and phase coefficients were calculated by the open source code NEMOH [102]. This code is a numerical solver for computation of firstorder hydrodynamic coefficients, such as added mass, radiation, damping and excitation forces, in the frequency domain. NEMOH is based on linear free surface potential flow theory with assumptions of an inviscid fluid and an incompressible and irrotational flow. The resulting linear Boundary Value Problem (BVP) for the free surface flow around a body is of first order with assumptions of small motions around the mean position and linearised free surface equations. We enable the simulation of large motions by calculating the hydrodynamic coefficients at numerous buoy submergence levels. The assumption of small motions mean that the calculated hydrodynamic parameters are only accurate in a small displacement range around the submergence level they were calculated. Therefore, we allow the parameters in the model to change based on the displacement, so they match the values corresponding to the instantaneous submergence level during the simulation.
In order to solve the linear BVP, panel methods are applied in NEMOH. For the present test case, the buoy geometry was discretised using 445 panels and 1120 nodes. The mesh created for NEMOH is illustrated in Figure 5. The mesh independence was checked by using two other meshes with 230 and 880 panels.  Figure 6 illustrates the parameters of the wave excitation force ( f e (z, ω), θ e (z, ω)) for the test case buoy. For this current test case, second-order polynomials are employed to approximate the wave excitation amplitude and phase as a function of the buoy position for each frequency, i.e., where the frequency-dependent coefficients a 0 (ω), a 1 (ω), a 2 (ω), b 0 (ω), b 1 (ω), b 2 (ω) are least-square fitted on the NEMOH data. The coefficients as the function of the wave frequency ω are illustrated in Figure 7. The accuracy of the fitted wave excitation amplitude f e (z, ω) and phase θ e (z, ω) are illustrated in Figure 8, where the continuous line denotes the fitted curve, while the dots are the data calculated by NEMOH.

Wave Radiation Force
The added mass m a (ω) and damping c a (ω) were also calculated using NEMOH for the z = 0 buoy position (as they are considered here as independent of the position). Figure 9 illustrates the added mass and radiation damping coefficients as a function of the wave frequency ω.
, Figure 9. The hydrodynamic coefficients for the (a) added mass m a (ω) and (b) damping c a (ω) coefficients.

The Nonlinear Hydrostatic Restoring Force Model
The nonlinear hydrostatic restoring force model extends the conventional linear hydrodynamic model (Equation (10)) with the nonlinear restoring force described in Section 2.2.1: (m + m a (ω))z + c a (ω)ż + ρg n ∑ k=1 c k z k = f e (ω)H cos(ωt − θ e (ω)). (23) This type of nonlinear hydrodynamic model has previously been utilised for WECs with position-varying waterplane areas [98-100] but does not enable parametric resonance to be captured. The parameters for the hydrostatic restoring force and the radiation force are the same as calculated in Section 3.2.1 and Section 3.2.3, respectively. The parameters for the exciation force are the same as calculated for the equilibrium postion (z = 0) in Section 3.2.2.

The Nonlinear Froude-Krylov Force Model
In this section, we describe the NLFK model used for the code-to-code verification of the analytical model developed in this paper. The hydrodynamic force in the NLFK is considered in the following way: where F FK z (z, t, ω) is the Froude-Krylov force acting in z-direction, and F r (ż,z, ω) is the linear radiation force (see Equation (8)). The resulting equation of motion is (m + m a (ω))z + c a (ω)ż = F FK z (z, t, ω).
Froude-Krylov (FK) forces correspond to the integral of the pressure (static and dynamic) of the undisturbed incident wave field over the wetted surface of the device. Such a pressure is defined, according to linear Airy's theory, with the application of Wheeler stretching [103], as: where p st and p dy are the static and dynamic pressure, respectively. The variable z is the change of coordinates required by Wheeler stretching in order to eliminate the free surface boundary condition error: z = h z+h η+h − h, whereη is an approximation of the free surface, computed at the center of gravity. The FK force can be divided into the static and dynamic FK forces The static FK force is defined as: where F g = (0, 0, −mg) T is the weight, S w the instantaneous wetted surface and n the unity vector normal to the surface. The dynamic FK force is similarly defined as the surface integral of the dynamic pressure: To calculate the Froude-Krylov force F FK z in the time domain simulations, the opensource NLFK4ALL MATLAB toolbox developed by Giorgi [101,104] was utilised. The Froude-Krylov force F FK z (z, t, ω) was calculated every time step using NLFK4ALL toolbox, and the equation of motion was integrated using the fourth order Runge-Kutta method. To choose the proper time step for the integration, we checked time sensitivity of the solution on the time step, see Figure 10. To reduce the simulation time and keep the proper accuracy, ∆t = 0.002 [s] was chosen as the simulation time step. Figure 11 illustrates a rendering of the test case geometry provided within the NLFK4ALL toolbox.

Comparison against Nonlinear Hydrostatic Restoring Force Model
Parametric resonance is expected to occur for input waves with a frequency twice the buoy's natural frequency and with amplitudes above a certain threshold. To examine the ability of the analytical model to capture parametric resonance, numerous simulations were performed at frequencies close to twice the buoy's natural frequency for a range of wave heights. The buoy's natural frequency was measured from free decay experiments (see Section 4.2) and equals ω 0 = 0.944 [rad/s]. Figure 12 shows the buoy steady-state oscillation amplitude, normalised against the wave amplitude, z max /H, as a function of the wave amplitude H. For a linear system, normalising the output by the input amplitude results in a constant as a function of the input amplitude. Therefore, due to the normalisation in Figure 12, any deviation from a constant response highlights the nonlinearity in the output and the threshold wave amplitude where this nonlinearity manifests. The results are shown for three different wave frequencies ω ∈ {1.87, 1.88, 1.89} [rad/s] in Figures 12a, 12b and 12c, respectively. The solid line denotes the oscillation amplitudes if we apply the new analytical parametric wave excitation model, while the dashed line is the oscillation amplitudes in the case of the simpler nonlinear hydrostatic restoring force model. The results show the analytical model is able to capture parametric resonance, signified by a large dramatic increase in the normalised oscillation amplitude when the wave amplitude reaches a certain threshold, whereas the nonlinear hydrostatic restoring force model is not able to capture this. Comparing Figure 12a-c shows that the threshold wave amplitude to trigger parametric resonance decreases as the wave frequency approaches twice the buoy's natural frequency (ω/ω 0 = 2). The results in Figure 13 show the normalised steady-state oscillation amplitude as a function of the input wave frequency, considering frequencies near, (a) half, (b) equal and (c) double the buoy's natural frequency. These frequencies correspond to the so-called (a) 1:2, (b) 1:1 and (c) 2:1 resonances, respectively. Note that the wave amplitude is constant for each frequency within a particular graph but differs between graphs due to the different threshold amplitudes required to trigger the 1:2 and 2:1 resonances (seen in more detail in Section 4.3). However, the oscillation amplitudes are normalised against the wave amplitude to allow a comparison between graphs. The analytical and nonlinear hydrostatic restoring force models are seen to give similar results in the vicinity of the 1:2 and 1:1 resonance frequencies. However, only the analytical model is able to capture the parametric resonance in the vicinity of the 2:1 resonance frequency.

Code-to-Code Verification against a NLFK Force Model
First, we compared the free response (no wave conditions) of the two models, see Figure 14a. We started the simulation from an initial displacement z(t = 0) = 4 [m], and the two models gave almost identical results. A small difference in the frequency of oscillations can be observed. This can be shown on the power spectral density plot of time series, see Figure 14b. Secondly, the response of the two models to the parametric was investigated, considering a range of input wave series in the frequency range near the 2:1 resonance. The major difference between the two models is the required computational time to run the simulation. Each input wave series is 3000 [s] in duration, requiring 2 [s] of computational time for the analytical model, compared to 14,000 [s] for the NLFK model (the simulation was run on a Core i7 7700HQ processor). This highlights a significant advantage of the analytical model, being 3-4 orders of magnitude faster than the NLFK model. Figure 15 shows an example of the predicted heave motions of the buoy for wave amplitude H = 2.2 [m] and frequency ω = 1.87 [rad/s]. The oscillation amplitude grows slightly faster and larger for the analytical model compared to the NLFK model. In addition, there appears to be a small phase difference between the response of the two models.  Figure 16a shows the normalised oscillation amplitudes, z max /H, for a range of wave amplitudes at a frequency of ω = 1.98ω 0 , comparing the results of the analytical, nonlinear hydrostatic and NLFK models. Due to the computational requirements, the NLFK results are obtained for less different wave amplitudes than the other two computationally efficient models. The results from the analytical model agree reasonably well with the NLFK model results. The main difference is in the oscillation amplitudes, which are approximately 15% higher in the analytical model compared to the NLFK model. Figure 16b illustrates the normalised oscillation amplitudes z max /H for varying wave frequencies in the vicinity of the 2:1 resonance for a fixed wave amplitude of 2.2 [m]. Once again, the oscillation amplitudes are slightly higher for the analytical model. In addition, the range of frequencies for which parametric resonance occurs is slightly larger for the analytical model.

Two-Parameter Analysis
The simulation speed of the new model with parametric excitation allows us to make multi-parameter investigations. In this subsection, the effect of the wave frequency, together with the amplitude on the buoy dynamics, is analysed. Figure 17 illustrates the buoy steady-state oscillation amplitudes z max as a function of the wave amplitude H and frequency ration ω/ω 0 . From Figure 17, we can conclude that the classical 1:1 resonance provides the largest steady-state oscillation amplitude for a given wave amplitude. To trigger the 2:1 resonance, large wave amplitudes are required. Below H = 1. 92 [m], the radiation damping mitigates the occurrence of parametric resonance. Indeed, the amount of damping in the system dictates the threshold at which parametric resonance occurs; therefore, an accurate measure of the total damping in the system is vital. The current model built upon the conventional linear hydrodynamic model that assumes an inviscid fluid and therefore neglects any type of viscous damping. However, the effect of viscosity becomes less negligible as the size of the buoy decreases, its velocity increases and/or its geometry become less streamlined. As such, a viscous damping term can be added to this model whose parameters can be identified from CFD simulations [105][106][107][108][109][110][111][112][113]. In addition, for the case of a WEC, the power take-off (PTO) will apply damping in order to extract energy. Therefore, compared to regular 1:1 resonance, the potential of parametric resonance for the present test case (whose geometry was chosen arbitrarily) seems limited. The efficient analytical model developed in this paper will allow future work to explore the vast space of possible geometries to examine the potential of exploiting parametric resonance and finding the appropriate relationships between the wave climates and the buoy dimensions.

Conclusions
In this paper, an analytical model for a heaving buoy with a position-varying waterplane area in monochromatic waves is introduced. The target application of the mathematical model is to evaluate the ability of heaving wave energy converters to leverage large amplitude oscillations triggered by the occurrence of parametric resonance due to a position-varying waterplane area. The introduced model extends the capability of current analytical models to be able to capture the occurrence of heave-to-heave Mathieu instabilities using an efficient nonlinear parametrically excited differential equation. The parametric excitation is captured by including position dependence to the excitation force parameters, which is achieved by fitting position-dependent functions (such as polynomials) to the excitation force coefficients calculated for various levels of buoy submergence. The performance of the model efficiently and accurately captures parametric resonance, as verified by a comparison against the results from a nonlinear Froude-Krylov force model, which numerically calculates the forces on the buoy at each time step based on the evolving wetted surface area. The introduced model is demonstrated to perform simulations in monochromatic waves nearly four orders of magnitude faster than the nonlinear Froude-Krylov force model simulations with comparable accuracy. The performance of the model is also compared against an analytical model that only considers the nonlinear effect of the position-varying waterplane area on the restoring force term rather than both the restoring and excitation force terms as the introduced model does. The results show that by only including the nonlinear restoring force term, the model is only able to capture the 1:2 and 1:1 resonances but is unable to capture the 2:1 resonance like the introduced analytical model and the nonlinear Froude-Krylov models can. In addition to the tremendous computational efficiency, by enabling many different WEC parameters to be evaluated using a brute force approach, the introduced model also allows analytical investigation techniques, such as perturbation or harmonic balance methods, to be utilised. These methods will be employed in future work to optimise the geometric design of a heaving WEC to exploit parametric resonance for increased power capture. Funding: The research reported in this paper and carried out at BME was supported by the ÚNKP-20-3 New National Excellence Program of the Ministry for Innovation and Technology of Hungary and the NRDI Funds (TKP2020 IES, Grant No. BME-IE-WAT; TKP2020 NC, Grant No. BME-NCS) based on the charter of bolster issued by the NRDI Office under the auspices of the Ministry for Innovation and Technology of Hungary. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 867453.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare that they have no conflict of interest.