Numerical simulation of parametric resonance in point absorbers using a simplified model

Parametric resonance is a non-linear phenomenon in which a system can oscillate at a frequency different from its exciting frequency. Some wave energy converters are prone to this phenomenon, which is usually detrimental to their performance. Here, a computationally efﬁcient way of simulating parametric resonance in point absorbers is presented. The model is based on linear potential theory, so the wave forces are evaluated at the mean position of the body. However, the ﬁrst-order variation of the body’s centres of gravity and buoyancy is taken into account. This gives essentially the same result as a more rigorous approach of keeping terms in the equation of motion up to second order in the body motions. The only difference from a linear model is the presence of non-zero off-diagonal elements in the mass matrix. The model is benchmarked against state-of-the-art non-linear Froude–Krylov and computational ﬂuid dynamics models for free decay, regular wave, and focused wave group cases. It is shown that the simpliﬁed model is able to simulate parametric resonance in pitch to a reasonable accuracy even though no non-linear wave forces are included. The simulation speed on a standard computer is up to two orders of magnitude faster than real time.


INTRODUCTION
Parametric resonance in an oscillating system is a resonance which is excited parametrically, that is, by a time variation of some parameter of the system, as opposed to some external excitation.A much-studied example is the parametric rolling of a ship moving in head or following seas.At certain wave encounter periods, the ship can experience amplified rolling motions with a period twice the encounter period.Because in head or following seas no roll exciting moment should be present due to the transverse symmetry of the ship, such motions cannot be a direct consequence of wave excitation, but are caused instead by periodic variations in the transverse stability of the ship [1].As the ship moves through the waves, the submerged portion of the hull varies with time, and this results in a periodic variation of the roll restoring stiffness.The problem of parametric rolling in ships is essentially that of a simple spring-mass system with a periodically varying spring constant.The equation of motion can be expressed as what is known as the Mathieu equation.This equation admits stable, unstable, or periodic solutions depending on the period and amplitude of the parametric excitation, that is, the variation of the spring constant (for the ship, this is the varying roll restoring stiffness due to the variation of the submerged portion of the hull).Given an initial disturbance, the response can be either stable, where we have decaying oscillatory response from the initial disturbance to rest; unstable, where we have oscillatory response with exponentially growing amplitudes with time; or periodic, where the amplitude of the oscillatory solution neither grows nor decays with time.In the absence of damping, solutions become unstable even with the slightest disturbance when the period of the parametric excitation equals some critical periods, the first of which is T = T n ∕2, that is, half the natural period of the spring-mass system.Away from these critical periods or in the presence of damping, solutions can still be unstable if the amplitude of the parametric excitation is sufficiently large.In a real physical system, however, the response will not grow unbounded due to the presence of non-linearities, which becomes non-negligible as the motions grow, and will limit the response amplitudes.
Floating wave energy converters (WECs) are in some respects similar to a ship without a forward speed, and we might expect that WECs are also prone to parametric resonance.Indeed, parametric resonance has been observed in model experiments of a variety of WECs, including floating oscillating water columns [2][3][4][5], self-reacting floating WECs of various configurations [6][7][8][9][10][11][12], and WECs tethered to the sea bed, either surface piercing [13] or completely submerged [14][15][16].A reduction in power capture usually accompanies the occurrence of parametric resonance, and thus the importance of avoiding it has been recognised for some time [17].
Numerical simulations of parametric resonance in WECs normally employ the so-called non-linear Froude-Krylov model, which is implemented as an option in state-of-the-art wave-body interaction software [18][19][20].This evaluates the hydrostatic and Froude-Krylov forces by integrating both the hydrostatic pressure and the pressure due to the undisturbed incident wave over the instantaneous wetted surface of the body [9,11,21].The other components of the hydrodynamic forces due to radiated and scattered waves are assumed to remain linear.Because the forces due to the radiated wave can be of the same order of magnitude as the forces due to the incident wave, especially when the body is in resonance, such an approach is adopted not for the reason of consistency, but rather for simplicity.It is difficult to judge the accuracy of such an approach because comparisons with model test measurements are relatively scarce.A rather poor agreement is reported in [9], whereas a better agreement is reported in [11].Going further than the non-linear Froude-Krylov model by including non-linear forces due to wave radiation would be a more consistent approach, but this is expected to be more computationally demanding.
The question is then whether a simpler approach is possible.For freely floating WECs, the mechanism that triggers parametric resonance is usually thought to be a time-varying metacentric height due to the heaving motion of the body relative to the free surface.This would modulate the roll or pitch restoring stiffness and excite the roll or pitch resonance.The variation of the metacentric height can be brought about by other modes of motion besides heave.In the SEAREV device [6,9], for example, the internal pendulum motion, which is induced by pitch of the body, results in a variation of the metacentric height of the body.This, in turn, induces parametric motion in roll.For tethered WECs, which are like inverted pendulums, it is the timevarying tether length due to the heaving motion of the body that is responsible for the occurrence of parametric resonance in roll or pitch.A model that includes these effects without considering any non-linear wave forces might be sufficiently accurate to simulate parametric resonance in WECs.Such a model would run much faster than real time and would be potentially useful for applications requiring real-time execution, such as control.
In the following, such a simplified approach is presented and illustrated by modelling a single freely floating axisymmetric body.The simplified model is based on the familiar linear model, but we consider the instantaneous, instead of the mean, positions of the centres of buoyancy and gravity up to the first order.This intuitive but less rigorous approach is shown to lead to essentially the same result as a more rigorous approach of keeping terms up to the second order in body motions in the non-linear equations of motion of the body.The final result is the presence of non-zero off-diagonal elements in the mass matrix such that heave and pitch are coupled.This is the only difference between the simplified model and a linear model.All wave forces are otherwise linear.Similar simplified approaches (i.e.without the use of non-linear wave forces) exist in the literature and have been used to model parametric resonances of ships [1], spars [22], and WECs [12,23,24], although we note that there are differences in these models.
Although hardly resembling a WEC, a freely floating axisymmetric body is chosen because of its simplicity to illustrate the method.The implication for practical WECs, such as selfreacting WECs consisting of two bodies moving in relative heave, is obvious.For such WECs, large pitch/roll angles will increase the friction between the two bodies, with a negative impact on power production [12].The present model can be easily adapted to model such WECs.
The accuracy of the model is assessed by comparing predictions of the body's response with those from state-of-the-art non-linear Froude-Krylov and computational fluid dynamics (CFD) models, for some selected cases, including free decay and excitation by regular wave trains and focused wave groups.Parametric resonance in moored floating buoys has been recently captured by CFD simulations [25,26].
We begin by outlining the properties of the floating body to be modeled.A detailed description of the simplified model is then given, followed by descriptions of the non-linear Froude-Krylov model and the CFD model used for comparison.The three simulation cases are then described and discussed in turn, including remarks on typical computation times.

