Stochastic analysis of nonlinear wave energy converters via statistical linearization

Applied Ocean Research 95 (2020) 102023

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

Stochastic analysis of nonlinear wave energy converters via statistical linearization

T

⁎,a,b

L.S.P. Silva

, N.Y. Sergiienkob, C.P. Pescec, B. Dingb, B. Cazzolatob, H.M. Morishitaa

a

University of São Paulo, Escola Politécnica, Dept. of Naval Architecture and Ocean Engineering, Brazil The University of Adelaide, Ocean Renewable Energy Research Group, School of Mechanical Engineering, Australia c University of São Paulo, Escola Politécnica, Offshore Mechanics Laboratory, Brazil b

ARTICLE INFO

ABSTRACT

Keywords: Statistical linearization Stochastic analysis Nonlinear wave energy converters Spectral domain,

In the initial stages of wave energy converter development, a wide range of design parameters should be tested and analyzed across all operational sea states. As a result, it is required to conduct a large number of time domain simulations that are computationally expensive. The following study presents the application of a fast and highfidelity alternative technique, named statistical linearization (SL), that can be used for the stochastic analysis of nonlinear wave energy converters. The method consists of an approximate solution that allows one to quickly estimate the contribution of nonlinear terms using a probabilistic model. The application of the statistical linearization is demonstrated using three conceptually different wave energy converters, a point absorber, an oscillating wave surge converter, and an oscillating water column, with distinct nonlinear dynamics. The performance of SL is assessed against nonlinear time-domain simulations in terms of spatial response distributions and power spectral densities of responses and excitation loads, mean offsets, mean absorbed power and equivalent linear terms. The results show a good agreement between the two types of models for all wave energy converters considered, while showing that SL is approximately 3 to 4 orders of magnitude faster than the time domain model.

1. Introduction Ocean wave energy has a great potential to contribute to the global renewable energy matrix due to its high power density and relatively predictable energy output [1]. Even though over a thousand prototypes have been proposed during the past decades [2,3] and novel concepts are being developed and patented every year, only a few full-scale devices have been tested in an open sea environment [4], and even fewer have been connected to the grid [5]. However, it is expected that in the near future, the wave technology will converge to one or several designs suitable for different installation locations (onshore, nearshore, and offshore). In order to enhance further development of wave energy converters (WECs), there is a need for a fast and high-fidelity modeling approach to be used at both preliminary design and optimization stages. Frequency-domain (FD) models based on the linear potential flow theory have been, and remain, the most widely used technique in the initial stages of WECs development [6]. It allows researchers to quickly estimate important parameters of the system such as the resonance frequency, average power output, steady-state response, etc. In these models, the system is described by a set of linear equations, where the

⁎

nonlinear forces are linearized usually by assuming small displacements and velocities around static equilibrium configurations. Therefore, the application of FD models is limited to systems that operate in the linear or close to linear regions. However, as WECs commonly operate close to resonant conditions to enhance power production, leading to large displacements, there is a range of nonlinear effects and forces that significantly influence their performance [7]. As a result, FD models might lead to an unreliable prediction of the system dynamics and mean power absorbed. More complex models capable of incorporating relevant nonlinear dynamics of the system are usually built in the time domain (TD) [8]. These models allow one to include nonlinear effects associated with the wave-device interaction [7], as well as the nonlinear dynamics of the power take-off (PTO) system [9,10], and thus can be used to estimate the WEC performance across a wide range of irregular sea states [11]. However, the computational cost associated with these TD simulations is considerably higher compared to a FD analysis. Especially, this issue becomes evident when there is a need to optimize the entire system with a large number of varying parameters [12] or to simulate the performance of the wave farm with a large number of WECs involved

Corresponding author at: The University of Adelaide, North Terrace Campus, University of Adelaide, Adelaide, SA 5005, Australia. E-mail address: [email protected] (L.S.P. Silva).

https://doi.org/10.1016/j.apor.2019.102023 Received 26 August 2019; Received in revised form 11 December 2019; Accepted 13 December 2019 0141-1187/ © 2019 Elsevier Ltd. All rights reserved.

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

[13]. Therefore, even though the TD model produces more reliable estimations compared to FD models, it might not be considered as an effective approach to be used during preliminary design stages or running optimization procedures. Based on the modeling requirements, this work deals with the nonlinear dynamics of WECs via the statistical linearization (SL) technique, which offers a systematic approach to estimate the contribution of the nonlinear terms in FD. The SL technique has been widely used in other engineering disciplines (e.g. earthquake, aerospace and offshore engineering) to study the stationary random response of nonlinear systems [14]. Despite its relative simplicity, reliable output and low computational, only a few research groups in the field of wave energy have used and adapted SL to their needs. The method was initially applied to wave energy devices in [15–18], which the authors referred to as spectral-domain modeling. However, the technique employed initially presented several limitations as highlighted by [19], which proposed a harmonic balance method to deal with the nonlinearities [19–22]. To overcome the limitations of some of the previous derivations, this work uses a general formulation of the SL technique as described in [14,23]. The SL has been applied to WECs such as an oscillating wave surge converter (OWSC) [15,24], point absorber (PA) [16,25–27], and an oscillating water column (OWC) [17,28,29]. The main nonlinear effects associated with a wave-structure interaction that have been modeled using SL are the viscous drag forces [15,24,27], water column dynamics [28,29] for OWCs, nonlinear hydrostatic stiffness [25,27], nonlinear behavior of the mooring lines and nonlinear PTO dynamics [16,26,27]. It should be noted that the majority of these studies are limited to the comparison of the predicted output from the equivalent linear model with (i) numerical simulations of the original nonlinear model, or with (ii) experimental wave-tank data, as in [17]. However, to the best of the authors’ knowledge, none of them have performed a comprehensive analysis with a general formulation including a detailed analysis of the equivalent linear terms in order to identify the relevance and compare different nonlinear effects acting on the WEC. The present study aims to demonstrate the potential of the SL technique using a general formulation, in which three conceptually different WECs with different nonlinear PTO systems and nonlinearities are investigated (see Fig. 1). Besides the diversity of WECs, the SL procedures differ in terms of the number of DOFs, and additional computational routines that are specified in Sections 3.1 and 3.2. Moreover, the analysis of each device focuses on different aspects of the SL technique, such as the analysis of the response probability distribution function (PDF), power spectral density (PSD) of the response,

equivalent linear terms, PDF and PSD of the force, force modulation, mean offsets and mean absorbed power. The paper is organized as follows. Initially, the basic principles of the FD analysis are presented in Section 2, followed by the mathematical formulation of SL in Section 3. Then, the description and numerical models of three WECs under consideration are presented and simulated in Section 4. Finally, the main outcomes are discussed in Section 5. 2. Frequency domain analysis Since the SL is formulated based on the FD model, this work commences with the basic concepts and terms of the FD analysis that will be widely utilized in the subsequent sections. In this regard, consider a general multi-degree-of-freedom (MDOF) linear system that can be described by the differential equation of the form: (1)

Mq¨ + Bq + Kq = Q (t ),

