Fast nonlinear Froude–Krylov force calculation for prismatic floating platforms: a wave energy conversion application case

Computationally fast and accurate mathematical models are essential for effective design, optimization, and control of wave energy converters. However, the energy-maximising control strategy, essential for reaching economic viability, inevitably leads to the violation of linearising assumptions, so the common linear models become unreliable and potentially unrealistic. Partially nonlinear models based on the computation of Froude–Krylov forces with respect to the instantaneous wetted surface are promising and popular alternatives, but they are still too slow when floaters of arbitrary complexity are considered; in fact, mesh-based spatial discretisation, required by such geometries, becomes the computational bottle-neck, leading to simulations 2 orders of magnitude slower than real-time, unaffordable for extensive iterative optimizations. This paper proposes an alternative analytical approach for the subset of prismatic floating platforms, common in the wave energy field, ensuring computations 2 orders of magnitude faster than real-time, hence 4 orders of magnitude faster than state-of-the-art mesh-based approaches. The nonlinear Froude–Krylov model is used to investigate the nonlinear hydrodynamics of the floater of a pitching wave energy converter, extracting energy either from pitch or from an inertially coupled internal degree of freedom, especially highlighting the impact of state constraints, controlled/uncontrolled conditions, and impact on control parameters’ optimization, sensitivity and effectiveness.


Introduction
Wave-structure interactions are complex physical phenomena, often difficult to model and predict with a degree of confidence satisfactory for several non-trivial ocean engineering applications. One field presenting major modelling challenges is wave energy conversion, for two main reasons: on the one hand, accurate mathematical models are essential for the design and development of economically performing wave energy converters (WECs); on the other hand, due to the inherent requirement of exaggerating the device motion to increase power capture, simple linear models are often inadequate, demanding consideration of an appropriate compromise between higher-complexity and computational time. The importance of accurate mathematical models for WECs has been widely stressed, even in recent years. Not only a reliable power production assessment is essential for correctly evaluating the performance of the device, but also the prediction of structural loads is fundamental for parsimonious but resilient design, for which high model accuracy is mandatory . Moreover, the control strategy, playing a fundamental role in modifying the response of the system, increasing power absorption while ensuring compliance with constraints , is very sensitive to modelling errors . Furthermore, particularly due to the action of the energymaximising controller, linear models may become not only inaccurate, but even physically absurd, with floating bodies clearing the water if no constraint is applied (Faedo et al. 2020b). In general, under controlled conditions, the inclusion of nonlinearities should always be at least considered.
Due to the established awareness of the need for nonlinear models, a plethora of modelling options is available (Wolgamot and Fitzgerald 2015;Penalba et al. 2017; Davidson and Costello 2020). However, what model for what device and what application condition remains an open question, so that recent great collaborative effort is leading to a wide crossmodelling comparison (Wendt et al. 2019;Ransley et al. 2020). Each model relies on a set of assumptions, more or less restrictive, with consequent impact on the complexity and computational burden. The resulting accuracy depends on how representative the assumptions are of the expected working conditions. In particular, while fully-nonlinear models may be of the highest accuracy, the time they require for computation is not compatible with control, optimization, or extensive power production assessment. On the contrary, partially-nonlinear models have the potential to define a convenient compromise between computation and fidelity, since only major nonlinearities are described. Among partially nonlinear models, nonlinear Froude-Krylov force (NLFK) models appear to be the most attractive for floaters relatively small compared to the characteristic wavelength, as WECs usually are: since radiation and diffraction effects are relatively small, only the static and dynamic undisturbed pressure field, known regardless of the floater motion, is integrated over the instantaneous wetted surface, resulting from the consideration of the floater displacement with respect to the actual free surface elevation.
The appeal of NLFK models is certified by the large number of independent numerical implementations witnessed in recent years. The vast majority of such models relies on a discretized spatial representation of the external surface of the floater, by means of planar mesh panels. The consideration of the instantaneous wetted surface requires either a very fine mesh or a re-meshing routine every time steps, both solutions being computationally demanding. The pioneer of NLFK models in the WEC field is LAMSWEC, developed for the nonlinear modelling of the SEAREV device , and applied to a heaving sphere (Merigaud and Gilloteaux 2012) and the WaveBob device (Tarrant and Meskell 2016). The currently widespread simulation software WEC-Sim (Lawson et al. 2014) also includes a mesh-based NLFK computation feature, which has been applied to a bottom-referenced point absorber (Tosdevin et al. 2019), to a self-referenced 2-body heaving point absorber (Van Rij et al. 2018), and a multi-body WEC (Chandrasekaran and Sricharan 2020). Further examples of software including NLFK computation via a mesh are SIM-DYN (Somayajula and Falzarano 2015;Wang et al. 2019), as well as various other in-house implementations (Guerinel et al. 2013;Kim 2019, 2020). Such NLFK models are already faster than weakly-or fully-nonlinear potential flow models (Penalba et al. 2017), where the potential problem is solved for each time-domain simulation at each time step on a continuously updating wetted surface (Letournel et al. 2014).
In all of such NLFK implementations, the computation bottleneck is the mesh discretization, leading to simula-tions about 2-orders of magnitude slower than real time (Gilloteaux 2007). Although for geometries of arbitrary complexity a mesh is likely the only possible solution, simpler shapes can offer symmetries that allow to reach an analytical representation of the wetted surface, hence tearing down the computational time. Extensive work has been recently done for axisymmetric floaters, which are very popular geometries for point absorbers WECs and spars. In case of purely heaving devices, an algebraic solution for the NLFK force can be obtained (Giorgi and Ringwood 2018a), with evident computational advantage (simulations in about two orders of magnitude faster than real time (Giorgi and Ringwood 2018b)). For multiple degrees of freedom (DoFs) simulations, a numerical solution is still required (Giorgi and Ringwood 2018c) but, since the wetted surface is represented analytically, the computational saving with respect to meshbased approaches is significant (about real time computation (Giorgi et al. 2020d)). This modelling approach is shown able to articulate parametric resonance (Giorgi et al. 2020c) and instabilities (Giorgi et al. 2020a), and has been validated via comparison with wave tank experiments (Giorgi et al. 2020b). An open source demonstration toolbox is also available (Giorgi 2019).
The computationally efficient analytical formulation of NLFK forces for axisymmetric floaters is particularly useful for the class of WECs exploiting mainly the heaving motion, absolute or relative, since such devices are usually designed to be insensitive to the wave direction, hence axisymmetric. This paper has the purpose to similarly exploit the geometric simplification of prismatic floaters in order to provide a computationally efficient alternative to general mesh-based approach. Prismatic floaters are often used for the class of WECs extracting energy from the pitching motion, as (Nemos GmbH 2020;WavePowerLab 2020;Heras et al. 2019;Ma et al. 2020;Bracco et al. 2014;Poguluri et al. 2019;Sirigu et al. 2020c;Kurniawan and Zhang 2018;Sirigu et al. 2019;Marcollo et al. 2017;Wu et al. 2017), that can benefit from an analytical formulation of NLFK forces.
The novel contribution of this paper is to provide a fast NLFK force computation framework, based on the analytical representation of the wetted surface of generic prismatic floaters. State-of-the-art mesh-based NLFK models, although able to handle geometries of arbitrary complexity, are too slow to be applied for computationally demanding tasks, such as power production assessment and control; moreover, mesh-based approaches are less transparent then fully-analytical models: based on a closed-form analytical description of the wetted surface, it is possible to apply advanced model-order reduction techniques (Faedo et al. 2020a) and geometry-based optimizations (Garcia-Teruel et al. 2020). The modelling approach has been validated against experimental data for an application based on an axisymmetric wave energy converter (Giorgi et al. 2020b), providing confidence that such a mathematical framework is appropriate to describe the dynamics also of a pitching wave energy converter under operational conditions, excluding extreme events.
The remainder of the paper is organized as follows: Sect. 2 presents the definition of nonlinear Froude-Krylov forces and the mathematical formulation of the computationally efficient computation for generic prismatic floaters; Section 3 presents the case study considered in this paper to quantify different nonlinear effects and their impact on the behaviour and performance of the WEC; in particular, two scenarios are considered, where energy is directly extracted from the pitching motion, in Sect. 3.1, or indirectly via inertial coupling, in Sect. 4.3. Results are discussed in detail in Sect. 4, while some concluding remarks are given in Sect. 4.