PROPERTIES OF THE FLOATING BODY
The method is illustrated by considering the motions of a floating body which has the shape of an upright circular cylinder with a hemispherical bottom, as sketched in Figure 1.A right-handed coordinate system is used with the origin on the mean free surface, the x-axis aligned with the incident wave direction, and the z-axis pointing vertically upwards.The body is constrained  to move in three degrees of freedom (surge, heave, and pitch), which are denoted by the index j = 1, 3, 5.The translational modes (surge and heave) are defined as translations along the xand z-axes, respectively, whereas the rotational mode (pitch) is defined as rotation about the y-axis.The dimensions and properties of the body are listed in Table 1.At equilibrium, there are no vertical forces acting on the body except buoyancy and gravity.That is, the mass of the body is m = V 0 , where  is the seawater density and V 0 is the mean displaced volume of the body.The location of the centre of buoyancy follows from geometry, whereas the centre of gravity and the radius of gyration about the xor y-axis are specified.

Equations of motion
The equations of motion in the three degrees of freedom can be written in a matrix form as where M is the mass matrix, u = ξ is the body velocity vector, with  being the displacement vector, and F is the force vector.The latter is the sum of all forces acting on the body: where F e is the wave excitation forces, which are the sum of the Froude-Krylov forces and the forces due to the scattered waves, F r is the forces due to the radiated waves, F d is the hydrodynamic drag forces, F k is the restoring forces due to hydrostatics and gravity, F s is the restoring forces due to moorings, and F u is the power take-off (PTO) forces.For convenience, we have assumed that it is possible to decompose the hydrodynamic forces into F e , F r , and F d .If there are no incident waves, F e = 0, and if there is no PTO, F u = 0.The rest of the forces are non-zero as long as the body oscillates in water.

