Derivation and Experimental Validation of a Parameterized Nonlinear Froude-Krylov Force Model for Heaving Point Absorber Wave Energy Converters

1. Introduction

The escalating global demand for clean and sustainable energy sources has fueled extensive research to harness the vast potential of ocean waves. With its consistent availability, wave energy is an attractive renewable resource. Site evaluations have confirmed the potential for ocean wave energy resources [1]. Among the various technologies developed, point absorber wave energy converters (PAWEC) have emerged as promising devices capable of efficiently capturing and converting wave energy into electrical power [2]. PAWECs, characterized by buoyant structures tethered to the seabed, offer adaptability to varying wave conditions. The wave-induced vertical motion of its buoyant platform is translated into useful work through a power take-off that can act as either a generator or actuator.

The successful implementation of PAWECs depends on carefully considering design aspects, including buoy size, shape, and mass distribution. Material selection must account for the harsh marine environment, and placement is critical, influenced by factors such as water depth and distance from the shore. Advancements in control systems enable adaptive strategies that maximize energy extraction under varying wave conditions. Improving PAWEC control mechanisms stands out as a crucial avenue with immense potential to improve the feasibility of wave energy utilization [3]. Recent advances highlight the importance of refined control methods, which demonstrate their ability to increase power absorption by up to 20% while simultaneously reducing structural loads [1].

Extensive research on energy-optimal control, particularly for buoy geometries, has been explored. Topics include: establishing upper bounds for energy extraction [4,5,6], impedance matching control [6,7,8,9,10], and exploring both closed-form [11,12] and numerical optimal control solutions [13] such as Model Predictive Control (MPC) [14,15,16,17]. A common theme is the importance of using accurate WEC math models for analysis and control design. Initially, research focused on buoy response regimes that could be approximated using linear differential equations, valid under the assumption of small motion amplitudes. However, a recent shift toward addressing the nonlinear response of buoys [18,19,20] is underway to improve energy extraction, particularly when the objective is to increase power production.

Computational fluid dynamics (CFD) models can accurately capture these nonlinearities, but their application in model-based real-time control laws is impractical [21]. A previous investigation [22] comparing various modeling approaches showed that an important nonlinearity comes from the Froude-Krylov (FK) force, particularly for axisymmetric buoys with varying cross-sectional areas. This insight underscores the importance of understanding and incorporating the nonlinear effects of FK forces in the modeling process to develop effective control strategies.

The FK forces arise from the pressure field around a submerged or partially submerged body as it moves through waves. Specifically, they are associated with the nonlinear effects induced by the varying shape and submersion of the structure. Together with the scattering forces, they constitute the entire non-viscous force exerted on a floating body subjected to regular waves [23]. Developing closed-form, computationally efficient FK force expressions is advantageous for implementing real-time, model-based control strategies. This is especially true for buoys whose shapes result in nonlinear FK forces.

Closed-form nonlinear FK force models have been derived in previous studies [24,25,26], aimed explicitly at using buoy shape effects in model-based control laws for large-motion, nonlinear operating regimes. Giorgi et al. [24] developed a method for generating closed-form FK force expressions using Airy’s wave theory to approximate the pressure on a buoy in regular waves. It has recently been used for model-based control solutions, including sliding model control [27], feedback linearization [28], and latching control [21,29,30].

These studies have taken advantage of various buoy shapes, ranging from conventional spheres and cylinders to more unconventional geometries such as the double cone-shaped buoy [18,28]. The FK force approach by Giorgi et al. has been extended in this paper [24] in several ways, including considering flat bottomed buoys, irregular waves, and a parametrization of a large class of buoy shapes using three parameters. Parameterization allows for rapid evaluation of different geometries, making shape optimization efficient. By generalizing the FK force form, this work provides a powerful tool for improving the design and control of wave energy converters. Furthermore, this paper compares the nonlinear FK model with the buoyancy force approach often used in model-based control [18]. Experiments are used to assess model performance using a spherical buoy subjected to regular waves at several frequencies and amplitudes.

