# Response Behavior of Bi-stable Point Wave Energy Absorbers under Harmonic
Wave Excitations

Mohammad A. Khasawneh
Department of Mechanical Engineering, Tandon School of Engineering, New York
University New York, 11201, USA.
Division of Engineering, NYU Abu Dhabi, Abu Dhabi, UAE.
<EMAIL_ADDRESS>
& Mohammed F. Daqaq
Department of Mechanical Engineering, Tandon School of Engineering, New York
University New York, 11201, USA.
Division of Engineering, NYU Abu Dhabi, Abu Dhabi, UAE.
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

To expand the narrow response bandwidth of linear point wave energy absorbers
(PWAs), a few research studies have recently proposed incorporating a bi-
stable restoring force in the design of the absorber. Such studies have relied
on numerical simulations to demonstrate the improved bandwidth of the bi-
stable absorbers. In this work, we aim to understand how the shape of the bi-
stable restoring force influences the effective bandwidth of the absorber. To
this end, we use perturbation methods to obtain an approximate analytical
solution of the nonlinear differential equations governing the complex motion
of the absorber under harmonic wave excitations. The approximate solution is
validated against a numerical solution obtained via direct integration of the
equations of motion. Using a local stability analysis of the equations
governing the slow modulation of the amplitude and phase of the response, the
loci of the different bifurcation points are determined as function of the
wave frequency and amplitude. Those bifurcation points are then used to define
an effective bandwidth of the absorber. The influence of the shape of the
restoring force on the effective bandwidth is also characterized by generating
design maps that can be used to predict the kind of response behavior (small
amplitude periodic, large amplitude periodic, or aperiodic) for any given
combination of wave amplitude and frequency. Such maps are critical towards
designing efficient bi-stable PWAs for known wave conditions.

_Keywords_ Wave energy, Point wave energy absorber, Bi-stability, Nonlinearity

##  1 Introduction

Wave energy constitutes one of the most promising and dense renewable energy
sources that is yet to be fully exploited. Since the beginning of human
civilization, several devices have been devised to harness energy from ocean
waves both at small and large scales. Today, methods used to exploit wave
energy can be categorized based on their working principle into three
different categories; namely, oscillating water columns, overtopping devices,
and point wave energy absorbers (PWAs) [1]. Among such approaches, PWAs
received the most attention due to their simple design and working principle.
In its simplest form, a PWA is composed of a partially-submerged body (buoy)
connected through a mooring mechanism to a linear electromagnetic generator
attached to the seabed. When waves set the buoy into motion, it pulls a cable
connecting it to a linear generator, which creates relative motions between
the translating part of the generator (translator) and stationary magnets
(stator). As per Faraday’s law of induction, this motion induces a current in
the generator coils.

Due to their fundamental principle of operation, traditional PWAs which employ
a linear restoring force can work efficiently only near resonance; i.e., when
the buoy’s velocity is in phase with the wave excitation force. Unfortunately,
for typical energetic marine sites, this condition cannot be easily satisfied
for reasonably sized systems. Because of the high stiffness of the hydrostatic
restoring force emanating from buoyancy, the resonance frequency of the
absorber is typically higher than the dominant frequencies in the spectrum of
the incoming ocean waves [2]. Furthermore, because linear PWAs has a narrow
bell-shaped frequency response with a distinct peak occurring at the resonance
frequency, they are incapable of efficiently extracting power from the wide
frequency content of the ocean waves; thereby leaving most of the wave energy
unexploited.

To overcome such issues, different ideas and solutions have been proposed [3].
These include the use of active control strategies to bring the natural
frequency of the absorber closer to the dominant frequency in the ocean wave
spectrum, and the introduction of a bi-stable restoring force to broaden the
frequency response bandwidth of the absorber [4, 5]. The idea of utilizing a
bi-stable restoring force in PWAs emanated from the field of vibration energy
harvesting, where it was shown that vibratory energy harvesters whose
potential energy function has two potential wells separated by a potential
energy barrier have a broader frequency bandwidth, and are, therefore, less
sensitive to changes in the excitation parameters [6].

A schematic diagram of a bi-stable PWA is shown in Figure 1. The only
difference between the linear and bi-stable PWAs is the addition of the bi-
stable spring attachment in parallel with the power take-off unit (PTO). This
attachment is specifically designed to create a bi-stable restoring force
behavior and can be created by using a set of pre-stretched springs [4] or by
using magnetic interactions [5, 7, 8, 9]. The shape of the potential energy
function associated with the bi-stable PWA is shown in Figure 2. The system
has two stable equilibria (nodes) separated by a potential barrier (saddle).
For some combination of the wave frequency and amplitude, the response of the
buoy remains confined to a single potential well (intra-well motion), while
for others, the dynamic trajectories overcome the potential barrier causing
the buoy to undergo large-amplitude inter-well motions that span the two
stable equilibria. This type of large-amplitude motion can extend over a wide
spectrum of frequencies which, depending on the shape of the restoring force,
can even extend to very low frequencies. These characteristics are key to
improving the energy capture from the lower frequency content of the ocean
waves.

In terms of performance, a comparison between linear and bi-stable PWAs has
revealed a superior bandwidth for the bi-stable absorbers under harmonic waves
conditions [4]. In addition, when considering irregular random waves, results
demonstrated superior robustness of the bi-stable absorber with less
sensitivity to variations in the frequency content of the waves. A numerical
analysis performed in Ref. [10] demonstrated that the performance of the bi-
stable PWA is dependent on the shape of its potential energy function and the
ability of the dynamic trajectories to escape the potential wells for any
given combination of wave frequency and amplitude. Thus, in order to improve
the ability of the absorber to perform large-amplitude inter-well motions for
a wide range of wave conditions, an adaptive bi-stable absorber, which can
adjust the depth of its potential barrier to match the waves excitation
intensity was proposed first in Ref. [11], followed by other studies [12, 13].

Figure 1: Schematic diagram of a bi-stable PWA. Figure 2: Typical potential
energy function of a symmetric bi-stable system showing the two types of
possible motions (inter- and intra-well).

We noticed that, despite the relatively large body of research focused on
studying the behaviour of bi-stable PWAs, all of the previous studies relied
on purely numerical means without attempting to investigate the complex
underlying dynamics of the system via analytical or semi-analytical
techniques. Utilizing approximate analytical solutions of the governing
nonlinear equations can provide key additional insights into the long-time
behavior and bandwidth characteristics of the PWA that cannot be otherwise
inferred by relying on numerical simulations alone [6].

Aiming to bridge this gap, we derive in this paper an approximate analytical
solution of the nonlinear differential equations governing the complex motion
of the absorber under harmonic wave excitations. Using a local stability
analysis of the equations governing the slow modulation of the amplitude and
phase of the response, we determine the loci of the different bifurcation
points as function of the waves’ excitation frequency and amplitude. Those
bifurcation points are then used to define an effective bandwidth of the
absorber. We generate design maps that characterize the influence of the shape
of the potential energy function of the PWA on its effective bandwidth. Those
maps can be used to predict the type of response behavior of the absorber;
e.g. small amplitude periodic, large amplitude periodic, or aperiodic, for any
given combination of wave amplitude and frequency. We believe that such maps
are valuable towards designing efficient bi-stable PWAs for known wave
conditions.

The rest of the paper is organized as follow: in Section 2, the mathematical
model governing the motion of the bi-stable PWA is presented and discussed. In
Section 3, an asymptotic analytical solution of the governing equations is
derived using the method of multiple scales for both intra-well and inter-well
oscillations. In Section 4, the different bifurcations of the asymptotic
solution are identified and analyzed using a stability analysis of the
equations governing the slow modulation of the response. In Section 5, design
maps that characterize the influence of the shape of the bi-stable restoring
force of the PWA on its effective bandwidth are generated and discussed.
Finally, in Section 6, the main conclusions of this work are presented.

##  2 Mathematical formulation

###  2.1 Governing equations

Figure 3: A lumped-parameter model of the bi-stable PWA.

Assuming that the buoy undergoes motions in the heave direction only, the
equations governing the motion of the absorber can be obtained by applying
Newton’s second law on the buoy, and Kirchhoff’s current law on the harvesting
circuit to obtain the following governing equations for the equivalent lumped
system shown in Figure 3:

$\displaystyle\begin{split}&(m+m_{\infty})y^{\prime\prime}+\int_{0}^{t}h(t-\tau)y^{\prime}d\tau+cy^{\prime}+(k_{hys}-k_{1})y+k_{3}y^{3}=f_{wave}\cos(\omega
t),\end{split}$ (1a) $\displaystyle
V_{L}^{{}^{\prime}}+\frac{R_{L}}{L}V_{L}=\alpha y^{\prime}.$ (1b)

Here, $y$ represents the displacement of the buoy in the heave direction and
the overprime represents a derivative with respect to time, $t$. In Equation
(LABEL:eq:Subeq1), $m$ and $m_{\infty}$ represent, respectively, the mass of
the buoy and the added mass of the fluid. The integral term is used to account
for the radiation damping effect; $c$ is a linear viscous damping coefficient;
$k_{hys}$ is the stiffness resulting from the hydrostatic buoancy force, which
is equal to $\rho gSy$. Here, $\rho$ is the density of water, $g$ is the
gravitational acceleration constant, and $S$ is the buoy’s wet surface area.
The coefficients $k_{1}>0$ and $k_{3}>0$ represent, respectively, the linear
and cubic coefficients of the nonlinear restoring force added to introduce the
bi-stable behavior, and $f_{wave}$ and $\omega$ represent, respectively, the
wave amplitude and frequency. Note that to induce a bi-stable potential energy
function, $k_{1}$ must be larger than $k_{hys}$.

In Equation (1b), $V_{L}$ represents the voltage induced by the generator
across a purely resistive load $R_{L}$; $L$ represents the inductance of the
harvesting coil, and $\alpha$ is the electromechanical coupling coefficient.

The added mass and the radiation damping depend on the fluid velocity field
around the buoy, and hence, are a function of the wave frequency. Following
the asymptotic analysis provided in the work of Holme [14] for a spherical
buoy of radius $R$, the curves governing the dependence of the normalized
added mass denoted here as $\overline{m}_{a}=\frac{m_{a}}{M}$
$(M=\frac{2}{3}\pi R^{3}\rho)$, and the normalized radiation damping
coefficient $\overline{B}=\frac{B}{M\omega}$ on the normalized wave frequency
$\Omega=\omega/\sqrt{g/R}$ are shown in Figure 4.