Wave forces
The wave excitation and radiation forces are assumed to be linear.That is, the wave excitation forces are assumed to be linearly proportional to the incident wave amplitudes, and the wave radiation forces to the body velocities.Furthermore, they are evaluated at the mean position of the body.
The wave excitation forces are therefore given as where (t ) is the incident wave elevation at a chosen location (here chosen to be the origin of the coordinate system), f(t ) is the vector of wave excitation impulse response functions, which are the inverse Fourier transforms of the complex excitation force coefficients f e (), and * denotes the convolution operation.Here, t is time and  is the angular frequency.
The wave radiation forces are given as [27] F r (t where A(∞) is the added mass matrix at infinite frequency and k(t ) is the radiation impulse response function matrix, whose elements are the inverse Fourier transforms of the elements of the radiation damping matrix B():

Hydrostatic and gravitational restoring forces
In a linear model, the hydrostatic and gravitational restoring moment in pitch for a floating body is (see, e.g.Section 6.16 in [28] or equation (B12) in the appendix) where  5 is the pitch displacement, S yy is the second area moment of the mean water plane about the y-axis, and z B0 and z G 0 are the mean vertical positions of the centres of buoyancy and gravity, respectively.The second equality follows from the fact that the body is freely floating, m = V 0 .
In the present model, instead of using the mean positions of the centres of gravity and buoyancy, we use the instantaneous positions of these centres.For the centre of gravity, this is, assuming small rotations (and as such, keeping terms up to the first order in pitch  5 ; see, e.g.[28,29] or Appendix A).Here, i and k are standard unit vectors in the direction of the xand z-axes.Thus, since x G 0 = y G 0 = 0. Furthermore, approximating the body as a uniform vertical cylinder, we can write Replacing S yy , V 0 , z B0 , and z G 0 in (6) with S yy + ΔS yy , V 0 + ΔV , etc. and neglecting terms containing the products of the variations give For small rotations, the term proportional to ΔS yy ∕S yy is much smaller than the other terms, and can therefore be neglected.However, the second area moment of the water plane increases as sec 3  5 , and ΔS yy ∕S yy > 1 for  5 > 37.47 degrees.This must be borne in mind as a possible limitation of the model.To be consistent with our assumption of small rotations, however, we neglect any term proportional to ΔS yy ∕S yy , and obtain the following expression for the pitch restoring moment: where we have made the substitutions ΔV = −S w  3 , Δz G =  3 (t ), and Δz B =  3 (t )∕2.Further simplification is possible by noting that V 0 ≈ −2z B0 S w , which upon substituting into (11) gives which is the same as the linear expression (6).
The hydrostatic restoring force in heave is modeled as Again, this is true for small rotation angles.The water plane area S w increases with pitch angle as sec  5 , and ΔS w ∕S w > 1 for  5 > 60 degrees.Figure 2 plots the functions sec  5 and sec 3  5 , which describe the growth of the water plane area and the second area moment of the water plane, respectively.The lowest-order approximations for ΔS w ∕S w and ΔS yy ∕S yy are both proportional to  2 5 .In Appendix B, we present a rigorous derivation of the hydrostatic and gravitational restoring forces for a freely floating body up to the second order in body motions (while neglecting second-order contributions of the free-surface elevation).This leads to the same results; see (B5), (B11)-(B13).Thus, (12) and ( 13) are valid up to second order in body motions if we do not include the contribution of the free-surface elevation.The exclusion of this contribution in the model is a potential source of inaccuracy.

Restoring forces due to moorings
As a simple mooring model, a horizontal spring with linear stiffness k s = 10 5 N/m is applied at the centre of gravity of the body.This induces restoring forces in surge and pitch: A rigorous derivation is given in the Appendix; see (B15) and (B17).The restoring force in heave due to the spring is assumed to be negligible.

Drag forces
A simplified lumped-parameter approach is used to model the hydrodynamic drag forces, which are assumed to be a quadratic function of the body velocities relative to the undisturbed fluid velocities: The drag coefficients, B d ,i are assumed to be constants, independent of the oscillation periods.The undisturbed fluid velocities, v i are calculated at the mean centre of gravity z G 0 .For deep water, the linear fluid velocity coefficients of proportionality in the horizontal and vertical directions are given as The area of the water plane increases with pitch angle,  5 as sec  5 , whereas the second area moment of the water plane grows as sec 3  5 .These variations are neglected in the simplified model where A is the (complex) incident wave amplitude and k is the wave number.

Power take-off forces
A linear PTO force is assumed to be applied to the heave mode.Thus, where B u is the PTO damping.The other elements of the PTO force vector, F u are assumed zero.

Mass matrix
The mass matrix is an important component of the present model.In a linear model, the non-zero elements of the mass matrix for the body are (see, e.g.Section 4.16 in [28] or Equation (A28) in the appendix) where r is the mean radius of gyration of the body about the xaxis.The M 35 = M 53 elements are equal to zero because x G 0 = 0 (see Table 1).Following the same approach as we use for the restoring forces, we replace the mean coordinates of the centre of gravity with their instantaneous values given in (7).However, since we express the equation of motion about the origin of the body system, the translations do not contribute.Thus, in the present model, Equations ( 20) and ( 21) are unchanged.The variation of the moment of inertia due to the variation of the radius of gyration is of a higher order, and is therefore neglected.Equation ( 22) thus also remains unmodified.However, instead of having M 35 = M 53 = 0, we have Notice that the heave displacement is now coupled to the pitch displacement through the non-zero M 35 = M 53 elements.
A rigorous derivation of the body inertia forces up to the second order in body motions leads to the same result, with an extra centripetal force term in heave (see (A29)), but this term is small and is neglected herein.

Summary and remarks
In sum, the simplified model is very much similar to a linear model.The only difference, apart from the introduction of drag, is the non-zero off-diagonal elements (23) in the mass matrix.
It should be noted that although the simplified approach, that is, taking into account the instantaneous positions of the centres of gravity and buoyancy (or the instantaneous length of the tether in the case of tethered devices) while keeping the wave forces linear is meant to be general in principle, the specific form of some of the expressions outlined above might differ from one WEC to another.The idea of the paper is more to verify the general principle by way of modeling a freely floating axisymmetric body as an example.

Computation of radiation impulse response functions, wave excitation forces, and undisturbed fluid velocities
A three-dimensional panel method [30] based on linear potential theory is used to compute the complex excitation force coefficients, f e (), the infinite-frequency added mass matrix, A(∞), and the radiation damping matrix, B().Due to the body symmetry, each of the A(∞) and B() matrix has only four unique non-zero elements.
The radiation impulse response functions k(t ) are precalculated according to (5).Because B i j () → 0 as  → ∞, a finite but sufficiently high upper integration limit,  u , may be chosen when calculating the integral.For a given panel size, the computed values of B i j () start to diverge for frequencies above a certain cut-off frequency,  c , which may be lower than  u .The values for  c <  ≤  u are therefore approximated using a fitting function in the form of c 1 exp(c 2 ) + c 3 exp(c 4 ), following the approach outlined in [31].The upper limit of integration is here chosen to be  u = 15 rad/s.The cut-off frequency is chosen to be  c = 5 rad/s for all elements of matrix B().
The wave excitation forces, because they do not depend on the body's response, are also precalculated.For a regular incident wave of frequency  and amplitude A, these are given as For an irregular wave train of significant wave height, H s , peak period, T p , and a given spectrum, S (), the incident wave elevation can be generated using the random coefficients method [32] as: where a n and b n are independent random values generated from a Gaussian distribution of zero mean and common variance,  2 n = S ( n )Δ.Here, N = t max ∕Δt is the number of values in the time series.Also,  n = nΔ, where Δ = 2∕t max .As an alternative to (3), the wave excitation forces corresponding to the incident wave elevation in (25) can be obtained as: where a n and b n are the same set of random values used in (25).Likewise, because the undisturbed fluid velocities v irequired for the calculation of the drag forces-depend only on the incident wave elevation, they can be precalculated in the same way as the wave excitation forces, that is, with the coefficients v x and v z in (17) and (18) replacing f e in (26).
The sums ( 25) and ( 26), as well as integral (5) (in the discrete form), can be efficiently evaluated using an inverse fast Fourier transform (FFT).The initial part of the wave excitation force and the fluid velocity time series is tapered using a cosine taper window.

NON-LINEAR FROUDE-KRYLOV MODEL
To assess the accuracy of the simplified model, the simulation results are compared with those of a non-linear Froude-Krylov model and a CFD model.
The non-linear Froude-Krylov model used for comparison in this study is the partially non-linear approach [33] implemented by the time-domain model WEC-Sim [18].The non-linear hydrostatic restoring and Froude-Krylov force components are calculated by integrating the static and dynamic pressures over each panel along the wetted body surface, which is discretised into triangular elements.Herein, the Wheeler stretching method is applied to correct the flow velocity and pressure due to the instantaneous wave elevation.The hydrodynamic coefficients including added mass, added damping, wave excitation, and restoring stiffness were obtained from a potential flow program [30].To model the PTO force in the heave mode only (19), a user-defined block is adopted for the current WEC-Sim model.Additionally, the mooring matrix is added to represent the linear surge stiffness.In order to take into account the drag forces, drag coefficients in surge, heave, and pitch are added.It should be noted, however, that all drag force components are calculated based on the body velocity instead of the relative velocity between the body velocity and the undisturbed fluid velocity as in the simplified model.
The origin of the coordinate system to define the body motions in WEC-Sim is at the centre of gravity instead of the mean free surface as employed in the simplified model.That is, a point with coordinates (x, y, z ) in the coordinate system defined in Figure 1 will have coordinates (x ′ , y ′ , z ′ ) = (x, y, z − z G 0 ) in WEC-Sim's coordinate system.The modes of motion are defined as translations along the x ′ -, y ′ -, and z ′ -axes, and rotations around the same axes.Thus, because pitch is defined in WEC-Sim as rotation about the y ′ -axis (as opposed to the yaxis), the surge and heave displacements in the two coordinate systems are related through the following transformations: where  ′ 1 ,  ′ 3 , and  ′ 5 are the surge, heave, and pitch displacements in WEC-Sim.To ease comparison, all responses from WEC-Sim are transformed to the same coordinate system defined in Figure 1.It should be mentioned, however, that because of the difference in the coordinate systems, the drag forces in the WEC-Sim model are different from those in the simplified model (the drag force model as defined in ( 16) is used in WEC-Sim without transforming the body velocities and without considering the fluid velocities).
As the non-linear hydrostatic restoring and Froude-Krylov force components are calculated from the discretised body surface, their accuracy relies on the body mesh resolution.Thus, a mesh convergence study is conducted to determine the proper mesh resolution.Figure 3 shows different mesh resolutions where the panel size ratio between two adjacent meshes is kept constant at √ 2. The panel sizes of the coarsest and finest meshes are 0.707 and 0.25 m, respectively.The proper mesh resolution is chosen to guarantee solution accuracy.In the present study, we conduct a mesh convergence study at a wave period that is close to the parametric resonance period, and the medium mesh size (0.5 m) is selected for further simulations.Unless stated otherwise, the result presented is based on the medium mesh.

CFD MODEL
High-fidelity CFD simulations are conducted using the opensource libraries provided by OpenFOAM [34], following an approach that was successfully utilised [35,36] as part of the CCP-WSI Blind Test Series [37,38] which considered a similar body to that in the present study.The basis of the approach is the interFoam solver [39], modified for wave generation and absorption [40], which solves the two-phase, incompressible, Reynolds-averaged Navier-Stokes (RANS) equations under the single fluid approximation: where p is the pressure, v is the fluid velocity, and g is acceleration due to gravity [39].The fluid density, , and dynamic viscosity, , are determined using the volume of fluid (VOF) interface capturing scheme [39].
The RANS equations ( 29)-(30) are solved using a first-order temporal scheme (Backwards Euler) and second-order spatial schemes (Central Differencing and Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) [41]).A variable time-stepping approach has been utilised based on a maximum Courant number of 0.5, with the pressurevelocity coupling achieved via the Pressure Implicit with Splitting of Operators (PISO) algorithm [42] using three correctors.The k- SST [43] turbulence closure model is used, which has been successfully demonstrated for other wavedriven flows [44].The motion of the device is captured using OpenFOAM's sixDoFRigidBodyMotion library, which uses a deforming mesh approach based on the Spherical Linear Interpolation (SLERP) algorithm [45], and the Newmark method [46] to determine the instantaneous position of a rigid body.The translational motion modes are defined as the displacement of the centre of mass (CoM) from its equilibrium position.The rotational modes are defined as the rotation about the CoM.For comparison with the simplified model, all responses from OpenFOAM are transformed according to (27) and (28).The deformation region is set to have an inner radius of 1 m, outer radius of 35 m and, for numerical stability, an acceleration relaxation of 0.9 is applied.
Wave generation is achieved via expression-based boundary conditions, which are supplied by the waves2Foam toolbox [40].For regular wave cases, linear wave theory is used, whereas the focused wave cases use a linear superposition of 100 first-order components.The relaxation zone technique [40] is utilised to absorb waves and reduce reflections within the numerical domain.
A linear stiffness restraint is applied along the x-axis based on the relative position of the CoM: where k s is the stiffness coefficient, and x G is the x-component of the CoM at the present time step.This force is used to apply an additional acceleration to the system at each time step using Newton's second law.A similar approach is used for the PTO damping (19), which is applied based on the vertical velocity of the body.

