Enhancement of a 3-DOF submerged wave energy device using bistability

The dynamic response of a submerged CETO shaped quasi-point absorbing wave energy converter coupled to a bistable power take off is presented in this study. Whilst the impact of bistability has been shown in a limited number of situations to improve the amount of power generated, many models have been restricted to a single degree of freedom and often ignore drag effects. To overcome these model limitations, a submerged single tether point absorber with a bistable power take off was modelled using both 1 and 3 degrees of freedom. The device was subjected to regular waves and included a simple model of viscous drag. The bistable mechanism was provided by a magnetic dipole model quantified by a dimensionless parameter applicable to any bistable system. The performance of the device was is assessed by the theoretical power generated. Over each model, the previously observed benefit of bistability was not consistently obtained. Simulations of regular waves demonstrated an increase in generated power for suboptimal conditions for some frequencies, while a reduction in generated power was observed in optimal conditions. The performance increase showed strong correlation to the phase relationship between the motion and exciting forces as a result of bistability.


I. INTRODUCTION
O CEAN wave energy has been the subject of over two centuries of research according to [1]. Over this time many individual designs have been proposed with differing modes of operation. A number of typical wave energy converters (WEC) are discussed in [2] and can be broadly classed as one of three types: an attenuator, a point absorber (PA), or a terminator. This paper focuses on a PA type WEC, which are systems in which the buoy is small relative to the wavelengths of the incident waves, and are subsequently relatively insensitive to wave direction. A representative diagram of a point absorber with a generic buoy is shown in Fig. 1. While the simple operation of a generic PA WEC is well known, there remains many challenges for wave energy, as identified in [3]. Areas in which further work is required before commercialisation are: materials and manufacture; fluid dynamics and hydrodynamics; survivability and reliability; environmental resources; devices and arrays; power conversion and control; infrastructure and grid connection; marine operations 1559 Wave device development and testing. This work was supported by the Australian Government Research Training Program Scholarship.
The authors are with the School of Mechanical Engineering, The University of Adelaide, Adelaide, South Australia, 5005 Australia.
The email of the corresponding author, B. W. Schubert, is benjamin.schubert@adelaide.edu.au.
Power Take O¡ Generic Buoy and maritime safety; socio-economic implications; and environmental impact. Due to the complexity of these unique challenges in the context of wave energy, the technology as a whole remains at a low technology readiness level, a measure used to quantify progress to commercialisation.
Within the area of power conversion and control, hybrid control, in which passively applied dynamics supplements the active controller, may provide some advantages in complexity over active control algorithms. References [4]- [7] propose different bistable mechanisms within the power take off unit (PTO) and show for simple single degree of freedom (DOF) models that a substantial benefit is seen. Methods to create a nonlinear stiffness without bistability have been preposed such as in [8] in which an assymetric mass distribution is used to change the effective stiffness. A bistable mechanism using magnets is given in [6] which, being an electromagnetic system, has the advantage of having relatively small reactive power losses compared to other reactive mechanisms, such as hydraulic systems. Furthermore, a multi-DOF device and the potential for bistability to improve the efficiency of WEC devices was investigated in [9]. Experimental work within [10] and [11] demonstrated on different devices that the realistic benefit of bistability applied to is significant and in some cases may be as much as a three times power increase. A method to optimize a nonlinear stiffness has been presented in [12]. Each device and mechanism is different and bistability is employed in various ways to either excite super harmonics to convert low frequency oscillations into high frequency oscillations, such as in [13], or to passively phase match floating buoys, such as in [11]. Of the work presented, only [9] is related to submerged WECs and does not operate in the vertical heave direction. In addition to the restricted DOF, many of the proposed models do not include even simplistic drag; how these model deficiencies impact the benefit of bistability, as well as the mechanism behind the benefit within the context of submerged WECs, remains a gap within current literature.
To advance the technology, a whole system approach is suggested in [14], where an understanding of all parts of the complex system are considered and an understanding of interaction between components is developed. With this guiding principle, the impact of bistable control on a submerged 3-DOF PA WEC was explored in the present study to quantify how simple hydrodynamics and the passive control system interact. Four models were constructed: an optimal 1-DOF model, a suboptimal 1-DOF model, an optimal 3-DOF model, and a suboptimal 3-DOF model, all subjected to a set of regular waves. The mathematical models are described in Section II, with the representation of the magnetic bistable mechanism given in Section II-A4, as well as a dimensionless parameter relevant to any bistable system as a means to generalize results. The considerations and limitations of the used modelling method are outlined in Section III, and the results of these simulations are presented in Section IV. These results are discussed in detail in Section V, with a summary of the findings in Section VI.