Nonlinear Froude-Krylov force formulation
Assuming homogeneous, ideal, incompressible fluid, and irrotational flow, a velocity potential ϕ can be defined (Newman 1977), as follows: where u is the velocity vector. Since the motion is assumed incompressible, ϕ satisfies Laplace's equation throughout the fluid described as follows: If a floating body is present in the fluid, the potential flow can be decomposed as the sum of the undisturbed incident flow potential (ϕ I ), the diffraction potential (ϕ D ), and the radiation potential (ϕ R )as follows: The pressure p can be derived from the potential flow applying Bernoulli's equation as follows: where ρ is the water density, g the acceleration of gravity, and p the total pressure, divided into static (−ρgz) and timevarying components. In a fully linear approach, the time-varying pressure is notionally integrated over the (constant) mean wetted surface of the floater, generating Froude-Krylov, diffraction, and radiation forces, respectively, for ϕ I , ϕ D , and ϕ R components of ϕ. Conversely, the balance between the static pressure and the weight of the floater is computed as if the cross-sectional area of the floater at the still water level (SWL) was constant, hence leading to a hydrostatic restoring force linearly proportional to the displacement. However, such a linearization is acceptable only under small steepness and small amplitude of motion assumptions, which are normally violated by WEC under control conditions, eventually leading to large inaccuracies or even unrealistic results (Faedo et al. 2020b).
Partially nonlinear models can preserve part of the computational convenience of linear models, while increasing the accuracy. If the characteristic dimension of the floating body is much smaller than the wavelength, ϕ D and ϕ R are much smaller than ϕ I (Clément and Ferrant 1988), so that the total wave field is sufficiently well approximated by the undisturbed pressure field. The NLFK force, accounting for the greater portion of the total excitation force, can be computed as the integral of the undisturbed pressure ( p u ) field over the instantaneous wetted surface (S w (t)).
Assuming two-dimensional waves in the (x, z) coordinate system, where x is the direction of propagation of the wave, and z is the vertical axis, positive upwards, with the origin at the (SWL), a the wave amplitude, ω the wave frequency, k the wave number, h the water depth, and z the vertical coordinate modified according to Wheeler's stretching (Giorgi and Ringwood 2018d), p u follows: where the second term at the right-hand-side is the dynamic pressure. Note that Wheeler's stretching is not a strict theory, since it implies that the potential function does not satisfy the Laplace equation in the fluid domain; however, it has been demonstrated in (Giorgi and Ringwood 2018d) that Wheeler's stretching is the best option for computationally efficient calculation of nonlinear Froude-Krylov forces, since it significantly improves the pressure field description accuracy without increasing its complexity. Froude-Krylov generalized forces (F F K ), divided into linear forces (f F K ) and torques (τ F K ), integrate the undisturbed pressure field ( p u ), shown in (5), as follows: where f g is the gravity force, τ g its contribution to the torque, n is the unity vector normal to the surface, r is the generic position vector, r R = (x R , y R , z R ) is the reference point around which the torque is computed, and likewise r g is the position vector of the centre of gravity. Finally, for future reference f p and τ p indicate the pressure integrals for the force and torque, respectively. While generic NLFK solvers for arbitrary complex floaters must rely on a meshed representation of S w (t), which becomes the computational bottleneck, a faster analytical representation, already available for axisymmetric floaters (Giorgi 2019), can be readily obtained also for prismatic floaters, as discussed hereafter.