The paper is organized as follows. The model form is introduced in Section 2, including the distinction between the buoyancy and Froude-Krylov terms that are derived in Section 3 and Section 4. Examples of these terms are compared in Section 5, followed by the experimental validation study in Section 6. Some concluding remarks are provided in Section 7.

2. System Description

Consider the axisymmetric buoy constrained to heave attached by a PTO to the sea floor in Figure 1. In the configuration shown, the buoy has a positive displacement $\zeta \left(t\right)$ of its draft line relative to the still-water line. The wave elevation, $\eta \left(t\right)$, is also positive in the figure and is assumed to have the general form of Eq. 1 that allows analysis of both regular and irregular waves.

$$\eta (x,t)=\sum _{i=1}^n{A}_{i}sin({\varphi }_i+{\chi }_{i}x-{\omega }_{i}t)=\sum _{i=1}^n{\eta }_i$$

(1)

where $\eta$’s ith component has an amplitude, wave number, angular frequency, and phase shift denoted $A_i$, ${\chi }_i$, ${\omega }_i$, and ${\varphi }_i$. The wavelengths, ${\lambda }_i=2\pi /{\chi }_i$, are assumed to be large compared to the maximum buoy radius, and therefore the free elevation is considered to be locally horizontal, as shown in Figure 1. A body-fixed coordinate system $(x,y,z)$ is referenced at the bottom of the buoy. The constant water depth, $h_w$, is the distance from the sea floor to the still water line and is assumed to be large, though it is shown as small in the figure.

Two model variants will be compared, both having the general form of Eq. 2

$$M\ddot{\zeta }=F_P+F_g+F_S+F_R+F_{PTO}$$

(2)

where M is the buoy’s mass. Assuming that the buoy does not affect the wave field, the pressure exerted on the buoy by the water is denoted $F_p$, and the buoy weight is $F_g$. The scattering force, $F_S$, when combined with the Froude-Krylov force, is called the diffraction force, $F_D$,31] and captures the effect of wave distortion around a fixed buoy [6]. The radiation force, $F_R$, is due to the energy transfer between the buoy and the water as the buoy moves. Finally, $F_{PTO}$ is the power take-off force that adds or removes buoy energy.

The two model variants shown in Eq. 3 are distinguished by how the $F_p+F_g$ term of Eq. 2 is treated. In one case, they are replaced by the Froude-Krylov force $F_{fk}=F_P+F_g$, and in the other by the buoyancy force $F_b=F_P+F_g$.

$$\begin{array}{cc}\hfill M\ddot{\zeta }& =F_{fk}+F_S+F_R+F_{PTO}\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill M\ddot{\zeta }& =F_b+F_S+F_R+F_{PTO}\hfill \end{array}$$

(3b)

In both cases, the $F_P+F_g$ term will be expressed in a parametric form where a wide range of buoy shapes can be represented using only three independent quantities. In addition, the advantages and disadvantages of the $F_{fk}$ and $F_b$ approaches are discussed. The most significant distinction between them is the assumed motion of the buoy. The $F_{fk}$ derivation assumes small motions, while the $F_b$ approach does not.

3. Derivation of ${\mathit{F}}_{\mathit{f}\mathit{k}}$

The derivation of $F_{fk}$ begins by developing the pressure field for irregular waves. The parameterization of the buoy shape is considered next, followed by a representation of $F_{fk}$ in compact form.

3.1. Pressure

The FK term of Eq. 3a, is the integration of the water pressure field, $\overrightarrow{P}(x,z,t)$, over the instantaneous wetted surface of the buoy, S, illustrated by the blue shaded region of Figure 2

$$F_{fk}=F_g-\int \int \overrightarrow{P}(x,z,t)·\widehat{n}dS$$

(4)