II. MATHEMATICAL MODELS
The governing equations for the 1-DOF and 3-DOF systems can be derived from the contributing forces represented by: where F represents a force and the subscripts e, r, h, PTO, D, and bi indicate the excitation, radiation, hydrostatic, PTO, drag, and bistable forces, respectively. The inertial mass of each DOF of the system in represented by the matrix M, and x is the position vector containing the surge (x), heave (z), and pitch (θ) coordinates. This equation forms the basis for the governing equations for all three models presented with constituting force components either simplified or represented differently.

A. Single DOF
By simplifying (1) to 1-DOF and ignoring the viscous drag force, the following governing equation can be constructed to represent a simple WEC model: where the subscript z represents only the heave component. A time domain representation of each force is provided in the following sections.

1) Excitation force
The excitation force for the 1-DOF system, F e,z , can be calculated using linear potential flow theory and is the superposition of the Froude-Krylov and diffraction forces. This force may be considered as the force on the buoy due to incoming waves. In the work presented here the boundary element method (BEM) solver NEMOH, described in [15], was used to calculate the amplitude (F e ) and phase (φ) of the excitation force (in all 3-DOF) and can be represented as where ω and t represent the excitation frequency and time, respectively.

2) Radiation force
The radiation force can be considered as the force on the buoy due to the waves radiated as the buoy moves in calm water. This force is represented in the time domain by the Cummins equation as [16], where A ∞ and K are termed the infinite frequency added mass and the memory function, respectively. The convolution integral addresses the influence of the previous state on the current state. In practice, this integral is difficult to properly quantify and an alternative approach is more commonly implemented. The frequency domain representation of the radiation force isF where B(ω) is the frequency-dependent radiation damping and A(ω) is the added mass. These hydrodynamic quantities are also able to be calculated using NEMOH. The method outlined in [17] uses these two hydrodynamic quantities to construct a transfer function to relate velocity to the value of the integral in (4). In the 1-DOF model, only the heave velocity contributes to the radiation force in the heave direction so a basic transfer function will suffice; whereas in 3-DOF, in general, the radiation force for any DOF may be influence by other DOFs, depending on geometry.

3) Hydrostatic and PTO forces
The hydrostatic (or buoyancy) force acts only in the heave direction and is the difference between the weight of displaced water and the weight of the buoy, given by where ρ, g, m, and V are the density of water, acceleration due to gravity, buoy mass, and buoy volume, respectively. For a fully submerged buoy, this force is constant. Accordingly, a pretension force equal in magnitude is included in the PTO force, in addition to a stiffness and damping term as follows: where b is the damping coefficient and k s is the spring constant. This defines the submerged nominal position to be z = 0.
Both the damping coefficient and spring constant may be optimized for a given frequency; the optimal values have been derived in [18] and for a heaving PA are defined to be The optimal PTO stiffness condition for the 3-DOF regular wave scenarios may be found approximately by using (8) as an initial estimate, and adjusting k s until the root mean square (RMS) of the instantaneous power is at a local maximum or until the buoy breaches the water surface. The PTO damping, b, was left unchanged in the 3-DOF model by assuming that the optimisation of the stiffness was sufficient. However, in some low frequency cases for the 1-DOF scenarios, the damping was artificially increased to ensure stability and prevent surface breaching. In the 3-DOF scenarios, heave was considered the most influential DOF, and therefore the same optimisation process may be used. These conditions were found to produce conditions in which the phsae of F e,z andż were matched. The instantaneous power generated by the WEC is approximated by