Numerical domain
The simulations are set up such that the waves propagate in the positive x-direction (Figure 4), with the coordinate and origin defined the same as in Figure 1.To increase computational efficiency, a symmetry condition is utilised at the transverse centre of the device (y = 0), and the simulation is constrained to three degrees of freedom (surge, heave, and pitch).
The body is truncated 10 m above the SWL.The water depth is set as 50 m to ensure that the generated waves are in the deepwater regime and to allow the mesh deformation algorithm to be accurate and robust at large deformations, whilst minimising the size of the computational domain.The numerical domain is 325 m long (−150 ≤ x ≤ 175 m), 50 m wide (0 ≤ y ≤ 50 m), and 100 m high (−50 ≤ z ≤ 50 m) (Figure 4).

Mesh discretisation
For the regular wave cases, the discretisation in the free-surface region (|z| ≤ 1.5 m) is set to Δz = 0.25 m (Figure 4a), which ensures that the wave is resolved within 1% of the theoretical solution [47].A cell aspect ratio of 1 is used in the working region (|x| ≤ 25 m), in the vicinity of the body (y ≤ 10 m) (Figure 4b).Mesh grading (indicated by double-headed arrows ↠ in Figure 4) is used to mitigate the computational cost by reducing the number of mesh cells: in the x-axis the mesh has constant discretisation (Δx = 0.25 m) for |x| ≤ 25 m, and it linearly increases up to Δx = 1.25 m for |x| > 25 m; in the yand zcoordinates (y > 10 m, |z| > 1.5 m), the mesh increases to Δy = Δz = 2.5 m (Figure 4).The focused wave cases are set up in a similar way but with additional vertical refinement around the free surface (Δz = 0.0625 m; Δx = 4Δz).The body is meshed using snappyHexMesh to one level of octree refinement finer than the free surface (Δz = 0.125 m for regular waves; Δz = 0.03125 m for focused waves), and based on preliminary simulations, additional layers are applied around the body to ensure that the dimensionless wall distance, y + , is less than 300, ensuring that the wall functions are applied within the appropriate region of the boundary layer.This leads to a total mesh size of 7.25M cells for the regular cases and 11.8M cells in the focused wave cases.Relaxation zones based on an exponential weighting function [40] are used at both the inlet and outlet boundaries with target solutions of linear wave theory, and still water conditions, respectively.These zones are 75 and 100 m in length, respectively (Figure 4).Therefore, depending on the case, the outlet relaxation zone will be 0.75-1.25 wavelengths, and is estimated to produce a reflection coefficient of less than 0.1% at the outlet boundary [40].