where $\widehat{n}$ is the normal unit vector pointing outwards from the buoy’s surface, $dS$ is the infinitesimal surface area of the buoy, and the red arrows show the pressure field at a few points. In this snapshot, the buoy’s draft line is below the free elevation, indicating that the buoy is moving relative to the wave.

To derive a closed-form expression for $F_{fk}$, it is necessary to obtain an expression for $\overrightarrow{P}(x,z,t)$ that can be integrated. The approach used here is the same as that introduced by Giorgi et al. [24] but with three modifications.

• Consideration of irregular waves of Eq. 1.
• Inclusion of flat-bottomed buoys.
• A three-parameter buoy shape.

The pressure applied to the buoy is derived using linear wave theory in the appendix and given by Eq. 5 where $\rho$ is the water density and g is the gravitational acceleration. Eq. 5 assumes that the fluid is incompressible, inviscid, and irrotational with small particle velocities. Furthermore, it’s valid only for the portion of the buoy below the SWL unrelated to the wave elevation.

$$P(x,z,t)=\rho \sum _{i=1}^n{\eta }_i(x,t){\displaystyle \frac{{\omega }_i^2}{{\chi }_i}}{\displaystyle \frac{cosh\left({\chi }_i(z+h_w)\right)}{sinh{\chi }_i{h}_w}}-\rho gz$$

(5)

The dispersion relationship between the wave number, frequency, and depth, ${\omega }^2=g\chi tanh\left(\chi h_w\right)$, is also derived in the Appendix and valid under the above mentioned assumptions. If we further assume deep water, $\chi h_w\gg 1,$ it can be approximated as ${\omega }^2=g\chi$. When applied to Eq. 5, it yields the compact form of Eq. 6.

$$P(x,z,t)=\rho g\sum _{i=1}^n{\eta }_i(x,t)e^{{\chi }_{i}z}-\rho gz$$

(6)

If the wave is regular with a single amplitude, frequency, and phase shift $\varphi ={\displaystyle \frac{\pi }{2}}$, Eq. 6 reduces to the equation used by Giorgi et al. in [24].

$$P(x,z,t)=\rho g{e}^{\chi z}Acos(\omega t-\chi x)-\rho gz$$

(7)

Finally, it is assumed that the pressure is uniform in x near the buoy, letting us set $x=0$ in Eq. 6. This is important since it allows for closed-form integration to obtain the FK forces and gives the final form of the pressure field of Eq. 8.

$$P(z,t)=P_d(z,t)+P_s\left(z\right)=\rho g\sum _{i=1}^n{\eta }_i\left(t\right)e^{{\chi }_{i}z}-\rho gz$$

(8)

where $P_d(z,t)$ and $P_s\left(z\right)$ are the dynamic and static pressure components respectively.

3.2. Pressure Integration and Buoy Shape Parameterization

The approach for generating a closed-form expression for $F_{fk}$ is the same as [24], namely, define the buoy shape by a radius $r\left(z\right)$ that varies vertically and then integrate the pressure field of Eq. 8 along the radius using cylindrical coordinates. The three additions to the procedure mentioned in Section 3.1 are discussed below. In preparation for a later comparison of FK and buoyancy models, we will write $F_{fk}$ with static and dynamic terms of Eq. 9.

$$F_{fk}=F_{fk,s}+F_{fk,d}$$

(9)

using the "s" and "d" subscripts of Eq. 8.

Figure 3 illustrates the notation used in the derivation where $r\left(z\right)$ is the revolution profile of the symmetry axis as a function of z that defines the shape of the buoy, $\theta$ is the angle of revolution, $\zeta \left(t\right)$ is the vertical displacement of the buoy’s draft line relative to the still water line (SWL) at $z=0$, and h defines the draft line relative to a fixed reference of the buoy, in this case, at the bottom of the buoy.