4) Bistable force
A bistable force is any force that induces a two-well potential energy barrier, as depicted on Fig. 2. Such systems exhibit unique dynamic features not seen in regular linear systems. These features, explored in [19], include multiple single period steady state responses for a given excitation frequency leading to bifurcation of responses, aperiodic or chaotic behaviour, a large dependence on initial conditions and excitation amplitudes, stochastic resonances (resonating using a combination of low amplitude inputs), and the excitation of input frequency harmonics. These features can lead to broader resonance bandwidths particularly at lower frequencies, frequency up-conversion which turns low frequency oscillations into higher frequency oscillations, and performance improvement in stochastic excitation contexts. The main benefit of bistability occurs during interwell motion rather than intrawell due to the snap through property. This snap through mechanism forms a more energetic system if the excitation is sufficient to overcome the dividing potential barrier. The escape from a potential well was shown to broaden the frequency range of a generator in [20]. While each of these characteristic features may be exploited for different circumstances, wave energy as an application may be well suited given the low frequency and stochastic nature of waves.
For the purposes of this study, a magnetic dipole model was used to create a bistable system. A depiction of the WEC PTO model with bistability is given in Fig. 3.
The force between two magnetic dipoles, in this case between the stationary outer dipole and the dipole within the translator, may be derived from [21]  where m d,1 is the dipole moment of the first dipole, and B d,2 is the magnetic field of the second dipole at the location of the first dipole. Assuming magnetisation in the z direction, the dipole moments are given by where V is the volume of the magnet, and M is the magnetisation field per unit volume. The magnetic field for a dipole is given by [22] where r is the displacement between the dipoles, and µ 0 is the permeability of free space. Using (10) and (12), the vertical force between the two dipoles may be derived as where r 0 is the horizontal distance between dipoles. Similarly, the potential energy of this force can be derived as using the relationship The derivation was generalized in both [22] and [23]. This force was incorrectly derived in [6] and [20], however, the end result was still a force which exhibits the intended bistable behaviour. The potential energy of the dipole representation superimposed with the potential energy from a linear stiffness, provides the intended bistable profile.
To generalize the degree of bistability between any mechanism, the dimensionless parameter γ is proposed, which is the ratio of the average potential energy of the wave over the buoy, U wave , and the potential energy of the bistable system at the nominal position, U bi,z=0 , that is For a regular wave, the time average potential energy per unit area,Û wave , may be expressed aŝ whereÂ is the amplitude of the incoming wave. The parameter γ may be varied to define the height of the central potential peak, represented in Fig. 4a. In addition, the location of the stable regions can also be varied, in this case by the changing the horizontal distance between dipole r 0 while keeping γ constant. The impact of changing r 0 on the potential wells is given in Fig. 4b. Varying these parameters gave a large range of possible force and potential profiles which were used to build an understanding around which bistable conditions are favourable for the ocean wave energy context.
The aforementioned forces were combined into (2) and the dynamic equation of a 1-DOF submerged PA WEC was formed. Further descriptions and model specifics may be found in Section III.

B. Three DOF regular waves
Expanding the model from a single degree of freedom is not trivial as the excitation, radiation, and PTO forces are significantly altered. A simple drag force is also included to improve the fidelity of the simulation. The alterations and additions are explained in the following sections.
1) Excitation force Previously the excitation force was represented as a scalar in the z direction. However, in general, the BEM solver NEMOH can provide the excitation amplitude and phase for all 3 DOFs. In 3 DOF, the phase of the excitation force shifts as the buoy moves in the surge direction relative to the wave. The forces may be included into (1) in a similar way as in the 1-DOF model, with the inclusion of the phase offset due to surge, φ s , as where ωt and φ s are applied to each component of φ.
The phase offset due to surge is where k is the wavenumber, which is the solution to the dispersion relationship where the depth of the water is represented by H.

2) Radiation force
The radiation force for the 1-DOF model was a single transfer function between velocity and the z component of the convolution integral in (4). This may be expanded to accommodate 3-DOF by constructing a transfer function between all velocity components and all convolution integral components then representing this as a state space model. This method accounts for the interaction between DOFs, for example how the surge velocity impacts the radiation torque experience in the pitch direction. The contribution of the convolution integral, µ, may be described aṡ where p is a state vector of non-physical variables, and A ss , B ss , and C ss are state space matrices constructed from the aforementioned transfer functions.

3) Drag force
The viscous drag force is a nonlinear force often neglected in simple models given by where C D , A D , andẋ r are the coefficients of drag, characteristic area, and velocity of the buoy relative to the surrounding fluid. As an approximation, the fluid velocity at the geometric center of the buoy was used. By assuming a linear wave, the fluid velocity in the heave and surge directions may be estimated aṡ where d s is the submergence depth of the geometric center of the buoy.