Boundary conditions
The 'Bottom', 'Outlet', and 'Side Wall' boundaries (see Figure 4) are considered as solid walls and hence no-slip conditions have been applied along with zero-gradient conditions for pressure and VOF.A zero-gradient condition is also applied on the 'Atmosphere' for the VOF but the velocity boundary condition varies according to the near boundary flux: using a zero-gradient condition for outflow and the internal cell value of the normal component to the patch face for the inflow.The 'Atmosphere' boundary condition for pressure is defined as the total pressure: where p 0 is the user defined reference value, and v is the velocity.For this case p 0 is set to zero because the solver uses the difference between total pressure and hydrostatic pressure.The 'Inlet' boundary uses expression-based wave generation boundary conditions for the VOF and velocity.The body boundary is modelled as a rigid, moving wall with zero-gradient applied for the VOF and a consistent no-slip condition for the velocity and pressure.Continuous wall functions are applied on the 'Side Wall', 'Bottom', and body boundaries, which switch between low and high Reynolds number flows, depending on whether the near wall cell centre lies in the log or laminar sublayer.The other boundaries use zero gradient for the turbulence model parameters, except at the 'Inlet' and 'Atmosphere' boundaries where the turbulent kinetic energy is zero gradient for outflow and fixed as for inflow, where I and c p are the turbulence intensity and phase speed of the wave, respectively [44,48].The dissipation rate, , is then adjusted such that the eddy viscosity is a fraction, , of the kinematic viscosity, that is,  t = .Following Lin and Liu [48], I and  are chosen as 0.0025 and 0.1, respectively, in this study.The initial conditions are set as still water and the turbulence parameters specified at the inlet.

RESULTS AND DISCUSSION
Only motions in surge, heave, and pitch are considered in this study.The parameters of the model are given in Table 1.
Simulations for the simplified model are performed using a fourth-order Runge-Kutta scheme with a constant time step of 0.02 s.The convolution part of the wave radiation force (4) is evaluated once for each time step.

Decay
For the first case, we study the response of the floating body to a given set of initial displacements, in otherwise still water (in the absence of any exciting forces), with the PTO damping set to zero, but with the mooring force applied.
From the Fourier spectra of the decay time histories, the natural periods of the body can be obtained.These are listed in Table 2.The pitch natural period, 19.0 s, is close to twice the heave natural period, 7.9 s.The natural periods obtained from eigenvalue analysis of the linear undamped equations of motion are also presented in the table.Surge and pitch as defined here are coupled, even in the linear case.Thus, depending on the initial condition, the surge decay may contain oscillations in the pitch natural period, and vice versa.Based on the instability condition of the Mathieu equation, we expect parametric resonance in surge/pitch to be triggered by heave oscillations at around 9.5 s.
To ensure that the simplified model and WEC-Sim have consistent setup, the linear free-decay displacements obtained from the two models (i.e. with non-linearities switched off) are first compared.They are shown to be in good agreement (Figure 5), confirming that the two models have consistent setup despite the different coordinate systems used.
Having ensured that the two models are consistent, we then run the models with non-linearities included and compare them also with OpenFOAM.In contrast to the linear case, heave is coupled to pitch in the non-linear case.Thus, free-decay oscillations in pitch can excite heave (as seen from Figure 6h), unlike in the linear case.In addition, free-decay oscillations with initial displacements in both heave and surge/pitch can excite surge and pitch oscillations with amplitudes greater than their initial displacements (not shown).As seen from Figure 6, non-linear responses obtained from the three models show more prominent differences as non-linear forces come into play.Because there is no incident wave, differences should arise only from restoring and drag forces and, with OpenFOAM, also radiation forces.The smaller amplitudes of the response predicted by the simplified model suggest that it has higher damping, presumably due to larger drag forces.Some phase differences are also evident.In Figure 6c), the pitch response predicted by WEC-Sim appears to oscillate at a shorter period than those predicted by the simplified model and OpenFOAM, which seem to agree in phase, despite having different decay rate.From Figure 6e, we see that the heave response predicted by OpenFOAM has higher amplitudes than the response predicted by WEC-Sim or the simplified model, presumably due to reflections from the boundaries of the numerical domain.