Unlike previous models [24], Eq. 10 is also supplemented by the pressure acting on the base of the buoy, calculated as the product of the base pressure and its area. The cylindrical coordinate rotation limits $[0,2\pi ]$ and the vertical limits, $[z_1,z_2]$, are determined by the instantaneous wetted surface below the SWL, or $z_1=-(h-\zeta )$ and $z_2=0$. The buoy’s non-uniform cross-sectional plane area is $S\left(z\right)=\pi r^2\left(z\right)$. Based on these definitions, and utilizing Eq. 4 with considering only the vertical component of $P(z,t)$ due to heave-only-motion, the $F_{fk}$ forces can be calculated as:

$$F_{fk}=F_g+\left[P_d(z_1,t)+P_s\left(z_1\right)\right]S\left(z_1\right)+{\int }_0^{2\pi }{\int }_{z_1}^{z_2}\left[P_d(z,t)+P_s\left(z\right)\right]r\left(z\right)r^{\prime }\left(z\right)dz$$

(10)

Eq. 10 could be rearranged to obtain the components of the Froude-Krylov (FK) forces as

$$\begin{array}{cc}\hfill F_{fk}=& F_g+S\left(z_1\right)P_s\left(z_1\right)+{\int }_0^{2\pi }{\int }_{z_1}^{z_2}{P}_s\left(z\right)r\left(z\right)r^{\prime }\left(z\right)dz+\hfill \\ \hfill & S\left(z_1\right)P_d(z_1,t)+{\int }_0^{2\pi }{\int }_{z_1}^{z_2}{P}_d(z,t)r\left(z\right)r^{\prime }\left(z\right)dz\hfill \end{array}$$

(11)

Substituting Eq. 8 into Eq. 11 and splitting it into static and dynamic components yields Eq. 12.

$$\begin{array}{cc}\hfill F_{fk,s}& =F_g-\rho g\left(S\left(z_1\right)z_1+{\int }_0^{2\pi }{\int }_{z_1}^{z_2}zr\left(z\right)r^{\prime }\left(z\right)dzd\theta \right)\hfill \\ \hfill F_{fk,d}& =\rho g\left(S\left(z_1\right)\sum _{i=1}^n{\eta }_i\left(t\right)e^{{\chi }_i{z}_1}+{\int }_0^{2\pi }{\int }_{z_1}^{z_2}\sum _{i=1}^n{\eta }_i\left(t\right)e^{{\chi }_{i}z}r\left(z\right)r^{\prime }\left(z\right)dzd\theta \right)\hfill \end{array}$$

(12)

Closed-form FK force expressions can be created for a wide range of buoy shapes by introducing a special case of a quadric surface given by Eq. 13 where R is the midline buoy radius and $\alpha$ is the inverse square of the profile slope’s asymptote in the vertical direction. Equation 13 limits the buoy shapes to being axisymmetric with circular cross sections.

$$x^2+y^2-\alpha z^2=R^2$$

(13)

Evaluation of Eq. 12 requires the profile equation $r\left(z\right)$ obtained from Eq. 13 by shifting the z axis by the buoy displacement $\zeta$ and solving for x or y while setting the other variable to zero as shown in Eq. 14.

$$r\left(z\right)=\sqrt{R^2+\alpha {(z-\zeta )}^2}$$

(14)

The area of the circular cross-section is

$$S=\pi r^2=\pi \left(R^2+\alpha {(z-\zeta )}^2\right)$$

(15)

and is used to describe the pressure effect on buoy shapes with flat bottoms in Eq. 12. Table 1 summarizes the classes of shapes possible with the parameterization of Eq. 14. The value of $\alpha$ lets the shape range from an oblate spheroid to the hourglass shape of a hyperboloid. The conical shape is a limiting case with zero radius at the midline. This has little practical application for a buoy but provides insight into the effect of the "necking" of the hyperboloid on nonlinear terms in the final differential equation model. Although h is not part of the parameterization of $r\left(z\right)$, it is an important parameter to define the vertical height of the buoy.