Figure 4: Variation of the normalized added mass $\overline{m}_{a}$, and the
normalized radiation damping coefficient $\overline{B}$ with the square of the
normalized wave frequency.

The radiation kernel can be further related to the radiation damping
coefficient $B(\omega)$ via the following equation [15]:

$h(t)=\frac{2}{\pi}\int_{0}^{\infty}B(\omega)\cos(\omega t)d\omega,$ (2)

which can be discretized as following:

$h(t)=\frac{2}{\pi}\lim_{\delta\omega\to
0}\sum_{i=1}^{\infty}B(\omega_{i})\cos(\omega_{i}t)\delta\omega.$ (3)

Here, the $B(\omega_{i})$ are the values of the radiation damping coefficient
calculated at the discrete value of $\omega_{i}$ shown in Figure 4.

The amplitude of the force $f_{wave}$ acting on the buoy due to the incident
waves can be related to the radiation damping coefficients by employing
Haskinds’s relation [16, 17], which states that

$f_{wave}(\omega)=A_{wave}\sqrt{\frac{2\rho g^{3}}{\omega^{3}}B(\omega)},$ (4)

where $A_{wave}$ is the regular wave amplitude. The physical interpretation of
Equation (4) is fairly intuitive, as it relates the tendency of the buoy to
radiate waves in a certain direction to the excitation forces the buoy
experiences from waves propagating in the same direction. For more insight,
the reader can refer to the seminal works of Haskind and Newman[16, 17].

###  2.2 Approximation of the convolution integral

As aforestated, the goal of our work is to obtain approximate analytical
solutions of the equations governing the motion of the buoy in order to gain
deeper insights into the influence of the shape of the restoring force on the
performance of the absorber. To achieve this goal, we will use the method of
multiple scales [18]. In order to facilitate the implementation of the method,
we obtain in this section an approximation of the convolution integral
governing the radiation damping in Equation (LABEL:eq:Subeq1).

To this end, we first use the fact that the output, $z(t)$, of a linear
dynamical system can be expressed as a convolution integral between its input
$u(t)$, and an impulse response function $h(t)$, in the form

$\textbf{z}(t)=\int_{0}^{t}h(t-\tau)\textbf{u}(\tau)d\tau\approx\textbf{C}_{r}\textbf{x}(t).$
(5)

Here, $\textbf{z(t)}\in\mathbb{R}^{q}$ is the output vector,
$\textbf{x(t)}\in\mathbb{R}^{n}$ is the state vector and $\textbf{C}_{r}$ is a
$q\times n$ output matrix governed by the following linear state-space
equation

$\displaystyle\textbf{x}^{\prime}(t)$
$\displaystyle=\textbf{A}_{r}\textbf{x}(t)+\textbf{B}_{r}\textbf{u}(t).$ (6)

Using the Eigensystem Realization Algorithm (ERA) detailed in Ref. [19], the
convolution integral can be expressed in terms of the realized state-space
matrices $\textbf{A}_{r},\textbf{B}_{r}$ and $\textbf{C}_{r}$ as:

$h(t)=\left(M\sqrt{\frac{g}{R}}\right)\textbf{C}_{r}e^{\textbf{A}_{r}t}\textbf{B}_{r},$
(7)

where the realized state-space matrices are :

$\textbf{A}_{r}=0.8\begin{pmatrix}-1&1&1\\\ -1&0&0\\\ -1&0&-2\end{pmatrix},$
$\textbf{B}_{r}=\begin{pmatrix}-0.48&-0.02&-0.22\end{pmatrix}^{T},$

and

$\textbf{C}_{r}=\begin{pmatrix}-0.46&0&0.18\end{pmatrix}.$

Note that the numerical values appearing in the realized state-space matrices
are general for any spherical buoy of radius, R. For more details on the ERA
procedure, the interested reader can refer to Appendix A.

The matrix exponential $e^{\textbf{A}_{r}t}$ in Equation (7) can be further
expressed in the following form

$e^{\textbf{A}_{r}t}=\mathcal{L}^{-1}\left(s\textbf{I}-\textbf{A}_{r}\right)^{-1}.$
(8)

Here, $\mathcal{L}^{-1}$ is the inverse Laplace transform, and I is the
identity matrix. Substituting Equation (8) into Equation (7), we obtain the
following analytical expression for $h(t)$:

$h(t)=\left(M\sqrt{\frac{g}{R}}\right)e^{-\mu
t}\left(\lambda_{1}+\lambda_{2}\cos(\mu t)+\lambda_{3}\sin(\mu t)\right),$ (9)

where $\mu,\lambda_{1},\lambda_{2}$ and $\lambda_{3}$ are constants listed in
Table 1, and are valid for any spherical buoy of radius, R.

Table 1: Numerical values of the constants appearing in Equation (9). $Parameter$ | $Value$
---|---
$\mu$ | 0.8
$\lambda_{1}$ | -0.44
$\lambda_{2}$ | 0.62
$\lambda_{3}$ | 0.24

Figure 5 depicts a comparsion between the analytical expression of $h(t)$ as
obtained using Equation (9) and that obtained using the original expression of
Equation (3) for a hemispherical buoy of radius $R=5$ [m]. It can be clearly
seen that the analytical expression for $h(t)$ is in an excellent agreement
with the original impulse response function obtained via Equation (3).

Figure 5: Impulse response function $h(t)$ for a hemispherical buoy of radius
5 [m]. Solid line represents the original function obtained by applying
Equation (3). Circles represent the values of $h(t)$ obtained from the
analytical expression in Equation (9) which is based on the third-order
realized state-space estimated through employing the eigensystem realization
algorithm.

##  3 Approximate Analytical Solution

In this section, we employ the method of multiple scales [18] to obtain an
approximate analytical solution of Equation (1). As aforementioned, and shown
in Figure 2, there are two possible steady-state motions. Those that are
confined to the potential well and known as intra-well oscillations and those
that span the two potential wells and known as the inter-well motions. We
first obtain analytical approximations of the motion trajectories within one
potential well, then we seek approximate solutions that govern the inter-well
motions.

###  3.1 Local intra-well response

In this subsection, we obtain an approximate solution for Equation (1) when
the buoy undergoes motions within a single potential well; i.e, the dynamics
around one stable equilibrium node. Upon using the radius of the buoy, $R$, as
a length scale and the $\sqrt{g/R}$ as a time scale, we obtain the following
dimensionless equivalent equation of motion:

$\displaystyle\ddot{Y}$
$\displaystyle+\delta_{1}\int_{0}^{t^{*}}\overline{h}(t^{*}-\tau)\dot{Y}(\tau)d\tau+\delta_{2}\dot{Y}-\omega^{2}_{n}Y+\gamma
Y^{3}=g_{wave}\cos(\Omega t^{*}),$ (10a) $\displaystyle\dot{v}+\theta
v=\dot{Y},$ (10b)

where the overdot represents the derivative with respect to the nondimensional
time, $t^{*}$, and the other parameters and nondimensional groups are defined
as:

$\delta_{1}=\frac{m+m_{\infty}}{M}$,
$\delta_{2}=\frac{c}{(m+m_{\infty})}\sqrt{\frac{R}{g}}$,
$\omega_{n}=\sqrt{\frac{(k_{1}-\rho gS)R}{(m+m_{\infty})g}}$,
$\gamma=\frac{R^{3}k_{3}}{(m+m_{\infty})g}$,
$g_{wave}=\frac{AM}{R(m+m_{\infty})}\Omega\sqrt{\frac{3\overline{B}}{\pi}}$ ,
$M=\frac{2}{3}\pi R^{3}\rho$, $Y=\frac{y}{R}$, $v=\frac{V}{\alpha R}$,
$\overline{h}=\frac{h}{M}\sqrt{\frac{R}{g}}$,
$\theta=\frac{R_{L}}{L}\sqrt{\frac{R}{g}}$, $t^{*}=t\sqrt{g/R}$.

We expand the dynamics governed by Equation (10) about the stable node,
$Y_{s}=\sqrt{\omega^{2}_{n}/\gamma}$, by introducing the transformation
$Y(t^{*})=z(t^{*})-Y_{s}$ into Equation (10). This yields

$\displaystyle\begin{split}\ddot{z}&+\delta_{1}\int_{0}^{t^{*}}\overline{h}(t^{*}-\tau)\dot{z}(\tau)d\tau+\delta_{2}\dot{z}+\omega^{2}_{0}z+\eta
z^{2}+\gamma z^{3}=g_{wave}\cos(\Omega t^{*}),\end{split}$ (11a)
$\displaystyle\dot{v}+\theta v=\dot{z}.$ (11b)

Here $z(t)$ represents the dynamic trajectories within the potential well,
$\omega_{o}=\sqrt{2}\omega_{n}$ is the corresponding local frequency of
oscillations. It is worth noting that the expansion about the stable node
introduces a new quadratic term which captures the asymmetric nature of the
motion trajectories about the stable node, $Y_{s}$.

To implement the method of multiple scales on Equation (11), we introduce
multiple time scales, $T_{n}=\epsilon^{n}t^{*}$, $n=0,1,2$. It follows that

$\begin{split}\dot{(.)}&=D_{0}+\epsilon
D_{1}+\epsilon^{2}D_{2}+O(\epsilon^{3}),\\\
\ddot{(.)}&=D_{0}^{2}+\epsilon^{2}D_{1}+2\epsilon
D_{0}D_{1}+2\epsilon^{2}D_{0}D_{2}+O(\epsilon^{3}),\end{split}$ (12)

where $\epsilon$ is a scaling parameter, and $D_{n}$ is the temporal
derivative operator with respect to the time scales, $T_{n}$. Next, we seek an
expansion for the response $z(t^{*})$ and voltage $v(t^{*})$ as

$\displaystyle z(t^{*},\epsilon)$
$\displaystyle=z_{0}(T_{0},T_{1},T_{2})+\epsilon
z_{1}(T_{0},T_{1},T_{2})+\epsilon^{2}z_{2}(T_{0},T_{1},T_{2})+O(\epsilon^{3}),$
(13a) $\displaystyle v(t^{*},\epsilon)$
$\displaystyle=v_{0}(T_{0},T_{1},T_{2})+\epsilon
v_{1}(T_{0},T_{1},T_{2})+\epsilon^{2}v_{2}(T_{0},T_{1},T_{2})+O(\epsilon^{3}).$
(13b)