Body-fixed frame and mapping
The free surface elevation and the pressure field are defined, according to the linear potential flow, in an inertial worldfixed frame of reference, with the origin at the still water level. Therefore, in order to evaluate the action of the pressure field onto the instantaneous wetted surface, it is necessary to displace and rotate the whole floater in the world frame. However, an alternative approach is to consider a body-fixed frame and rotate the wave and pressure fields around it. The great advantage of the body-fixed approach is to preserve all geometric parameters invariant during the simulation; such a simplification is paid with a more complex pressure and free surface formulation, which has to be mapped from the worldframe to the body-frame. In (Giorgi and Ringwood 2018b) it is shown that the body-fixed formulation exhibits the best compromise in terms of both simplicity and computational time.
The body-fixed frame x,ŷ,ẑ is defined such that it is coinciding with the world-frame when the floater is at rest, i.e. right-handed, with the origin at the still water level,x parallel and concordant with the positive direction of wave propagation, andẑ upwards. Hereafter, every quantity decomposed in the body-fixed directions is shown with aˆabove.
The mapping function from the world-to the body-frame is based on the displacement of the floater, here assumed planar (just one rotation), due to the unidirectional wave. Let us define the pure rotation matrix (R θ ) as The rotation matrix R θ , if directly applied, generates a rotation around the origin. Conversely, it should be applied in order to induce a rotation around the reference point r R . Therefore, notionally starting at rest, i.e. when the two frames are coincident, the following steps are taken: a translation is applied to have the origin at the reference point; then the rotation is applied; then a translation opposite to the first one takes the origin back; finally the linear displacement is applied.
All such operations can be executed via a single composite homogeneous transformation by expanding the dimension of coordinate vector and using the following mapping matrix: It follows that the transformation between frames is given by the following:

Geometry formulation
In the body-fixed frame, a generic pointr belonging to the external surface of a prismatic body invariant in theŷ direction can be mapped by a change of coordinates R 3 → R 2 as as follows: r : where γ x and γ z are generic parametric curves, and α is the sweep parameter. Typically, but not necessarily, α goes from 0 to 1 as the directional curve moves from one to the other end. The parametric formulation enables the use of arbitrary complex shapes. A degenerate case is when γ x (α) = α, equivalent toẑ = γ z (α), which has the intrinsic limiting constraint of one-to-one relationship betweenẑ andx. According to the formulation in (10), hereafter avoiding to report the explicit dependence on α for simplicity, the consequent directional vectors are as follows: that can be used to compute the Jabobian (J ), the normal unity vector (n), and the ( r −r R ×n) product as follows: Note that the denominator of the multiplying coefficient of bothn and (r ×n) is exactly the Jacobian, which will simplify the analytical formulation of the NLFK integrals hereafter.

Integral formulation
Pressure integrals are computed in the body-fixed frame in order to take advantage of the invariant geometry formulation, and consequently mapped onto the world-frame: f p = R θf p , while τ p =τ p , due to the planar motion.
Due to the change of coordinates in (10), the pressure integrals in (6) become as follows: where W is the width of the floater (along the y-direction). Note that (13b) simplifies because the incoming wave is unidirectional ( p u invariant withŷ) and the geometry, moving in the plane, is prismatic alongŷ; therefore, the components of the torque in roll and yaw directions integrate to zero.
The last required step is to compute the instantaneous intersection between the floater and the free surface elevation (η(x, t)), namely the bounds α 1 and α 2 of the sweep parameter range. Since the integral formulation is in the body-frame, while the free surface is in the world-frame and depends on the world-frame coordinates, appropriate mapping is required. Let us define the instantaneous free surface in the x − z world-frame frame as follows: The mapped free surface elevation, in homogeneous coordinates, is defined asΓ = R −1 Γ . However, in order to compute the intersection, the explicit dependence on x should be substituted with the respective α, i.e The intersection is finally found numerically as the zero of I (α) =η(α) − γ z (α).