Regular wave excitation
We now compare the response of the body to regular wave excitation of various periods (T = 7, 8, 9.5 s) and wave heights (H = 0.5, 1 m).The PTO damping is set to 2 × 10 4 kg/s.At T = 7 s, no parametric resonance is observed.The displacements from the simplified model agree closely with those from WEC-Sim for both incident wave heights.For H = 0.5 m, OpenFOAM gives a similar response, but the motion envelopes display small oscillations, presumably arising from some reflections from the boundaries of the numerical domain.
Parametric resonance is observed for OpenFOAM at T = 8 s and H = 0.5 m (see Figure 8).Here, the surge and pitch responses start off small and oscillate with a period equal to the incident wave period, but, after some time, they increase gradually to much larger amplitudes, oscillating at twice the wave period.At this wave height, parametric resonance is not predicted by WEC-Sim and the simplified model, at least for the simulated duration.When the incident wave height is increased to 1 m, parametric resonance is observed for both WEC-Sim and the simplified model.The simplified model predicts similar response amplitudes to those predicted by WEC-Sim, for all the three modes.
At T = 9.5 s, H = 0.5 m, parametric resonance is not observed in the response predicted by WEC-Sim within the duration of the simulation.Both OpenFOAM and the simplified model predicts parametric resonance towards the  (d-f) initial heave displacement of 1.5 m, and (g-i) initial pitch displacements of 10 degrees about the centre of gravity: simplified model (blue), WEC-Sim (red), OpenFOAM (yellow) FIGURE 7 Surge, heave, and pitch displacements of the body in response to incident regular waves of two different wave heights (H = 0.5 m and H = 1 m) and wave period T = 7 s: simplified model (blue), WEC-Sim (red), OpenFOAM (yellow).The PTO damping is 2 × 10 4 kg/s.The incident wave amplitude is gradually increased from zero for the first 200 s using a cosine tapering window FIGURE 8 As in Figure 7, for wave period T 8 s FIGURE 9 As in Figure 7, for wave period T = 9.5 s end of the 800-s simulation period, with the simplified model predicting it slightly later than OpenFOAM.For H = 1 m, parametric resonance is predicted by both WEC-Sim and the simplified model, with the simplified model predicting it slightly earlier and with response amplitudes that are slightly higher than those predicted by WEC-Sim.
In all cases, the response predicted by the simplified model practically agrees with that predicted by the non-linear Froude-Krylov model (WEC-Sim) up to the time when parametric resonance takes place.Some differences are observed as soon as parametric resonance takes place, but the fact that the simplified model predicts parametric resonance at wave periods where OpenFoam or WEC-Sim predicts it, with similar response amplitudes to those predicted by WEC-Sim or OpenFOAM, is encouraging, given its simplicity.
The simulation results confirm that the degree of parametric resonance in surge/pitch is tied with the vicinity of the incident wave period to half the pitch natural period, as well as the amplitude of parametric excitation available at this wave period.The closer the wave period is to half the pitch natural period and the larger the heave oscillations are at this wave period, the more likely and the more immediate parametric resonance is to happen.
The normalised displacement amplitudes obtained from the simplified model are plotted in Figure 10 as functions of the incident wave period.These are obtained from regular wave simulations with different discrete periods (with an increment of 0.3 s around T = 9.5 s, an increment of 0.5 s around T = 19 and 27.5 s, and an increment of approximately 1 s elsewhere) and three different wave heights, making up a total of 159 regular wave simulations.The displacement amplitudes are nor-malised by the incident wave amplitude A for surge and heave, and by A∕a for pitch.The subharmonic and superharmonic amplitudes are also plotted alongside the wave-frequency amplitude.These amplitudes are obtained by least-squares fitting sinusoidal functions with frequency equal to the wave frequency, half the wave frequency, and twice the wave frequency, to relevant steady-state portions of the displacement response time series.For the heave displacement, the wave-frequency, subharmonic, and superharmonic amplitudes are all obtained by fitting sinusoidal functions to a portion towards the end of the heave displacement time series.For the surge and pitch displacements, the subharmonic and superharmonic amplitudes are obtained from the end of the time series, whereas the wave-frequency amplitudes are obtained from the end if no parametric resonance occurs; otherwise, they are obtained from the beginning of the time series.
Three sets of curves are plotted for each of the subharmonic and superharmonic responses, corresponding to different incident wave heights, H = 0.5, 1, and 2 m.For the wave-frequency response, a curve corresponding to the response obtained from the linear model (but with drag) for H = 0.5 m is also plotted.The wave-frequency responses peak as expected around the natural periods identified earlier (see Table 2).Significant subharmonic response, corresponding to parametric resonance, is evident for surge and pitch at around T = 7.9 s (the heave natural period) and T = 9.5 s (half the surge/pitch natural period).
When parametric resonance happens, we see that as the subharmonic response in surge/pitch amplifies, the wavefrequency response in heave reduces relative to the linear response.It is clear how this would be detrimental to power FIGURE 10 Normalised displacement amplitudes of the body: wave-frequency response (blue), subharmonic response (red), superharmonic response (yellow); for three wave heights: H = 0.5, 1, and 2 m (thinnest to thickest).The PTO damping is 2 × 10 4 kg/s.Also plotted is the response from a linear model with drag (dashed black).Vertical dotted lines indicate (from left to right) T = 7.9, 9.5, 19, and 27.5 s, which are related to the system's natural periods (see Table 2) absorption for devices whose absorbing mode is heave.An additional observation is that parametric resonance is relatively narrow banded, with increasing bandwidth as the wave amplitude is increased, in agreement with the Mathieu equation's stability diagrams (see [15] or [16]).We also observe that the normalised amplitudes, including those of the subharmonic surge and pitch motions, are reduced as the incident wave height is increased, which is expected due to the greater drag forces, which increase quadratically with motion amplitudes.