Based on the typically small values of the damping coefficients, quadratic and
cubic nonlinearities, electromechanical coupling, the time constant of the
harvesting circuit, and the excitation amplitude, we scale them to be at order
$\epsilon^{2}$; that is

$\delta_{1}=\epsilon^{2}\delta_{1}$, $\delta_{2}=\epsilon^{2}\delta_{2}$,
$\eta=\epsilon^{2}\eta$, $\gamma=\epsilon^{2}\gamma$ ,
$g_{wave}=\epsilon^{2}g_{wave}$.

Since large-amplitude intra-well motions occur near the primary resonance of
the system; i.e. when $\Omega$ is near $\omega_{o}$, we limit the intra-well
analysis to wave frequencies that are close to $\omega_{o}$ by introducing the
detuning parameter $\sigma$ such that

$\Omega=\omega_{o}+\epsilon^{2}\sigma.$ (14)

Upon substituting Equations (12 \- 14) into Equation (11), then collecting
terms of equal powers of $\epsilon$, we obtain the following perturbation
problems at the different scales:

$O(\epsilon^{0})$:

$\displaystyle D_{0}^{2}z_{0}+\omega_{o}^{2}z_{0}$ $\displaystyle=0,$ (15a)
$\displaystyle D_{0}v_{0}+\theta v_{0}$ $\displaystyle=D_{0}z_{0},$ (15b)

$O(\epsilon^{1})$:

$\displaystyle D_{0}^{2}z_{1}+\omega_{o}^{2}z_{1}$
$\displaystyle=-2D_{0}D_{1}z_{0}-\eta z_{0}^{2},$ (16a) $\displaystyle
D_{0}v_{1}+\theta z_{1}$ $\displaystyle=D_{0}z_{1}+D_{1}z_{0}-D_{0}v_{0},$
(16b)

$O(\epsilon^{2})$:

$\displaystyle\begin{split}D_{0}^{2}z_{2}+\omega^{2}_{o}z_{2}=&-2D_{0}D_{1}z_{1}-2D_{0}D_{2}z_{0}-D_{1}^{2}z_{0}-\delta_{1}\int_{0}^{T_{0}}\overline{h}(T_{0}-\tau)D_{0}z_{0}(\tau)d\tau\\\
&-\delta_{2}D_{0}z_{0}-2\eta z_{0}z_{1}-\gamma
z^{3}_{0}+g_{wave}\cos\left((\omega_{o}+\epsilon^{2}\sigma)T_{0}\right),\end{split}$
(17a) $\displaystyle\begin{split}D_{0}v_{2}+\theta
v_{2}=&D_{0}z_{2}+D_{2}z_{0}+D_{1}z_{1}-D_{2}v_{0}-D_{1}v_{1}.\end{split}$
(17b)

Upon solving Equations (15) and (16), we obtain the following expressions for
$z_{0}$, $v_{0}$, $z_{1}$ and $v_{1}$:

$\displaystyle z_{0}$ $\displaystyle=A(T_{1},T_{2})e^{i\omega_{o}T_{0}}+cc,$
(18a) $\displaystyle v_{0}$
$\displaystyle=\Gamma_{0}A(T_{1},T_{2})e^{i\omega_{o}T_{0}}+cc,$ (18b)

and

$\displaystyle z_{1}$
$\displaystyle=\frac{\eta}{\omega^{2}_{o}}\left(\frac{A^{2}(T_{1},T_{2})}{3}e^{2i\omega_{o}T_{0}}-2A(T_{1},T_{2})\overline{A}(T_{1},T_{2})\right)+cc,$
(19a) $\displaystyle v_{1}$
$\displaystyle=\eta\Gamma_{1}\frac{A^{2}(T_{1},T_{2})}{3\omega_{o}^{2}}e^{2i\omega_{0}T_{0}}+cc,$
(19b)

where

$\Gamma_{0}=\frac{\omega^{2}_{o}+i\theta\omega_{o}}{\omega^{2}_{o}+\theta^{2}},\qquad\Gamma_{1}=\frac{4\omega^{2}_{o}+2i\theta\omega_{o}}{4\omega^{2}_{o}+\theta^{2}},$

and $cc$ stands for the complex conjugate of the preceding terms.

Elimination of the secular terms from the second-order perturbation problem,
yields $D_{1}A(T_{1},T_{2})=0$, which implies that the complex valued function
$A$ is only dependent on the third time scale, $T_{2}$.

Elimination of the secular terms from the third-order problem associated with
Equation (17) requires evaluating the convolution integral
$\int_{0}^{T_{0}}\overline{h}(T_{0}-\tau)D_{0}z_{0}(\tau)d\tau$. To this end,
we use Equation (9) to write:

$\overline{h}(T_{0}-\tau)=\frac{e^{-\mu
T_{0}}}{M}\sqrt{\frac{R}{g}}\Big{(}\lambda_{1}e^{\mu\tau}+\lambda_{2}e^{\mu\tau}\cos(\mu\tau-\mu
T_{0})-\lambda_{3}e^{\mu\tau}\sin(\mu\tau-\mu T_{0})\Big{)}.$ (20)

Also, from Equation (18), we have

$D_{0}z_{0}(\tau)=i\omega_{o}Ae^{i\omega_{o}\tau}+cc.$ (21)

Thus, using Equations (20) and (21), we can express the convolution integral
as

$\begin{split}\int_{0}^{T_{0}}\overline{h}(T_{0}-\tau)D_{0}z_{0}(\tau)d\tau=i\omega_{o}Ae^{-\mu
T_{0}}\int_{0}^{T_{0}}\Big{(}&\lambda_{1}e^{(\mu+i\omega_{o})\tau}+\lambda_{2}e^{(\mu+i\omega_{o})\tau}\cos(\mu\tau-\mu
T_{0})\\\ &-\lambda_{3}e^{(\mu+i\omega_{o})\tau}\sin(\mu\tau-\mu
T_{0})\Big{)}d\tau,\end{split}$ (22)

which upon integration by parts yields

$\int_{0}^{T_{0}}\overline{h}(T_{0}-\tau)D_{0}z_{0}(\tau)d\tau=A\omega_{o}(\xi_{0}+i\bar{\xi}_{0})e^{i\omega_{o}T_{0}}+NST+cc,$
(23)

where $NST$ stands for non-secular terms, and

$\begin{split}\xi_{0}&=\left(\frac{\lambda_{1}\omega_{o}}{\mu^{2}+\omega_{o}^{2}}+\frac{2\lambda_{2}\mu^{2}\omega_{o}}{4\mu^{4}+\omega_{o}^{4}}-\frac{\lambda_{3}\omega_{o}^{3}}{4\mu^{4}+\omega_{o}^{4}}\right),\\\
\bar{\xi}_{0}&=\left(\frac{\lambda_{1}\mu}{\mu^{2}+\omega_{o}^{2}}+\frac{\lambda_{2}(2\mu^{3}-\mu\omega_{o}^{2})}{4\mu^{4}+\omega_{o}^{4}}-\frac{\lambda_{3}(2\mu^{3}+\mu\omega_{o}^{2})}{4\mu^{4}+\omega_{o}^{4}}\right).\end{split}$
(24)

Substituting Equations (23) and (24) into the third perturbation problem
presented in Equation (17), and using the following polar transformation for
the complex valued function $A(T_{2})$

$\begin{split}A(T_{2})&=\frac{a(T_{2})}{2}e^{i\beta(T_{2})},\\\
\overline{A}(T_{2})&=\frac{a(T_{2})}{2}e^{-i\beta(T_{2})},\end{split}$ (25)

then eliminating the secular terms, yields the following modulation equations:

$D_{2}a=-\left(\frac{\delta_{1}\bar{\xi}_{0}}{2}+\frac{\delta_{2}}{2}\right)a+\frac{g_{wave}}{2\omega_{o}}\sin\psi$
(26)
$aD_{2}\psi=\left(\sigma-\frac{\delta_{1}\xi_{0}}{2}\right)a+\left(\frac{5\eta^{2}}{12\omega_{o}^{3}}-\frac{3\gamma}{8\omega_{o}}\right)a^{3}+\frac{g_{wave}}{2\omega_{o}}\cos\psi,$
(27)

where $a$ and $\beta$ represent, respectively, the amplitude and phase of
oscillations, and $\psi=\sigma T_{2}-\beta$.

###  3.2 Global inter-well response

In this section, we seek an approximate solution for the high energy orbits;
i.e., when the system undergoes global inter-well oscillations. Since Equation
(10) has a negative linear stiffness, it is difficult to employ the method of
multiple scales in its straightforward fashion. To overcome this, we first
expand the natural frequency of symmetric oscillations to be in the following
form

$\omega_{N}^{2}=-\omega_{n}^{2}+\sigma_{1},$ (28)

where $\sigma_{1}$ is a detuning parameter. Since we are seeking an
approximate solution in the vicinity of the primary resonance of the inter-
well motions, we express the nearness of the wave frequency to the natural
frequency, $\omega_{N}$, by introducing the detuning parameter $\sigma_{2}$,
such that

$\Omega^{2}=\omega_{N}^{2}+\epsilon^{2}\sigma_{2}.$ (29)

Upon adding Equations (28 \- 29), we arrive at

$-\omega_{n}^{2}=\Omega^{2}-(\omega_{N}^{2}+\omega_{n}^{2})-\epsilon^{2}(\Omega^{2}-\omega_{N}^{2}).$
(30)

Substituting Equation (30) into Equation (10), we obtain

$\begin{split}\ddot{Y}+&\epsilon^{2}\delta_{1}\int_{0}^{t^{*}}h(t^{*}-\tau)\dot{Y}(\tau)d\tau+\epsilon^{2}\delta_{2}\dot{Y}+\Omega^{2}Y\\\
+&\epsilon(-(\omega_{N}^{2}+\omega_{n}^{2})Y+\gamma
Y^{3})-\epsilon^{2}(\Omega^{2}-\omega_{N}^{2})Y=\epsilon^{2}g_{wave}\cos(\Omega
t^{*}).\end{split}$ (31)

Next, we employ the method of multiple scales on Equation (31) by introducing
slow and fast time scales $T_{n}=\epsilon^{n}t^{*}$, $n=0,1,2$, and seek an
expansion of the response in the form