Three $n\times 1$ arrays, $\eta$, $ϵ$ and $\kappa$, whose elements are given in Eq. 16 are introduced to make the closed form $F_{fk}$ expression more compact. Note that $\eta$ contains the components ${\eta }_i$, while $\eta$ is the sum of the components or the free elevation.

$${\eta }_i=A_{i}sin({\varphi }_i-{\omega }_{i}t),{ϵ}_i=e^{z_1{\chi }_i},{\kappa }_i={\int }_{z_1}^{z_2}{e}^{z{\chi }_i}r\left(z\right)r^{\prime }\left(z\right)dz$$

(16)

Substituting Eq. 14 into ${\kappa }_i$ and carrying out the integration gives the closed-form expression for ${\kappa }_i$ of Eq. 17.

$${\kappa }_i={\displaystyle \frac{\alpha }{{\chi }_i^2}}\left[{ϵ}_i(1+h{\chi }_i)-{\chi }_i\zeta -1\right]$$

(17)

Substituting $\eta$ and $ϵ$ from Eq. 16 and $\kappa$ of Eq. 17 into Eq. 12 provides compact expressions for the static and dynamic FK force, which are combined in Eq. 19.

$$\begin{array}{cc}\hfill F_{fk,s}& =-{\displaystyle \frac{1}{3}}\pi \rho g\zeta \left[3R^2+\alpha {\zeta }^2\right]\hfill \\ \hfill F_{fk,d}& =\pi \rho g{\eta }^T\left[(R^2+\alpha h^2)ϵ+2\kappa \right]\hfill \end{array}$$

(18)

$$F_{fk}={\displaystyle \frac{1}{3}}\pi \rho g\left[3{\eta }^T\left((R^2+\alpha h^2)ϵ+2\kappa \right)-\zeta (3R^2+\alpha {\zeta }^2)\right]$$

(19)

While this work focuses on closed-form expressions, it is important to note that $r\left(z\right)$ could be any set of points describing an axisymmetric buoy shape. The integrals of Eq. 12 could be computed numerically in real time as the buoy moves. The form of the FK force would be the same except for a modification of the ${\kappa }_i$ term.

4. Derivation of ${\mathit{F}}_{\mathit{b}}$

Before developing the $F_b$ model, note that $F_{fk,s}$ of Eq. 18 is indeed a buoyancy force, but it neglects submersion due to wave elevation. This effect is included in the $F_b$ model below.

The total buoyancy force acting on the buoy is

$$F_b=F_g+\rho g{V}_{sub}$$

(20)

where $V_{sub}$ is its submerged volume, calculated as the integral of its varying cross-sectional area $S\left(r\right(z\left)\right)$ over the wetted surface.

$$V_{sub}={\int }_{z_1}^{z_2}S\left(z\right)dz$$

(21)

To account for wave elevation, the area function $S\left(z\right)$ and the lower integration bound $z_1$ are modified from those of Section 3 to include the instantaneous wave elevation and are shown in Eq. 22

$$\begin{array}{cc}\hfill S\left(z\right)& =\pi \left(R^2+\alpha {(z+\eta -\zeta )}^2\right)\hfill \\ \hfill z_1& =-h+\zeta -\eta \hfill \end{array}$$

(22)

Substituting Eq. 22 into the volume integral of Eq. 21, the buoyancy $F_b$ is given in Eq. 23.

$$F_b={\displaystyle \frac{1}{3}}\pi \rho g(\eta -\zeta )\left[3R^2+\alpha {(\eta -\zeta )}^2\right]$$

(23)

5. Model Form Comparison

The difference between Eq. 23 and $F_{fk}$ of Eq. 19 is how wave elevation, $\eta$, is used. The $F_{fk}$ expression uses the pressure field containing ${\eta }_i$. However, the boundary conditions of linear wave theory, used to solve Laplace’s equation, are applied to the still water line instead of the free elevation, $\eta$. This means that the pressure field of Eq. 8 is only valid up to $z_2=0$, which is fine as long as both $\eta$ and $\zeta$ are small. In contrast, the $F_b$ model can be applied to scenarios where both waves and buoy movement are large, though it omits dynamic pressure effects.