4) Hydrostatic and PTO forces
The 1-DOF PTO force acts only in the z direction, whereas in 3-DOF the force acts in the direction of the tether. Therefore, (7) becomes where the extension of the tether, ∆l, is used instead of the heave coordinate z. The coordinate transform between (x, z, θ) → (∆l, α, φ) is represented by T with coordinates defined on Fig. 5. The inverse transformation is where a and l are the distance between buoy center and tether attachment point, and the total tether length, respectively.
It should be noted that the proposed bistable mechanism is contained within the PTO. Therefore the bistable force undergoes the same transformation T to convert to (x, z, θ) coordinates. Combining these new descriptions of forces into (1) gives the model of a 3-DOF submerged PA WEC.

III. SIMULATION CONSIDERATIONS
Within each scenario, an identical WEC geometry and setup was developed. In both the 1-DOF and 3-DOF simulations, the WEC system was excited by a series of monochromatic waves and the resulting dynamic behaviour was recorded for each frequency. Within the scenarios, to show the sensitivity to tuning, the system was simulated with the PTO settings tuned to a single frequency found by (8), and then optimal conditions at each frequency. The bistable parameter γ was also varied and the resulting root mean square (RMS) motion, time-averaged power P avg , and maximum steady state forces were recorded. The timeaveraged power is given by where t i , t f , ∆t j , and T are the initial time, final time, the j th time step, and the total time interval, respectively. These are chosen to discard any transient effects. In addition to the time-averaged power metric, which is the power absorbed by the PTO, the power absorbed by the buoy from the wave can be quantified by the time-averaged excitation power, P e , given by The phase between the excitation force and the tether velocity was also explored to further understand the impact of the bistable mechanism. The model dimensions and scenario conditions are outlined below.

A. Model dimensions
The geometry of the buoy was selected to be similar to the CETO design presented by [24]. The diagram in Fig. 6, details various relevant physical quantities. These physical quantities as well as some general quantities used in the simulations are listed below in Table I, and the drag properties based on [25] for a buoy of this shape is given in Table II.

B. Scenario conditions
The parameters pertinent to the 1-DOF and 3-DOF scenarios are detailed in Table III

IV. RESULTS
The results for the simulation scenarios described in section III are displayed below.

A. Single DOF
The 1-DOF scenarios were simulated with the PTO configuration given by (8) at each frequency and the resulting RMS of the motion as well as the timeaveraged power for varying γ are given in Figs 7a and 7b, respectively. The same results for the 1-DOF buoy simulated with PTO settings optimized for a frequency of ω = 0.77 rad/s are given in Figs 8a and 8b. This frequency was chosen to avoid peaks whilst having a reasonable response.
Given the significant difference between the optimal and non-optimal scenarios, the time domain relationship between the excitation force and the velocity for the non-optimal scenario at a frequency which showed benefit (ω = 1 rad/s) with varying γ is given in Fig. 9. The phase difference between the excitation force and velocity, as well as the excitation power for this case is presented in Fig. 10 and a phase plot relative to the potential energy profiles of these scenarios are given in Fig. 11.

B. Three DOF regular waves
The 3-DOF scenarios were also simulated with locally optimized and then suboptimal PTO conditions and the resulting motion RMS and time-averaged power for the optimal condition with varying γ are shown in Figs 12a and 12b, respectively. The same variables for the suboptimal conditions tuned to frequency ω = 0.77 rad/s is given in Figs 13a and 13b, respectively. These results will be analysed and discussed in Section V.