$\displaystyle Y(t^{*},\epsilon)$
$\displaystyle=Y_{0}(T_{0},T_{1},T_{2})+\epsilon
Y_{1}(T_{0},T_{1},T_{2})+\epsilon^{2}Y_{2}(T_{0},T_{1},T_{2})+O(\epsilon^{3}),$
(32a) $\displaystyle v(t^{*},\epsilon)$
$\displaystyle=v_{0}(T_{0},T_{1},T_{2})+\epsilon
v_{1}(T_{0},T_{1},T_{2})+\epsilon^{2}v_{2}(T_{0},T_{1},T_{2})+O(\epsilon^{3}).$
(32b)

Implementing the method of multiple scales on the scaled equations as
described in the previous subsection yields the following solution for $Y$ and
$v$.

$\displaystyle Y$ $\displaystyle=a\cos(\Omega t^{*}+\psi)+O(\epsilon),$ (33a)
$\displaystyle v$ $\displaystyle=\Gamma_{2}a\cos(\Omega
t^{*}+\psi)+O(\epsilon),$ (33b)

where

$\Gamma_{2}=\frac{\Omega^{2}+i\theta\Omega}{\Omega^{2}+\theta^{2}}$

and the amplitude, $a$, and phase, $\psi$, of the response are governed by the
following modulation equations:

$\omega_{n}D_{2}a=-\left(\frac{\omega_{n}\delta_{1}\bar{\xi_{1}}+\omega_{n}\delta_{2}}{2}\right)a-\frac{f_{wave}}{2}\sin\psi$
(34)
$\omega_{n}aD_{2}\psi=-\left(\frac{\Omega^{2}+\omega_{n}^{2}}{2}-\omega_{n}\delta_{1}\xi_{1}\right)a+\frac{3\gamma}{8}a^{3}+\frac{3\gamma^{2}}{256\omega_{n}^{2}}a^{5}-\frac{f_{wave}}{2}\cos\psi,$
(35)

where

$\begin{split}\xi_{1}&=\left(\frac{\lambda_{1}\omega_{n}}{\mu^{2}+\omega_{n}^{2}}+\frac{2\lambda_{2}\mu^{2}\omega_{n}}{4\mu^{4}+\omega_{n}^{4}}-\frac{\lambda_{3}\omega_{n}^{3}}{4\mu^{4}+\omega_{n}^{4}}\right)\\\
\bar{\xi}_{1}&=\left(\frac{\lambda_{1}\mu}{\mu^{2}+\omega_{n}^{2}}+\frac{\lambda_{2}(2\mu^{3}-\mu\omega_{n}^{2})}{4\mu^{4}+\omega_{n}^{4}}-\frac{\lambda_{3}(2\mu^{3}+\mu\omega_{n}^{2})}{4\mu^{4}+\omega_{n}^{4}}\right)\end{split}$
(36)

###  3.3 Steady-state response

The long-time steady-state behavior of the absorber is of particular interest
to assess its performance. Thus, the time derivatives in Equations (26 \- 27)
and (34 \- 35) are set to zero and the resulting algebraic equations are then
solved numerically for the steady-state amplitude $a_{o}$, and phase
$\psi_{o}$. For intra-well oscillations, we obtain:

$\displaystyle\begin{split}Y=a_{o}\cos(\Omega
t^{*}-\psi_{o})+\frac{\eta}{2\omega_{0}^{2}}\left(-a_{o}^{2}+\frac{a_{o}^{2}}{3}\cos(2\Omega
t^{*}-2\psi_{o})\right)+O(\epsilon^{2}),\end{split}$ (37a) $\displaystyle
v=\frac{\omega_{0}^{2}}{\omega_{0}^{2}+\theta^{2}}a_{o}\cos(\Omega
t^{*}-\psi_{o})-\frac{\omega_{0}\theta}{\omega_{0}^{2}+\theta^{2}}a_{o}\sin(\Omega
t^{*}-\psi_{o})+O(\epsilon^{2}),$ (37b)

and for the inter-well oscillations, we obtain

$\displaystyle\begin{split}Y&=a_{o}\cos(\Omega
t^{*}-\psi_{o})+\left(\frac{\gamma}{32\omega_{n}^{2}}a_{o}^{3}+\frac{3\gamma^{2}}{1024\omega_{n}^{4}}a_{o}^{5}\right)\cos(3\Omega
t^{*}-3\psi_{o})\\\
&+\frac{\gamma^{2}}{1024\omega_{n}^{4}}a_{o}^{5}\cos(5\Omega
t^{*}-5\psi_{o})+O(\epsilon^{2}),\end{split}$ (38a) $\displaystyle v$
$\displaystyle=\frac{\omega_{n}^{2}}{\omega_{n}^{2}+\theta^{2}}a_{o}\cos(\Omega
t^{*}-\psi_{o})-\frac{\omega_{n}\theta}{\omega_{n}^{2}+\theta^{2}}a_{o}\sin(\Omega
t^{*}-\psi_{o})+O(\epsilon^{2}).$ (38b)

Using the steady-state response, the average power available at the buoy can
be expressed as

$P_{avg}=\frac{1}{T^{*}}\int_{0}^{T^{*}}\delta_{2}\dot{Y}^{2}dt^{*},$ (39)

where $T^{*}$ is the period of oscillations. For intra-well oscillations
Equation (39) reduces to

$P_{avg}=\delta_{2}\left(\frac{\Omega^{2}a_{o}^{2}}{2}+\frac{\eta^{2}\Omega^{2}a_{o}^{4}}{18\omega_{0}^{4}}\right)+O(\epsilon^{2}),$
(40)

while for inter-well oscillations, we get

$P_{avg}=\delta_{2}\left(\frac{\Omega^{2}a_{o}^{2}}{2}+\frac{9\gamma^{2}\Omega^{2}a_{o}^{6}}{2048\omega_{n}^{4}}\right)+\
O(\epsilon^{2}).$ (41)

In addition to the averaged power, we are also interested in evaluating the
capture width ratio (CWR), which is a common parameter used to evaluate a
PWA’s performance. CWR or absorption width as commonly coined in some
literature is defined as the ratio between the average absorbed power,
$P_{avg}$, and the power available at the wave front, $P_{wave}$. The latter
can be obtained by multiplying the wave energy flux per unit crest length with
the buoy’s characteristic length, which is the diameter for a hemispherical
buoy. This yields [20]:

$CWR=\frac{6(m+m_{\infty})\Omega}{\rho RA^{2}}P_{avg}.$ (42)

##  4 Stability Analysis and Numerical Simulations

Bi-stable PWAs are known to produce large-amplitude responses over certain
frequency ranges under harmonic waves excitations. Unfortunately, such desired
motions can be incited only when the excitation is capable of channelling
enough energy to the PWA to overcome the potential energy barrier and perform
periodic inter-well motions. Even when excited, those desired motions can only
be uniquely realized over a specific bandwidth of the wave frequency, which we
coin here as the effective bandwidth of the PWA. Outside the effective
bandwidth, the large-amplitude motions are often accompanied with other, less
desirable responses; e.g., small periodic or aperiodic motions. To define this
effective bandwidth, it is important to study the stability of the steady-
state periodic solutions $a_{o}$ and $\psi_{o}$ as the excitation frequency is
varied. At the first level, this can be realized by finding the eigenvalues of
the Jacobian matrices associated with Equations (26) - (27) and (34)-(35).

Figure 6 depicts variation of the steady-state amplitude of the absorber with
the excitation frequency as obtained using equations (26) - (27) and
(34)-(35). The parameters used in the simulation are listed in Table 2 and the
associated symmetric potential energy function of the PWA is depicted in
Figure 7.

Table 2: WEC parameters. $Parameter$ | $Value$
---|---
$\omega_{n}$ | 0.78
$\gamma$ | 50
$\delta_{2}$ | 0.13

On the figure, the solid lines represent stable steady-state single-period
periodic solutions of amplitude $a_{o}$, and dashed lines represent unstable
unrealizable periodic solutions. It is evident that there are three branches
of stable periodic solutions based on the Jacobian-based stability analysis.
The branches $B_{r}$ and $B_{n}$ which represent, respectively, the resonant
and non-resonant branches of the small-amplitude intra-well motions, and the
branch $B_{L}$, which represents the large amplitude inter-well motions. The
stability analysis reveals two cyclic-fold bifurcations, $Cf_{1}$, and
$Cf_{2}$, which result from the stable and unstable periodic orbits colliding
and destructing each other.

Figure 6: Stroboscopic and analytical bifurcation diagrams for the bi-stable
PWA under regular wave excitation at nondimensional wave amplitude of 0.1.
Solid lines: stable solution. Dashed lines: unstable solutions. Figure 7:
Potential energy function associated with the system’s parameters listed in
Table 2.

The superimposed stroboscopic bifurcation map in Figure 6 obtained by a
numerical integration of the original equations of motion, Equation (1),
reveals a much more complex behaviour than what can be seen by relying on the
Jacobian-based stability analysis. In particular, The bifurcation map reveals
regions of aperiodic motions that extend over a wide range of frequencies.
While we notice good agreement between the analytical solution and the
numerical solution on the $B_{r}$ branch in the higher range of frequencies
down to about $\Omega=1.2$. The Jacobian-based stability analysis does not
reveal that the intra-well solutions on the branch, $B_{r}$, actually undergo
a cascade of period-doubling, $pd$, bifurcations starting near $\Omega=1.2$.
These bifurcations ultimately lead to a window of chaos, $CH$, which extends
down to about $\Omega=0.9$. This window ultimately disappears in a boundary
crisis. Phase portraits and the Fast Fourier Transform (FFT) of the period-
doubled and chaotic solutions are shown in Figure 8 clearly demonstrating the
period-doubling route to chaos.

$(a)$ | $(b)$ | $(c)$
---|---|---
|  |
|  |
Figure 8: Phase portraits and FFT spectra showing the instability route of the
intra-well branch. $(a):\Omega_{wave}=1$, $(b):\Omega_{wave}=1.2$, and
$(c):\Omega_{wave}=1.4$.