Another way to compare $F_{fk}$ and $F_b$ is to write $F_b$ in two parts where the first term is identical to $F_{fk,s}$ of Eq. 18

$$F_b=-{\displaystyle \frac{1}{3}}\pi \rho g\zeta \left[3R^2+\alpha {\zeta }^2\right]+{\displaystyle \frac{1}{3}}\pi \rho g\eta \left[3R^2+\alpha ({\eta }^2-3\eta \zeta +3{\zeta }^2)\right]$$

(24)

and the second term is an alternate representation of $F_{fk,d}$.

To illustrate the differences in the models, Eq. 19 and Eq. 23 are applied to four sample buoy geometries in Table 2 for regular waves, where $\eta$ has one component. The $F_b$ expressions illustrate the effect of the buoy shape on the nature of the nonlinear contribution to the differential equation model of Eq. 3b. For example, the cylinder is linear in $\eta -\zeta$, whereas the other shapes have a cubic effect. Compared to the linear term, the relative effect of the nonlinearity depends on R and, to some extent, $\alpha$. The limiting case of the cone-cone shape has no linear term, which means that it has zero stiffness at its equilibrium position with a draft line at its apex. The hourglass shape approximates this behavior for small R.

6. Experimental Validation

The two models of Eq. 3 are evaluated below using regular wave experiments. The objective was to compare $F_b$ with $F_{fk}$ without buoy motion, $\zeta =0$, and then examine the effect of the scattering force $F_S$.

6.1. Experimental Setup

The tests were carried out at MTUWave, the wave tank laboratory at Michigan Technological University, shown in Figure 4. A spherical buoy with a 10 cm radius was mounted to a dynamometer with a Sensing Systems load cell to measure the vertical force exerted by incoming waves. The range of the load cell was $\pm 222N$ with a maximum error of $\pm 0.556N$ and sampled at 100 Hz. The buoy was positioned to have a midline draft and the load cell output was adjusted to zero Newtons.

Two Edinburgh Designs resistance wave gauges, sampled at 128 Hz, were placed on either side of the buoy, laterally aligned with the center of the buoy. The wave gauge measurements were averaged to estimate the wave elevation at the buoy. The wave gauges were calibrated at five vertical positions, so their uncertainty was approximately 0.5 mm.

A dSPACE MicroLabBox was used to log the load cell and provide a synchronization signal to start the wave maker paddles and wave gauge logging. This allowed for the alignment of the data collected across multiple acquisition platforms.

The test conditions are shown in Table 3 and were chosen to exercise the model while not exceeding the capabilities of the wave tank. The nine tests included three different frequencies and amplitudes with wave steepness, the ratio of wave amplitude to wavelength, provided for reference.

6.2. Data Comparison

Wave elevation measurements, $\eta$, were used to calculate $F_{fk}$ and $F_b$ using Eq. 25, taken from the sphere entry of Table 2, where $\zeta =0$, $\rho =1000\mathrm{kg}/{\mathrm{m}}^3$, $g=9.81\mathrm{m}/{\mathrm{s}}^2$, and $R=0.1\mathrm{m}$.

$$\begin{array}{cc}\hfill F_{fk}& ={\displaystyle \frac{2\pi \rho g}{{\chi }^2}}\eta \left[1-(1+R\chi )e^{-R\chi }\right]\hfill \\ \hfill F_b& ={\displaystyle \frac{1}{3}}\pi \rho g\eta (3R^2-{\eta }^2)\hfill \end{array}$$

(25)