Case study
Prismatic wave energy converters typically exploit the pitching motion response to extract power, either directly or indirectly. Direct extraction requires a power take-off (PTO) mechanism directly acting on the pitch motion, either via a fixed/semi-fixed reference (sea-floor or external structure such as a breakwater or a floating multi-purpose platform), or via a multi-body floating system (such as a hinge-raft device). Indirect extraction is based upon inertial coupling to transfer part of the incoming energy from hydrodynamically-excited DoFs (pitch) to internal (dry) DoFs. A few examples of indirect extraction from pitching motion are the ISWEC , the PeWEC (Sirigu et al. 2020c), WELLO (Wello 2020), and SEAREV (Clément et al. 2005).
For sake of generality, direct extraction is first and mainly analysed in this paper, in order to investigate the direct effect of NLFK on power production performance. However, also indirect extraction is considered, since the inertial coupling is likely to attenuate the impact of hydrodynamic nonlinearities, due to the fact that energy is transferred to a DoF not in contact with the waves. Moreover, inertial coupling acts as a stabilizer of the pitching motion, potentially further decreasing the significance of hydrodynamic nonlinearity. Therefore, the hull of an inertially activated WEC is considered, first in a direct-extraction configuration, in Sect. 3.1, then in an indirect-extraction configuration, in Sect. 3.2.

Heave-pitch device: direct pitch extraction
The floater is inspired at the Inertial Sea Wave Energy Converter (ISWEC), whose geometry is shown in Fig. 1. The cross-sectional area of the prismatic floater is defined by three tangent circumferences. However, in order to define a single curve (instead of a piece-wise composition of three curves), a Bezier curve has been fitted to the profile. Bezier curves are generic parametric curves and a flexible tool to represent complex geometries (Garcia-Teruel et al. 2020).
Note that, in principle, the floater shown in Fig. 1 could sink part of the deck under the water. Such instances could be represented geometrically, simply considering the equation of the line joining bow and stern upper limits. Such conditions are taken into account in the static analysis of the submerged portion of the floater for different displacements (Sect. 4.1). However, the formulation of the dynamic pressure from the undisturbed pressure field would be hardly applicable to waves climbing above the deck, so that its integration on the deck surface would be inaccurate. Moreover, normal working operations should avoid such green-water effects, which fall under survivability conditions. Note that green-water effects can be meaningfully described only by fully-nonlinear approaches such as CFD (Computational Fluid Dynamics) (Rosetti et al. 2019) or SPH (Smoothed-Particle Hydrodynamics) (Le Touzé et al. 2010). Therefore, both for simplicity and because the NLFK model should be used in the powerproduction region and not under extreme conditions, the total pressure is integrated over the bottom portion of the hull only when response to incoming waves is considered, even if the deck becomes partially submerged.
Although unidirectional waves provide a 3-DoF planar excitation (surge, heave, and pitch), a simpler 2-DoF (heave and pitch) is used. In fact, in order to clearly highlight nonlinear effects, a fully-uncoupled model under linear conditions is preferred, so that any coupling found in the nonlinear model is due to NLFK forces. Therefore, surge is neglected, since it is linearly coupled with pitch. Moreover, NLFK in surge would generate mean drift forces that would make the distinction between different nonlinear effect more cumbersome. Similarly, for a clear inference of causality, only NLFK effects are included, hence the effect of mooring stiffness (Paduano et al. 2020) and viscous drag has been linearized (Fontana et al. 2020). However, note that, depending on the device geometry, working principle and mooring configuration, surge and drift forces may have a significant impact on the device behaviour and performance. Nevertheless, for the particular case of the ISWEC-like device, it is found that the surge DoF has a negligible impact on power production (Sirigu et al. 2020b).
The linear equation of motion about the center of gravity, including linear uncoupled viscous drag (B v ) and linear uncoupled mooring stiffness (K m ), written in the frequency domain for compactness, follows: where ξ 2 is the 2×1 state vector, composed of heave (z) and pitch (θ ), M the diagonal inertia matrix, A(ω) and B(ω) the diagonal frequency-dependent added mass and radiation damping, K h the diagonal linear hydrostatic stiffness, F d and F F K d are the diffraction and linear dynamic FK forces, and B PT O and K PT O the diagonal PTO damping and stiffness. Note that the linear hydrodynamic curves are computed just once for a given mean wetted surface of a floater; the computation is performed in frequency domain using a linear Boundary Element Method (BEM) software, such as Nemoh (Babarit and Delhommeau 2015) or WAMIT (WAMIT 2019). It is assumed that only the pitching motion is used to extract energy, so that the PTO matrices are defined as follows: where 1 2 is the 2×2 identity matrix. The NLFK version of (16) alternatively computes F F K d − K h ξ 2 in a nonlinear way, namely through equations (13).
Finally, note that both linear and nonlinear equations of motion are solved in the time domain through a constant time step Runge-Kutta integration scheme, substituting the radiation frequency domain components by their time-domain state-space approximation, identified via the FOAMM toolbox Peña-Sanchez et al. 2019).