Furthermore, as shown in Figure 6, the period one solutions on branch $B_{L}$,
also lose stability via symmetry breaking bifurcations, $SB_{1}$ and $SB_{2}$
resulting in asymmetric periodic orbits as shown in the phase portraits and
FFTs of Figure 9 . In the region between $SB_{1}$ and $SB_{2}$, the desired
large orbit period one periodic solution is unique and, unlike the rest of the
frequency bandwidth considered, does not coexist with other less desirable
solutions. We coin this desired bandwidth as the effective bandwidth of the
PWA.

$(a)$ | $(b)$ | $(c)$
---|---|---
|  |
|  |
Figure 9: Phase portraits and FFT spectra showing the instability route of the
inter-well branch. $(a):\Omega_{wave}=0.42$, $(b):\Omega_{wave}=0.62$, and
$(c):\Omega_{wave}=0.8$.

Based on the numerical analysis, the response bandwidth of the absorber can be
divided into different regions as shown in Figure 10. The first region ($I$)
occurs in the low frequency range and contains inter-well chaotic motions
coexisting with small magnitude intra-well motions. The second region ($II$)
is the effective bandwidth, which contains the unique large-magnitude period
one inter-well motions. The third region ($III$) contains chaotic motions
coexisting with asymmetric periodic motions. Region four ($IV$) contains
chaotic motions only. Finally, region five ($V$) contains unique small-
amplitude period one intra-well periodic motions.

Figure 10: Stroboscopic diagram for the Bi-stable PWA showing regions with
different types of motion. Simulation performed at nondimensional wave
amplitude $A_{wave}/R$ = 0.1.

At lower wave amplitudes, the bi-stable PWA reveals a different frequency
response behavior as shown in Figure 11. The most important feature is the
birth of a new cyclic-fold bifurcation point denoted by $Cf_{3}$ and located
at the tip of the resonant branch, $B_{r}$. In addition, we notice that the
entire $B_{r}$ branch remains stable and the $pd$ bifurcation is shifted
towards the nonresonant intra-well branch, $B_{n}$. This arrangement creates a
region where two stable branches coexist, $B_{r}$ and $B_{n}$, which leads to
jumps between them as depicted in the enlarged window in Figure 11. In
addition, the stroboscopic bifurcation map points shows a chaotic region that
spans the entire frequency spectrum to the left of $Cf_{1}$. This Implies that
the large-amplitude branch $B_{L}$ breaks symmetry at $Cf_{1}$, and no
periodic inter-well motion can be realized.

Figure 11: Stroboscopic and analytical bifurcation diagrams for the bi-stable
PWA under regular wave excitation at nondimensional wave amplitude of 0.034.
Solid lines: stable solution. Dashed lines: unstable solutions.

Understanding the behavior of the absorber requires characterizing these
regions as function of the frequency and amplitude of the incident waves. This
can be realized by approximating the loci of the different bifurcations:
$Cf_{1}$, $Cf_{2}$, $Cf_{3}$ $pd$, $SB_{1}$ and $SB_{2}$. The first three
bifurcations can be approximated using the Jacobian stability analysis, while
the last three need the implementation of the Floquet theory. In what follows,
we obtain analytical approximation of all of these bifurcations as function of
the design parameters of the PWA.

###  4.1 Cyclic-fold bifurcation $Cf_{1}$ on the inter-well branch

The cyclic-fold point $C_{f1}$ satisfies the relation $d\Omega/da_{o}=0$ on
the resonant branch, $B_{L}$. Thus, its location can be obtained by
differentiating Equations (34-35) with respect to $a_{o}$ and setting
$d\Omega/da_{o}$ to zero to arrive at the following 8th order polynomial:

$\displaystyle\frac{45\gamma^{4}}{4096\omega_{n}^{4}}a_{b}^{8}+\frac{9\gamma^{3}}{128\omega_{n}^{2}}a_{b}^{6}+\left(\frac{99\gamma^{2}}{64}-\frac{3\gamma^{2}\Omega_{b}^{2}}{128\omega_{n}^{2}}+\frac{3\gamma^{2}\delta_{1}\xi_{1}}{64\omega_{n}}\right)a_{b}^{4}-3\gamma\left(\Omega_{b}^{2}+\omega_{n}^{2}-\delta_{1}\omega_{n}\xi_{1}\right)a_{b}^{2}$
(43)
$\displaystyle+\left(\Omega_{b}^{2}+\omega_{n}^{2}-\delta_{1}\omega_{n}\xi_{1}\right)^{2}+\left(\delta_{2}\omega_{n}+\delta_{1}\omega_{n}\bar{\xi}_{1}\\}\right)^{2}=0.$

where $(\Omega_{b},a_{b})$ represent the wave frequency and response amplitude
at which the bifurcation occurs. For any set of design parameters, the
polynomial in Equation (43) can be solved to obtain the locus of the $Cf_{1}$
bifurcation in the parameter space of wave amplitude and frequency.

###  4.2 Cyclic-fold bifurcations $Cf_{2}$ and $Cf_{3}$ on the intra-well
branch

These cyclic fold bifurcations also satisfy the relation $d\Omega/da_{o}=0$
but on the intra-well branch, $B_{n}$; i.e. Equations (26-27). Differentiating
Equations (26-27) with respect to $a_{o}$ and setting $d\Omega/da_{o}=0$, we
obtain

$a_{b}^{2}=\frac{-2\left(\Omega_{b}-\omega_{o}-\frac{\delta_{1}\xi_{0}}{2}\right)\pm\sqrt{\left(\Omega_{b}-\omega_{o}-\frac{\delta_{1}\xi_{0}}{2}\right)^{2}-3\left(\frac{\delta_{1}\bar{\xi}_{0}}{2}+\frac{\delta_{2}}{2}\right)^{2}}}{3\left(\frac{5\eta^{2}}{12\omega_{o}^{3}}-\frac{3\gamma}{8\omega_{o}}\right)}.$
(44)

Here, $a_{b}$ and $\Omega_{b}$ are, respectively, the amplitude of
oscillations and the wave frequency at which the cyclic-fold bifurcations,
$Cf_{2}$ and $Cf_{3}$, occurs. Upon solving Equation (44), we can obtain the
loci of the $Cf_{2}$ and $Cf_{3}$ bifurcations in the parameter space of wave
amplitude and frequency.

###  4.3 Period doubling bifurcation

One of the main dynamical characteristics of the intra-well oscillations is
the appearance of a period-doubling bifurcation on the $B_{r}$ branch as the
wave frequency is decreased toward lower values. To estimate this bifurcation
point as function of the wave frequency and amplitude, we perturb the periodic
orbit ${Y}(t^{*})$ by introducing the perturbation $\mathcal{p(t^{*})}$ such
that

$\tilde{Y}(t^{*})=Y(t^{*})+\mathcal{p(t^{*})}.$ (45)

Here, $\tilde{Y}(t^{*})$ is a stable periodic orbit if $\mathcal{p(t^{*})\to
0}$ as $t^{*}\to\infty$; otherwise it is unstable. Substituting Equation (45)
into Equation (11) and retaining only linear terms in $\mathcal{p(t^{*})}$, we
obtain, after simplifications, the following scaled equations:

$\begin{split}\mathcal{p}^{\prime\prime}+\epsilon\delta_{1}\int_{0}^{t^{*}}\overline{h}(t^{*}-\tau)\mathcal{p}^{\prime}(\tau)d\tau+\epsilon\delta_{2}\mathcal{p}^{\prime}&+[G_{0}+\epsilon
G_{1}\cos(\Omega t^{*}-\psi)+\epsilon^{2}G_{2}\cos(2\Omega t^{*}-2\psi)\\\
&+\epsilon^{2}G_{3}\cos(3\Omega t^{*}-3\psi)+\epsilon^{2}G_{4}\cos(4\Omega
t^{*}-4\psi)]\mathcal{p}=0,\end{split}$ (46)

where Equation (46) represents a Mathieu’s-type differential equation with
four parametric excitation terms. Each of these terms produces a principle
subharmonic parametric instability when their frequency is half the natural
frequency $\sqrt{G_{0}}$ of the perturbation [21, 22]. The interested reader
can refer to Appendix B for the full expressions of the parametric terms
constants, $G_{i}$.

In order to examine the dynamics of the perturbation $\mathcal{p(t^{*})}$, we
derive an asymptotic approximation for the evolution of $\mathcal{p(t^{*})}$
using the method of multiple scales assuming the following first-order
expansion:

$\mathcal{p}(t^{*},\epsilon)=\mathcal{p}_{0}(t^{*})+\epsilon\mathcal{p}_{1}(t^{*})+O(\epsilon^{2}).$
(47)

Using the time scales and time derivatives defined in Equation (12), and
substituting Equation (47) into Equation (46), we obtain the following
differential equations at the different orders of $\epsilon$:
$O(\epsilon^{0})$:

$D_{0}^{2}\mathcal{p}_{0}+G_{0}\mathcal{p}_{0}=0,$ (48)

which admits the following homogeneous solution

$\mathcal{p}_{0}=Q(T_{1})e^{i\sqrt{G_{0}}T_{0}}+\overline{Q}(T_{1})e^{-i\sqrt{G_{0}}T_{0}},$
(49)

and
$O(\epsilon^{1})$:

$D_{0}^{2}\mathcal{p}_{1}+G_{0}\mathcal{p}_{1}=-2D_{1}D_{0}\mathcal{p}_{0}-\delta_{1}\int_{0}^{T_{0}}\overline{h}(T_{0}-\tau)D_{0}\mathcal{p}_{0}d\tau-\delta_{2}D_{0}\mathcal{p}_{0}-G_{1}\mathcal{p}_{0}\cos(\Omega
t^{*}-\psi).$ (50)

Here, $Q(T_{1})$ and its complex conjugate $\overline{Q}(T_{1})$ are unknowns
that can be expressed in the following polar form:

$\begin{split}Q(T_{1})&=\frac{q(T_{1})}{2}e^{i\nu(T_{1})},\\\
\overline{Q}(T_{1})&=\frac{q(T_{1})}{2}e^{-i\nu(T_{1})},\end{split}$ (51)

where $q$ and $\nu$ are, respectively, the amplitude and phase of the
perturbation $\mathcal{p}$. Since we are interested in obtaining the point at
which the first period-doubling bifurcation occurs, we seek to approximate the
solution when $\Omega$ is near $2\sqrt{G_{0}}$. Thus, we express the proximity
of the excitation frequency to twice the natural frequency by introducing

$\Omega=2\sqrt{G_{0}}+\epsilon\sigma.$ (52)