The calculated and measured forces are compared in Figure 5, Figure 6 and Figure 7. Both models perform well at low frequencies, the uppermost plot in Figure 5 through Figure 7. As the frequency of the wave increases, both models overpredict the force; however, the $F_{fk}$ approach outperforms the $F_b$ model.

The comparison of $F_{fk}$ and $F_b$ with the measured buoy forces ignores the contribution of the scattering force, $F_S$, of Eq. 3. This force is generally considered negligible when the buoy diameter is significantly smaller than the wavelength [32]. Neglecting $F_S$ aligns well with the 0.2 Hz frequency data, where $\lambda /R\approx 150$ and the difference between the measured and calculated forces is less than 0.5 N for the 20 mm amplitude case. However, neglecting $F_S$ for the 0.6 Hz and 1.0 Hz cases where $\lambda /R$ is about 40 and 16, respectively, is a possible explanation for the increase in the deviation between the measured and computed forces in the second and third plots in Figure 5 through Figure 7. The maximum deviation for $F_b$ is 1.7 N, and for $F_{fk}$, it is 0.8 N, both occurring in the 1.0 Hz case of Figure 7.

The total diffraction force, $F_D$, is defined as the sum of the scattering and Froude-Krylov forces [33].

$$F_D=F_{fk}+F_S$$

(26)

While the Froude-Krylov component can be obtained as described in Section 3, the scattering force can be extracted from the numerical solution of the diffraction problem, [32]. The complex scattering force coefficients, $F_S\left(\omega \right)$, were obtained using the boundary element method solver, WAMIT, applied to the spherical buoy using a 972-panel mesh from $\omega =0.01-15$ rad/s in increments of 0.01 rad/s. The real part of $F_s\left(\omega \right)$ is shown in Figure 8 and the time domain force used in the model was created using Eq. 27

$$F_S\left(t\right)=Re\left[F_S\left(\omega \right)\right]\eta \left(t\right)$$

(27)

for the nine cases of frequencies and amplitudes considered. Eq. 27 is the special case of the irregular wave form, Eq. 28, for a regular wave with a single amplitude, frequency, and phase shift.

$$F_S\left(t\right)=\sum _{i=1}^{n}Re\left[F_{\mathrm{S}}\left({\omega }_i\right)\right]A_{i}sin({\varphi }_i-{\omega }_{i}t)=Re\left[F_S^T\left({\omega }_i\right)\right]\eta \left(t\right)$$

(28)

Figure 9 compares the effect on the maximum force error of using the $F_s$ expression of Eq. 27 in the model. Each $3\times 3$ matrix shows the maximum force error between the model and the experiment for three frequencies and three amplitudes. The top two matrices are generated from the data of Figure 5 through Figure 7, and the bottom two matrices show maximum errors between the experiments and the model with $F_s$. As expected, as the frequency of the wave increases, $F_s$ has a larger effect. It is also readily observed that the $F_{fk}$ model estimates force better than the $F_b$ model as frequency increases. The $F_{fk}$ force is about twice as accurate as the $F_b$ force at 0.6 Hz and four times as accurate at 1.0 Hz.

As an example, Figure 10 shows the time history comparison of the measured and modeled forces when including $F_s$, and is a direct comparison to Figure 7.

7. Conclusion and Future Work

Both the $F_{fk}$ and $F_b$ models contain geometry-induced nonlinearity. The $F_{fk}$ model is based on linear wave theory and is applicable to small wave elevations and motion. In contrast, the $F_b$ model, while not as accurate as the $F_{fk}$ model in small-wave conditions, does not assume small waves and motion as long as the waves are far from breaking. Including a scattering force term improves force prediction even for large $\lambda /R$ ratios, e.g., 40.

The ability to express $F_{fk}$ using three parameters that span a variety of buoy shapes is convenient for the design of control laws based on the model and possibly the optimization of the shapes.

Irregular wave model validation needs to be performed to further bound the model applicability. Experiments with different buoy shapes and nonlinear control law implementation are additional areas to explore.