n × 1 are the vectors of generalized displacements, where q, q, and q¨ n× n are mass/ velocity and acceleration respectively; M, B, and K inertia, damping and stiffness matrices respectively; and Q(t) is a vector of generalized excitation forces acing on the system. The steady-state response of the system under stochastic loads can be generally determined based on the transfer function of the system. In such formulations, it is assumed that the system is excited by a harmonic load and follows a harmonic response described by a complex ^ ei t } . Based on this assumption, the ^ , where q (t ) = Re {q amplitude q transfer function matrix of the system can be written as: 2M

( )=[

(2)

+ i B + K] 1 .

If the stochastic nature of the force Q(t) can be described by its PSD matrix Sf, then the PSD matrix of the response can be obtained as:

Sq ( ) =

( ) Sf ( )

H(

(3)

),

where Sq is the response spectrum matrix, and denotes the Hermitian transpose. For a single-input single-output system, Eq. (3) reduces to:

()H

(4)

Sq ( ) = | ( )|2 Sf ( ). 3. Statistical linearization formulation

Now consider a MDOF nonlinear system described by the equation of motion:

Mq¨ + Bq + Kq +

(q, q , q¨) = Q (t ),

(5)

where (q, q , q¨) is a vector of nonlinear terms, functions of the generalized coordinate vector q and its first and second time derivatives (q and q¨ ). In order to obtain an approximate solution to Eq. (5), it is useful to introduce an equivalent linear system as: n ×1

(M + Meq) q¨ + (B + Beq) q + (K + K eq) q = Q (t ),

(6)

where Meq, Beq and Keq represent equivalent inertia, damping and stiffness matrices respectively coming from the nonlinearities. The SL procedure implies that the difference between the actual nonlinear system (Eq. (5)) and the linearized equivalent system (Eq. (6)) is minimized:

min

T

(7)

,

where ⟨ · ⟩ denotes mathematical expectation, and the difference ε is given by:

(q, q, q¨)

Meq q¨

Beq q

K eqq.

(8)

Under the assumption of a Gaussian distribution of the response, the elements of the equivalent linear system can be calculated as [14]:

Fig. 1. Types of WECs and PTO mechanisms investigated. 2

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

i

Meqi, j =

q¨j i

Beqi, j =

(10)

i

,

qj

(11) th

2 (M

+ Meq) + i (B + Beq) + (K + K eq)] 1 ,

(12)

and the response spectrum obtained as:

Sq ( ) =

eq (

H eq (

) Sf ( )

(13)

).

3.1. Parametric forces For some cases when the excitation to the system has a nonlinear character, Eq. (5) can be modified to accommodate parametric excitation [30] in the form:

Mq¨ + Bq + Kq +

(q, q , q¨) =

(q, q , q¨) Qp (t ),

(14)

(q) Qp (t )

Q p (t ) ,

2 q1

(15)

V=

The equivalent linear excitation force spectrum is given by:

S f, eq ( ) =

Sf ( )

which is used in Eq. (13) to account the parametric force contribution.

2 qk

3.2. Treatments of asymmetric nonlinearities

qi qj

=

(25)

(26)

(27)

q1q2

q1q1

q1q2

2 q2

q2q1

q2q2

2 q1

q1q2 2 q2

, (28)

where :

(20)

qk qk

(21)

2 qk

Subtracting Eq. (21) from Eq. (20), an equivalent form of Eq. (5) is obtained:

Mq¨ + Bq + Kq + G (q, q, q¨ ) = Q ,

,

2 qn

sym

and the mean offset is estimated taking the expectation of the response:

(q + q, q, q¨ ) = Q .

q2qn

Sqi qj ( ) d .

0

V=

¯ denotes the mean component such as wave drift force. For where Q such systems, Eq. (5) can be rewritten as:

Kq +

q1qn

2 q2

Sqk ( ) d ,

0

2 q1

(19)

(q + q, q, q¨ ) = Q + Q,

(24)

For a two DOF system, the covariance matrix of the body response and their first time derivatives can be written as:

Similarly, the force vector can be also represented by two components:

Mq¨ + Bq + Kq + Kq +

=

(18)

˜ (t ), ¯ +Q Q (t ) = Q

q¯ ) ,

and the non-diagonal terms (i ≠ j) are based on the cross-spectrum and obtained by:

In general, the nonlinear force vector Q(t) may have a non-symmetric character leading to a constant offset in the system response (the solution of Eq. (5) has a non-zero mean value). Therefore, the system response can be approximately expressed in terms of a mean value ( ¯. ) and a random zero mean component ( ˜. ) [14]:

q (t ) = q¯ + q˜ (t ).