Upon substituting (49) and (52) into Equation (50), then eliminating the
secular terms, we obtain the following equation which governs the locus of the
period-doubling bifurcation as function of the design parameters of the
absorber:

$(\sqrt{G_{0}}\Omega-2G_{0}-\delta_{1}\sqrt{G_{0}}\xi_{2})^{2}+(\delta_{1}\sqrt{G_{0}}\bar{\xi}_{2}+\delta_{2}\sqrt{G_{0}})^{2}=\frac{G^{2}_{1}}{4},$
(53)

where

$\begin{split}\xi_{2}&=\left(\frac{\lambda_{1}\sqrt{G_{0}}}{\mu^{2}+G_{0}}+\frac{2\lambda_{2}\mu^{2}G_{0}}{4\mu^{4}+G_{0}^{2}}-\frac{\lambda_{3}G_{0}^{\frac{2}{3}}}{4\mu^{4}+G_{0}^{2}}\right),\\\
\bar{\xi}_{2}&=\left(\frac{\lambda_{1}\mu}{\mu^{2}+G_{0}}+\frac{\lambda_{2}(2\mu^{3}-\mu
G_{0})}{4\mu^{4}+G_{0}^{2}}-\frac{\lambda_{3}(2\mu^{3}+\mu
G_{0})}{4\mu^{4}+G_{0}^{2}}\right).\end{split}$ (54)

###  4.4 Symmetry-break bifurcations of the inter-well branch

The other dynamical features that we have a particular interest in estimating
are the points of symmetry-break bifurcation on the symmetric inter-well
solution branch. To this end, we examine the stability of the approximate
expansion given in Equation (38) by introducing an infinitesimal perturbation
$\mathcal{p}(t^{*})$ to the periodic solution $Y(t^{*})$ as

$\tilde{Y}(t^{*})=Y(t^{*})+\mathcal{p}(t^{*}).$ (55)

Upon substituting Equation (55) into Equation (10), we obtain the following
Mathieu’s-type differential equation:

$\mathcal{p}^{\prime\prime}(t^{*})+\delta_{1}\int_{0}^{t^{*}}h(t^{*}-\tau)\mathcal{p}^{\prime}(\tau)d\tau+\delta_{2}\mathcal{p}^{\prime}(t^{*})+\left(K_{0}+\sum_{n=1}^{N}K_{n}\cos(2n\Omega
t^{*})\right)\mathcal{p}(t^{*})=0.$ (56)

Because of the presence of the term, $K_{n}\cos(2n\Omega t^{*})$ in Equation
(56), $\mathcal{p}(t^{*})$ must admit solutions that have even frequencies,
which breaks the symmetry of the original solution, $Y(t^{*})$. Since the
parametric terms are periodic with period $T=\pi/\Omega$, it follows by the
virtue of Floquet theory that $\mathcal{p}(t^{*})$ satisfies the following
equation:

$\mathcal{p}(t^{*}+T)=\lambda\mathcal{p}(t^{*}),$ (57)

where $\lambda$ is an eigenvalue called the Floquet multiplier. This
multiplier is an eigenvalue of a matrix $C$ associated with the solution
$\Phi(t^{*})$ of the Equation (56), which satisfies
$\Phi(t^{*}+T)=\Phi(t^{*})C$. The matrix $C$, also known as the monodromy
matrix, can be thought of as a transformation that maps $\Phi(0)$ to
$\Phi(T)$. Specifying the initial condition $\Phi(0)=I$ yields

$C=\Phi(T).$ (58)

In order to assess the stability of $\mathcal{p}(t^{*})$, we solve for
$\Phi(t^{*})$ starting at $t^{*}=0$ and ending at $t^{*}=T$. We then examine
the eigenvalues of the matrix $C=\Phi(T)$. The set of differential equations
governing the fundamental matrix solution can be constructed using the
following equation:

$\frac{d}{dt^{*}}\Phi(t^{*})=\left[\includegraphics[width=108.405pt,valign={c}]{matrix.pdf}\right]\Phi(t^{*}),$
(59)

where $A_{r}$, $B_{r}$ and $C_{r}$ are the realized state-space accounting for
radiation damping, and $f(t^{*})$ is the parametric excitation term given as:

$f(t^{*})=K_{0}+\sum_{n=1}^{N}K_{n}\cos(2n\Omega t^{*}).$ (60)

The reader can refer to Appendix B for the constants $K_{i}$. Equation (59)
results in a set of 25 linear differential equations subjected to the initial
conditions $\Phi(0)=I$. Here, $I$ is an identity matrix of dimension 5.

Integrating Equation (59) numerically in [0, $T$] using the initial conditions
of $\Phi(0)=I$, finding the Floquet multipliers, $\lambda$, of the resulting
numerical matrix, then inspecting their location with respect the unit circle,
we can find the loci of the symmetry breaking bifurcations in the wave
amplitude versus frequency parameter space. In particular, a symmetry breaking
bifurcation occurs when one of the Floquet multipliers exits the unit circle
through $\lambda=1$.

###  4.5 Bifurcation diagram

The bifurcation diagram based on the aforedescribed stability analysis is
shown in Figure 12 (a). It is evident that the bifurcations $SB_{1}$,
$SB_{2}$, and $pd$ are in excellent agreement with the points where the
stroboscopic numerical bifurcation maps undergo qualitative changes. The
average power and the CWR curves shown in Figures 12 (b-c) also illustrate
good agreement with the analytical solution in the range when the response
exhibits unique period one periodic motion. It is evident that maximum average
power and CWR are realized within the effective bandwidth of the absorber and
that outside that bandwidth the average power drops whether the response is of
the chaotic or intra-well type.

$(a)$
---

$(b)$ | $(c)$
---|---
|

Figure 12: Bi-stable PWA response under regular wave excitation at
nondimensional wave amplitude $A_{wave}/R=0.1$. (a): Analytical and
stroboscopic bifurcation map. (b): Nondimensional averaged power. (c): CWR.
(Solid lines represent stable solutions. Dashed lines represent unstable
solutions based on the Jacobian matrix criterion. Red dashed-dotted lines
represent unstable solutions based on period-doubling instability analysis.
Blue dashed-dotted lines represent unstable solutions based on Floquet
analysis. Circles represent the numerical results).

##  5 The Effective Bandwidth

In this section, we use the stability analysis furnished in the previous
sections to define an effective bandwidth for the bi-stable wave energy
absorber by marking the boundaries at which the PWA switches to different
types of motion in the wave amplitude versus frequency $(A_{wave}/R,\Omega)$
parameters space. We also draw a clearer picture of how the shape of the bi-
stable potential influences the effective bandwidth.

Using the loci of the different bifurcations, we create the design map shown
in Figure 13 which characterizes the type of motion realized for every
combination of wave amplitude and frequency. The largest region on the figure
is the one denoted by $B_{r}$, which represents $(A_{wave}/R,\Omega)$
combinations that incite small-amplitude inter-well motions. We can see that
this type of response occurs for any frequency when the wave amplitude is
small and for any wave amplitude when the wave frequency is large. In terms of
size, the second region on the figure is that which results in coexisting
chaotic, $CH$, and $n$-period, $nT$, periodic solutions. This region occurs
slightly to the left of $\Omega=1$ and increases in bandwidth as $A_{wave}/R$
is increased. The basin of attraction of the different coexisting solutions in
this region are shown in Figure 14 for the two points indicated on the map. We
can see that there are three different basins: the black one represents
symmetric period one motions, the white one represents asymmetric period one
motions, and the one where the black and white colors blend together
represents chaotic motions.

The most important region on the map is that denoted by $B_{L}$, which
corresponds to the combination of wave parameters leading to a unique large-
amplitude inter-well motion. We see that this region exists around $\Omega=1$
for wave amplitudes larger than $A_{wave}/R=0.05$. The size of the effective
bandwidth of the absorber increases as the wave amplitude is increased up to a
value of $A_{wave}/R=0.125$ beyond which the size of the bandwidth remains
almost constant. The $B_{L}$ region is bordered from the right by the region,
$B_{L}+CH$, where large-amplitude periodic motions coexist with chaotic
motions, and from below by the region, $CH$, where the chaotic attractor is
unique. The region denoted by $CH+B_{L}+B_{n}$ corresponds to wave parameters
that results in a chaotic attractor coexisting with two types of periodic
orbits: $B_{L}$ and $B_{n}$ each with competing basins of attraction.

Figure 13: A map demarcating regions of quantitatively different PWA
responses. Figure 14: Basins of attraction showing different coexisting
orbits.

###  5.1 Influence of the shape of the potential energy function

In this section, we use the understanding developed in the previous section to
examine the influence of the potential function on the effective bandwidth of
the bi-stable PWA. In particular, we want to understand how the depth of the
potential wells and their separation influences the size of the effective
bandwidth and the other regions in the bifurcation map. To this end, we
consider the three different potential energy functions shown in Figure 15,
which were obtained using $\gamma=$30, 50, and 90. The potential energy
function associated with $\gamma=$30 is deeper with a larger separation
between the stable equilibria, while the potential energy function associated
with $\gamma=$90 has shallower potential wells and smaller separation.

The bifurcation maps showing the effective bandwidth associated with each of
the potential energy functions considered are depicted in Figure 16. In order
to make a quantitative comparison, three critical wave amplitudes were marked
on the figures, and labeled as $(A_{wave}/R)_{cr1}$, $(A_{wave}/R)_{cr2}$, and
$(A_{wave}/R)_{cr3}$. The first critical amplitude occurs at the intersection
between the $Cf_{1}$ and $SB_{1}$ curves and can be used to define the wave
amplitude at which the effective bandwidth of the absorber approaches its
maximum size. We notice that this critical level of wave amplitude increases
as the potential wells become deeper; that is, larger excitation levels become
necessary to attain the unique periodic large-amplitude motions when
increasing the depth and the separation distance between the potential energy
wells.

It is interesting to note that, for any wave amplitude above
$(A_{wave}/R)_{cr1}$, the size of the effective bandwidth remains almost
constant regardless of the shape of the potential energy function. However,
while the effective bandwidth remains unchanged, the power levels within the
effective bandwidth change considerably with the shape of the potential energy
function as shown in Figure 17. The bi-stable PWA with the deeper potential
energy wells produces higher average power levels within the effective
bandwidth (the region bounded between $SB_{1}$ and $SB_{2}$).