Heave-pitch-gyro device: indirect pitch extraction
The actual ISWEC-like device, shown in Fig. 2, exploits the inertial coupling between the floater and a spinning flywheel sealed inside the hull (Sirigu et al. 2020a), with the major advantage of no moving parts in contact with salty water. Fig. 1 Cross-sectional area of the prismatic floater of the ISWEC-like device, symmetric about the vertical z-axis, defined by a small circumference for the bow (and symmetric stern) sections (radius R (1) ) and center C (1) ), tangent with the keel circumference (radius R (2) ) and center C (2) ), determining the deck length (L), draft (D), and height (H ). Each circumference is defined between the boundaries [x 1 , x 2 ] Fig. 2 Representation of a possible structural layout of the ISWEClike device, with two counter-rotating spinning flywheels enclosed in a sealed hull. As the wave propagates in the x direction, the floater moves in heave (z) and pitch (θ). Through the gyroscopic coupling, the pitching motion induces the inner gyroscopic unit to oscillate around the precession axis (ε) Therefore, energy is extracted from the oscillation of the gyroscopic unit (ε). The state vector is expanded to 3 DoFs: ξ 3 = (z, θ, ε) T , with the (linearized) gyroscopic effect providing the coupling between the θ and ε: where n g is the number of gyroscopic units (here 2, with counter-spinning flywheels), J f the axial inertia of the flywheel, andφ f the rotational speed of the flywheel. For a comprehensive description of the gyroscopic effect and the resulting mechanical coupling between the flywheel and the hull, please refer to (Battezzato et al. 2015). Three control variables are available as follows:φ f , B PT O , and K PT O , now acting only on ε: where 1 3 is the 3×3 identity matrix.
Finally, (16) is further expanded to include the inertia and a linear stiffness of the gyroscopic unit around the ε axis ).

Results
This section progressively investigates the NLFK model, first considering the consistency of supporting hypothesis and the correctness of the numerical implementation, then isolating specific nonlinear effects under simplified scenarios, and finally evaluating the impact of nonlinearities in response to incoming waves, considering first a simple uncontrolled and unconstrained case, and then applying constrained reactive control. In particular, Sect. 4.1 focuses on the hydrodynamics of the floater, Sect. 4.2 studies the response to incoming waves of a WEC extracting energy directly from the pitching motion (as discussed in Sect. 3.1), while Sect. 4.3 considers the same floater with an internal PTO excited through inertial coupling (as discussed in Sect. 3.2).