V. DISCUSSION
The results for the 1-DOF simulation in the nonoptimized condition show a small improvement in average power due to bistability. The motion seems to indicate that there are bands of frequencies in which there is large motion reductions. In the optimal PTO condition cases, bistability seems to greatly reduce the effectiveness of the PTO system at absorbing power. It is well known that the optimal case is when the excitation force and the velocity are in phase and any change to this caused by the bistability is detrimental. However, for the non-optimized case, bistability seems to improve the tendency of the velocity and the excitation (a) (b) Fig. 7. The resulting peak motion (7a) and the time-averaged power (7b) for the 1-DOF scenario and varying γ, using PTO parameters found from (8). The damping of the PTO was increased for frequencies below 0.54 rad/s to prevent breaching; therefore, these results may be considered only near optimal (region shaded).  force to be in phase which helps in harvesting more power. The relationship between γ, velocity-excitation lag, and excitation power shows there is a strong trend indicating that the most beneficial scenario is one in which the lag time is minimized. Even when restricted to one side of the potential well, seen in Fig. 11 at γ = 12, the force due to the nonlinear stiffness seems to passively adjust the motion leading to near optimal results. The phase plots shown in Fig. 11 show the transition between an interwell motion and intrawell motion. The corresponding excitation power demonstrates that intrawell motion leads to an improvement. One explanation of this may be that the asymmetry in the potential well, when restricted to one side, allows the device to be partially tuned to a range of frequencies. This mechanism would lead to a benefit for higher frequencies at high levels of bistability due to the higher stiffnesses in those regions (seen in Fig.8b). Similarly, bistability would benefit lower frequencies at low levels of bistability due to regions of lower stiffnesses (seen in Fig. 13b). Such a property may be useful to passively partially tune a device making a system more robust. Another property of bistability seen in the γ = 11 response in Fig. 9 is the excitation of multiple frequencies. In this case, the multiple frequencies cause significant deviation from a simple sinusoid and are detrimental for energy generation. For the optimal case, in which the excitation force and velocity are already in phase, the resulting multiple harmonic motion invariably lead to a suboptimal scenario in which only part of the motion is matched. For a perfectly tuned 1-DOF submerged WEC, bistability is detrimental, otherwise, bistability seems to provide a benefit.
In practice, it is difficult to find the exact PTO tuned conditions for more complex scenarios such as 3-DOF or for irregular waves. The 3-DOF simulations were all tuned using the aforementioned local optimization method. Accordingly, the inclusion of bistability also showed a decrease in average power, similar to the 1-DOF scenario. The resulting motion RMS exhibited a large reduction in heave motion and an amplification of the surge and pitch motions at some frequencies. This is likely due to the phase altering property and the geometric nonlinear coupling between DOFs. However, overall this did not improve the time-averaged power. This finding further supports the assumption that heave is the DOF most closely associated with power production for a single tether PA WEC. For a monochromatic wave, it is feasible to run multiple simulations to determine the optimal PTO stiffness as is done in this study. In irregular waves representing real ocean waves, optimal conditions are challenging to predict and require foreknowledge of the excitation force. Therefore, exploring the characteristics of bistability in sub-optimally tuned scenarios is useful to gain insight into the performance impact to real conditions. When the 3-DOF model uses optimal conditions for frequency ω = 0.77 rad/s, a notable improvement is seen at lower frequencies, even with viscous drag forces. In some cases, over double the power is absorbed by the PTO. The highest peaks occur when γ ≤ 2 indicating that the potential energy of the bistable mechanism relative to average potential energy of the wave is a contributing factor. However, for all these simulations the magnetic dipole spacing was held constant at r 0 = 1.5 m. To further extend this work, the influence of the magnetic dipole spacing should be investigated to characterize favourable nonlinear potential profiles for submerged PA. Bistable mechanics within the context of irregular waves should also be considered as the passive phase matching property seen in the 1-DOF sub-optimal condition could be effective in irregular conditions. Furthermore, the impact of nonlinear hydrodynamics on bistable mechanics should quantified either by including a nonlinear representation of Froude-Krylov forces, or by simulating using computational fluid dynamics. In principle, if an optimal nonlinear stiffness potential profile is identified, a passive system can be constructed by springs or magnets to realize the benefit. One potential advantage of the nonlinear magnetic spring over conventional hydraulic or electrodynamic PTOs is that the reactive energy of such magnetic mechanisms is highly efficient, allowing the natural frequency to be altered with almost no parasitic losses.

VI. CONCLUSION
The impact of a bistable force acting on a submerged PA WEC was explored by simulating a CETO-like device in 1-DOF and 3-DOF regular wave scenarios, where the 3-DOF simulations included geometric nonlinearity and drag effects. The bistable force was represented as a magnetic dipole and combined with a linear stiffness to give an overall nonlinear stiffness in the system. The results showed that for a device that is not optimally tuned, bistability can provide a substantial benefit by exciting multiple frequencies and passively matching the phase of the excitation force and velocity. For optimally tuned situations, the excitation force and velocity are already in phase and the addition of nonlinear stiffness reduces the effectiveness of the device. In practice tuning a device for a given wave condition is challenging, but bistability was shown in  where the PTO is tuned for ω = 0.77 rad/s. regular waves to provide a passive benefit. This work builds upon existing work in bistability by applying such a mechanism to the submerged WEC context and develops an understanding around what properties of bistable mechanics is providing the enhanced power performance.