# Mathematical Framework for Hydromechanical Time-Domain Simulation of Wave Energy Converters

Authors: J. Seixas de Medeiros; S. Brizzolara
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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.,∇2ϕ=0), and the resulting linear Boundary Value Problem (BVP) to be solved for is summarized in Table 1.Table 1
Boundary Value Problem (BVP) for a floating body in finite and infinite depth.
Boundary condition
Finite depth (-h<z<0)
Infinite depth (h→∞)
Bottom
∂
ϕ
∂
z
=
0
at
z
=
-
h
∇
ϕ
→
0
as
z
→
-
∞
Free-surface kinematic
∂
ϕ
∂
z
=
∂
η
∂
t
=
0
at
z
=
0
Free-surface dynamic
∂
ϕ
∂
t
+
g
η
=
0
at
z
=
0
Floating body
∂
ϕ
∂
n
=
v
n
at
S
b
Radiation condition
ϕ
,
∇
ϕ
→
0
as
|
x
,
y
|
→
∞The most common approach to study the body’s motion is to assume harmonic inputs (i.e., forceXicos⁡(k→x→-ωt+ϵf) due to an incident regular wave Acos⁡(k→x→-ωt)), for which the linear equation of motion, with frequency ω, becomes(1)∑j=16ξj-ω2Mij+aijω+iωbijω+cij=AXi.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].(2)ϕ=ϕR+ϕD=ϕR+ϕI+ϕS,where∂ϕS∂n=0.The radiation potential is(3)ϕR=iω∑j=16ξjϕj,∂ϕj∂n=nj,where nj is the normal vector n→ pointing out of the body and into the fluid for j=1,2,3 and x→×n→ for j=4,5,6. 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].(4)aij-iωbij=ρ∬Sbϕj∂ϕi∂ndS,Xi=-iωρ∬Sb∂ϕi∂nϕI-ϕi∂ϕI∂ndS.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(5)∑i=16Mij+aijoξ¨j+bijoξ˙j+cijξj+∫-∞tKijt-τξ˙jτdτ=fkt.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).(6)aijω=aijo-1ω∫0∞Ktsin⁡ωt,aijo=limω→∞⁡aijω,bijω=∫0∞Ktcos⁡ωtdω,Kijτ=2π∫0∞bijωcos⁡ωtdω.Finally, to solve the BVP stated in Table1 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].(7)αx→ϕx→+∬Sbϕξ′→∂Gξ′→;x→∂nξ′→dξ′→=∬Sb∂ϕξ′→∂nξ′→Gξ′→;x→dξ′→,(8)ξ′→=ξi^+ηj^+ζk^,(9)x→=xi^+yj^+zk^.The solid angle can be calculated by noting that a uniform potential applied over a closed body produces no flux (i.e.,ϕn=0) [10]:(10)αx→=-∬∂V∂Gξ′→;x→∂nξ′→dξ′→.The Green functionG(ξ′→;x→) used by WAMIT is the wave source potential, where Jo(x) is the zero-order Bessel function and x→ and ξ′→ are the position vectors of the point considered and the source, respectively [9].(11)Gx→;ξ′→=1r+1r′′+2∫0∞dkk+Kcosh⁡kz+hcosh⁡kζ+hksinh⁡kh-Kcosh⁡khe-khJokR,r2=x-ξ2+y-η2+z-ζ2,r′′2=x-ξ2+y-η2+z+ζ+h2.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 S(ω) 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.Figure 1
(a) ISWEC hull, with its taper from bottom to top. (b) IOwec hull, with additional tapper at the stern and bow.
(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 thez-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.Figure 2
Pair of gyroscopes spinning counter to each other. The yaw torque is canceled, while the pitch is doubled.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. Table2 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.Table 2
IOwec main dimensions.
Part
Dimension
Value
Unit
Hull
L
45
m
B
20
m
D
10
m
T
d
7
m
V
C
G
-
1.485
m
Δ
5217.87
ton
R
x
x
8.45
m
R
y
y
14.65
m
R
z
z
3.84
m
Gyroscope
d
6.8
m
m
f
140
ton
I
x
x
5.17
⋅
10
5
kg⋅m2
I
y
y
5.17
⋅
10
5
kg⋅m2
J
1.02
⋅
10
6
kg⋅m2
I
x
f
3.82
⋅
10
3
kg⋅m2
I
y
f
3.32
⋅
10
4
kg⋅m2
I
z
f
3
⋅
10
4
m


### 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 referencexyz0 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].(12)Mij+aijoξ¨jt+cijξj=Xjt-∫-∞tKt-τξ˙jτdτ,where K(t) is the velocity impulse function and Xj(t) is the excitation force in the jth mode. The hydrostatic restoring coefficients cij are easily calculated following classical naval architecture theory [24]. Figure 4 shows the body-fixed xyz1 reference frame, where (12) is applicable, as well as the gyroscope’s local frame xyz2.Figure 3
Domain boundary surfacesSF, S∞, and SB, along with the inertial frame of reference xyz0 [8].Figure 4
Hull and gyroscope Frames of Reference.θ is the rotation along the PTO axis and only unconstrained DoF of the gyroscope.The fully nonlinear gyroscopes equations of motion are [25](13)Ixθ¨=Iy-Izξ˙5ξ˙6cos2⁡θ-sin2⁡θ+ξ˙62-ξ˙52sin⁡θcos⁡θ-Jφ˙ξ˙5cos⁡θ+ξ˙6sin⁡θ-klθ-clθ˙,(14)Ty1=Iycos2⁡θ+Izsin2⁡θξ¨5+Iy-Izsin⁡θcos⁡θξ¨6+2Iz-Ix-Iyξ˙5θ˙cos⁡θsin⁡θ+sin2⁡θ-cos2⁡θJξ˙6θ˙+Ix+Iycos2⁡θξ˙6θ˙-Jφ˙θ˙cos⁡θ,(15)Tz1=Iy-Izsin⁡θcos⁡θξ¨5+Iysin2⁡θ+Izcos2⁡θξ¨6-Ix+Iy-Izξ˙5θ˙sin2⁡θ-Ix+Iy-2Izξ˙6θ˙sin⁡θcos⁡θ-Izξ˙5θ˙cos2⁡θ-Jφ˙θ˙sin⁡θ.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 vectory→ written as a combination of displacements and velocities [8].(16)dy→dt=f→twithy→t=y→1ty→2t=ξ˙jtξjt.Gyroscopic termsT→ia(t), proportional to the hull’s acceleration, will go alongside mass terms, just like the added-mass. The forcing vector becomes(17)ft→=-M+ao+∑i=1nTiat-1Cy→2t+∫0tKt-τy→1τdτ+∑i=1nTivy→1,ty→1t,where n 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.(18)Tiay→˙1,t=0⋯0⋮⋱⋮T55,iaT56,ia0⋯T65,iaT66,ia,Tivy→1,t=0⋮T5,ivT6,iv,where, from (13), (14), and (15),(19)T55,ia=Iycos2⁡θ+Izsin2⁡θ,T56,ia=Iy-Izsin⁡θcos⁡θ,T65,ia=Iy-Izsin⁡θcos⁡θ,T66,ia=Iysin2⁡θ+Izcos2⁡θ,T5v=-2Iz-Ix-Iyξ˙5θ˙cos⁡θsin⁡θ-sin2⁡θ-cos2⁡θJξ˙6θ˙-Ix+Iycos2θξ˙6θ˙+-1iJφ˙θ˙cos⁡θ,T6v=Ix+Iy-Izξ˙5θ˙sin2⁡θ+Ix+Iy-2Izξ˙6θ˙sin⁡θcos⁡θ+Izξ˙5θ˙cos2⁡θ+-1iJφ˙θ˙sin⁡θ.Assuming (13) can be linearized by getting rid of all higher-order terms in ξj 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(20)Θ=JφiωIxω2-iωcl-klΞ5.Introducing the PTO natural frequency asωPTO2=kl/J and substituting it in (20) [22],(21)Θ=JφiωIxω2-JωPTO2-iωclΞ5.It is clear now that we must tuneωPTO to match the relevant incident wave frequency. The ideal PTO spring constant becomes(22)kl=Ix·ω2.The hull is designed to have a pitch resonance periodTp = 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(23)kl=Ix·ωn2≃604,921.50kg·m2s2.The power is extracted by the PTO through the linear damping coefficientcl. The power metric utilized for this work is the average power extracted over one wave period, which can be written as [22, 25](24)PE¯ω=12clθ˙2=12clω2Θ2.Substituting (21) into (24) and considering an optimum reactive control, the power extracted becomes(25)PE¯ω=Jφω2Ξ5Θ2.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(26)φ=φISWEC0.5645≃240RPM.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.(27)φ=120RPM.


## 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 referencexyz0 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].(12)Mij+aijoξ¨jt+cijξj=Xjt-∫-∞tKt-τξ˙jτdτ,where K(t) is the velocity impulse function and Xj(t) is the excitation force in the jth mode. The hydrostatic restoring coefficients cij are easily calculated following classical naval architecture theory [24]. Figure 4 shows the body-fixed xyz1 reference frame, where (12) is applicable, as well as the gyroscope’s local frame xyz2.Figure 3
Domain boundary surfacesSF, S∞, and SB, along with the inertial frame of reference xyz0 [8].Figure 4
Hull and gyroscope Frames of Reference.θ is the rotation along the PTO axis and only unconstrained DoF of the gyroscope.The fully nonlinear gyroscopes equations of motion are [25](13)Ixθ¨=Iy-Izξ˙5ξ˙6cos2⁡θ-sin2⁡θ+ξ˙62-ξ˙52sin⁡θcos⁡θ-Jφ˙ξ˙5cos⁡θ+ξ˙6sin⁡θ-klθ-clθ˙,(14)Ty1=Iycos2⁡θ+Izsin2⁡θξ¨5+Iy-Izsin⁡θcos⁡θξ¨6+2Iz-Ix-Iyξ˙5θ˙cos⁡θsin⁡θ+sin2⁡θ-cos2⁡θJξ˙6θ˙+Ix+Iycos2⁡θξ˙6θ˙-Jφ˙θ˙cos⁡θ,(15)Tz1=Iy-Izsin⁡θcos⁡θξ¨5+Iysin2⁡θ+Izcos2⁡θξ¨6-Ix+Iy-Izξ˙5θ˙sin2⁡θ-Ix+Iy-2Izξ˙6θ˙sin⁡θcos⁡θ-Izξ˙5θ˙cos2⁡θ-Jφ˙θ˙sin⁡θ.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 vectory→ written as a combination of displacements and velocities [8].(16)dy→dt=f→twithy→t=y→1ty→2t=ξ˙jtξjt.Gyroscopic termsT→ia(t), proportional to the hull’s acceleration, will go alongside mass terms, just like the added-mass. The forcing vector becomes(17)ft→=-M+ao+∑i=1nTiat-1Cy→2t+∫0tKt-τy→1τdτ+∑i=1nTivy→1,ty→1t,where n 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.(18)Tiay→˙1,t=0⋯0⋮⋱⋮T55,iaT56,ia0⋯T65,iaT66,ia,Tivy→1,t=0⋮T5,ivT6,iv,where, from (13), (14), and (15),(19)T55,ia=Iycos2⁡θ+Izsin2⁡θ,T56,ia=Iy-Izsin⁡θcos⁡θ,T65,ia=Iy-Izsin⁡θcos⁡θ,T66,ia=Iysin2⁡θ+Izcos2⁡θ,T5v=-2Iz-Ix-Iyξ˙5θ˙cos⁡θsin⁡θ-sin2⁡θ-cos2⁡θJξ˙6θ˙-Ix+Iycos2θξ˙6θ˙+-1iJφ˙θ˙cos⁡θ,T6v=Ix+Iy-Izξ˙5θ˙sin2⁡θ+Ix+Iy-2Izξ˙6θ˙sin⁡θcos⁡θ+Izξ˙5θ˙cos2⁡θ+-1iJφ˙θ˙sin⁡θ.Assuming (13) can be linearized by getting rid of all higher-order terms in ξj 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(20)Θ=JφiωIxω2-iωcl-klΞ5.Introducing the PTO natural frequency asωPTO2=kl/J and substituting it in (20) [22],(21)Θ=JφiωIxω2-JωPTO2-iωclΞ5.It is clear now that we must tuneωPTO to match the relevant incident wave frequency. The ideal PTO spring constant becomes(22)kl=Ix·ω2.The hull is designed to have a pitch resonance periodTp = 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(23)kl=Ix·ωn2≃604,921.50kg·m2s2.The power is extracted by the PTO through the linear damping coefficientcl. The power metric utilized for this work is the average power extracted over one wave period, which can be written as [22, 25](24)PE¯ω=12clθ˙2=12clω2Θ2.Substituting (21) into (24) and considering an optimum reactive control, the power extracted becomes(25)PE¯ω=Jφω2Ξ5Θ2.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(26)φ=φISWEC0.5645≃240RPM.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.(27)φ=120RPM.


## 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Ψ(x→,t), which must satisfy Laplace’s equation (i.e., ∇2Ψ=0). The flow pressure can also be calculated in the xyzo frame Bernoulli’s equation [2]. The fully nonlinear kinematic and dynamic free-surface, body boundary, and radiation conditions are, respectively [8],(28)ddt+∇Ψx→,t·∇z-ηx,y,t=0onz=ηx,y,t,(29)Ψdt=-gη-12∇Ψ·∇Ψonz=ηx,y,t,(30)∂Ψx→,t∂n=V→B·n→onSB,(31)∇Ψ⟶0atS∞.The total disturbance potentialΨ(x→,t) can be written as a superposition of two flows, the local ϕl (i.e., instantaneous fluid response due to the impulsive body motion) and the memory ψ (i.e., wave) [5, 8].(32)Ψx→,t=ϕlx→,t+ψx→,t.To arrive at the linear boundary conditions we apply a Taylor expansion about the mean body position for the body condition and aboutz=0 for the free surface, keeping only its linear terms [8].(1)
Kinematic free-surface condition is(33)∂η∂t=∂ϕl∂z+∂ψ∂zapplied  onz=0.(2)
Dynamic free-surface condition is(34)∂∂tϕl+ψ=-gη-12∇ϕl·∇ϕlapplied  onz=0.(3)
Body boundary condition is(35)∂ϕl∂n=∂δ∂t·n→-∇ψ·n→applied  onSB.Equations (33) and (34) will be our evolution equations for the free-surface deformation η and wave potential ψ. Before, however, we must calculate the unknown ϕl and the derivatives of ψ, as we are only given the value of the wave potential at SF and the derivative of the local potential at SB.The displacementδ→(x→,t) of any point within the rigid body can be described by combining translation ξ→T(t) and rotation ξ→R(t) [8].(36)δ→x→,t=ξ→Tt+ξ→Rt×x→.Substituting (36) into (35) yields the body boundary condition in Ogilvie and Tuck’s notation, minus the m-terms. These omitted terms would be fundamental if our body had a translation velocity [8, 32].(37)∂ϕl∂n=∑j=16∂ξj∂t.Particularly for the pressure calculation, the linearized Bernoulli equation can be decomposed into localpl, memory pm, and hydrostatic ph components. The total pressure is p=pl+pm+ph [8].(38)pl=-ρ∂ϕl∂t,pm=-ρ∂ψ∂t,ph=-ρgz.Finally, the force acting on the body is the integration of the pressure onSB.(39)Fj=∬SBp·njdSforj=1,…,6.


### 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 thea0 term written in [8](40)ajk0=ρ∬SBNknjdS,where(41)Nk=0onz=0Nk∂n=nkonSBfork=1,…,6n1,n2,n3=n→n4,n5,n6=x→×n→.


### 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 (ψ=ψ+ΨI and η=η+ηI). The body boundary condition becomes [8](42)∂ψ∂n=-ΨI∂nonSB.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 atSF and the derivative of the local potential at SB. 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]:(43)2πΨx→-∬SF∪SB∂Ψx→′∂nGx→′;x→dx′+∬SF∪SBΨx→′∂Gx→′;x→∂ndx′=0.Finally, AEGIR is called a Rankine panel method due to its choice on the Rankine source potential as the known function in the integral:(44)Gx→′;x→=1x→-x→′.


## 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 thea0 term written in [8](40)ajk0=ρ∬SBNknjdS,where(41)Nk=0onz=0Nk∂n=nkonSBfork=1,…,6n1,n2,n3=n→n4,n5,n6=x→×n→.


## 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 (ψ=ψ+ΨI and η=η+ηI). The body boundary condition becomes [8](42)∂ψ∂n=-ΨI∂nonSB.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 atSF and the derivative of the local potential at SB. 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]:(43)2πΨx→-∬SF∪SB∂Ψx→′∂nGx→′;x→dx′+∬SF∪SBΨx→′∂Gx→′;x→∂ndx′=0.Finally, AEGIR is called a Rankine panel method due to its choice on the Rankine source potential as the known function in the integral:(44)Gx→′;x→=1x→-x→′.


## 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 become(45)ψx→,t≃∑jψjtBjx→,ηx→,t≃∑jηjtBjx→,∂ψ∂nx→,t≃∑jψzjtBjx→with(46)Bx,y=bxby.The basis function within a panel is a second-order B-spline, meaning first and second derivatives can be obtained analytically [8]:(47)bx=12lp2x+3lp22,on-3lp2<x<-lp21lp2-x2+3lp24,on-lp2<x<lp212lp2-x+3lp22,onlp2<x<3lp2,where lp 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 (Ax=b) whose size depends on the number of panels in the problem [8].(48)∂ηj∂tBij=∂ϕl∂z+ψj∂zBij,∂ψj∂tBij=-ηjgBij,2πψjBij+ψjDij-∂ψj∂zSij=0,where Dij and Sij are the dipole and source influence coefficients, respectively:(49)Bij=Bjx→i=Bi-j,Dij=∬-∞∞Bjx→′dGx→i;x→′dndx→′=Di-j,Sij=∬-∞∞Bjx→′Gx→i;x→′dx→′=Si-j.


### 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:(50)ηjn+1-ηjnΔtBij=∂ϕl∂z+∂ψjn∂zBij,ψjn+1-ψjnΔtBij=-ηjn+1gBij,2πψjn+1Bij+ψjn+1Dij-∂ψjn+1∂zSij=0,where n is the variable value at the current (nth) 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(51)y→n+1/2∗=y→n+Δt2dy→ndt.(2)
Implicit Euler corrector is(52)y→n+1/2∗∗=y→n+Δt2dy→n+1/2∗dt.(3)
Midpoint rule predictor is(53)y→n+1∗=y→n+Δtdy→n+1/2∗∗dt.(4)
Simpson’s rule corrector is(54)y→n+1=y→n+Δt6dy→ndt+dy→n+1/2∗dt+dy→n+1/2∗∗dt+dy→n+1∗dt.Figure5 better illustrates the information exchange between AEGIR and the MATLAB code. The latter calculates the true state vector from AEGIR’s bare hull estimation.Figure 5
Schematics describing the information change between AEGIR and the MATLAB code for each time step, whereF→I(t) and F→m(t) are the incident and memory force vectors, respectively.


#### 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].(55)dψdt=-gη,dηdt=dψdz-2νη+ν2gψ,where ν is an artificial viscosity imposed on the waves for their decay near the edges.


## 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 become(45)ψx→,t≃∑jψjtBjx→,ηx→,t≃∑jηjtBjx→,∂ψ∂nx→,t≃∑jψzjtBjx→with(46)Bx,y=bxby.The basis function within a panel is a second-order B-spline, meaning first and second derivatives can be obtained analytically [8]:(47)bx=12lp2x+3lp22,on-3lp2<x<-lp21lp2-x2+3lp24,on-lp2<x<lp212lp2-x+3lp22,onlp2<x<3lp2,where lp 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 (Ax=b) whose size depends on the number of panels in the problem [8].(48)∂ηj∂tBij=∂ϕl∂z+ψj∂zBij,∂ψj∂tBij=-ηjgBij,2πψjBij+ψjDij-∂ψj∂zSij=0,where Dij and Sij are the dipole and source influence coefficients, respectively:(49)Bij=Bjx→i=Bi-j,Dij=∬-∞∞Bjx→′dGx→i;x→′dndx→′=Di-j,Sij=∬-∞∞Bjx→′Gx→i;x→′dx→′=Si-j.


## 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:(50)ηjn+1-ηjnΔtBij=∂ϕl∂z+∂ψjn∂zBij,ψjn+1-ψjnΔtBij=-ηjn+1gBij,2πψjn+1Bij+ψjn+1Dij-∂ψjn+1∂zSij=0,where n is the variable value at the current (nth) 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(51)y→n+1/2∗=y→n+Δt2dy→ndt.(2)
Implicit Euler corrector is(52)y→n+1/2∗∗=y→n+Δt2dy→n+1/2∗dt.(3)
Midpoint rule predictor is(53)y→n+1∗=y→n+Δtdy→n+1/2∗∗dt.(4)
Simpson’s rule corrector is(54)y→n+1=y→n+Δt6dy→ndt+dy→n+1/2∗dt+dy→n+1/2∗∗dt+dy→n+1∗dt.Figure5 better illustrates the information exchange between AEGIR and the MATLAB code. The latter calculates the true state vector from AEGIR’s bare hull estimation.Figure 5
Schematics describing the information change between AEGIR and the MATLAB code for each time step, whereF→I(t) and F→m(t) are the incident and memory force vectors, respectively.


### 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].(55)dψdt=-gη,dηdt=dψdz-2νη+ν2gψ,where ν is an artificial viscosity imposed on the waves for their decay near the edges.


## 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:(50)ηjn+1-ηjnΔtBij=∂ϕl∂z+∂ψjn∂zBij,ψjn+1-ψjnΔtBij=-ηjn+1gBij,2πψjn+1Bij+ψjn+1Dij-∂ψjn+1∂zSij=0,where n is the variable value at the current (nth) 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(51)y→n+1/2∗=y→n+Δt2dy→ndt.(2)
Implicit Euler corrector is(52)y→n+1/2∗∗=y→n+Δt2dy→n+1/2∗dt.(3)
Midpoint rule predictor is(53)y→n+1∗=y→n+Δtdy→n+1/2∗∗dt.(4)
Simpson’s rule corrector is(54)y→n+1=y→n+Δt6dy→ndt+dy→n+1/2∗dt+dy→n+1/2∗∗dt+dy→n+1∗dt.Figure5 better illustrates the information exchange between AEGIR and the MATLAB code. The latter calculates the true state vector from AEGIR’s bare hull estimation.Figure 5
Schematics describing the information change between AEGIR and the MATLAB code for each time step, whereF→I(t) and F→m(t) are the incident and memory force vectors, respectively.


## 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].(55)dψdt=-gη,dηdt=dψdz-2νη+ν2gψ,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ω2=gk, where ω and k 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., T = 5 s).


#### 5.1.1. Domain Size

Figure6 shows the final geometry utilized for the sensibility analyses.Figure 6
NURBS representation of IOwec’s hull.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.Table 3
Initial domain size.
Extent
Direction
Reference
Value
Unit
Domain
Upstream
λ
112.5
m
Downstream
-
3
⋅
λ
−337.5
m
Side
3
⋅
L
135
m
Beach
Upstream
-
45
m
Downstream
-
112.5
m
Side
-
54
mTo be certain the domain was sized correctly, the longest wave (T = 12 s) with amplitude A = 1 m was simulated for three domain sizes of 1.0, 1.5, and 2.0 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 0.2%. The domain size corresponding to the 1.5 multiplier was selected.Figure 7
Domain sensitivity analysis.


#### 5.1.2. Mesh Density

Three mesh density factors of0.5, 1.0, and 1.5, 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.Figure 8
Free-surface mesh sensitivity analysis.


#### 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 theT = 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 Δt = 0.2 s, which was chosen for all subsequent runs.Figure 9
Time-step sensitivity analysis.Finally, Table4 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.Table 4
Final domain size.
Extent
Direction
Reference
Value
Unit
Domain
Upstream
1.5
⋅
λ
168.75
m
Downstream
-
4.5
⋅
λ
−506.25
m
Side
4.5
⋅
L
202.5
m
Beach
Upstream
-
67.5
m
Downstream
-
202.5
m
Side
-
81
m
Time
Ramp time
-
60
s
Time step
-
0.2
sFigure 10
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. Figures11 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.Figure 11
Heave RAO comparison between WAMIT and AEGIR for the bare hull case (i.e., no gyroscope).Figure 12
Pitch RAO comparison between WAMIT and AEGIR for the bare hull case (i.e., no gyroscope).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.Figure 13
Heave scattering force coefficient.Figure 14
Pitch scattering moment coefficient.We now oscillate the hull under a prescribed motion of amplitudeA 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(56)aii=Xicos⁡ϵf-ϵmΞiω2,bii=-Xisin⁡ϵf-ϵmΞiω.Figures15(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 λ/L starting at 3. 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.Figure 15
(a) Heave added-mass coefficient. (b) Heave damping coefficient. (c) Pitch added-mass coefficient. (d) Pitch damping coefficient.
(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.Figure 16
Average power extracted for different multiples of the base spin rate. This is estimated through (25), since we still do not know the appropriate PTO damping. The waves have amplitude of 0.1 m.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 Figure17.Figure 17
Pitch RAO for different multipliers of the basis spin rate. We can see a shift in the natural frequency towards longer waves as the spin is increased and the system gets stiffer.The best spin found, considering our 8 s resonant period, is1.25φ˙ (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.(57)cl=2P¯Eωω2Θ2≃455,667.00kg·m2/s.The actual instantaneous power extracted can now be integrated to calculate the average power:(58)PE¯t=1T∫0Tclθ˙2dt.We can observe, from Figure18, 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 5cl is attained, after which the device quickly loses energy extraction efficiency, especially for longer waves.Figure 18
Average power extracted over one wave period. The optimum damping identified is2cl, after which the device quickly loses efficiency for longer lengths. All waves have amplitude of 0.1 m.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. Figure19 shows how the pitch RAO changes for all the damping variation.Figure 19
Pitch RAO for all the damping values considered. The optimal damping of2cl tries to minimize the motion throughout all wave lengths, while still retaining the natural frequency at 8 s.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.(59)CW=PE¯tF¯=32πρg2T2H2∫tt+TPtdt.Figure20 shows the IOwec is able to, at resonance, extract a monochromatic wave 25 m wide, 1.25 times its beam.Figure 20
Capture width of the IOwec. The device is able to, at resonance, extract a wave front equivalent to 1.25 times its beam.In all cases, the yaw torques induced by both counterrotating flywheels are seen to cancel one another. Figure21 shows the torques from both gyroscopes for the T = 8 s wave, with frequency exactly twice that of the wave.Figure 21
Yaw torque from both counterrotating gyroscopes for theT = 8 s wave. Their summation goes perfectly to zero, and the frequency is exactly twice that of the wave.


## 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ω2=gk, where ω and k 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., T = 5 s).


### 5.1.1. Domain Size

Figure6 shows the final geometry utilized for the sensibility analyses.Figure 6
NURBS representation of IOwec’s hull.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.Table 3
Initial domain size.
Extent
Direction
Reference
Value
Unit
Domain
Upstream
λ
112.5
m
Downstream
-
3
⋅
λ
−337.5
m
Side
3
⋅
L
135
m
Beach
Upstream
-
45
m
Downstream
-
112.5
m
Side
-
54
mTo be certain the domain was sized correctly, the longest wave (T = 12 s) with amplitude A = 1 m was simulated for three domain sizes of 1.0, 1.5, and 2.0 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 0.2%. The domain size corresponding to the 1.5 multiplier was selected.Figure 7
Domain sensitivity analysis.


### 5.1.2. Mesh Density

Three mesh density factors of0.5, 1.0, and 1.5, 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.Figure 8
Free-surface mesh sensitivity analysis.


### 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 theT = 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 Δt = 0.2 s, which was chosen for all subsequent runs.Figure 9
Time-step sensitivity analysis.Finally, Table4 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.Table 4
Final domain size.
Extent
Direction
Reference
Value
Unit
Domain
Upstream
1.5
⋅
λ
168.75
m
Downstream
-
4.5
⋅
λ
−506.25
m
Side
4.5
⋅
L
202.5
m
Beach
Upstream
-
67.5
m
Downstream
-
202.5
m
Side
-
81
m
Time
Ramp time
-
60
s
Time step
-
0.2
sFigure 10
Free-surface and body surfaces after meshing.


## 5.1.1. Domain Size

Figure6 shows the final geometry utilized for the sensibility analyses.Figure 6
NURBS representation of IOwec’s hull.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.Table 3
Initial domain size.
Extent
Direction
Reference
Value
Unit
Domain
Upstream
λ
112.5
m
Downstream
-
3
⋅
λ
−337.5
m
Side
3
⋅
L
135
m
Beach
Upstream
-
45
m
Downstream
-
112.5
m
Side
-
54
mTo be certain the domain was sized correctly, the longest wave (T = 12 s) with amplitude A = 1 m was simulated for three domain sizes of 1.0, 1.5, and 2.0 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 0.2%. The domain size corresponding to the 1.5 multiplier was selected.Figure 7
Domain sensitivity analysis.


## 5.1.2. Mesh Density

Three mesh density factors of0.5, 1.0, and 1.5, 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.Figure 8
Free-surface mesh sensitivity analysis.


## 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 theT = 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 Δt = 0.2 s, which was chosen for all subsequent runs.Figure 9
Time-step sensitivity analysis.Finally, Table4 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.Table 4
Final domain size.
Extent
Direction
Reference
Value
Unit
Domain
Upstream
1.5
⋅
λ
168.75
m
Downstream
-
4.5
⋅
λ
−506.25
m
Side
4.5
⋅
L
202.5
m
Beach
Upstream
-
67.5
m
Downstream
-
202.5
m
Side
-
81
m
Time
Ramp time
-
60
s
Time step
-
0.2
sFigure 10
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. Figures11 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.Figure 11
Heave RAO comparison between WAMIT and AEGIR for the bare hull case (i.e., no gyroscope).Figure 12
Pitch RAO comparison between WAMIT and AEGIR for the bare hull case (i.e., no gyroscope).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.Figure 13
Heave scattering force coefficient.Figure 14
Pitch scattering moment coefficient.We now oscillate the hull under a prescribed motion of amplitudeA 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(56)aii=Xicos⁡ϵf-ϵmΞiω2,bii=-Xisin⁡ϵf-ϵmΞiω.Figures15(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 λ/L starting at 3. 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.Figure 15
(a) Heave added-mass coefficient. (b) Heave damping coefficient. (c) Pitch added-mass coefficient. (d) Pitch damping coefficient.
(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.Figure 16
Average power extracted for different multiples of the base spin rate. This is estimated through (25), since we still do not know the appropriate PTO damping. The waves have amplitude of 0.1 m.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 Figure17.Figure 17
Pitch RAO for different multipliers of the basis spin rate. We can see a shift in the natural frequency towards longer waves as the spin is increased and the system gets stiffer.The best spin found, considering our 8 s resonant period, is1.25φ˙ (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.(57)cl=2P¯Eωω2Θ2≃455,667.00kg·m2/s.The actual instantaneous power extracted can now be integrated to calculate the average power:(58)PE¯t=1T∫0Tclθ˙2dt.We can observe, from Figure18, 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 5cl is attained, after which the device quickly loses energy extraction efficiency, especially for longer waves.Figure 18
Average power extracted over one wave period. The optimum damping identified is2cl, after which the device quickly loses efficiency for longer lengths. All waves have amplitude of 0.1 m.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. Figure19 shows how the pitch RAO changes for all the damping variation.Figure 19
Pitch RAO for all the damping values considered. The optimal damping of2cl tries to minimize the motion throughout all wave lengths, while still retaining the natural frequency at 8 s.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.(59)CW=PE¯tF¯=32πρg2T2H2∫tt+TPtdt.Figure20 shows the IOwec is able to, at resonance, extract a monochromatic wave 25 m wide, 1.25 times its beam.Figure 20
Capture width of the IOwec. The device is able to, at resonance, extract a wave front equivalent to 1.25 times its beam.In all cases, the yaw torques induced by both counterrotating flywheels are seen to cancel one another. Figure21 shows the torques from both gyroscopes for the T = 8 s wave, with frequency exactly twice that of the wave.Figure 21
Yaw torque from both counterrotating gyroscopes for theT = 8 s wave. Their summation goes perfectly to zero, and the frequency is 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 Section2.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.



---

*Source: 1710253-2018-01-17.xml*