Results: Floater hydrodynamics
The correctness of the mathematical formulation and the numerical implementation can be probed by considering the dynamic FK forces under linear conditions, which can be evaluated in two independent ways, and compared. On the one hand, linear boundary element software, such as NEMOH (Babarit and Delhommeau 2015), return, by definition, linear hydrodynamic data, which refer to the mean wetted surface and flat free surface elevation. On the other hand, the NLFK approach can be used under significantly linear conditions (zero displacement and very small η, wave height of 1cm) to obtain equivalent results. Figure 3 shows amplitude and phase of such frequency-dependent dynamic FK force, with perfect agreement between linear (LFK) and nonlinear (NLFK) computation. Furthermore, the diffraction force is also shown, which is the complementary part of the dynamic FK force to build up the total excitation force, according to linear potential flow theory. In the regions of operational wave periods (T w > 5.5s) of this WEC, and especially in the working DoF (pitch), the diffraction component is considerably smaller than the Froude-Krylov part, which is a promising hint for the effectiveness of the NLFK approach.
The impact of static Froude-Krylov force, or restoring force, is considered in a pitch free decay in Fig. 4, which clearly show a coupling effect with heave. As extensively discussed in Sect. 3.1, the linear model is fully-uncoupled but for NLFK force, which is hence the cause of the coupling. Although small in amplitude, the heave response shows a pronounced nonlinear behaviour, with a main low frequency response (equal to the pitch natural frequency, about 5.5 s), and a higher frequency component creating doubling peaks, especially evident for larger initial displacements. This can be explained by analysing the variation of submerged volume with respect to the equilibrium (ΔV , positive if the volume increases) for different heave and pitch displacements, shown in Fig. 5.
Vertical equilibrium, hence zero heave displacement, is ensured when ΔV is zero, along the dashed black line in Fig. 5; if also θ is zero, the static equilibrium is reached (black dot in Fig. 5). With reference to the free decay in Fig. 4, Fig. 5 shows that, for both positive and negative pitching angles, the submerged volume increases, so that a positive heave is required to restore the equilibrium. In particular, the initial pitch of 12 • generates a ΔV of 61 m 3 , corresponding to a required displacement of 0.16 m; this is consistent with the maximum displacement in Fig. 4, which is slightly smaller since pitch is concurrently decaying. However, as pitch becomes negative and increases in absolute value, there is a consequent increase in submerged volume, hence heave starts to increase again, generating a second peak before pitch completes a full period motion.
The map in Fig. 5, generated for given displacements in calm water, also highlights a change in convexity of the isolines, marked by the red dotted lines, which correspond to the instances when the deck of the floater peaks under water: as the floater sinks, the maximum allowable pitching angle before the deck sinks under water decreases. Such conditions are to be avoided during normal operation.
A final overview of nonlinearities in the static FK force can be achieved by comparing several free decay tests from different initial conditions, normalized by maximum displacement, as shown in Fig. 6. By definition, a linear model would produce normalized time series perfectly overlapping; conversely, any difference is ascribable to nonlinear effects, the more important the larger the displacement.
While pitch is remarkably linear in the range of initial displacements of ±12 • , the heave response is both nonlinear and asymmetric, with normalized amplitude of motion decreasing as the initial condition goes from negative to positive values; when the initial condition is very small, the normalized curves overlap (as to be expected). The asymmetry is also evident in the half-natural period: although the first half is shorter for negative z 0 , the full natural period is longer for the positive z 0 . This means that the half-cycle from negative to positive is faster than the half-cycle from positive to negative displacements, which is consistent with the asymmetric rate of change of the submerged volume. In particular, the rate of change of the submerged volume is monotone for monotone displacement, due to the convex shape of the hull, and not symmetric with respect to the still water level.
Finally, a further tool to investigate the actual impact of nonlinearities is the prescribed-displacement test, as suggested in , i.e. to impose a displacement (computed with the linear model) to the floater and compute forces with different methods, namely linear and nonlinear. In this way it is easier to compare the effects of different inherent structure of force calculations, since the same displacement   The red dots show the boundaries of the three circular sections, as shown in Fig. 1 is considered; conversely, in a dynamic response to incoming waves, different forces generate different displacements, making effects mutually indiscernible. Therefore, considering the free (no PTO) response, computed with the linear model, to an incoming regular wave, with wave period (T w ) of 6.5 s and wave height (H w ) of 3 m, Fig. 7 shows the total (static and dynamic) FK force, according to the LFK and NLFK models. A compelling graphical representation of the difference between linear and nonlinear approaches is provided in Fig. 8, showing a snapshot of the prescribed motion test of Fig. 7. On the one hand, the linear approach considers the mean wetted surface (intersection of floater at rest with the still water level), regardless of the actual position of the hull or the incoming wave; on the other hand, the NLFK approach analytically computes the instantaneous wetted sur-face and submerged volume of the hull with respect to the instantaneous relative position between floater and free surface elevation, remarkably different from mean values, with consequent different centre of buoyancy and hydrodynamic forces.
The difference between η and z gives a notional estimate of the variations of the wetted surface. In Fig. 7, NLFK in heave is above LFK for all t, i.e. with a magnitude smaller for negative values and larger for positive values. While LFK is a perfect sum of two sinusoids (one is −K h¸2 , the other is F F K d ), NLFK convexity differs from the simple sinusoid function. In particular, in the region where z − η is positive, forces are smaller and with a lower curvature. This is explained by the fact that when z − η > 0, the hull is less submerged than in equilibrium; moreover, due to the curvature of the hull, the rate of change of the wetted surface is higher closer to the flat deck, which is reached when z − η is positive.
Since the importance of nonlinear effects depends on the relative motion between the hull and the free surface, a more systematic analysis of prescribed displacement test results is shown in Fig. 9 with a complete map of the difference between LFK and NLFK forces (ΔF F K ), using an error metric defined as the root mean square (rms) of ΔF F K , normalized by the maximum (equal to the amplitude) of the LFK force. The absolute value is used to avoid partial errors to cancel out, since NLFK is always on one side of LFK, as shown in Fig. 7. Figure 9 presents errors up to 35% for heave, up to 18% for pitch, localized in the region of T w between 5 s and 8 s, which is where the pitch resonance lays. Clearly, such errors increase with increasing H w . However, it is worth remarking that, although meaningful, such differences are obtained for prescribed displacement. Section 4.2 considers real coupled simulations, where the displacement depends on the forces, and vice versa.
Note that, hereafter, from the complete grid of regular waves analysed, and from all colour maps in forthcoming figures, waves with excessive steepness (over 6%) are excluded. Fig. 9 Map of the normalized root mean square (rms) of the absolute value of the difference between linear and nonlinear FK forces in heave (left) and pitch (right) for the prescribed-displacement test

Results: Response of the heave-pitch device
A second-order Runge-Kutta time progressing scheme is adopted to solve both linear and nonlinear models, using 75 time steps per wave period, simulating 50 wave periods and applying a smoothing sigma ramp to the first two wave periods. In order to verify the correctness of the timedomain implementation and the accuracy of the state space approximation of the radiation forces, the response amplitude operator (R AO) in heave and pitch is computed under linear conditions (very small wave height), both with linear and nonlinear time domain models, and compared with independent frequency domain calculation. Figure 10 shows good agreement in both amplitude and phase.
Since the major perk of the NLFK approach herein proposed is the low computational burden, a representative wave is used for evaluating the computational time, i.e. a wave at the pitch resonant frequency: T w of 6 s and H w of 2 m. Calculations are performed on a single core of a standard laptop, with processor Intel(R) Core (TM) i7-7500U CPU @ 2.70GHz and 8GB RAM. In the NLFK model, a slightly non-zero initial condition is provided ([0.001 m, 0.001 rad]) in order to favour convergence of the numerical integration (Giorgi et al. 2020b). Finally, 300 s of simulation with 3750 time steps requires just 1.4 s, which is 0.45% real-time computation and about 0.3 ms per time step. Therefore, the simulation time is about two orders of magnitude faster than real-time, hence 3 to 4 orders of magnitude faster than stateof-the-art mesh-based approaches.