This critical level of excitation marked by $(A_{wave}/R)_{cr2}$ occurs at the
intersection between the $Cf_{1}$ and $pd$ lines. It represents the minimum
value of wave amplitude necessary to generate unique inter-well motions. It is
evident that this amplitude decreases as the potential wells become shallower.
Finally, the amplitude level $(A_{wave}/R)_{cr3}$, represents the wave
amplitude below which no bifurcations occur as the wave frequency is varied.
At such a low level of wave amplitude, the response of the PWA resembles the
bell-shaped response of the traditional linear PWA.

Figure 15: Potential energy function of the bi-stable wave energy absorber at
different values of $\gamma$ and $\omega_{n}=0.78$. (1): $\gamma=30$, (2):
$\gamma=50$, (3): $\gamma=90$.

(a) $\gamma=30$

(b) $\gamma=50$

(c) $\gamma=90$

Figure 16: Bifurcation maps in the wave amplitude versus frequency parameter
space for the different potential functions shown in Figure 15.

---
$(a)$: $\gamma=30$ |  $(b)$: $\gamma=50$
---|---
|

$(c)$: $\gamma=90$
---

Figure 17: Comparison of the average generated power in the wave amplitude -
frequency parameter space for the different potential functions shown in
Figure 15. The numerical simulations were performed on Equation (10) at
initial conditions $(Y_{0},\dot{Y}_{0})$ of $(0,0)$.

##  6 Conclusion

This paper presented a theoretical analytical analysis of the response of bi-
stable PWAs to harmonic wave excitations. To this end, approximate asymptotic
solutions of the governing equations of motion were derived by implementing
the method of multiple scales. A stability analysis of the attained solutions
revealed the presence of key bifurcations that can be used to define an
effective bandwidth of the generator. This effective bandwidth is
characterized by the presence of a unique large-orbit inter-well motion for a
set of wave amplitudes and frequencies. This effective bandwidth occurs
slightly below the resonant frequency of the absorber and exists only above a
certain threshold in the wave amplitude. This threshold increases as the depth
of the potential well is increased. The size of the effective bandwidth
increases as the wave amplitude is increased up to a certain threshold above
which the effective bandwidth remains almost constant even when the wave
amplitude is substantially increased. The size of the effective bandwidth is
observed to be insensitive to variations in the depth of the potential well of
the absorber. However, the power levels within the effective bandwidth change
considerably with the shape of the potential energy function. In particular, a
bi-stable PWA with deeper potential energy wells produces higher average power
levels within its effective bandwidth. It is our belief that this
comprehensive analytical treatment is key to designing effective bi-stable
PWAs for known wave conditions and provide backbone results for future studies
addressing more realistic regular non-harmonic wave excitations.

## Funding

This research was funded by Abu Dhabi Education and Knowledge Council (ADEK)
under grant number AARE2019-161: Exploiting Bi-stability to Develop a Novel
Broadband Point Wave Energy Absorber.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data availability

The data that support the findings of this study will be made available upon
reasonable request.

##  Appendix A Eigensystem realization algorithm

Consider the following single-input single-output discrete-time dynamical
system:

$\begin{split}\textbf{x}_{k+1}&=\textbf{A}\textbf{x}_{k}+\textbf{B}u_{k},\\\
h_{k}&=\textbf{C}\textbf{x}_{k}+\textbf{D}u_{k},\end{split}$ (A.1)

and a discrete-time scalar input $u$:

$u_{k}^{\delta}\equiv u^{\delta}(k\Delta t)=\begin{cases}1,&\text{if $k=0$}\\\
0,&\text{if $1,\dots\infty$}\end{cases}$ (A.2)

According to linear system theory, the discrete-time impulse response data
$h_{k}^{\delta}$ could be expressed as:

$h_{k}^{\delta}\equiv h^{\delta}(k\Delta
t)=\textbf{C}\textbf{A}^{k}\textbf{B},\hskip
14.22636pt\left(k=0,1,\dots,\infty\right)$ (A.3)

In our analysis, the discrete-time impulse response data are obtained through
substituting the radiation damping coefficients $B(\omega_{i})$ into Equation
(3). The next step is to proceed by constructing the generalized Hankel matrix
$\textbf{H}_{r\times s}(g)$ for $g=0,1$ which consists of $r$ rows and $s$
columns, generated by stacking time-shifted impulse response data in the
following order:

$\textbf{H}_{r\times
s}(g)=\begin{pmatrix}h_{g}^{\delta}&h_{g+1}^{\delta}&\dots&h_{g+s-1}^{\delta}\\\
h_{g+1}^{\delta}&h_{g+2}^{\delta}&\dots&h_{g+s}^{\delta}\\\
\vdots&\vdots&\ddots&\vdots\\\
h_{g+r-1}^{\delta}&h_{g+r}^{\delta}&\dots&h_{g+r+s-2}^{\delta}\end{pmatrix}$
(A.4)

Using Equation (A.3), we can express the generalized Hankel matrix in Equation
(A.4) in terms of the realized state-space matrices $\textbf{A}_{r}$,
$\textbf{B}_{r}$ and $\textbf{C}_{r}$ as:

$\textbf{H}_{r\times
s}(g)=\begin{pmatrix}\textbf{C}_{r}\textbf{A}_{r}^{g}\textbf{B}_{r}&\textbf{C}_{r}\textbf{A}_{r}^{g+1}\textbf{B}_{r}&\dots&\textbf{C}_{r}\textbf{A}_{r}^{g+s-1}\textbf{B}_{r}\\\
\textbf{C}_{r}\textbf{A}_{r}^{g+1}\textbf{B}_{r}&\textbf{C}_{r}\textbf{A}_{r}^{g+2}\textbf{B}_{r}&\dots&\textbf{C}_{r}\textbf{A}_{r}^{g+s}\textbf{B}_{r}\\\
\vdots&\vdots&\ddots&\vdots\\\
\textbf{C}_{r}\textbf{A}_{r}^{g+r-1}\textbf{B}_{r}&\textbf{C}_{r}\textbf{A}_{r}^{g+r}\textbf{B}_{r}&\dots&\textbf{C}_{r}\textbf{A}_{r}^{g+r+s-2}\textbf{B}_{r}\end{pmatrix}$
(A.5)

which could be reduced as:

$\textbf{H}_{r\times s}(g)=\mathcal{O}\textbf{A}_{r}^{g}\mathcal{C}$ (A.6)

where

$\displaystyle\mathcal{O}$ $\displaystyle=\left(\textbf{C}_{r}\hskip
14.22636pt\textbf{C}_{r}\textbf{A}_{r}\hskip 14.22636pt\dots\hskip
14.22636pt\textbf{C}_{r}\textbf{A}_{r}^{r-1}\right)^{T}$
$\displaystyle\mathcal{C}$ $\displaystyle=\left(\textbf{B}_{r}\hskip
14.22636pt\textbf{A}_{r}\textbf{B}_{r}\hskip 14.22636pt\dots\hskip
14.22636pt\textbf{A}_{r}^{r-1}\textbf{B}_{r}\right)^{T}$

are respectively the generalized observability and controllability matrices,
with observability and controllability indices of $r$ and $s$. Upon taking the
singular value decomposition (SVD) for the first Hankel matrix
$\textbf{H}_{r\times s}(0)$ we get the following definition:

$\displaystyle\textbf{H}_{r\times s}(0)$
$\displaystyle=\textbf{U}\Sigma\textbf{V}^{T}$ (A.7)
$\displaystyle=\begin{pmatrix}\tilde{\textbf{U}}&\textbf{U}_{t}\\\
\end{pmatrix}\begin{pmatrix}\tilde{\Sigma}&0\\\ 0&\Sigma_{t}\\\
\end{pmatrix}\begin{pmatrix}\tilde{\textbf{V}}^{T}\\\
\tilde{\textbf{V}}_{t}^{T}\\\ \end{pmatrix}$
$\displaystyle\approx\tilde{\textbf{U}}\tilde{\Sigma}\tilde{\textbf{V}}^{T}$

where,

$\displaystyle\tilde{\textbf{U}}^{T}\tilde{\textbf{U}}$
$\displaystyle=\textbf{I}$
$\displaystyle\tilde{\textbf{V}}^{T}\tilde{\textbf{V}}$
$\displaystyle=\textbf{I}$ $\displaystyle\tilde{\Sigma}$
$\displaystyle=\begin{pmatrix}\sigma_{1}&&&\\\ &\sigma_{2}&&\\\ &&\ddots&\\\
&&&\sigma_{N}\\\ \end{pmatrix}$

The diagonal matrix $\tilde{\Sigma}$ which is constructed from the first
$N\times N$ block of $\Sigma$ contains the dominant singular values
$\sigma_{i}$ in the following order
$(\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{N}\geq 0)$. While $\Sigma_{t}$
contains the small truncated singular values. This truncation step is very
vital for reducing the order of the realized state-space model; such that the
realized dynamics matrix $A_{r}$ has size $N$. Also, vectors in
$\tilde{\textbf{U}}$ and $\tilde{\textbf{V}}^{T}$ contain the dominant modes
associated with the singular values retained in $\tilde{\Sigma}$. As a result,
the product $\tilde{\textbf{U}}\tilde{\Sigma}\tilde{\textbf{V}}^{T}$ is
considered to be a faithful representation of the original Hankel matrix H for
the smallest size of $\tilde{\Sigma}$.
It follows from Equation (A.6) that

$\displaystyle\textbf{H}_{r\times s}(0)$
$\displaystyle=\tilde{\textbf{U}}\tilde{\Sigma}\tilde{\textbf{V}}^{T}$ (A.8)
$\displaystyle=\left(\tilde{\textbf{U}}\tilde{\Sigma}^{\frac{1}{2}}\right)\left(\tilde{\Sigma}^{\frac{1}{2}}\tilde{\textbf{V}}^{T}\right)=\mathcal{O}\mathcal{C}$

Using the above balanced decomposition of $\textbf{H}_{r\times s}(0)$ we can
write:

$\mathcal{O}=\tilde{\textbf{U}}\tilde{\Sigma}^{\frac{1}{2}}\hskip
28.45274ptand\hskip
28.45274pt\mathcal{C}=\tilde{\Sigma}^{\frac{1}{2}}\tilde{\textbf{V}}^{T}$

Also from Equation (A.6) we can express the second Hankel matrix
$\textbf{H}_{r\times s}(1)$ as:

$\displaystyle\textbf{H}_{r\times s}(1)$
$\displaystyle=\mathcal{O}\textbf{A}_{r}\mathcal{C}$ (A.9)
$\displaystyle=\left(\tilde{\textbf{U}}\tilde{\Sigma}^{\frac{1}{2}}\right)\textbf{A}_{r}\left(\tilde{\Sigma}^{\frac{1}{2}}\tilde{\textbf{V}}^{T}\right)$

Using the properties of U and V we can write:

$\tilde{\Sigma}^{\frac{1}{2}}\textbf{A}_{r}\tilde{\Sigma}^{\frac{1}{2}}=\tilde{\textbf{U}}^{T}\textbf{H}_{r\times
s}(1)\tilde{\textbf{V}}$ (A.10)

It follows that the matrix $\textbf{A}_{r}$ could be obtained through:

$\textbf{A}_{r}=\tilde{\Sigma}^{-\frac{1}{2}}\tilde{\textbf{U}}^{T}\textbf{H}_{r\times
s}(1)\tilde{\textbf{V}}\tilde{\Sigma}^{-\frac{1}{2}}$ (A.11)

Let:

$\textbf{E}^{T}_{1}=\begin{pmatrix}1&0&\dots&0\\\ \end{pmatrix}\hskip
28.45274pt\textbf{E}^{T}_{2}=\begin{pmatrix}1&0&\dots&0\\\ \end{pmatrix}$

where, $\textbf{E}^{T}_{1}$ and $\textbf{E}^{T}_{2}$ are respectively $1\times
r$ and $1\times s$ vectors. We use that along with Equation (A.5) to write the
following balanced expression for $h_{k}^{\delta}$:

$\displaystyle h_{k}^{\delta}$
$\displaystyle=\textbf{E}^{T}_{1}\textbf{H}_{r\times s}(g)\textbf{E}^{T}_{2}$
(A.12)
$\displaystyle=\textbf{E}^{T}_{1}(\mathcal{O}\textbf{A}_{r}^{g}\mathcal{C})\textbf{E}_{2}$
$\displaystyle=(\textbf{E}^{T}_{1}\tilde{\textbf{U}}\tilde{\Sigma}^{\frac{1}{2}})(\tilde{\Sigma}^{-\frac{1}{2}}\tilde{\textbf{U}}^{T}\textbf{H}_{r\times
s}(1)\tilde{\textbf{V}}\tilde{\Sigma}^{-\frac{1}{2}})^{g}(\tilde{\Sigma}^{\frac{1}{2}}\tilde{\textbf{V}}^{T}\textbf{E}_{2})$
$\displaystyle\equiv\textbf{C}_{r}\textbf{A}_{r}^{g}\textbf{B}_{r}$

We use the decomposed expression above to obtain the reduced input and output
matrices $\textbf{B}_{r}$ and $\textbf{C}_{r}$ as:

$\textbf{B}_{r}=\tilde{\Sigma}^{\frac{1}{2}}\tilde{\textbf{V}}^{T}\textbf{E}_{2}$
(A.13)
$\textbf{C}_{r}=\textbf{E}^{T}_{1}\tilde{\textbf{U}}\tilde{\Sigma}^{\frac{1}{2}}$
(A.14)

##  Appendix B Parametric terms constants

$\displaystyle G_{0}$
$\displaystyle=\omega_{o}^{2}-\frac{\eta^{2}}{\omega_{o}^{2}}a_{o}^{2}+\frac{3\gamma}{2}a_{o}^{2}+\frac{3\gamma\eta^{2}}{4\omega_{o}^{4}}a_{o}^{4}+\frac{\gamma\eta^{2}}{24\omega_{o}^{4}}a_{o}^{4}$
(B.1) $\displaystyle G_{1}$ $\displaystyle=2\eta
a_{o}-\frac{5\gamma\eta}{2\omega_{o}^{2}}a_{o}^{3}$ (B.2) $\displaystyle
G_{2}$
$\displaystyle=\frac{\eta^{2}}{3\omega_{o}^{2}}a_{o}^{2}+\frac{3\gamma}{2}a_{o}^{2}-\frac{\gamma\eta^{2}}{2\omega_{o}^{4}}a_{o}^{4}$
(B.3) $\displaystyle G_{3}$
$\displaystyle=\frac{\gamma\eta}{2\omega_{o}^{2}}a_{o}^{3}$ (B.4)
$\displaystyle G_{4}$
$\displaystyle=\frac{\gamma\eta^{2}}{24\omega_{o}^{4}}a_{o}^{4}$ (B.5)
$\displaystyle K_{0}$
$\displaystyle=-\omega_{n}^{2}+3\gamma\left(\frac{R_{1}}{2}+\frac{R_{3}^{2}}{2}+\frac{R_{5}^{2}}{2}\right)$
(B.6) $\displaystyle K_{2}$
$\displaystyle=3\gamma\left(\frac{R_{1}^{2}}{2}+R_{1}R_{3}+R_{3}R_{5}\right)$
(B.7) $\displaystyle K_{4}$
$\displaystyle=3\gamma\left(R_{1}R_{3}+R_{1}R_{5}\right)$ (B.8) $\displaystyle
K_{6}$ $\displaystyle=3\gamma\left(\frac{R_{3}^{2}}{2}+R_{1}R_{5}\right)$
(B.9) $\displaystyle K_{8}$ $\displaystyle=3\gamma\left(R_{3}R_{5}\right)$
(B.10) $\displaystyle K_{10}$
$\displaystyle=3\gamma\left(\frac{R_{5}^{2}}{2}\right)$ (B.11)

where:

$\displaystyle R_{1}$ $\displaystyle=a_{o}$ (B.12) $\displaystyle R_{3}$
$\displaystyle=\frac{\gamma}{32\Omega^{2}}a_{o}^{3}+\frac{3\gamma^{2}}{1024\Omega^{4}}a_{o}^{5}$
(B.13) $\displaystyle R_{5}$
$\displaystyle=\frac{\gamma^{2}}{1024\Omega^{4}}a_{o}^{5}$ (B.14)

## References

  * [1] Elie Al Shami, Ran Zhang, and Xu Wang.  Point absorber wave energy harvesters: A review of recent developments.  Energies, 12(1):47, 2019.
  * [2] Johannes Falnes and Jørgen Hals.  Heaving buoys, point absorbers and arrays.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1959):246–277, January 2012.
  * [3] Benjamin Drew, Andrew R Plummer, and M Necip Sahinkaya.  A review of wave energy converter technology.  Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 223(8):887–902, 2009.
  * [4] Davood Younesian and Mohammad-Reza Alam.  Multi-stable mechanisms for high-efficiency and broadband ocean wave energy harvesting.  Applied energy, 197:292–302, 2017.
  * [5] Benjamin W Schubert, William SP Robertson, Benjamin S Cazzolato, Mergen H Ghayesh, and Nataliia Y Sergiienko.  Performance enhancement of submerged wave energy device using bistability.  Ocean Engineering, 213:107816, 2020.
  * [6] Mohammed F Daqaq, Ravindra Masana, Alper Erturk, and D Dane Quinn.  On the role of nonlinearities in vibratory energy harvesting: a critical review and discussion.  Applied Mechanics Reviews, 66(4), 2014.
  * [7] Ru Xi, Haicheng Zhang, Huai Zhao, Ramnarayan Mondal, et al.  High-performance and robust bistable point absorber wave energy converter.  Ocean Engineering, 229:108767, 2021.
  * [8] Xiaolong Xiao, Longfei Xiao, and Tao Peng.  Comparative study on power capture performance of oscillating-body wave energy converters with three novel power take-off systems.  Renewable Energy, 103:94–105, 2017.
  * [9] Haicheng Zhang, Ru Xi, Daolin Xu, Kai Wang, Qijia Shi, Huai Zhao, and Bo Wu.  Efficiency enhancement of a point wave energy converter with a magnetic bistable mechanism.  Energy, 181:1152–1165, 2019.
  * [10] Xian-tao Zhang, Jian-min Yang, and Long-fei Xiao.  An oscillating wave energy converter with nonlinear snap-through power-take-off systems in regular waves.  China Ocean Engineering, 30(4):565–580, 2016.
  * [11] Xiantao Zhang, Xinliang Tian, Longfei Xiao, Xin Li, and Lifen Chen.  Application of an adaptive bistable power capture mechanism to a point absorber wave energy converter.  Applied Energy, 228:450–467, 2018.
  * [12] Xiantao Zhang, XinLiang Tian, Longfei Xiao, Xin Li, and Wenyue Lu.  Mechanism and sensitivity for broadband energy harvesting of an adaptive bistable point absorber wave energy converter.  Energy, 188:115984, 2019.
  * [13] Yang Song, Xiaoxian Guo, Hongchao Wang, Xinliang Tian, Handi Wei, and Xiantao Zhang.  Performance analysis of an adaptive bistable point absorber wave energy converter under white noise wave excitation.  IEEE Transactions on Sustainable Energy, 2020.
  * [14] A Hulme.  The wave forces acting on a floating hemisphere undergoing forced periodic oscillations.  Journal of Fluid Mechanics, 121:443–463, 1982.
  * [15] T Francis Ogilvie and Ernest O Tuck.  A rational strip theory of ship motions: part i.  Technical report, University of Michigan, 1969.
  * [16] MD Haskind.  The exciting forces and wetting of ships in waves.  report, 2010.
  * [17] John Nicholas Newman.  The exciting forces on fixed bodies in waves.  Journal of ship research, 6(04):10–17, 1962.
  * [18] Ali H Nayfeh.  Perturbation methods.  John Wiley & Sons, 2008.
  * [19] Steven L Brunton and J Nathan Kutz.  Data-driven science and engineering: Machine learning, dynamical systems, and control.  Cambridge University Press, 2019.
  * [20] Liguo Wang.  Modelling and advanced control of fully coupled wave energy converters subject to constraints: the wave-to-wire approach.  PhD thesis, Acta Universitatis Upsaliensis, 2017.
  * [21] Robert A Meyers.  Mathematics of complexity and dynamical systems.  Springer Science & Business Media, 2011.
  * [22] Ivana Kovacic, Richard Rand, and Si Mohamed Sah.  Mathieu’s equation and its generalizations: overview of stability charts and their features.  Applied Mechanics Reviews, 70(2), 2018.