Wave group excitation
Rather than simulating the response of the body to a long irregular wave record, for the third case we study the response of the body to wave group excitation, in the form of a scaled autocorrelation function of some chosen wave spectrum S (), representing the average (linear) shape of large wave events in the record [49].That is, where  is the surface elevation at t = t 0 and The integral (34) as well as the corresponding exciting force vector may be written in discrete forms as in (25) and (26), FIGURE 11 Wave spectra corresponding to the focused wave groups (with T p = 7, 8, 9, 10, and 20 s).Vertical dotted lines (from right to left) indicate T = 7.9, 9.5, 19, and 27.5 s, which are related to the system's natural periods (see Table 2) with S ( n )Δ n ∕m 0 replacing (a n − ib n ), and similarly evaluated using an inverse FFT.
The use of a wave group as input facilitates comparison with CFD, for which a long irregular wave simulation would be computationally prohibitive.At the same time, it allows responses under broad-banded excitation to be studied, since a wave group can be thought of as a compact representation of an irregular sea state (as it is generated from the same underlying spectrum).As our interest is in comparing the response of the body to a deterministic train of irregular waves, only isolated wave groups are used as input to the simulation.If one is interested in obtaining the correct statistical properties of the response in a true random sea state, a wave group embedded into a random train of background waves could be used [50].
The focused wave groups are generated using a Pierson-Moskowitz (PM) spectrum with different peak periods T p .The spectra are shown in Figure 11 and the time series of the incident surface elevation at the focused location x = 0 are shown in Figure 12.The amplitude at focus is  = 1 m and the focus time is t 0 = 200 s.Each simulation is run for 400 s.
The simulated responses of the body to the focused wave groups are shown in Figure 13a-c.We can see that although the wave groups are very localised (in time), the responses are long lived.The surge and pitch amplitudes predicted by the simplified model after the wave group passage are generally much smaller than those predicted by WEC-Sim and OpenFOAM.This could be due to the larger drag forces in the simplified model and/or the exclusion of the contribution of the freesurface elevation in the hydrostatic restoring forces in the simplified model.An interesting behaviour is observed for the heave response predicted by OpenFOAM, where the response increases after an initial decay.This behaviour is not seen with the simplified model and WEC-Sim, and hence is thought to be due to imperfect absorption within the CFD numerical domain, particularly on the side walls because there is no absorption applied in the transverse direction.
For comparison, the response of the device from a linear model with drag is also shown in Figure 13d-f (the response from a linear model without drag is similar).The linear heave response is similar to the non-linear one.However, the linear surge and pitch responses after the wave group passage are smaller than the non-linear ones.In particular, there is very little response in surge and pitch after the wave group passage for T p = 7 to 10 s, in contrast to the non-linear surge and pitch responses, which persist after the wave passage.The wave spectra for T p = 7 to 10 s are well away from the pitch natural period (as indicated by the second vertical line from the left in Figure 11), and so there is no mechanism, linearly, to excite these motions.There is, however, quite significant energy content at the heave natural period for these spectra, which excite the heave resonance.The significant motions in heave in turn excite the pitch/surge motions in a non-linear, parametric, fashion.

Computation time
For the simplified model, simulations are run with an Intel Core i7-7700 CPU @ 3.60 GHz, 64 GB RAM.A 800-s long simulation with a time step of 0.02 s takes approximately 8 s to complete.This is equivalent to approximately 5000 time steps per second.
For the non-linear Froude-Krylov model (WEC-Sim), which is based on body surface discretisation, a computer specification with Intel Core i7-7820 CPU @ 2.9 GHz, 32 GB RAM is utilised to run the simulations.A 800-s long simulation with a time step of 0.05 s takes approximately 33 s to complete, equivalent to approximately 500 time steps per second.
For the CFD model (OpenFOAM), all simulations were run on 64 cores of the in-house high performance computing service at the University of Plymouth.This facility consists of 52 2U Twin Sq. (4 nodes) networked with Intel Omni-Path cabling, and equipped with dual Intel E5-2683v4 8 core 2.5 GHz processors with 128 GB of memory per motherboard.For the regular wave cases, the simulations had an execution time of 60-85 h, whereas each focused wave case took 65-70 h.

CONCLUSION
Floating bodies of certain mass distributions are prone to parametric resonance in pitch/roll.WECs are no exception.Here, a simplified model to capture parametric resonance in pitch has been presented and applied to model a hypothetical WEC in the form of a floating axisymmetric body, where power is taken from heave only.We hypothesised that a simple model which includes a heavemodulated pitch/roll restoring stiffness without including any non-linear wave forces is able to simulate parametric resonance in WECs with sufficient accuracy.Contrary to our initial hypothesis, our derivation shows that linear expressions for the hydrostatic and gravitational restoring forces remain valid at second order in body motions, if we ignore variations of the free-surface elevation.Compared to a linear model, the only modification to the equation of motions arising from keeping terms up to the second order in body motions is the introduction of non-zero off-diagonal elements in the mass matrix.
We confirm that our model is capable of predicting parametric resonance in pitch at the wave periods where it is expected to occur, as well as predicting the associated reduction in heave response at these periods.The accuracy of the model has been assessed by comparing the simulated response of the body in surge, heave, and pitch with those obtained from state-ofthe-art non-linear Froude-Krylov and CFD models, generally showing satisfactory agreement.
With a typical simulation speed of two orders of magnitude faster than real time, the proposed model offers a computationally efficient means to simulate parametric resonance in the time domain.
The presented numerical simulations are yet to be validated by physical model tests.It will then be necessary to extend the analysis to six degrees of freedom.This is expected to introduce additional coupling between modes of motion which do not exist in the present model.Our purpose here is to rigorously derive the expressions of the body inertia forces for different orders of approximations in the body motions.We will show that the expressions given in ( 20)-( 23) agree with the second-order expression of the body inertia forces.We start with the general form of equations for rigid body motion, as given, for example, in [51].For completeness, the expression is derived again here.As in the main body of the paper, we restrict our treatment to surge, heave, and pitch motions.

