Abstract
Efficient design of wave energy converters based on floating body motion heavily depends on the capacity of the designer to accurately predict the device’s dynamics, which ultimately leads to the power extraction. We present a (quasi-nonlinear) time-domain hydromechanical dynamic model to simulate a particular type of pitch-resonant WEC which uses gyroscopes for power extraction. The dynamic model consists of a time-domain three-dimensional Rankine panel method coupled, during time integration, with a MATLAB algorithm that solves for the equations of the gyroscope and Power Take-Off (PTO). The former acts as a force block, calculating the forces due to the waves on the hull, which is then sent to the latter through TCP/IP, which couples the external dynamics and performs the time integration using a 4th-order Runge-Kutta method. The panel method, accounting for the gyroscope and PTO dynamics, is then used for the calculation of the optimal flywheel spin, PTO damping, and average power extracted, completing the basic design cycle of the WEC. The proposed numerical method framework is capable of considering virtually any type of nonlinear force (e.g., nonlinear wave loads) and it is applied and verified in the paper against the traditional frequency domain linear model. It proved to be a versatile tool to verify performance in resonant conditions.
1. Introduction
We will first introduce the standard mathematical approach followed to evaluate performance of wave energy converters (WECs), whose design requires early quantification of the new device’s power extraction capabilities. Concepts which rely on a floating body’s motion in waves for power extraction have taken advantage of methods created to simulate motion of ships and offshore structures. The very first objective is to tune the converter’s natural frequency (of its translation and rotation considered motions) with the predominant wave period from the operating site spectrum. Traditionally, the waves and equations of motion are linearized (i.e., small wave amplitude and body motion compared to the device’s nominal size), potential flow is assumed with Laplace as the governing equation (i.e., ), and the resulting linear Boundary Value Problem (BVP) to be solved for is summarized in Table 1.
The most common approach to study the body’s motion is to assume harmonic inputs (i.e., force due to an incident regular wave ), for which the linear equation of motion, with frequency , becomes
The linearization of the motion allows a decomposition of the total fluid potential as a summation of the classical radiation and diffraction wave potentials [1]. The former models the waves generated by the moving body (with the same frequency as the incident wave) in calm water. It can be represented by the summation of each degree of freedom (DoF) displacement, subject to its own body BC [2].
The radiation potential iswhere is the normal vector pointing out of the body and into the fluid for and for . The added-mass and damping come from the real and imaginary components of the radiation forces [2]. The exciting forces may be calculated straight from Haskind relations [3]. All of these integral equations come from Green’s second identity [4].
To bridge frequency and time domain, Cummins developed general equations of motion in time domain [5]. These equations are based on impulsive motion, decomposing the flow around the body into two components, one due to an impulsive displacement on the body and the other one representing the decaying wave motion generated by such displacement. For a body translating on the water surface, the equation of motion becomes
The kernel of the convolution integral is related to the frequency dependent quantities through Fourier Transform [5–8]. It is the radiation impulse response function and tells us how the history of radiation forces influences the body motion for some time (i.e., fluid memory effects).
Finally, to solve the BVP stated in Table 1 and find the potential on the body, the second Green identity is used. WAMIT, an industry standard frequency domain Boundary Element Method (BEM) code, presents the equation in the form below, solving it for discretized panels of the body wetted surface with center at [9].
The solid angle can be calculated by noting that a uniform potential applied over a closed body produces no flux (i.e., ) [10]:
The Green function used by WAMIT is the wave source potential, where is the zero-order Bessel function and and are the position vectors of the point considered and the source, respectively [9].
With the computation of the hydrodynamic coefficients of (1) from the solution of the BVP through (7), the transfer function, or Response Amplitude Operator (RAO), can be calculated and the body’s response to a given sea state calculated through a simple multiplication [11]. This is what justifies the predominance of frequency versus time-domain approaches.
The frequency domain approach has seen consistent use for WEC dynamics analysis. More recently, Rhinefrank et al. have applied such methodology to evaluate the motion of their respective devices [12–14]. The first, in particular, improved the model of their cylinder-shaped WEC by introducing viscous effects in the form of Morrison’s equation.
Equation (5) is widely used, together with frequency domain BEM, to describe the motion of WECs in time domain. This is a very powerful approach for most situations, yielding accurate results in many cases with small computational time. Works such as the two-body heave converter from Andres et al., the SEAREV from Babarit et al., and heaving WEC are all examples of recent successful use of linear time-domain approach for WEC design [15–19]. However, the description of the problem through Fourier Transform makes it impossible to incorporate nonlinear (potential flow wave-dependent) effects into these models.
In this work, we approach the design of a gyroscope powered WEC using Kring’s time-domain Rankine panel method AEGIR, which solves for both the body and free-surface motions directly [8]. To model this, the state-space vector and external forces calculated by AEGIR at every time step will have to be corrected to account for the spinning flywheel dynamics. The change in angular momentum of the gyroscope will induce torques in all 3 directions, coupling all the rotation degrees of freedom (DoF). The corrections and time marching were implemented in a MATLAB code, which constantly shares information with AEGIR through TCP/IP protocol. This numerical model is the initial (linear) framework for future inclusion of effects such as fully nonlinear boundary conditions, pressure integration up to the deformed free-surface [17], and retention of high-order terms in the Bernoulli equation for pressure calculation [20].
2. The IOwec Case Study
The WEC analyzed in this paper is the Inertial Ocean Wave Energy Converter (IOwec), which originated from a collaboration project between the University of Torino, Wave4Energy, and iShip lab at MIT (now moved to Virginia Tech). This concept was developed and first presented as a contender in the Department of Energy Wave Energy Prize (WEP) [21]. It is an evolution of the ISWEC, another gyroscope driven WEC designed in Italy [22], by modifying the simple hull shape with the intent of increasing the resonance period of the ISWEC up to the predominant wave period of 8 s given by the WEP [23]. Figures 1(a) and 1(b) shows the hull change, from the ISWEC to the IOwec shape.
(a)
(b)
Inside the IOwec, due to the precession motion of the spinning masses, the gyroscope roll will induce torques in both pitch and yaw directions. The former activates the angular motion of the PTO shaft and it is used to directly produce electric power from a variable frequency alternator. The latter, however, is extremely undesirable, as it excites the yaw motion of the hull. In this case an offset between the prevalent sea direction and the hull’s longitudinal symmetry plane will always exist, and part of the wave energy will be transmitted to undesirable modes in the horizontal diametral plane (i.e., sway, roll, and yaw).
To eliminate the yaw torque, a pair of counterrotating gyroscopes is used. Both flywheels will generate the same pitch torque, with equal magnitude and direction, but for yaw the torque will be opposite, negating the -axis induced rotation completely (Figure 2). In fact, the counterrotating gyroscope pair is also desirable to create redundancy in the power extraction system. This means the device will continue to extract a reduced amount of energy, even if one of the gyroscopes fails.
The gyroscopes were sized to have the largest diameter possible, while leaving sufficient clearance on the sides, bottom, and top to allow for a complete 360° turn. Table 2 lists IOwec’s principal characteristics. The VCG is measured from the free surface, with positive values above the waterline, while the rest of the vertical quantities are measured from the keel.
2.1. Floating Body’s Mathematical Model
Consider an unrestricted, stationary, floating body, free to move in its six, rigid-body degrees of freedom about a noninertial frame of reference fixed to its equilibrium position in calm water (Figure 3) [8]. We assume small body motions with respect to the equilibrium position and small wave disturbance with respect to the body’s nominal size. The problem is linearized, allowing the application of potential flow together with impulse theory and the decomposition of the forcing into impulsive (local) and wave (memory) forces, as proposed by Cummins [5].where is the velocity impulse function and is the excitation force in the th mode. The hydrostatic restoring coefficients are easily calculated following classical naval architecture theory [24]. Figure 4 shows the body-fixed reference frame, where (12) is applicable, as well as the gyroscope’s local frame .
The fully nonlinear gyroscopes equations of motion are [25]
It is important to notice how, in our present model, we chose a linear PTO and disregarded possible control dynamics. Nonetheless, such extra complications of a particular system may be easily included in our framework through careful consideration of the dynamics. Recent works have considered, for example, mooring forces, nonlinear PTOs, and control techniques such as latching [26–31].
The mathematical formulations regarding the hull are written in a state-space format in AEGIR, so the MATLAB code which implements the gyroscope and PTO mechanics must also follow the same convention. Going down the same path as outlined by Kring, the six-DoF, second-order equations of motion can be modified into twelve of first degree, with the state vector written as a combination of displacements and velocities [8].
Gyroscopic terms , proportional to the hull’s acceleration, will go alongside mass terms, just like the added-mass. The forcing vector becomeswhere is the total number of gyroscopes inside the hull. Both acceleration and velocity proportional gyroscope torque matrices, coming from (14) and (15), are nonzero only for pitch and yaw.where, from (13), (14), and (15),
Assuming (13) can be linearized by getting rid of all higher-order terms in and while, at the same time, assuming small roll angles for the gyroscope, we may apply the Fourier Transform and move into frequency domain to express
Introducing the PTO natural frequency as and substituting it in (20) [22],
It is clear now that we must tune to match the relevant incident wave frequency. The ideal PTO spring constant becomes
The hull is designed to have a pitch resonance period = 8 s in order to match the predominant sea given by the WEP committee. Equation (22) can then be used to size our PTO spring. At this stage, we will consider only the hull resonance frequency as reference for the PTO tuning, since the motion response is greatly magnified around that frequency. The reactive control constant becomes
The power is extracted by the PTO through the linear damping coefficient . The power metric utilized for this work is the average power extracted over one wave period, which can be written as [22, 25]
Substituting (21) into (24) and considering an optimum reactive control, the power extracted becomes
We now have derived a direct way to estimate what would be the power prediction of the device by just knowing the hull and gyroscope motions. Also, the optimum reactive control can be designed for every incoming wave frequency, guaranteeing resonance. It is important to note that the damping would always be optimized for such condition.
Finally, we must choose a first guess for the nominal spin rate of the gyroscopes. Scaling the ISWEC’s spin of 2150 RPM from its 0.56 m length we find
Our numerical model later showed us this value is too high, taking us too long to converge to the optimal quantity. We then define our first guess as half of the ISWEC’s scaled value.
3. The Time-Domain Boundary Value Problem
Here, we also consider the fluid flow incompressible, irrotational, and inviscid so the flow can be represented by a velocity potential , which must satisfy Laplace’s equation (i.e., ). The flow pressure can also be calculated in the frame Bernoulli’s equation [2]. The fully nonlinear kinematic and dynamic free-surface, body boundary, and radiation conditions are, respectively [8],
The total disturbance potential can be written as a superposition of two flows, the local (i.e., instantaneous fluid response due to the impulsive body motion) and the memory (i.e., wave) [5, 8].
To arrive at the linear boundary conditions we apply a Taylor expansion about the mean body position for the body condition and about for the free surface, keeping only its linear terms [8].(1)Kinematic free-surface condition is(2)Dynamic free-surface condition is(3)Body boundary condition is
Equations (33) and (34) will be our evolution equations for the free-surface deformation and wave potential . Before, however, we must calculate the unknown and the derivatives of , as we are only given the value of the wave potential at and the derivative of the local potential at .
The displacement of any point within the rigid body can be described by combining translation and rotation [8].
Substituting (36) into (35) yields the body boundary condition in Ogilvie and Tuck’s notation, minus the -terms. These omitted terms would be fundamental if our body had a translation velocity [8, 32].
Particularly for the pressure calculation, the linearized Bernoulli equation can be decomposed into local , memory , and hydrostatic components. The total pressure is [8].
Finally, the force acting on the body is the integration of the pressure on .
3.1. The Local Flow Contribution
The local flow can be decomposed into terms proportional to the acceleration, velocity, and displacement of the body. Only the ones due to acceleration are nonzero for a stationary body (i.e., the basis flow is absent), yielding the term written in [8]where
3.2. The Memory Flow Contribution
If an incident wave is imposed in the domain, it is represented as an extra component of the memory flow ( and ). The body boundary condition becomes [8]
The free-surface conditions outlined in (33) and (34) have the unknown velocity potentials. They will be solved by AEGIR at every time step, calculating the free-surface deformation [8].
3.3. The Boundary Integral Formulation
As stated before, we are only given the value of the wave potential at and the derivative of the local potential at . To calculate for the derivatives of , Green’s second identity is applied to the Boundary Value Problem, yielding a very similar equation to (7) [8]:
Finally, AEGIR is called a Rankine panel method due to its choice on the Rankine source potential as the known function in the integral:
4. Numerical Model
4.1. Spatial Discretization
AEGIR is a high-order panel method, meaning the surfaces are broken down into the desired number of discrete quadrilateral facets and represented mathematically as nonuniform basic spline surfaces (NURBS). The variables of interest becomewith
The basis function within a panel is a second-order B-spline, meaning first and second derivatives can be obtained analytically [8]:where is the panel width. Applying such basis function on (33), (34), and (43) yields the discrete formulation of the kinematic and dynamic boundary conditions and boundary integral formulation, respectively. The first two are time dependent, while the last is a linear system of equations whose size depends on the number of panels in the problem [8].where and are the dipole and source influence coefficients, respectively:
4.2. Temporal Discretization
4.2.1. Free-Surface Boundary Conditions
The free-surface temporal discretization uses an “implicit” method, which is a mix of explicit and implicit Euler schemes [8]. The following discretization is used to integrate equations (48) in time:where is the variable value at the current (th) time step.
4.2.2. Body Motion
The time marching scheme for the body motion is the Runge-Kutta 4th order. It consists of four steps, the first is a forward Euler of half step, and the second is an implicit Euler corrector, still with half step, but now using the previous prediction as current state vector. The third is a full-step midpoint rule, using the predicted value. Finally, the last consists of a Simpson’s rule which uses all previous steps to build the true state vector [33]. Using a multipoint method allows the use of larger time-steps while still maintaining numerical stability, while, in particular for the Runge-Kutta, not needing extra points for calculation.(1)Forward Euler predictor is(2)Implicit Euler corrector is(3)Midpoint rule predictor is(4)Simpson’s rule corrector is
Figure 5 better illustrates the information exchange between AEGIR and the MATLAB code. The latter calculates the true state vector from AEGIR’s bare hull estimation.
4.2.3. Radiation Condition
Finally, in a discrete scenario, our domain is bounded. This means the radiation condition (see (31)) must be imposed on the edges of the domain. This is achieved in AEGIR by placing “numerical beaches” which artificially damp the outgoing waves far from the body [8, 34].where is an artificial viscosity imposed on the waves for their decay near the edges.
5. Results
5.1. Sensibility Analysis
All the analyses carried out are done with monochromatic waves in head seas, with periods ranging from 5 s to 12 s. The numerical simulation has to prove its stability and convergence within the wave range. Considering infinite depth (i.e., linear dispersion relation , where and are the wave frequency and number, resp.), the wave lengths range from 39 m to 225 m, the higher one driving the domain size. Time step and mesh density must then be fine enough to perceive the one with the lowest period (i.e., = 5 s).
5.1.1. Domain Size
Figure 6 shows the final geometry utilized for the sensibility analyses.
The domain initial sizing was based on the model length and the longest wave to be simulated, which has a period of 12 s and length = 225 m. Table 3 presents the values used.
To be certain the domain was sized correctly, the longest wave ( = 12 s) with amplitude = 1 m was simulated for three domain sizes of , , and times the dimension of the initial one. The pitch motion is plotted against the simulation time in Figure 7. For this run a 60 s of ramp time and 0.1 s time step was used, which should be more than enough for the considered wave. Also, a simulation time of 150 s was enough to achieve steady state motion. No noticeable difference was observed between the first two multipliers, while the last diverged in less than . The domain size corresponding to the multiplier was selected.
5.1.2. Mesh Density
Three mesh density factors of , , and , corresponding to panel sizes of 3.75 m and 2.50 m and 1.25 m, respectively, were tested against the smallest wave (39 m). For this run, the same ramp time and time-step size presented in Section 5.1.1 were used. Since no noticeable difference was observed between the two finest meshes, and the coarsest ones diverged about 2.36% (Figure 8), the 2.5 m panel width was chosen for subsequent runs.
5.1.3. Time Step
The simulation time step has to be chosen to accurately solve the smallest wave period. To select a proper step, four values were tested while simulating the = 5 s wave: 0.05 s, 0.1 s, 0.2 s, and 0.5 s. Figure 9 shows the pitch motion versus simulation time for each step, with the exception of 0.5 s which diverged. The largest one (i.e., 0.2 s) had a deviation of only 0.35% from the smallest one. Therefore, a good compromise between accuracy and speed for our simulations is achieved with = 0.2 s, which was chosen for all subsequent runs.
Finally, Table 4 summarizes the domain size, mesh density, and time step chosen based on the aforementioned sensibility analyses. Figure 10 shows the free surface and body surfaces after meshing.
5.2. Bare Hull Motion
Before assessing the gyroscope dynamics and power extraction we will take a look into the hull’s bare motion. Figures 11 and 12 show the heave and pitch RAO, respectively. The former has an excellent agreement with WAMIT results, while the latter is just slightly off at the resonant peak.
WAMIT solves the radiation potential and, through (4), calculate the added-mass, damping, and exciting forces [9]. Figures 13 and 14 show the comparison between AEGIR and WAMIT exciting force and moment, which are in good agreement, with minimal discrepancy between both methods.
We now oscillate the hull under a prescribed motion of amplitude and frequency on an undisturbed free surface. It will generate outgoing waves whose pressure, when integrated on the hull, describes a force out of phase with the motion. Breaking the force into a component in phase with the acceleration and another in phase with the velocity yields the added-mass and damping components, respectively. The diagonal coefficients are
Figures 15(a)–15(d) show both AEGIR and WAMIT added-mass and damping coefficients for heave and pitch. We can see an excellent agreement for the former, while pitch presents higher damping at resonance and a slight offset on the added-mass for starting at . This means the memory potential in AEGIR, which accounts for the radiation problem, yields a slightly different solution than WAMIT. The extra damping may explain why the pitch RAO peak is a bit smaller for AEGIR.
(a)
(b)
(c)
(d)
5.3. Gyroscope Spin Sensibility
The base spin rate was defined by (27). We would now like to determine the optimal PTO damping. Equation (25), which only uses hull pitch and gyroscope roll amplitudes, will be utilized to estimate the averaged power extraction over one wave period. Figure 16 shows such power by different multipliers of the base spin rate.
We can notice that, starting from low spin rates, the power extracted increases until an optimum value. After that the peak extraction moves to longer waves, until it is not noticeable within the selected wave bandwidth. Increasing the flywheel spin makes the system stiffer, moving the hull’s natural frequency towards longer waves, as can be seen in Figure 17.
The best spin found, considering our 8 s resonant period, is (i.e., 150 RPM). This rate will be used for the subsequent runs.
5.4. PTO Damping Sensibility
Maintaining the same PTO spring, optimized for the hull natural frequency, and with the chosen spin rate, we can do a sensibility analysis with the PTO damping. Equation (24) can be utilized, using the power predicted during the spin analysis, to estimate the required damping.
The actual instantaneous power extracted can now be integrated to calculate the average power:
We can observe, from Figure 18, a shift of the peak energy extracted towards smaller waves as the damping is increased, due to the decrease of gyroscopic roll. An optimum damping of is attained, after which the device quickly loses energy extraction efficiency, especially for longer waves.
The effects of damping are evident on the hull’s pitch as well. The optimum damping guarantees that its natural frequency remains at 8 s, while minimizing the motion as much as possible. Figure 19 shows how the pitch RAO changes for all the damping variation.
Reporting power extracted by itself is not a good representation of IOwec’s capabilities, especially since it varies with the wave amplitude. A better way to quantify efficiency is to plot the capture width, which is the ratio between the average powers extracted and incoming from a 2D section of the wave profile.
Figure 20 shows the IOwec is able to, at resonance, extract a monochromatic wave 25 m wide, 1.25 times its beam.
In all cases, the yaw torques induced by both counterrotating flywheels are seen to cancel one another. Figure 21 shows the torques from both gyroscopes for the = 8 s wave, with frequency exactly twice that of the wave.
6. Conclusions
We presented the mathematical formulation of the hydromechanics time-domain simulation model for a gyroscopic, pitch-resonant, floating WEC. First the initial design of the WEC to be simulated was outlined, the hull was sized to guarantee resonance for the most energetic wave period given by the WEP committee, while the gyroscopes’ diameter was maximized while guaranteeing clearance. A maximization of the flywheel size allowed the use of a smaller nominal spin of 120 RPM, which is easier to maintain.
Next, the Boundary Value Problem was outlined under a potential flow assumption. Green’s second identity provided a boundary integral equation to be solved for the unknown derivatives of on the free surface, which, after being discretized into quadratic panels, allowed for the calculation of the potential on the surface of the hull. The free-surface boundary conditions and body equation of motion can then be integrated in time, giving values for , and the body’s state vector for the next time step.
The time-domain numerical method, after validation against a state-of-the-art frequency domain panel method, was fitted with the gyroscope and PTO dynamics, allowing for a quick evaluation of the optimal flywheel spin, PTO damping, and average power extracted. In fact, the gyroscope mathematical model presented in Section 2.1 could be reformulated for any different type of external mechanics (e.g., a pendulum), making this a very versatile framework for designers and engineers.
One could also envision the extension of the present method to account for nonlinear effects (while also abandoning the state-space format), particularly important around resonance. The higher-order terms of the free-surface boundary conditions could be retained for a more accurate representation of steeper wave profiles. Integration of pressures on the hull could also be done up to the deformed free-surface, while also retaining the nonlinear Bernoulli terms. An even better extension would be a Mixed Eulerian-Lagrangian approach, in which collocation points are translated at each time step to the regions of higher gradient of the free-surface, guaranteeing better resolution at highly nonlinear areas [35–38].
Next steps of the IOwec design must include thorough and consistent experimental testing of a scale model of the WEC, including flywheel and PTO systems. This will provide much needed data for further validation of the present linear numerical model and future implementation of nonlinear mechanics.
Symbols
| : | Added-mass matrix |
| : | Added-mass matrix due to impulsive motion of the body |
| : | Wave amplitude |
| : | Hull’s maximum beam |
| : | Damping matrix |
| : | Damping matrix due to body impulsive motion |
| : | Restoring forces matrix |
| : | Power Take-Off linear damping constant |
| : | Hull’s depth |
| : | Gyroscope diameter |
| : | Wave force signal exciting the body motion |
| : | Gravitational acceleration |
| : | Water depth |
| : | Wave height |
| : | Gyroscope’s flywheel roll inertia |
| : | Gyroscope’s flywheel pitch inertia |
| : | Gyroscope’s flywheel spin inertia |
| : | Gyroscope frame roll inertia |
| : | Gyroscope frame pitch inertia |
| : | Gyroscope frame yaw inertia |
| : | Wave number |
| : | Power Take-Off spring constant |
| : | Length Over All, the hull’s maximum length |
| : | Body mass matrix |
| : | Gyroscope’s flywheel mass |
| : | IOwec’s roll radius of gyration |
| : | IOwec’s pitch radius of gyration |
| : | IOwec’s yaw radius of gyration |
| : | Submerged body wetted surface |
| : | Wave period |
| : | Hull’s design draft (submergence) |
| VCG: | Vertical position of the center of gravity, measured from the undisturbed free surface, with positive values above the waterline |
| : | Normal velocity on the surface of the body. The normal vector points into the fluid |
| : | Wave excitation force |
| : | Device’s displaced mass at the design draft |
| : | Wave elevation |
| : | Wave radian frequency |
| : | Force signal phase angle |
| : | Motion signal phase angle |
| : | Total disturbance velocity potential of the Boundary Value Problem |
| : | Velocity potential of the fluid |
| : | Diffracted wave potential |
| : | Flywheel spin rotation |
| : | Incident wave potential |
| : | Local flow potential. One of the three components of the total disturbance potential |
| : | Radiated wave potential |
| : | Scattered wave potential |
| : | Wave flow potential. One of the three components of the total disturbance potential |
| : | Fluid density |
| : | Gyroscope roll rotation |
| : | Body displacement for the th degree of freedom due to wave excitation |
| : | Body velocity for the th degree of freedom due to wave excitation |
| : | Body acceleration for the th degree of freedom due to wave excitation. |
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors would like to thank Navatek, especially Dr. David Kring, for generously offering AEGIR, as well as instructions on how to couple external dynamics to it. The present work was made possible through the financial support of the Brazilian Council for Scientific and Technological Development (CNPq).