Unconstrained uncontrolled conditions
The response of a free floater without the action of an external controller modifying its dynamic, also without energy extraction, is first analysed. Figure 11 shows the mean (left) and amplitude (right) of heave (top) and pitch (bottom) response, obtained according to the NLFK model. A major evidence of nonlinear response is the non-zero mean obtained for zeromean excitation, present mainly in the range of T w between 5 s and 7 s, which is the area where the floater is pitching the most, suggesting that pitch is the main cause of nonlinearity. The non-zero mean in heave is due to the nonlinear restoring force and is consistent with the free-decay and volume variations discussed in Sect. 4.1; the non-zero mean in pitch is likely due to a combination of second-order effects of the dynamic FK force and nonlinear coupling with heave. In fact, in the nonlinear range of T w , a clear coupling with heave is evident, which presents a local drop from the otherwise T windependent response.
Further insight on the importance of nonlinearities is achievable by comparison of the predicted amplitude with the LFK and NLFK models, shown in Fig. 12  where nonlinear coupling appears. The pitch nonlinear response has a higher peak, which moves to larger T w as H w (and the nonlinearity) increases. However, the free response seems overall more narrow-banded.
While the analysis of the free response is valuable, it is just preparatory for the investigation of the controlled conditions. In fact, the energy-maximising control, which should always be implemented in real applications yearning for economic viability, substantially modifies the dynamic response of the system and exaggerates the motion, hence potentially escalating the importance of nonlinear effects.

Controlled conditions
This section purports to highlight how dangerous can be to rely on a linear model when an energy-maximising control strategy is implemented, returning absurd unphysical results when the states are unconstrained, and suboptimal control parameters when the states are constrained. Since regular waves are considered, the simplest constant-coefficients reactive control is applied, namely with a PTO torque defined by a stiffness term (K PT O ) and damping term (B PT O ). In the unconstrained case, such coefficients are optimally defined by the linear complex-conjugate control (CCC) (Falnes 2002). The left side of Fig. 13 shows the resulting pitch amplitude (top) and mean power extracted (P, bottom) according to the unconstrained CCC, clearly unrealistic with peak pitching angles up to 600 • . Note that the CCC changes the dynamics of the device in order to extract as much energy as possible, so that the gradient of the power curve naturally follows the direction of increase of energy content in the incoming wave (proportional to T w H 2 w ). A common measure to limit unrealistic response is to include a constraint on one or more states of the system, the peak pitching angle in this case, and perform a numerical constrained optimization. The right side of Fig. 13 shows the resulting pitching angle, which saturates at the prescribed constraint (θ max of 60 • ), and the consequent mean extracted power, which is inherently much lower than in the linear unconstrained CCC. In fact, note that the peak of the power map shifted to lower periods, at the lower boundary of the θ max plateau, where the constraint is naturally respected. Constraint violation is prevented by modifying the control parameters when the pitching angle would exceed the given boundary. Figure 14 shows that K PT O remains substantially the same in unconstrained and constrained cases. In fact, K PT O modifies the natural period in order to bring the system into resonance (it is consistently zero at the natural period, negative at larger periods, positive at smaller periods); since the model is linear, the natural period is insensitive to H w , so that the K PT O map is H w -invariant. Similarly, also B PT O is insensitive to H w in the unconstrained CCC. Conversely, the constrained optimization increases B PT O to damp out excessive pitching angles. However, since the unconstrained model is so unrealistic, the required B PT O to verify the constraint must increase significantly, hence being overestimated. Therefore, using the control parameters in Fig. 14 in a higher-fidelity model or eventually in a real physical device deployed in sea is likely to be greatly suboptimal. It is evident that a more realistic numerical model is required, returning physically plausible results already when an unconstrained control is implemented, so that the results of the constrained optimization can be accepted with higher confidence. Figure 15 shows the resulting pitching angles and mean power extracted according to the NLFK model. On the one hand, note that already in the unconstrained case the pitching angles remain way below the prescribed constraint, with a peak pitching angle of about 30 • . On the other hand, since the linear CCC is optimal for a linear model, it is suboptimal for a nonlinear model; in fact, P increases in the constrained optimization. Furthermore, since the peak pitching angle is much lower than θ max , the constrained optimum is also the global (unconstrained) optimum.
The pair of control parameters for the NLFK model are mapped in Fig. 15. On the one hand, the K PT O in the constrained optimization show some dependence on H w , certifying some change in the pitch resonance frequency due to nonlinear effects under controlled conditions. On the other hand, the optimal B PT O shows great differences from the CCC B PT O , especially at longer T w .
It is important to highlight that the optimal control parameters, especially B PT O , are substantially different according to the LFK model and NLFK model, as shown in Fig. 14 and Fig. 16, respectively. Accordingly, the optimal predicted pitching response and power production assessment are significantly different if a linear or nonlinear model is used, as shown in Fig. 17. In particular, note that the linear model