How to cite this article:
It is useful to define a body-fixed coordinate system which moves with the body, but coincides with the inertial coordinate system (Figure 1) when the body is at rest.Thus, we define r b G as the coordinates of the centre of gravity in the body frame, that is, We make no distinction between the centre of gravity and the centre of mass.
The coordinates of a point on the moving body, expressed in the inertial frame, are given as where is the translation of the origin of the body system relative to the inertial system, r b is the coordinates of the point in the body frame, and R is the rotation matrix defined as with If pitch is the only rotational mode, then  4 =  6 = 0, and both R x and R z become identity matrices, and so R = R y .
The velocity of a point on the body, expressed in the inertial frame, is obtained by differentiating (A2) with respect to time, giving This is the general form of equation for an arbitrary rigid body.We can expand (A25) in terms of the translations and rotations and their derivatives and obtain expressions for different orders of approximations.The expressions up to the first order are and from which the familiar expression for the mass matrix in linear theory (see, e.g.Section 4.16 in [28]) is obtained: For a body with x G 0 = y G 0 = 0 and with x b = 0 and y b = 0 planes of symmetry, some of the off-diagonal elements in this matrix vanish.If we also exclude sway, roll, and yaw as we have done throughout, we end up with (20)-( 22) as the only non-zero elements.
The general second-order expressions corresponding to (A26) and (A27) are more involved, but they simplify considerably if the body has x b = 0 and y b = 0 planes of symmetry and x G 0 = y G 0 = 0, and if we also exclude sway, roll, and yaw.The result can be written as where the elements of the mass matrix M are as in the linear case, but with instead of zeros.The mz G 0 ξ2 5 term is the centripetal force due to pitch and can be moved to the left-hand side.This term is generally much smaller than the other force contributions.

B.1
Hydrostatic and gravitational restoring forces The restoring force due to hydrostatic pressures and gravity is given as where S b is the wetted surface of the body and n is a unit normal vector pointing outwards of the body.
In the following derivation, we assume that the surface integral is taken up to the mean free surface, z = 0, instead of the actual free surface z = .At second order, the free-surface elevation  will contribute to the integral, but we neglect this contribution for simplicity.We can therefore add a surface bounded by the intersection of the body and the z = 0 plane to the integral in (B1), so it becomes an integral over a closed surface.We can then express this surface integral as a volume integral by means of Gauss's theorem, such that The integral over the displaced volume of the body, ∀, can be expressed as a sum of an integral over the mean displaced volume, ∀ 0 , and an integral over a volume bounded by the z b = 0 and z = 0 planes and the surface of the body.For a cylindrical body, we can write where S w is the mean water plane area of the body and z b w is the z b -coordinate of a point on z = 0.This is obtained by setting z in (A2) to zero and solving for z b in terms of x b , y b , and the body displacements.Different orders of approximation for F k can then be obtained based on Taylor series expansion of z b w .The first-and second-order expressions turn out to be equal and are given as The restoring moment is given as where we have again expressed the surface integral as a volume integral and made use of (A2).Splitting the integral as we did before, we arrive at We can expand (B9) to any order of approximation.If we write where S x , S xx , etc. are the area moments of the mean water plane.For a body with x b = 0 and y b = 0 planes of symmetry, only the first term of each of these expressions remains.Equations (B5) and (B11)-(B13) are the familiar equations for the linear hydrostatic and gravitational restoring forces for an arbitrary freely floating body; see, for example, Section 6.16 in [28].The corresponding second-order expressions are more involved, but they reduce to the first-order expressions if the body is freely floating and if we also exclude yaw.

B.2
Restoring forces due to mooring The restoring force due to the horizontal spring, applied at the centre of gravity, is given as The first-order approximation, excluding yaw, is The second-order approximation, for a body with x G 0 = y G 0 = 0, and excluding yaw, is equal to the first-order one.
The restoring moment due to the horizontal mooring is given as The first-order approximation, for a body with x G 0 = y G 0 = 0, is The second-order approximation, if we also exclude yaw, is equal to the first-order one.

FIGURE 1
FIGURE 1 Sketch of the floating body

FIGURE 3
FIGURE 3 Four different mesh resolutions of floating body used for non-linear WEC-Sim simulation (left to right: coarse to extra fine)

FIGURE 4
FIGURE 4 Scale diagram of the numerical setup used in the OpenFOAM simulations: (a) side view, (b) top view.Information in red denotes mesh properties with double-headed arrows (↠) showing the direction of increasing mesh grading

Figures 7 to 9 FIGURE 5
FIGURE 5Linear free-decay displacements of the body without any PTO damping and with (a-c) initial surge displacement of 0.5 m, (d-f) initial heave displacement of 1.5 m, (g-i) initial pitch displacements of 10 degrees about the centre of gravity: simplified model (blue), WEC-Sim (red)

FIGURE 6
FIGURE 6Non-linear free-decay displacements of the body without any PTO damping and with (a-c) initial surge displacement of 0.5 m, (d-f) initial heave displacement of 1.5 m, and (g-i) initial pitch displacements of 10 degrees about the centre of gravity: simplified model (blue), WEC-Sim (red), OpenFOAM (yellow)

FIGURE 12 FIGURE 13
FIGURE12 Focused wave groups corresponding to the spectra in Figure11.The peak periods T p are as indicated and the focus amplitude is 1 m

TABLE 1
Dimensions and properties of the floating body

TABLE 2
Natural periods of the floating body identified from non-linear decay simulations and from linear undamped model (in parentheses)