q¯ )TV 1 (˜q

where the diagonal terms (i = j ) are based on the auto-spectrum and obtained from:

(17)

H.

q1q2

sym

(16)

(q) .

1 (q˜ 2

where V is the covariance matrix of q, and |V| denotes the determinant of the covariance matrix. The symmetric covariance matrix of the body response is given by:

resulting in:

=

1 exp (2 )n/2|V|1/2

f (q) =

where the vector of parametric forces Qp is modulated by the matrix (q, q, q¨), which can depend on the system response and its time derivatives. A modulation based on the displacement is investigated in this paper, in which the difference between the parametric excitation force and its equivalent constant modulation Λ is here minimized by taking their expectation based on the Gaussian distribution as:

min

(23)

The excitation force is considered to be a linear transformation of the wave elevation, which is assumed to be described by a Gaussian distribution. Therefore, the excitation force is considered Gaussian. In addition, despite the presence of nonlinearities in the system dynamics, the assumption of a Gaussian distribution of the WEC response is considered valid as long as the nonlinear forces are non-dominant [18], and the system has a unique equilibrium condition. However, ocean waves are non-Gaussian for cases of severe sea states traveling to transitional and shallow water depth [31]; and during severe sea and transient meteorological conditions in deep waters [32]. Similarly, strongly nonlinearities into the system dynamics lead to non-Gaussian response distributions [14], and the modulation of the excitation force is also non-Gaussian. Therefore, the WECs simulated are considered to be operating within a power production mode, where sea states are moderate. Further, the formulation proposed is an approximated solution that should not be used for the analysis of non-stationary random vibration, or for excitation from regular waves. For a MDOF system, the PDF of the displacement based on the Gaussian distribution is given by:

where the subscripts i,j refer to the i row and j column of the corresponding matrices. As a result, the equivalent linear transfer function of the nonlinear system is expressed as:

)=[

(q + q, q, q¨ ) .

3.3. Gaussian assumption

th

eq (

(q + q, q, q¨ )

(9)

,

qj

K eqi, j =

G (q, q, q¨ ) =

,

=

qi qj

=i

qi qj

=

(22)

where: 3

(29)

= 0, 2S ( qk

0

0

0

)d ,

Sqi qj ( ) d , 2S qi qj (

)d .

(30) (31) (32)

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

3.4. Statistical linearization procedure The main idea of the SL technique is to calculate the equivalent linear terms Meq, Beq and Keq in order to obtain the approximate dynamic response of the nonlinear system, such as in Eq. (5). Since there is no analytical solution to this problem, the SL technique utilizes an iterative procedure starting from an initial guess and converging towards approximate results that meet predetermined criteria. The standard step-by-step procedure for the iterative procedure is given by: Step 1: Define the linearized FD model as an initial guess (Meq = Beq = K eq = 0). Set the PSD matrix of the linear force Sf and evaluate the PSD of the body response Sq using Eq. (3) or (4). Step 2: Calculate the PDF of the body response, f(q), using Eq. (24) based on the covariance matrix V, and calculate the mean offsets using Eq. (21). Step 3: Calculate the equivalent linear terms Meq, Beq, Keq using Eqs. (9)–(11) and the equivalent linear excitation force using . Eq. (17). The equivalent linear terms are based on the response characteristics given in Eq. (24). Step 4: Update the equivalent transfer function using Eq. (12) and evaluate the spectral response of the body based on the equivalent system using Eqs. (13) and (17). Step 5: Check the convergence of the new values found in Step 4. If convergence is not achieved, return to Step 2.

Fig. 2. Illustration of a Point Absorber with nonlinear forces acting upon it.

displacement in the heave direction; Fr(t), Fhs,l(t), Fs(t), and Fexc(t) are the radiation force, the linearized hydrostatic force, the linear mechanical restoring force, and the excitation force in the vertical direction respectively (see in Appendix B), which are linear terms in the governing equation. Note that the negative sign occurs as the nonlinear terms are placed on the left side in Eq. (5), while in Eq. (33) the nonlinear terms are placed on the right side. The nonlinear terms, (z, z ), are given by:

The convergence criteria is verified comparing the relative error between the response used to estimate the equivalent coefficients (and mean offsets) with the response generated using those equivalent coefficients, namely comparing each term of the covariance matrix and the mean offset. In this regard, the results are considered converged when the relative error is less than 1%. Typically, simple solvers converge in 3 to 15 iterations depending on the source of nonlinearity and the ratio between the linear and nonlinear terms, and number of DOF. In some cases, an iteration relaxation method is valuable to guarantee the convergence [18], especially when the solution of mean values, such as in Eq. (21), is extremely sensitive to the nonlinear forces.

(z , z ) = Fhs, nl (t ) + Fls (t ) + Fes (t ) + Fvd (t ) + FEG (t ),

where Fhs,nl(t), Fls(t), Fes(t), Fvd, and FEG(t) are the nonlinear hydrostatic restoring force, lateral restoring force, the end-stop force, the viscous drag force, and electrical generator force respectively. The nonlinear hydrostatic restoring force of a semi-submerged sphere is given by a cubic term as [25,27]:

4. Case studies

3

z 3,

(35)

where ρ is the water density, and g denotes the gravitational acceleration. The lateral stiffness is composed of two linear stiffness connected in the transverse direction. For such mooring configuration, the force is described as [27]:

Fls =

2Kls 1

lls z 2 + dls2

z,

(36)

where the length of the unstretched lateral spring is given by lls, the distance between the connection points and the center is given by dls, and Kls is the magnitude of the lateral stiffness (see in [27]). In this work, lls is equal to dls, therefore, the system has only one stable condition. An additional stiffness, named end-stop, is employed in the PTO system to preserve the structure from damage due to the contact between the generator and the PTO extremities. A standard governing equation to describe the end-stop force is given by [27]:

4.1. Point absorber The PA considered in this study is a small bottom-referenced heaving buoy with a spherical shape connected to the electric power generator as illustrated in Fig. 2. The buoy is constrained to move in heave only. 4.1.1. Governing equations Based on the system description and the illustration given in Fig. 2, the dynamics of the heaving PA in the TD can be modeled as [33]:

(z, z ),

g

Fhs, nl =

To demonstrate the capability of the SL technique, three case studies are investigated as shown in Fig. 1: (a) a heaving PA with a nonlinear electric generator, (b) an OWSC connected to the hydraulic PTO system, and (c) an OWC with an air-turbine. Each device contains different sources of nonlinearities, that in some cases, require additional treatments as described in Sections 3.1 and 3.2. It should be noted that this work focuses on the SL procedure, and idealized nonlinear devices are considered. In addition, the magnitudes of nonlinear forces are presented for illustration purposes only, and the devices are not tuned to operate in optimal conditions.

Mz¨ (t ) = Fr (t ) + Fhs, l (t ) + Fs (t ) + Fexc (t )

(34)

0, Fes =

(33)

where M contains the mass of the PA and the PTO system, z denotes the

for les > |z| K es (z les ), for z > les K es (z + les ), for les > z ,

(37)

where Kes is the end-stop stiffness, and les is the distance from the 4

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

generator to the spring contact. The viscous drag force contribution can be written based on Morison’s equation [34]:

1 CD S z |z |, 2

Fvd (t ) =

(38)

where CD denotes the drag coefficient, and S⊥ is the cross-sectional area of the body perpendicular to the z direction. The electric generator exerts a damping associated with the electromagnetic force, which can be represented as a function of the active area of the stator as [35]:

FEG =

(39)

Cpto Afac (z ) z ,

where Cpto denotes the damping magnitude, and Afac is the active area of the stator given by:

0, Afac = 1, 1/ls [1/2(lp + ls )

for |z|

1/2(lp + ls )

for |z|

1/2(lp

Fig. 3. Time series of the buoy displacement in heave (Hs= 2 m and Tp= 5.75 s).

ls )

|z|], else,

(40)

where ls and lp denotes the length of the stator and piston respectively. Based on the electromagnetic force, the mean absorbed power can be obtained by: T

1 P¯EG = T

Cpto Afac z 2dt .

(41)

0

4.1.2. Frequency domain and statistical linearization As stated previously, FD models represent a set of linear equations where some forces are linearized around the mean operating position. Based on Eq. (33), the transfer function in FD model of the PA is simply given by: 2 (M

( )=[

+ A ( )) + i (B ( ) + Cpto) + Ks + Khs, l] 1 ,

Fig. 4. Comparison between the response distribution of the displacements using TD, and its respective PDF using FD and SL (Hs= 2 m and Tp= 5.75 s).

(42)

where A(ω) and B(ω) are hydrodynamic added mass and radiation damping respectively; Ks is the linear mechanical stiffness, and Khs,l is the linear hydrostatic stiffness at the mean position. The natural frequency of the buoy in heave can be formulated as: n

Ks + Khs, l . M + A ( n)

=

(43)

The modal analysis is generally applied as a tool to optimize the dynamics of the system [36]. It is valuable to remember that the FD model is used as an initial approximate solution (Step 1) of the iterative procedure. The equivalent transfer function based on the SL technique described in Section 3 has the form: eq (

)=[

2 (M

Fig. 5. Comparison between TD, SL and FD results of the PSD of the body displacement response (Hs= 2 m and Tp= 5.75 s). The shaded area is the standard deviation of the TD simulations.

+ A ( )) + i (B ( ) + Beq) + Ks + Khs, l + K eq] 1 , (44)

where Beq and Keq account for the equivalent linear terms:

Beq = Beq, EG + Beq, VD =

FEG z

+

K eq = K eq, HSNL + K eq, LS + K eq, ES =

Fvd z

,

Fhs, nl + z

4.1.3. Simulations The parameters of the point absorber WEC used in this study are specified in Table B.2 and relevant details of the simulation including generation of the irregular waves are described in Appendix E. It can be observed from Eq. (33) and (34) that the numerical model does not contain any parametric excitation forces and all nonlinear effects have zero mean. Therefore, the main objective of this case study is to analyze the equivalent linear terms and to demonstrate the reliability of the SL technique by comparing its output with FD (Eq. (42)) and TD simulations. As a result, the standard SL procedure described in Section 3.4 is used to estimate the WEC performance over the range of sea states. The time series of the first 500 s of the buoy displacement in heave in irregular sea of Hs= 2 m and Tp= 5.75 s is shown in Fig. 3. In order to demonstrate the applicability of the SL to the PA WEC, it must be shown that the buoy displacement follows a Gaussian

(45)

Fls + z

Fes z

(46)

The power absorbed by the electric generator based on Eq. (41) is given by:

PEG = Cpto Afac z 2 .

(47)

However, the SL allows writing the mean power based on the equivalent coefficient as:

PEG = Beq, EG where

z

2 z.

(48)

denotes the standard deviation of the velocity. 5

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

allows the identification of the resonance frequency of the spectrum that might be shifted due to the nonlinearities. Although only a single case is illustrated for the response distribution and PSD, these similarities are preserved to some extent for all conditions simulated in this case study. Once the Gaussian assumption is verified, the analysis of the equivalent linear terms is assessed. For comparison purposes, equivalent coefficients in TD are estimated to compare with the SL results. In this regard, the amount of work done by the equivalent linear system is set to be equivalent to the nonlinear one, namely:

K eq T

T 0

q2dt =

1 T

T 0

(q, q, q¨) qdt ,

(49)

where T is the period of time simulated. Analogously, the amount of power dissipated by the equivalent linear system is set to be equivalent to the nonlinear system:

Beq T

T 0

q2dt =

1 T

T 0

(q, q, q¨) qdt .

(50)

The equivalent stiffness and damping coefficients calculated over the range of peak wave periods (Hs= 2 m) are shown in Fig. 6 a and b. Such plots demonstrate a clear comparison between nonlinear forces and their effect on the system dynamics. It can be seen from Fig. 6 a that for the PA selected as an example in this case study, the lateral stiffness associated with the nonlinear mooring force Keq,LS has the greatest effect on the system dynamics as compared to other nonlinear forces (hydrostatic restoring force and force associated with an end-stop). Moreover, the value of this equivalent term Keq,LS has the same order of magnitude as the linear mechanical stiffness Ks = 100 kN/m. Also, it is interesting to note that the nonlinear component of the hydrostatic restoring force imposes a negative stiffness effect as Keq,HSNL < 0. Fig. 7 shows the mean power estimated using TD, SL, and FD models over the range of peak wave periods (Hs= 2 m). The level of absorbed power estimated using SL coincide with results obtained from the TD model showing the effectiveness of the SL technique. Interestingly, the FD model predicts twice as much output power from the system. The significant discrepancy between the models is mainly because the viscous damping force is not considered in the traditional FD models, and the system is not optimized. Moreover, the resonance frequency of the fully nonlinear system (refer to TD and SL) is shifted towards lower peak wave periods (higher frequencies) due to the additional contribution of all nonlinear stiffnesses. As mentioned previously, this change in the resonance frequency can be observed in the PSD of the response, such as in Fig. 5. Fig. 8 a and b demonstrate how the equivalent stiffness and damping coefficients change with increasing wave heights assuming the same peak wave period of Tp= 5.75 s. The equivalent terms follow the expected behavior in a qualitative and quantitative sense. Thus, the contribution of the end-stop mechanism Keq,ES becomes more pronounced with increasing wave height as the PA displacement reaches a certain limit; and the equivalent lateral stiffness Keq,LS approaches a certain level equal to the sum of both stiffnesses. In terms of the equivalent damping coefficients, the viscous damping Beq,VD increases with wave height, while the damping associated with the electric generator Beq,EG stays relatively constant across all wave heights due to the nonlinearity and the parameters utilized. The effect of the wave height on the mean absorbed power estimated using FD, TD and SL models is shown in Fig. 9. For relatively small values of the wave height (Hs < 1 m), the system operates in a linear region and the power predicted using FD and TD models are similar. However, the difference between models become more obvious with an increasing wave height. Moreover, for the linear system (FD model) the power output depends on the wave height quadratically, while the dependency between power and wave height behaves more linearly when the nonlinear dynamics of the WEC are considered (TD and SL).

Fig. 6. Comparison of the equivalent coefficients of the Point Absorber for a range of Tp (Hs= 2 m). The dashed line refers to SL, the solid line to TD, and the shaded area is the standard deviation of the TD simulations.

Fig. 7. Comparison of the mean power absorbed by the Point Absorber for a range of Tp (Hs= 2 m). The shaded area is the standard deviation of the TD simulations.

distribution. Based on that, the PDF of the buoy response calculated using FD and SL models are compared with the response distribution of TD simulations in Fig. 4. It can be seen that the response distribution of the nonlinear TD simulation recovers the theoretical Gaussian distribution using the SL technique, by preserving the standard deviation. The Gaussian distribution shown in Fig. 4 is based on the covariance matrix, whose terms are obtained using the PSD of the response as in Eqs. (24) to (32). Fig. 5 illustrates, the PSD of the PA response in heave. As can be observed, the spectral energy of the displacement using TD and SL are comparable in the entire range of frequency. In addition, the PSD 6

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

Fig. 10. Illustration of an Oscillating Wave Surge Converter with the nonlinear torques acting upon it.

where, I is the inertia of the paddle and ballast at the connection, Mr(t) denotes the radiation torque, and Ms(t) includes the linear restoring torque from the PTO and hydrostatic restoring torque (see in Appendix C), which are the linear terms in the governing equation; and Φ(θ) is the torque modulation (explained below). The nonlinear terms, ( , ), considered in this case are given by: (52)

( , ) = Mgb (t ) + Mvd (t ) + MHpto (t ),

where Mgb is the torque due to the balance between the weight and buoyancy force, Mvd(t) is the torque due to the viscous drag loads, and MHpto(t) is the torque of the hydraulic PTO system. Assuming that the center of gravity coincides with the center of buoyancy, where their distance from the hinged connection is given by d, the restoring torque due to the balance between the weight and buoyancy force of the OWSC can be expressed as [8]:

Fig. 8. Comparison of the equivalent coefficients of the Point Absorber for a range of Hs (Tp= 5.75 s). The dashed line refers to SL, the solid line to TD, and the shaded area is the standard deviation of the TD simulations.

Mgb =

g( V

(53)

M ) d sin( ),

where V is the volume displaced, and M accounts for the mass of the flap and ballast. The moment caused by the viscous drag can be modeled as a summation of each area section contribution, assuming an infinitesimal area, the moment can be expressed in the integral form as [38]: L

Mvd (t ) =

MHpto (t ) =

4.2. Oscillating wave surge converter

1 P¯Hpto = T

4.2.1. Governing equations Based on the system description and the illustration given in Fig. 10, the governing equation in TD of the OWSC around the hinged connection can be written as:

( ) Mexc (t ),

1 Cv wL4 8

,

(54)

Hptosgn( (t )),

(55)

where Hpto is the Coulomb coefficient, for which its magnitude depends on the cross-sectional area of the hydraulic cylinder and pressure difference during a cycle of stroke motion. For this type of a hydraulic PTO system, the mean absorbed power is given by:

The OWSC considered in this study is a buoyant flap type device connected to the hydraulic PTO system as illustrated in Fig. 10. The device is constrained to move in pitch only due to the hinged connection. The flap has a thickness of 2 m, and the water depth is 13 m [37].

( , )+

1 Cv l2 l wdl = 2

where w and L denotes the width and the submerged length of the OWSC flap respectively, and Cv denotes the viscous drag coefficient. A simplified model is used to describe the hydraulic system, in which the hydraulic moment is modeled by a Coulomb-like torque regulated by the angular speed as [8,16]:

Fig. 9. Comparison of the mean power absorbed by the Point Absorber for a range of Hs (Tp= 5.75 s). The shaded area is the standard deviation of the TD simulations.

I ¨ (t ) = Mr (t ) + Ms (t )

0

T

Hpto | |dt . 0

(56)

The wave excitation torque is highly dependent on the angular position of the flap due to its large motion amplitudes. Therefore, it would be a reasonable approximation to assume that the wave torque depends on the instantaneous flap angle [39] and can be modeled as a parametric excitation modulated by the cosine of the angle of rotation.

(51) 7

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

( )=[

2 (I

+ A ( )) + i B ( ) + Ks + g ( V

M ) d] 1 .

(58)

Note that the FD model is not able to represent the hydraulic PTO. Therefore, its damping contribution and mean power absorbed cannot be estimated. In the FD model, the modulation of the parametric excitation torque is neglected. Therefore, only the frequency dependent term is included; the spectrum of excitation torque can be expressed as: (59)

Sf ( ) = |HM, ( )|2 S ( ),

^ ( ) is the torque amplitude operator in the angular direcwhere H M, tion, and Sξ(ω) is the wave spectrum. Based on the free undamped FD formulation, the natural frequency of the system based on the inverted pendulum concept is given by:

Fig. 11. Time series of the exciting torque acting on the OWSC (Hs= 3 m and Tp= 12 s).

n

=

Ks + g ( V I + A(

M)d n)

.

(60)

The equivalent transfer function between the wave torque and body angular displacement based on the SL will have the form: eq (

)=[

2 (I

+ A ( )) + i (B ( ) + Beq) + Ks + K eq] 1 ,

(61)

with:

MHpto

Beq = Beq, Hpto + Beq, vd =

K eq = K eq, gb =

Mgb

.

+

Mvd

,

(62)

(63)

Due to the presence of the parametric torque, the SL procedure should be modified according to Section 3.1. As a result, the torque modulation in the SL technique is given by:

Fig. 12. Comparison between TD, SL and FD results of the torque distribution and its respective PDF estimation (Hs= 3 m and Tp= 12 s).

=

( ) = cos( ) ,

(64)

which is included into the spectrum of the excitation torque, Sf,eq as in Eq. (17). For the hydraulic PTO, the mean power produced is based on Eq. (56), which in the SL can be obtained as:

PHpto = Hpto | | .

(65)

As mentioned earlier, unlike the previous nonlinear PTO (electric generator with active area factor), the mean power produced by the hydraulic PTO system cannot be obtained in FD due to its source of nonlinearity. 4.2.3. Simulations The OWSC dynamics is characterized by several nonlinearities of zero-mean torque, and with parametric excitation torque (see in Section 3.1). Therefore, instead of comparing the equivalent terms like the PA, this case study focuses on the parametric torque and power produced. However, the analysis of the equivalent terms still applies to the OWSC and OWC devices. Once again, the reliability of the method is verified against TD simulations, where the main parameters of the OWSC are described in Table C.3 and relevant details of the simulation are described in Appendix E. For the case of the hydraulic PTO system, the Coulomb-like torque is assumed to not impede the body from moving, which would modify the response distribution. Fig. 11 illustrates a time series of a TD simulation used to compare against the SL and FD results. The assumption of a Gaussian distribution of the body response occurs because of the characteristics of the excitation torque, which is also assumed to be Gaussian. In this regard, the hypothesis of a Gaussian distribution of the exciting torque must be valid in order to apply the SL. Fig. 12 shows the torque distribution of the TD simulations and compares with the estimations using SL and FD. As observed, the torque distribution using TD is comparable with the estimations using SL based on the theoretical Gaussian Distribution, preserving the

Fig. 13. Comparison between TD, SL and FD results of the PSD of the exciting torque acting on the OWSC (Hs= 3 m and Tp= 12 s). The shaded area is the standard deviation of the TD simulations.

This type of a parametric torque was previously investigated in [15], and can be represented as:

( ) Mexc (t ) = cos( ) Mexc (t ).

(57)

This modulation tends to decrease the torque and mean power absorbed for high amplitudes of motion, which is relevant at resonant conditions and severe sea states. Therefore, the influence of the modulation must be included to obtain realistic estimations of the system dynamics as the device usually operates under these conditions. 4.2.2. Frequency domain and statistical linearization Based on the system dynamics described by Eq. (51), the linearized equation form in the FD of the OSWC device is described by: 8

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

Fig. 16. Comparison between TD, SL, and FD models for the torque modulation, Λ, for a range of Hs (Tp= 12 s). The shaded area is the standard deviation of the TD simulations.

Fig. 14. Comparison between TD, SL, and FD models for the torque modulation, Λ, for a range of Tp (Hs= 2 m). The shaded area is the standard deviation of the TD simulations.

Fig. 17. Comparison between TD, and SL models for the mean power absorbed, P¯Hpto, for a range of Hs (Tp= 12 s). The shaded area is the standard deviation of the TD simulations.

Fig. 15. Comparison between TD, and SL models for the mean power absorbed, P¯Hpto, for a range of Tp (Hs= 2 m). The shaded area is the standard deviation of the TD simulations.

Figs. 16 and 17 show the variation of the torque modulation and the mean power estimation respectively for several significant wave heights and a Tp of 12 s. The torque modulation reduced the torque acting on the OWSC, especially for higher significant wave heights. As a result, the mean power produced by the OWSC sees a diminishing slope in the power absorbed as the wave height is increased.

standard deviation. Fig. 13 shows the PSD of the exciting torque via TD simulations and compares against the SL and FD results. The spectral energy of the excitation torque is conserved over the entire frequency range. The spectral density depends on the system response, which includes the nonlinear terms present in the governing equation. Even though only one condition is shown, the torque distribution and PSD were similar for all conditions simulated. Based on that, the assumption of a Gaussian distribution for the nonlinear excitation torque is reinforced. Once the assumption of Gaussian distribution of the torque is verified, the equivalent torque modulation is assessed. For comparison purposes, the torque modulation in the TD is estimated to compare with the SL results. In this regard, the modulation is simply calculated based on Eq. (59) as:

=

1 T

4.3. Oscillating water column The OWC wave energy device considered in this study is composed of a fixed hollow structure that pierces the sea surface. The pressure variation caused by the waves leads the inner surface to oscillate, which compresses and decompresses the air inside and moves a bidirectional air turbine connected to a generator as illustrated in Fig. 18.

T

cos( ) dt . 0

(66)

Fig. 14 illustrates the variation of the torque modulation for several peak periods considering a Hs of 2 m. The FD model neglects the torque reduction caused by the modulation, while the TD model and SL results have a good agreement. The modulation depends on the body response. Therefore, this effect can be accentuated or attenuated depending on the operating condition. Fig. 15 shows the mean power produced for several peak periods for the same Hs of 2 m. The FD domain calculation is not able to predict the mean power produced by the hydraulic PTO system, therefore, it is not shown. Besides the expected reduction of the power absorbed caused by viscous losses, the torque modulation also presents a considerable influence reducing the mean power.

Fig. 18. Illustration of an Oscillating Water Column wave energy device. 9

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

4.3.1. Governing equations The OWC is characterized by a variable mass system in which the water column is coupled with the air chamber pressure. The water column dynamics is here described based on the plug-flow representation via the generalized coordinate zeta, ζ, as [33]:

( + H) ¨ +

1 Cv 2

+g +

1

p = Fexc ,

system to non zero mean values. Therefore, the SL procedure is performed based on Section 3.2. In this regard, the mean pressure is obtained taking the expectation of Eq. (68):

KT D p¯ NS

¨ + g ¯ + 1 p¯ = 0.

with

KT D 1 p + 2 (L NS ca

Meq =

air

= 0,

The mean values are used to estimate the equivalent coefficients and power produced. The equivalent transfer function based on the SL will have the matrix form: eq (

(68)

(70)

=

Pt , 3D5 N air

(71)

air N

3D5 T

T

fp

0

p (t ) dt . 2D 2

air N

P¯t =

2M

+ i B + K] 1 .

(72)

(73)

with

M=

B=

K=

H 0 , 0 0

0 air

0 , L / ca2

g 1/ w , 0 KT D/ NS

and the natural frequency of the oscillating water column in the absence of the impedance from the air pressure, may be approximated by: n

=

g . H

(Cv | |/2)

0

.

¯/ ca2

air N

3D5

fp

p 2D 2 N air

.

(78)

4.3.3. Simulations The OWC dynamics are characterized by a system with asymmetric nonlinearities and without parametric excitation. The asymmetric nonlinearities lead to nonzero mean force, consequently, offsets of the generalized coordinates are expected (see in Section 3.2), which can change the response of the system. For example, the surface elevation can change the resonance frequency of the water column, and the mean pressure affects the mean power absorbed by the air turbine. In this regard, this study case focuses on the analysis of the mean values of the air chamber pressure and water column displacement, and the mean power absorbed. However, the mean water column displacement can be associated with the equivalent mass and damping terms in this case. The mean values are compared against TD simulations. The main parameters of the OWC used during the simulation is described in Table D.4, relevant details of the simulation are described in Appendix E, and the air turbine coefficients were extracted from [44]. Figs. 19 a, 19 b and 19 c show the mean manometric pressure inside the chamber, mean water column displacement, and mean power absorbed by the air turbine respectively, using TD, SL, and FD models for several peak periods and considering a Hs of 2 m. While the FD model is not capable of estimating the mean values of manometric pressure and water displacement, the SL estimation obtained comparable results with the TD simulations. The higher magnitudes of mean pressure, mean displacement and mean power produced occurred in the resonance regime as expected. Regarding the mean power absorbed, the prediction based on FD results overestimates the mean power produced due to the large losses of power caused by the viscous loads. Even though the mean power absorbed of the OWC is a nonlinear function, the turbine contribution in the system dynamics can be described linearly. Based on that, the mean power absorbed can be estimated using FD model. Figs. 20 a, 20 b and 20 c show the mean manometric pressure inside the chamber, mean water column displacement, and mean power absorbed by the air turbine respectively, using TD, SL, and FD models for several wave heights and considering a Tp of 5.75 s. The mean values of

4.3.2. Frequency domain and statistical linearization Based on the linearized form of Eqs. (67) and (68), the transfer function between the wave force and the water elevation and pressure inside the chamber of the OWC device in FD model is given by:

( )=[

(77)

The mean power estimation for the FD and SL analysis are determined as:

where Ψ and Π are the pressure head and power coefficient respectively. Therefore, the mean power absorbed by the air turbine is given by:

P¯t =

+ Meq) + i (B + Beq) + K] 1 ,

¯ 0 , 0 0

0

with:

p , 2D 2 N air

2 (M

)=[

Beq =

(69)

=

(76)

w

where KT represents the linear relationship between the mass flow rate and pressure fluctuation of the air turbine [43], D denotes the rotor diameter, N is the rotational speed of the turbine (rad/s), S denotes the cross-sectional area of the chamber, ca is the speed of sound in atmospheric conditions, L is the air chamber height at the nominal position, and ρair is the air density. The power absorbed by the air turbine is obtained based on the power coefficient in dimensionless form [44]:

= f p ( ),

(75)

and for the water column taking the expectation of Eq. (67):

(67)

where H is the nominal draft, p is the chamber pressure, Cv is the viscous coefficient. In Eq. (67), the oscillating water column dynamic follows the derivations in [40], where the equation was obtained via two distinct approaches: by the extended Lagrange equations applied to systems of variable mass, and by classical hydrodynamic techniques; see in [41,42]. The governing equation of the air-chamber pressure, p, is expressed as [29]:

)p

1 p = 0, ca2

(74)

The source of nonlinearities described in Eqs. (67) and (68) lead the 10

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

Fig. 19. Comparison between TD, SL and FD models of the OWC mean results for a range of Tp (Hs = 2 m). The shaded area is the standard deviation of the TD simulations.

Fig. 20. Comparison between TD, SL and FD models of the OWC mean results for a range of Hs (Tp = 5.75 s). The shaded area is the standard deviation of the TD simulations.

manometric pressure, water displacement, and mean power absorbed obtained using SL were in good agreement with the nonlinear TD simulations. The mean water displacement for higher Hs was about 5–6% of the OWC draft, which shifts the resonance frequency to lower values. Although, this change in resonance frequency is relatively small, the magnitude of the mass and equivalent mass can become significant for cases where the OWC has a small draft. In addition, the reduction of OWC draft shifts the resonance frequency to higher values. As a result, the nonlinearities dependent on the acceleration and velocities, which is the case of OWC, becomes more important especially under severe wave conditions.

fidelity modeling approaches with low computational cost are desired at preliminary design and optimization stages. In this regard, FD models produce an unreliable prediction of the system dynamics, while TD simulations have a high computational cost. Therefore, either the FD or TD approaches might not be effective for the initial stages of development. On the other hand, these requirements can be easily achieved by using the SL method. However, it is important to notice that the SL technique is limited to stochastic analysis. The SL offers a systematic approach that can be applied to several nonlinearities, WECs, MDOF systems (described in Section 3), including the analysis of wave farms. However, like all models, the validation of the Gaussian assumption must be verified prior to the application of SL which can be achieved using TD simulations or experimental testing. The technique consists of an approximated solution founded on the assumption of a Gaussian distribution. This assumption is assumed to be

5. Discussion Wave energy technology remains under development and high11

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

valid due to the random phase relationship between the frequency components in the excitation force and the non-dominant nonlinear forces. In this paper, the assumption of a Gaussian distribution is verified during the simulations in the case of the body response and the torque with parametric excitation. However, as the nonlinearities increase their contribution in the system dynamics, the Gaussian distribution might not be achieved, and discrepancies between the SL and TD results might occur. This can be demonstrated comparing the simulations over a range of significant wave heights. However, assuming that for high significant wave heights the WEC might operate in the survival mode, the method is assumed to be effective for the operating zone and applied over an extensive range. The SL technique also recovers the spectral energy of the response and parametric excitation over the entire range of frequency. The PSD of the body response is used to estimate the equivalent terms and mean offsets. Consequently, the estimation of the mean offsets and mean power absorbed were comparable with the TD simulations. It is important to notice that the equivalent transfer function in the SL is sensitive to the excitation force. Therefore, besides the wave parameters such as the peak period and significant wave height, the type of spectrum changes the response of the system. The analysis of the equivalent linear terms gives a better understanding of the contribution of nonlinear forces with a simple physical

interpretation, such as the average work done by the restoring force and the average power dissipated or absorbed by the damping force. Although the SL technique is an approximate solution, the analysis of the equivalent terms demonstrated a good agreement with TD simulations for all nonlinearities reported in this work. In general, the SL results converged in approximately 4 to 8 iterations for the PA and OSWC, while the OWC ranged from 5 to 18. Higher numbers of iterations occurred for severe sea states and close to the natural frequency, where the nonlinearities are more relevant. Ten TD simulations containing 300 cycles were performed per sea state in order to make more plausible the hypothesis of ergodicity. The duration of the TD simulations was around 500–700 s per sea state, while for the SL the duration was around 0.04-0.09 s depending on the WEC. Therefore, the SL technique was in general 3 to 4 orders of magnitude faster than the TD simulations. Acknowledgment Leandro S. P. da Silva acknowledges CAPES, the Coordination for the Improvement of Higher Education Personnel, for his Master grant Process number 1754875. C.P. Pesce acknowledges a CNPq Research Grant, nr. 308230/2018-3.

Appendix A. Mathematical expectation In this work, the equivalent linear coefficients of mass, damping and stiffness were calculated numerically, where in the continuous form can be expressed as:

d dq

=

d f (q) dq, dq

where f(q) is the probability density function, see Eq. (24); and Θ denotes the nonlinear function. However, some analytic solutions are available in the literature. Table A.1 presents some analytic solutions of the expectations for the case with a zero mean response, where σq is the standard deviation of the variable q. Table A1 Evaluation of the expectations for a Gaussian variable q [14]. Θ(q)

d dq

sgn(q)

1 2 2

() ()

q|q|

1 8 2

q3

3

q

1

q

2 q

Appendix B. Point absorber The radiation force in the TD is represented here by means of Cummins’ equation [45]:

Fr (t ) =

A

, z z¨ (t )

t

h z (t

)z ( )d ,

(B.1)

where A∞,z denotes the added mass at infinite frequency in heave, and hz is the radiation impulse response function (RIRF), also known as memory function. The hydrostatic force is given by the increment and decrement of the differential volume immersed, which for a sphere submerged around the center and is given by:

Fhs (t ) =

Sw (z ) gdz ,

Fhs = Fhs, l + Fhs, nl =

gR2z

(B.2)

+

(B.3)

g /3z 3,

where Sw is the cross-sectional waterplane area which is a function of the vertical displacement, and R is the radius of the sphere. The nonlinear hydrostatic contribution, Fhs,nl, is used in Eq. (35). The linear mechanical restoring force is given by:

12

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

Table B1 Simulation parameters of the PA used in Section 4.1. Property M R Ks Cpto Kes les Kls dls lls lp ls CD

Fs (t ) =

Value

Unit 5

6 × 10 5 100 25 250 1.5 100 1 1 3.5 3 0.5

[kg] [m] [kN/m] [kN s/m] [kN/m] [m] [kN/m] [m] [m] [m] [m] [ ]

(B.4)

Ks z (t ),

where Ks denotes the stiffness of the vertical spring. The excitation force is represented by a summation of regular components expressed as [18]: N

Fexc (t ) =

^ ( Re H F ,z

^

k) a ( k ) e

i kt

, (B.5)

k=1

^ is the complex force response amplitude operator in the z direction, and ^ is the complex wave amplitude of a random angle distribution, where H F ,z a ℜ is the real operator, and the magnitude of the wave amplitude is given by a wave spectrum. Appendix C. Oscillating wave surge converter The radiation moment in the TD can be written based on Cummins’ equation:

Mr (t ) =

A

,

t

¨ (t )

h (t

) ( )d .

(C.1)

The linear restoring torque is given by: (C.2)

Ms (t ) = Ks ,

where Ks represents the linear rotational stiffness from the PTO and hydrostatic restoring torque. The excitation torque in the TD is given by a summation of regular components as: N

^ ( Re H M,

Mexc (t ) =

^

k) a ( k) e

i kt

, (C.3)

k=1

^ where H M, is the complex force response amplitude operator in the angular direction θ. Table C1 Simulation parameters of the OWSC used in Section 4.2. Property

Value

Unit

I V M w L d Cv Ks Hpto

5.3 × 106 442 1.8 × 105 26 9 4.5 1.2 21 × 103 1 × 103

[kg m2] [m3] [kg] [m] [m] [m] [ ] [kN m/rad] [kN m]

Appendix D. Oscillating water column The excitation force in the TD is given by: N

Fexc (t ) =

^ ( Re H F,

^

k) a ( k ) e

i kt

, (D.1)

k=1

13

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al.

Table D1 Simulation parameters of the OWC used in Section 4.3. Property

Value

Unit

ρw ρair ca g Cv KT H D L h S N

1025 1.25 344 9.81 0.5 0.28 8.21 1.75 9 50 10 30

[kg/m3] [kg/m3] [m/s] [m/s2] [ ] [ ] [m] [m] [m] [m] [m2] [rad/s]

with:

^ ( ) = 2 g cosh(k (h H )) , H F, cosh(kh)

(D.2)

where h denotes the water depth, and k is the wavenumber which is a function of the frequency. Appendix E. Simulations The hydrodynamic and excitation force coefficients of the PA and OWSC were obtained using the open-source software NEMOH [46]. Based on the hydrodynamic coefficients, the radiation impulse response function were replaced by an equivalent state space representation of order 4 via system identification described in [47]. The irregular waves simulated were characterized by the JONSWAP spectrum ( = 3.3), in which the wave spectrum was discretized into 300 frequencies from 0.25 to 2.5 rad/s. As the TD simulations depend on the phase distribution of the excitation force components, ten simulations containing a minimum of 300 cycles were calculated using different randomly generated phase angles to make more plausible the hypothesis of ergodicity. Based on the TD simulations, the PSD were calculated using the pwelch function in MATLAB. As ten simulations per condition were performed, the TD results are presented in terms of the mean value (solid line) and standard deviation (shaded area).

References [1] A. Clément, P. McCullen, A. Falcão, A. Fiorentino, F. Gardner, K. Hammarlund, G. Lemonis, T. Lewis, K. Nielsen, S. Petroncini, et al., Wave energy in europe: current status and perspectives, Renew. Sustain. Energy Rev. 6 (5) (2002) 405–431. [2] F.O. Antonio, Wave energy utilization: a review of the technologies, Renew. Sustain. Energy Rev. 14 (3) (2010) 899–918. [3] J. Falnes, A review of wave-energy extraction, Marine Struct. 20 (4) (2007) 185–201. [4] I. López, J. Andreu, S. Ceballos, I.M. de Alegría, I. Kortabarria, Review of wave energy technologies and the necessary power-equipment, Renew. Sustain. Energy Rev. 27 (2013) 413–434. [5] M. Karimirad, Offshore Energy Structures: for Wind Power, Wave Energy and Hybrid Marine Platforms, Springer, 2014. [6] J. Falnes, Ocean Waves and Oscillating Systems: Linear Interactions Including Wave-Energy Extraction, Cambridge University Press, 2002. [7] G. Giorgi, J.V. Ringwood, Comparing nonlinear hydrodynamic forces in heaving point absorbers and oscillating wave surge converters, J. Ocean Eng. Marine Energy 4 (1) (2018) 25–35. [8] A. Babarit, J. Hals, M. Muliawan, A. Kurniawan, T. Moan, J. Krokstad, The numWEC project. numerical estimation of energy delivery from a selection of wave energy converters – final report, Starkraft, Centrale Nantes, NTNU (2011) 1–317. [9] M. Penalba, N. Sell, A. Hillis, J. Ringwood, Validating a wave-to-wire model for a wave energy converter â- part i: the hydraulic transmission system, Energies 10 (7) (2017) 977. [10] M. Penalba, J.-A. Cortajarena, J. Ringwood, Validating a wave-to-wire model for a wave energy converter â- part ii: the electrical system, Energies 10 (7) (2017) 1002. [11] N. Sergiienko, A. Rafiee, B. Cazzolato, B. Ding, M. Arjomandi, Feasibility study of the three-tether axisymmetric wave energy converter, Ocean Eng. 150 (2018) 221–233. [12] A. Scialò, G. Malara, F. Arena, Geometrical optimization of u-oscillating water columns in random waves, ASME 2019 38th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers, 2019. [13] M. Göteman, J. Engström, M. Eriksson, J. Isberg, M. Leijon, Methods of reducing power fluctuations in wave energy parks, J. Renew. Sustain. Energy 6 (4) (2014) 043103. [14] J.B. Roberts, P.D. Spanos, Random Vibration and Statistical Linearization, Courier Corporation, 2003. [15] M. Folley, T. Whittaker, Spectral modelling of wave energy converters, Coastal Eng. 57 (10) (2010) 892–897. [16] M. Folley, T. Whittaker, Preliminary cross-validation of wave energy converter

[17] [18] [19] [20] [21]

[22] [23] [24] [25] [26] [27]

[28]

[29]

14

array interactions, ASME 2013 32nd International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers, 2013. V008T09A055–V008T09A055 M. Folley, T. Whittaker, Validating a spectral-domain model of an owc using physical model data, Int. J. Marine Energy 2 (2013) 1–11. M. Folley, Numerical Modelling of Wave Energy Converters: State-of-the-Art Techniques for Single Devices and Arrays, Academic Press, 2016. A. Mérigaud, J.V. Ringwood, A nonlinear frequency-domain approach for numerical simulation of wave energy converters, IEEE Trans. Sustain. Energy 9 (1) (2017) 86–94. A. Mérigaud, J.V. Ringwood, Power production assessment for wave energy converters: overcoming the perils of the power matrix, Proc. Inst. Mech. Eng. Part M 232 (1) (2018) 50–70. R. Novo, G. Bracco, S.A. Sirigu, G. Mattiazzo, A. Mérigaud, J.V. Ringwood, Non linear simulation of a wave energy converter with multiple degrees of freedom using a harmonic balance method, ASME 2018 37th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers Digital Collection, 2018. A. Mérigaud, A harmonic balance framework for the numerical simulation of nonlinear wave energy converter models in random seas, National University of Ireland Maynooth, 2018 Ph.D. thesis. T.S. Atalik, S. Utku, Stochastic linearization of multi-degree-of-freedom non-linear systems, Earthquake Eng. Struct. Dyn. 4 (4) (1976) 411–420. N. Tom, M. Lawson, Y.-H. Yu, A. Wright, Spectral modeling of an oscillating surge wave energy converter with control surfaces, Appl. Ocean Res. 56 (2016) 143–156. P. Gunawardane, M. Folley, S. Sanjaya, Spectral-domain modelling of the non-linear hydrostatic stiffness of a heaving-sphere wave energy converter, The 28th International Symposium on Transport Phenomena, (2017). P.D. Spanos, F. Arena, A. Richichi, G. Malara, Efficient dynamic analysis of a nonlinear wave energy harvester model, J. Offshore Mech. Arctic Eng. 138 (4) (2016) 041901. L.S.P. Silva, H.M. Morishita, C.P. Pesce, R.T. Gonçalves, Nonlinear analysis of a heaving point absorber in frequency domain via statistical linearization, ASME 2019 38th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers, 2019. P.D. Spanos, F.M. Strati, G. Malara, F. Arena, Stochastic dynamic analysis of u-owc wave energy converters, ASME 2017 36th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers, 2017, pp. 1–8. L.S.P. Silva, C.P. Pesce, H.M. Morishita, R.T. Gonçalves, Nonlinear analysis of an oscillating water column wave energy device in frequency domain via statistical linearization, ASME 2019 38th International Conference on Ocean, Offshore and

Applied Ocean Research 95 (2020) 102023

L.S.P. Silva, et al. Arctic Engineering, American Society of Mechanical Engineers, 2019. [30] I. Elishakoff, C. Soize, Nondeterministic Mechanics, 539 Springer Science & Business Media, 2013. [31] M.K. Ochi, Ocean Waves: The Stochastic Approach, 6 Cambridge University Press, 2005. [32] F. Fedele, J. Brennan, S.P. De León, J. Dudley, F. Dias, Real world ocean rogue waves explained without the modulational instability, Sci. Rep. 6 (2016) 27715. [33] L.S.P. da Silva, Nonlinear Stochastic Analysis of Wave Energy Converters via Statistical Linearization, University of São Paulo, Brazil, 2019 Master’s thesis. [34] J. Morison, J. Johnson, S. Schaaf, et al., The force exerted by surface waves on piles, J. Petroleum Technol. 2 (05) (1950) 149–154. [35] M. Eriksson, R. Waters, O. Svensson, J. Isberg, M. Leijon, Wave power absorption: experiments in open sea and simulation, J. Appl. Phys. 102 (8) (2007) 84910. [36] B. Ding, N. Sergiienko, F. Meng, B. Cazzolato, P. Hardy, M. Arjomandi, The application of modal analysis to the design of multi-mode point absorber wave energy converters, Ocean Eng. 171 (2019) 603–618. [37] A. Babarit, J. Hals, M.J. Muliawan, A. Kurniawan, T. Moan, J. Krokstad, Numerical benchmarking study of a selection of wave energy converters, Renew. Energy 41 (2012) 44–63. [38] G. Bacelli, J.V. Ringwood, Nonlinear optimal wave energy converter control with application to a flap-type device, IFAC Proc. Vol. 47 (3) (2014) 7696–7701. [39] A. McCabe, G.A. Aggidis, T. Stallard, A time-varying parameter model of a body

oscillating in pitch, Appl. Ocean Res. 28 (6) (2006) 359–370. [40] C.P. Pesce, L. Casetta, Systems with mass explicitly dependent on position, Dynamics of Mechanical Systems with Variable Mass, Springer, 2014, pp. 51–106. [41] C. Pesce, The application of lagrange equations to mechanical systems with mass explicitly dependent on position, J. Appl. Mech. 70 (5) (2003) 751–756. [42] C.P. Pesce, E.A. Tannuri, L. Casetta, The lagrange equations for systems with mass varying explicitly with position: some applications to offshore engineering, J. Brazil. Soc. Mech. Sci. Eng. 28 (4) (2006) 496–504. [43] A. de O Falcão, R. Rodrigues, Stochastic modelling of OWC wave power plant performance, Appl. Ocean Res. 24 (2) (2002) 59–71. [44] J. Henriques, J. Portillo, L. Gato, R. Gomes, D. Ferreira, A. Falcão, Design of oscillating-water-column wave energy converters with an application to self-powered sensor buoys, Energy 112 (2016) 852–867. [45] W. Cummins, The Impulse Response Function and Ship Motions, Technical Report, David Taylor Model Basin Washington DC, 1962. [46] A. Babarit, G. Delhommeau, Theoretical and numerical aspects of the open source bem solver nemoh, 11th European Wave and Tidal Energy Conference (EWTEC2015), (2015). [47] T. Perez, T.I. Fossen, A Matlab toolbox for parametric identification of radiationforce models of ships and offshore structures, Modeling, Identification and Control 30 (2019) 1–15.

15