Results: Response of the heave-pitch-gyro device
In this section, the indirect energy extraction is considered, presented in Sect. 3.2, operated by means of a spinning flywheel sealed inside an enclosed hull, which transfers energy from the pitching DoF to the precession DoF (ε) by means of the gyroscopic effect. It is expected an overall decrease of the importance of the hydrodynamic nonlinearity, quantified in this section, since part of the energy is channelled away from wetted DoF into a dry DoF. A constrained numerical optimization is used to determine the control parameters, limiting the maximum pitching angle (as in Sect. 4.2) and the maximum precession angle (ε max = 85 • ), avoiding overturning of the gyroscopic unit. The objective function is the gross extracted power, since no power loss term is included in the model. Consequently, the flywheel speedφ f , despite being one potential control variable, is kept fixed, since the higher theφ f the higher the inertial coupling (beneficial for power extraction) but the higher the power losses (Sirigu et al. 2020a), which are not modelled in this case study. Therefore, as in Sect. 4.2, only K PT O and B PT O are tuned. No constraint on the maximum PTO torque is included, in order to provide the controller with due freedom to operate, and making the analysis of the impact of nonlinearities easier. Figure 18 shows the resulting heave and pitch response, according to the LFK model (optimized using the LFK model) and the NLFK model (optimized using the NLFK model). The heave response and nonlinear coupling is similar to the 2-DoF case, shown in Fig. 12. However, the linear pitch response does not saturate to the provided constraint, as in the direct pitch extraction, shown in Fig. 17, since energy is transferred from the pitching to the precession angle by means of the inertial coupling. The only remarkable difference between LFK and NLFK is the shift of the peak of the pitch response to larger T w as H w increases. Figure 19 shows the resulting precession angle and mean produced power, according to the LFK model (optimized using the LFK model) and the NLFK model (optimized using the NLFK model). In the LFK, ε always reaches the constraints: this is because the linear hydrodynamics reached unreasonable pitch response in the case of the device without gyroscope; therefore, there is abundant excitation in pitch, which is optimally transferred by the control parameters in order to maximize power while not violating the constraint on ε. Conversely, in the NLFK model there are regions where constraints are not reached, showing the influence of hydrodynamic nonlinearity. Furthermore, in such regions, control parameters are such to optimize the power in absolute terms, since no constraints are naturally violated. Finally, while the peak power in the LFK model is at the same period, it shifts to larger periods for larger waves in the NLFK model, which is consistent with the nonlinear behaviour found in the gyroscope-free case in Fig. 12. Moreover, in the NLFK model there is more power extracted between 6 s and 7 s, while less between 5 s and 6 s.
The optimized control parameters, used in Fig. 19, are shown in Fig. 20. As expected, in the linear model the K PT O is mainly insensitive to H w , while B PT O is in charge of limiting ε in order not to violate ε max . On the other hand, the NLFK model shows clear nonlinear behaviour in the variations of K PT O that accommodate variations of the resonance period. The B PT O has the role of respecting the constraints on ε and to maximise the power (where the constraints are not reached). Overall, since there are different control maps between LFK and NLFK models, although not as much as in Sect. 4.2, it can be asserted that nonlinearities are still important for modifying the system dynamics and, ultimately for power extraction. Therefore, it is essential to use nonlinear models as a base for energy-maximising control strategies.

Conclusions
All wave energy converters should implement an energymaximising control strategy to enhance power extraction, essential contribution towards reaching economic viability. However, especially under controlled conditions, the accuracy of the mathematical model has major consequences on the design, control, and operation of the device. In fact, while linear models are often sufficiently accurate under uncontrolled conditions, nonlinearities become important when the control exaggerates the motion of the floater, so that linear models become unreliable and potentially unrealistic.
Although the usual practice of applying constraints on the states of the system may seem enough to obtain plausible results, they still inherit the structural inaccuracies of the linear model, potentially leading to misleading and suboptimal results.
Therefore, it is especially important that the optimization of the control parameters is based upon reliable nonlinear models, so that the control strategy is effective when applied on the real system. However, such models must be computationally quick enough to run in iterative optimization algorithms. This paper proposes and implements a computationally performing nonlinear Froude-Krylov approach for prismatic floating platforms, which is 4 orders of magnitude faster than generic state-of-the-art mesh-based approaches. Thanks to about 0.5 ms per force evaluation on a standard laptop, computation in about 1% of real-time is achieved, which is compatible with numerical optimization. In the case study analysed, nonlinear behaviour is found in non-zero mean free-decay, asymmetric response, and nonlinear coupling between heave and pitch degrees of freedom. Furthermore, the dynamic excitation takes into account the actual instantaneous wetted surface, while being constant in a linear approach, with major impact on the unconstrained controlled response and on the optimized control parameters of the constrained controlled response.
In general, since nonlinearities are device-dependent, the proposed nonlinear computation approach can be an effective tool for investigating and understanding the nonlinear behaviour of a floater, since it is at the same time straightforward to implement and computationally efficient.