Modelling of linear and non-linear two-body wave energy converters under regular and irregular wave conditions

Renewable Energy 147 (2020) 487e501

Contents lists available at ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

Modelling of linear and non-linear two-body wave energy converters under regular and irregular wave conditions Xueyu Ji a, Elie Al Shami a, Jason Monty b, Xu Wang a, * a b

School of Engineering, Royal Melbourne Institute of Technology University, Melbourne, Victoria, Australia Department of Mechanical Engineering, The University of Melbourne, Parkville, Victoria, Australia

a r t i c l e i n f o

a b s t r a c t

Article history: Received 12 July 2018 Received in revised form 23 January 2019 Accepted 3 September 2019 Available online 6 September 2019

This paper studies a two-body wave energy converter oscillating in heave with a floating body of variable geometry connected to a submerged body, which is designed for the specific ocean wave condition with wave excitation frequencies ranging from 0.08 Hz to 0.12 Hz. The study focuses on the parameter and non-linear effects of the converter. Taguchi method has been applied to investigate the system model parameters' influences on the maximum average power output. ANSYS AQWA is employed to obtain the hydrodynamic parameters for calculating the output power. Both linear and non-linear dynamic models of the two-body wave energy converter will be analysed and simulated in both the time and frequency domains. The new scientific contribution is identifying the best parameter combination of the two-body wave energy converter for irregular wave excitations for the Australian ocean wave conditions, evaluating and comparing the output power of both linear and non-linear dynamic models of the two-body wave energy converter. Another contribution is having established a non-linear model for those buoys of non-uniform shape in vertical direction which has largely improved the simulation accuracy in comparison with a linear model where the buoy hydrostatic forces were simulated using a linear approximation model. The most important results are having found that the power-take-off stiffness coefficient and submerged body geometry are the most important parameters to change for the largest output power. It is found that the significant wave height has very little influence on the efficiency whilst the variations of the peak period of the irregular wave can change the efficiency largely. The non-linear hydrostatic force is found to have considerable effects on the buoy's dynamic performance. The best combination of the system parameters for the harvesting performance has been identified using the Taguchi method, which produces a peak efficiency of 51% in irregular waves. This manuscript aims to disclose the parameter effect of a two body WEC on its performance using Taguchi method and to provide a guide for industry and research community for designing a good two body WEC. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Two-body wave energy converter Point absorber Taguchi method Optimisation Non-linear

1. Introduction Wave energy is a prominently renewable energy that is both periodical and predictable. Because of the energy crisis and greenhouse effects, wave energy has been an interesting topic in last few decades. The power contained in wave propagation along the ocean's surface can reach a few Terra-Watts. Average power flow intensity is typically 2e3 kW/m2 perpendicular to direction of wave propagation, which is denser than the power available from solar (0.1e0.3 kW/m2) or wind (0.5 kW/m2) energy [1]. Many

* Corresponding author. E-mail address: [email protected] (X. Wang). https://doi.org/10.1016/j.renene.2019.09.010 0960-1481/© 2019 Elsevier Ltd. All rights reserved.

proposed devices such as surface attenuator [2], water oscillating columns [3], water overtopping devices [4], and so on are being developed to efficiently harvest the energy contained in ocean waves. However, none of them has really stood out as a definitive choice of a wave energy-harvesting device to be developed commercially and connected to power grids around the world. Apart from that, the point absorber, which makes use of the resonant effect when coupled with a linear generator, is considered as an attractive option to extract energy from waves. It presents low complexity and maintenance and potentially high efficiency. Also, it can harvest energy from all wave directions, which is suitable for high ocean wave energy density places like Australia. However, due to the low frequency of the energetic ocean waves in the regions such as Australia (0.08 Hze0.12 Hz near the south) [5], huge mass

488

X. Ji et al. / Renewable Energy 147 (2020) 487e501

and volume of a buoy are typically needed in order to push the resonant frequency down and absorb as much energy as possible, which is impractical and costly to realize. In this work, a two-body wave energy converter with a submerged body and a floating buoy is adopted for the Australian sea state as it has two degrees of freedom (DOF) which can increase the equivalent mass of the system greatly due to the increase of the hydro-dynamically added mass. This increase in the mass results in a reduction of the natural frequency, thus resulting in the ability to operate in the ocean regions of low wave frequencies. Wavebob [6] and Powerbuoy [7] are currently available two body wave energy converters and shown in Fig. 1. The two-body wave energy converter's (WEC) numerical model has been discussed by Ref. [8] and is shown in Fig. 2, but this model did not take into considerations the non-linear effects which may exist in the buoy's hydrostatic stiffness (non-uniform shape). The influences of the floating buoy geometry, diameter, draft, damping coefficient, wave height and wave period on the power output have been researched by Ref. [9], but the viscous damping of the submerged bodies was neglected. Even though a non-linear spherical buoy was proposed by Ref. [10], its non-linear hydrostatic stiffness was not elaborated properly. Using a two-body system coupled with a linear driven generator has been proposed by Ref. [11], but the main issue is neglecting the viscous damping in the numerical model based on the assumption that the effects of the viscous damping on the power are very small for a submerged oscillating sphere. Based on the conclusion of [12], the phase difference between the induced voltage and current leads to the variation of the Power-Take-OFF (PTO) forces: one is linked with the velocity which can be modelled as a damping force and the other one is related to displacement and is modelled as a spring stiffness force. In this paper, the proposed model includes the viscous damping and non-linear hydrostatic force. The PTO is modelled as a spring damper device. Even though there are many other types of PTO, like hydraulic actuators [13], rack pinion systems [14,15] and hydraulic turbines [16], the linear generator [11] is one of the outstanding options with high efficiency, simplicity and low maintenance cost. Therefore, a generic linear generator is utilized in this research. The operational conditions in the ocean are complex. The influence of each parameter is somehow dependent on that of the other parameters. For example, the sea state determines the optimal size and geometry of the buoy, which alters the hydrodynamic coefficient. The hydrodynamic coefficient affects the control algorithm, thus affects the PTO coefficients. Therefore, Taguchi method which has been used to conduct parameter sensitivity studies and optimizations in the chemistry process industry [17] can be applied to evaluate the influence of the different parameters, as it can analyse the response of the system to the variation of each parameter and detect the best combination of parameters to produce the highest power. Additionally, the real sea state cannot only be regarded as regular waves but can have the irregular waves with random elevation. The performance of WEC,

Fig. 2. Simulation model of Wavebob and Powerbuoy [8].

like the efficiency, is one of the most fundamental criterions to judge if such device is worthwhile being commercially implemented in real world. The influences of the significant wave height and peak period on the efficiency of the two-body wave energy converter have not been researched comprehensively [11]. Thus, these two parameters' effects will be further studied for the twobody wave energy converter under irregular waves. To study the influence of non-linear hydrostatic force of the non-uniform shape of buoy, hydrostatic forces of different non-uniform buoy will be investigated in this research where the hydrostatic force is modelled as a linear and non-linear one. Moreover, the influences of the non-linear and linear hydrostatic forces on the time domain buoy displacement and velocity will be analysed. A flow chart for illustration of the whole simulation procedure has been added in Fig 13 in Appendix to showcase how the data come out and how the calculations are done either analytically or via simulation. 2. Numerical simulation model Fig. 3 (a) shows the outline of a two-body wave energy converter for this research and Fig. 3 (b) demonstrates the details of the PTO system where a cylinder with strong magnet arrays such as the Halbach array inside is connected to the submerged body and a coil bar is connected to the buoy, this is for demonstration purposes only, the PTO system is simulated as a spring damper system. The regular waves are assumed to have a wave height of 2 m for all the cases. For the two-body wave energy converter, the submerged body is set to have neutral buoyancy when placed in a stationary wave. In this section, both linear and non-linear models will be introduced. The added mass, radiation damping, viscous damping and excitation forces will be presented in the dynamic equations. A buoy with non-linear hydrostatic force is shown in Fig. 4 (a). For a linear model, the buoy is set to have a constant cross section with a zero-inclination angle and its shape profile is shown in Fig. 4 (b). The equations of the 2 DOF mass spring dashpot system can be developed based on the model in Fig. 2. The wave excitation forces exerted on the buoy and submerged body are fe1 and fe2 respectively. Radiation damping coefficients of the buoy and submerged body are b11 and b22 respectively, but b22 is usually ignored since it is very small compared to b11 [8]. The added masses for the buoy and submerged body are a11 and a22 respectively. Once the geometry features of the buoy are set, the hydrostatic force of the buoy fb and the hydrostatic stiffness ks can be calculated by

fb ¼ rg pR2 x1 Fig. 1. Wavebob and powerbuoy [8].

(1)

X. Ji et al. / Renewable Energy 147 (2020) 487e501

489

Fig. 3. (a) Illustration of a two-body wave energy converter device; (b) Illustration of a power take-off system.

Fig. 4. (a) Non-linear floating buoy; (b) Linear two body wave energy converter; (c) Non-linear two body wave energy converter.

8 < ks ¼ rgpR2

(2)

where r is the density of the ocean water; g is the gravity acceleration of 9.81 m/s2; p is circle rate which is 3.1415928; R is the radius of the buoy; x1 is the instantaneous displacement of the buoy. The draft of the buoy is set to be half its height. Based on the submerged volume of the buoy and submerged body in water, Archimedes’ principle can be utilized to calculate the physical masses of these two bodies. The dynamic heaving motion equations of the linear model under the excitation of regular waves have been implemented for the Fourier transformation and analysed in the frequency domain and are given by

2

3 m þ a ð u Þ a ð u Þ 11 12 5  u2 4 1 m1 þ a22 ðuÞ a21 ðuÞ : 2 3 b þ b ð u Þ þ b b ð u Þ pto 11 12 vis1 5 þ i u4 bpto þ bvis2 þ b22 ðuÞ b21 ðuÞ 2 39 2 3 = k þ k 0 s 5 , 4 X1 5 þ 4 pto 0 kpto ; X2 3 2 f ð uÞ 5 ¼ 4 e1 fe2 ðuÞ

(3)

490

X. Ji et al. / Renewable Energy 147 (2020) 487e501

where X1, X2 are the maximum displacements amplitude of the buoy and submerged body respectively, m1, m2 are the masses of the buoy and submerged body respectively; u is the excitation frequency; a11(u) and a22(u) are the added masses of the buoy and submerged body; kpto and bpto are the stiffness and damping coefficients of the power take off unit; b11(u) is the radiation damping coefficient of the buoy; b12(u), b22(u) and b22(u) are the radiation damping and a12(u) and a21(u) are the coupling added masses which typically ignored if buoy and submerged body are far from each other based on [8,18,19]; bvis1 and bvis2 are the linearized viscous damping coefficients but only for the submerged body will be accounted [18,20]; fe1(u) and fe2(u) are the wave excitation forces exerted on the buoy and submerged body respectively. A non-linear simulation numerical model is presented which has been proven to be a valid tool to broaden the harvesting frequency bandwidth [22]. It is found that if a cylinder has a certain inclination angle, it will lead to the non-linear hydrostatic force. In this research, a non-linear model seen in Fig. 4 (c) induced by the non-linear hydrostatic force is proposed and the governing equations of the non-linear system for the oscillation in the heaving direction are given by

ðm1 þ a11 ð∞ÞÞx€1 þ ðt þ

rg p 3

! 3Rx21 x31 2 þ 3R x1 þ tanð90  qÞ tan2 ð90  qÞ

hðt  tÞ,x_1 dt þ kpto ðx1  x2 Þ þ bpto ðx_1  x_2 Þ

¼ fe1 (4) ðm2 þ a22 ð∞ÞÞx€2 þ kpto ðx2  x1 Þ þ bpto ðx_2  x_1 Þ þ bvis x_2 ¼ fe2 (5)

h@t A ¼

2

∞ ð

p

1

0

1

bðuA , cos@ut Adu

4. MATLAB simulation To execute the Matlab simulation in the frequency domain, the viscous damping force is linearized for simplification and given by

Fd ¼

1  r  Cd  As  V 2 2

Fd ¼ bvis  V bvis ¼

∞

01

2 m, height of 2 m and draft of 1 m in the frequency domain will be presented in Fig. 5 (a)-(c). However, due to the limitations of the ANSYS AQWA, the viscous damping coefficient cannot be solved using the software. In order to calculate the viscous damping forces, the software based on the finite volume method such as ANSYS CFX where fluid is set as viscous is required.

(6)

0

where a11(∞) and a22(∞) are the added masses of the buoy and submerged body; q is the inclination angle; x1, ẋ1 and ẍ1 are the instantaneous displacement, velocity and acceleration of the buoy respectively; x2, ẋ2 and ẍ2 are the instantaneous displacement, velocity and acceleratrion of the submerged body respectively; h(t) is the impulse response for buoy radiation damping; t is the integral variable; t is the time variable; and b(u) is the radiation damping frequency impulse response.

3. ANSYS AQWA All these parameters such as the added mass, radiation damping coefficient, and excitation force can be obtained via simulations done in ANSYS AQWA. The geometries of the buoy and submerged bodies are modelled using a CAD software and the models are deployed in ANSYS AQWA where the hydrodynamics are simulated under the excitation of regular waves. ANSYS AQWA is a type of boundary element method software based on the linear potential wave theory to simulate the wave interaction with the WEC. Both the radiation and diffraction theories are utilized in this software. The potential flow theory is used with the assumption that the velocity potential is irrotational and incompressible. The viscous damping coefficient data can be estimated from experimental data in the literature. The added mass and radiation damping coefficient as well as the excitation force for a cylindrical buoy with a radius of

(7)

(8)

1  r  Cd  As  V 2

(9)

where r is the density of the ocean water of 1025 kg/m3, Cd is the drag coefficient set as 0.16 for a sphere and 1 for a cylinder (typically the submerged body works in very high Reynolds number more than 106) [23], and As is the characteristic area which is 28.274 m2 for all the submerged bodies in this research. In the linearized model, V in Equation (9) is set to be 0.6283 m/s which is the maximum oscillating velocity in the heaving direction for a wave height of 2 m. Since the radius is the same in the heaving direction, the characteristic radius is 3 m. From Equation (9), it can be found that the maximum linearized viscous damping coefficients bvis for the submerged body of both spherical and cylindrical shapes are 1457 kg/s and 9107 kg/s respectively. All the parameters simulated in ANSYS AQWA are imported into the Matlab codes. For the linear model, the average power output in the frequency range of 0.08e0.12 Hz can be simulated and obtained via the Fourier transform by Ref. [8].

Pave ¼

1 T

ðT Psim ðtÞdt ¼

bpto u2 ,jX1  X2 j 2

(10)

0

where Pave is the average power output in the frequency range; the Psim(t) is the simultaneous power generated by the PTO, u is angular frequency; T is the period of the wave excitation. Since the Fourier transform cannot be implemented for nonlinear model and therefore the non-linear model simulation has to be conducted in time domain. For this research, a Matlab Simulink code (with a fixed time step of 0.01 s through the RungeKutta solver) is programmed to conduct the time domain simulation to calculate the average power output for a specific frequency. For the time domain analysis, 9 frequency samples with an equal frequency interval of 0.0025 Hz are selected for the frequency range from 0.08 Hz to 0.12 Hz where the corresponding hydrodynamic parameters calculated using ANSYS AQWA at each of individual frequencies will be used as inputs of the Matlab Simulink code to calculate the average power output at that frequency. Then the seventeen average power outputs at the seventeen frequencies can be utilized to produce an average power vs. frequency curve in the frequency domain to find the maximum average power output in the desired frequency range for the non-linear model.

X. Ji et al. / Renewable Energy 147 (2020) 487e501

491

Fig. 5. (a) Added mass in the frequency domain a11(u); (b) Radiation damping coefficient in the frequency domain b11(u); (c) Excitation force in the frequency domain fe1(u).

1 Pave ¼ T

ðT

1 Psim ðtÞdt ¼ T

0

ðT

2

bpto ,ðx_1  x_2 Þ dt

(11)

0

5. Irregular wave simulation

hðx; tÞ ¼ z

N X i¼1

Ai cosðki x  ui t þ εi Þ

i¼1

Ai cosðki x  ui t þ εi Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2,Sðui Þ,Du

(13)

u2i ¼ g,ki ,tanhðki , dh Þ

(14)

Ai ¼

Irregular waves are more realistic for simulating real ocean wave conditions. Irregular waves can be envisaged as the superposition of a series of regular waves with random wave amplitudes, phases and frequencies and the irregular wave elevation h(x,t) is given by Ref. [24]. ∞ X

where εi is the random phase of the irregular waves and the x is the position of the WEC and selected as 0 at the original point. Ai is the wave amplitude at a frequency ui and the wave number ki can be given by,

(12)

where S(ui) is the value of the Jonswap spectrum function and Du is the frequency bandwidth (0.01 rad/s in this case). dh is the ocean depth. In order to calculate the wave amplitude Ai, the Jonswap spectrum function is selected to describe the wave excitation spectrum and given by,

492

X. Ji et al. / Renewable Energy 147 (2020) 487e501

(  u 4  exp ag2 p ,g SðuÞ ¼ 5 exp  1:25

0:5

u u s ¼ 0:07 up > u; s ¼ 0:09 up  u a ¼ 5:061

u 4 p

2p

H2s ½1  0:287 log g

2 )

 uup sup

(15)

(16)

where up is the peak frequency, Hs is the significant wave height and g is a constant which is typically selected as 3.3 [25]; u is the wave excitation frequency. For regular waves, the wave excitation force is given by,

fen ðtÞ ¼ fen ðuÞ  A cosðkx  ut þ εÞ

(17)

Where fen ðuÞ is the excitation force coefficient in the frequency domain seen in Fig. 5 (c) which is obtained from ANSYS AQWA; k is the wave number; A is the excitation wave displacement amplitude; ε is the phase. For the excitation of the irregular waves, the superposition of the generated regular waves can be applied [26] and thereby the irregular wave excitation force can be calculated by

fen ðtÞz

N X

fen ðui Þ  Ai cosðki x  ui t þ εi Þ

(18)

i¼1

where Ai is calculated from Equations 13 and 14. In order to uncover the energy flux for the irregular wave in a propagating direction, the group velocity cg and phase velocity cp are utilized. Based on dispersion relation, the definitions of cg and cp in deep water [27] are given by

cp ¼

rffiffiffi g g ¼ k u

cg ¼

vu cp ¼ ¼ vk 2

(19) rffiffiffiffiffiffi g 4k

(20)

where k is the wave number. In this study, the efficiency also known as capture width ratio (CWR) is involved to evaluate the performance of WEC in irregular wave. The energy flux and efficiency can be calculated by utilizing Equations (11), (21) and (22).

pwave ¼ rg

∞ ð

cg ðuÞ,SðuÞdu ¼

rg2 H 2s Te 64p

0

(21)

pave h¼ pwave ,D (22)where Te is the energy period which is typically 0.86 of the peak period Tp and h is the efficiency and D is the surface diameter of the buoy. 6. Taguchi method and the variables After the numerical simulation models are established, Taguchi method is applied to conduct a parameter study and optimisation of the WEC. Taguchi method is a practical method in quality engineering to maximize the robustness of a product and recognize the most essential controllable design inputs in new product development [17]. In this study, seven controllable parameters are selected for the WEC and listed in Table 1. The eight different two-

body WEC groups are formed and listed in Table 2. Each of them has seven identical controllable parameters as listed in Table 1 but the values of these parameters in each of WEC groups are different (either Level 1 or Level 2). Taguchi method can be utilized to identify the parameters’ influences of the WEC design on the power output. In this research, no control algorithm was implemented, as the goal of this work is to study the effect of the mechanical and hydrodynamic parameters on the system performance at low frequencies. In Taguchi method, the first parameter is the buoy's inclination angle which generates the linear and non-linear influences on the power output. The second and seventh parameters are the draft and diameter of the buoy which are mainly set for studying the influence of the buoy's dimension. The submerged body geometry affects the viscous damping coefficient and added mass. A shape with high added mass increases the system's inertia reducing the resonance frequency but also tends to increase the system's viscous damping, which lowers the captured power. The third parameter is the submerged body geometry which is set to find the effect of the spherical and cylindrical shapes of the submerged body on the power output through the difference of the viscous damping coefficient and hydrodynamic forces acting on the oscillating bodies. The large added mass of the submerged body can reduce the natural frequency and increase the harvested power in low frequency ranges but the large viscous damping coefficient of a submerged cylinder would also reduce the captured power. To maintain a consistent size in the heaving direction, radiuses of the submerged body of cylindrical and spherical shapes are set to be 3 m. According to Ref. [11], the submerged depth can affect the final power output outcome. If the submerged body is placed deep enough, the hydrodynamic interactions between the two bodies can be reduced. The PTO parameters like the damping and stiffness coefficients can have a large effect on the conversion frequency bandwidth [21] and natural resonant frequency. The PTO damping and stiffness coefficient are presented by Parameters five and six in Table 1. To better study the sensitivity of power output under the considerable variation of damping and stiffness coefficients, 50000 kg/s and 100000 kg/s are set as the extreme values of the PTO damping coefficients for bpto and 50000 N/m and 100000 N/m are set as the extreme values of the PTO stiffness coefficients for kpto in the simulation tests. The outcomes of the simulations at these nine frequencies for Groups 1e4 are presented in Fig. 6 (a). It is seen from Fig. 6(a) that Groups 1 and 4 have the output increasing up to a peak value then decreasing with the wave excitation frequency, while Groups 2 and 3 have the output power increasing with the wave excitation frequency. Group 4 have the largest peak output power value. In Fig. 6 (b)e(e), the x axis of the graph implies the selected ocean wave frequency simulated in time domain with Simulink and the y axis indicates the average power output at the corresponding frequency. The fitting average power curves of the Taguchi simulation groups using the non-linear model are shown in Fig. 6 (b)e(e). It is seen from Fig. 6 (b) that Group 5 has the

Table 1 Taguchi settings. Parameter Variable Setting

Level.1

Level.2

1.Inclination angle of floating cylinder 2.Cylinder draft 3.Submerged body geometry 4.Submerged depth 5.bpto 6.kpto 7.Diameter of the floater

0 deg 1m Cylinder 20 m 100000 kg/s 100000 N/m 4m

25 deg 2m Sphere 30 m 50000 kg/s 50000 N/m 6m

X. Ji et al. / Renewable Energy 147 (2020) 487e501

493

Table 2 Taguchi method data sheet. Input variable

Angle

Draft

Submerged geometry

Submerged Depth

bpto

kpto

Diameter of floater

Output

Settings Group 1 Group 2 Group 3 Group 4 Group 5 Group 6 Group 7 Group 8 Average 1 (kW) Average 2 (kW) Effect (kW)

1 1 1 1 1 2 2 2 2 79.55 78.69 0.86

2 1 1 2 2 1 1 2 2 89.8 68.45 21.35

3 1 1 2 2 2 2 1 1 45.76 112.48 66.72

4 1 2 1 2 1 2 1 2 90.85 67.4 23.45

5 1 2 1 2 2 1 2 1 58.88 99.36 40.48

6 1 2 2 1 1 2 2 1 124.5 33.75 90.75

7 1 2 2 1 2 1 1 2 62.31 95.93 33.62

Average (kW) 76.92 36.81 65.15 139.33 216.87 28.58 4.44 64.88

output increasing up to a peak value of 216.87 kW then decreasing with the wave excitation frequency. The peak output power of Group 5 is larger than that of Group 4. It is seen from Fig.6 (c) and Fig.6 (d) that Groups 6 and 7 have the output power decreasing with the wave excitation frequency. It is seen from Fig. 6 (d) that the output power frequency curve of Group 8 initially increases with frequency increments and then it drops down with larger frequency. The maximum power output values of the fitting average power frequency spectrum curves of the non-linear model and the peak average power output of the linear model in the frequency domain are listed in Table 2. From Table 1, every parameter variable is set as either the Level 1 or Level 2 values. For each group in Table 2, it is comprised of seven parameters. Number 1 or 2 in a group corresponds to the parameter variable setting of Level 1 or Level 2. Apart from that, it can be noticed for each input parameter variable (in a column), the numbers of Level 1 and Level 2 are identical and equal to 4. The average power output of a specific level of an input variable (Average 1 or Average 2) is calculated by averaging the output power of all the simulation groups in the same column of the input parameter variable at the same specific level, which is given by 4 P

APn ¼ i¼1

Pavg;i 4

(23)

where Pavg,i indicates one maximum time-average power output of a simulation group for a specific level of a specific input variable. AP1 and AP2 are the average of the maximum averaged power output of all the groups for Level 1 and Level 2 of a specific input variable. The effect of the variable in the Taguchi method can be found through subtraction of the average maximum power value of Level 2 from that of Level 1 (Average 2 minus Average 1) using the following equation,

Effect ¼ AP2  AP1

(24)

Typically, the large difference means that parameter has a large effect on the average power output. The absolute value of the effect is presented in Fig. 7. It is seen from Fig. 7 that the PTO stiffness has the largest effect on the output power. Decreasing the PTO stiffness will decrease the output power. The submerged body geometry has the second largest effect and the PTO damping coefficient has the third largest effect on the output power. Changing the geometry of the submerged body from the cylinder to the sphere will increase the output power. Decreasing the PTO damping coefficient will

increase the output power. The inclination angle of the buoy has the least effect on the output power.

7. Two-body WEC study in irregular waves From Table 2, it can be seen that the WEC with the Group 5 parameter combination has the maximum average power output of 216.87 kW in regular wave. However, the performance of the Group 5 WEC device in irregular waves is yet to be known. Thereby, the Group 5 WEC device is implemented in irregular waves to determine its power output to detect if the optimized WEC from Taguchi method can have good performance well in irregular wave. The Jonswap spectrum function (eq. (15)) coupled with Australian sea state (frequency range from 0.08 to 0.12 Hz with peak frequency as 0.1 Hz) of the significant wave height of 2 m and peak period of 10 s is utilized to generate the irregular wave in time domain using eqs.(12)e(16). The time domain traces of the wave elevation and excitation forces for floating and submerged bodies (by solving eqs. (17) and (18)) are presented in Fig. 8. It is seen from Fig. 8 that the wave elevation fluctuates from 1.2 m to 1.2 m. The buoy excitation force fluctuates from 200 kN to 200 kN. The submerged body excitation force fluctuates from 35 kN to 35 kN which is much less than that of the buoy excitation force. Moreover, the time domain traces of the radiation damping force for the Group 5 WEC buoy and output power generated in irregular waves obtained from eqs. (4)e(6) and (11) are shown in Fig. 9. It is seen from Fig. 9 that the buoy radiation damping force fluctuates from 20 kN to 35 kN. The peak value of output power can reach up to 900 kW. The influences of the significant wave heights (2e4 m for Australian sea state) on the performance and efficiencies of the Group 5 non-linear WEC under different peak periods Tp are presented in Fig. 10 (a). Additionally, the effect of kpto and Tp on the efficiency for the Group 5 WEC is studied in this research. kpto changes with three values of 25000, 50000 and 100000 N/m in this case. The peak period Tp varies from 8 to 12 s corresponding with the Australian sea states. In Fig. 10 (b), (c) and (d), the influences of kpto and Tp on the efficiency of the Group 5 WEC under the same significant wave height are presented. It is shown from Fig. 10 (a) that the wave energy conversion efficiency does not change much with the significant wave height for those irregular waves where the peak periods are 8, 9, 10s. As the wave excitation period increases to 11 and 12s, the wave energy conversion efficiency increases with higher significant wave heights. It is seen from Fig. 10 (b) that given the constant PTO damping coefficient, the wave energy conversion efficiency increases with the PTO stiffness. It is seen from Fig. 10 (b), 10 (c) and 10 (d) that the peak wave energy conversion efficiency increases with the significant wave height in the cases of kpto ¼ 50000 N/m,

494

X. Ji et al. / Renewable Energy 147 (2020) 487e501

Fig. 6. (a) Group1-4 average power; (b) Group 5 average power fitting curve; (c) Group 6 average power fitting curve; (d) Group 7 average power fitting curve; (e) Group 8 average power fitting curve.

X. Ji et al. / Renewable Energy 147 (2020) 487e501

495

Fig. 7. Absolute values of the effect of the seven selected parameter variables on the maximum average power output for a regular wave excitation.

kpto ¼ 25000 N/m. 8. The comparison between the linear and non-linear hydrostatic force The most commonly used geometry for the buoy is the cylinder and its hydrostatic force is generally modelled as a linear spring force [28], whereas it is not suitable to apply this model to other geometries with a non-uniform shape in heaving direction such as the sphere. Based on Archimedes' principle, the buoyancy force of an object is equal to the gravitational weight force of the water replaced by the object's submerged part. The volume of the water replaced by a sphere or the inclined cylinder in this research is no longer proportional to its displacement in heaving direction. Therefore, it is not suitable to model the hydrostatic force of the non-uniform shaped buoy as a linear spring force. The hydrostatic forces of the Group 5 WEC buoy (25 inclined cylinder), a sphere and a cylinder at the static water surface are shown in Fig. 11. It is seen from Fig. 11 that the nonlinear hydrostatic force of the Group 5 WEC buoy is larger than nonlinear hydrostatic force of a sphere and a linear hydrostatic force of a cylinder over the positive range of the oscillating displacement of the buoy and it is smaller when the floating buoy is dragged down. The nonlinear hydrostatic force of a sphere is close to the linear hydrostatic force of a cylinder. The values of the hydrostatic forces are calculated by

fgfbuoy ¼

rgp 3

! 3Rx2 x3 þ 3R2 x þ tanð90  qÞ tan2 ð90  qÞ

(25)

  1 fsphere ¼ rg p xR2  x3 3

(26)

flinear ¼ rg pR2 x

(27)

where x is the instantaneous buoy displacement, and R is the radius of the buoy of 2.5337m; fgfbuoy is the nonlinear hydrostatic force of

the Group 5 WEC buoy; fsphere is the nonlinear hydrostatic force of a sphere; flinear is the linear hydrostatic force of a cylinder. It can be observed that if the displacement of buoy is small then the hydrostatic force can be modelled as a linear spring force but if the displacement is further increased, the deviation between the non-linear buoyancy force and the linear buoyancy force becomes large. Moreover, the influences of linear and non-linear hydrostatic forces on the displacement of the Group 5 WEC buoy in irregular waves are shown in Fig. 12 (a). The impacts of the linear and the non-linear models on the buoy velocity are shown in Fig. 12 (b). It is seen from Fig. 12 (a) that the displacement and velocity curves of the Group five WEC buoy in irregular waves simulated by the linear and non-linear models are very close to each other, but the non-linear model has larger displacement in the negative direction. The peak amplitudes of the buoy displacement and velocity curves simulated by the non-linear model are larger than those simulated by the linear model at some instances.

9. Discussions From Table 2, it can be noticed that the submerged body geometry and the PTO's stiffness coefficient have the effects of 66.72 kW and 90.75 kW respectively and are most prominent two parameters in influencing the power output. The positive 66.72 kW means from Level 1 of the submerged body of the cylinder shape to Level 2 of the submerged body of the sphere shape, the output power increases by 66.72 kW. The negative 90.75 kW means from Level 1 of the PTO stiffness of 100,000 N/m being reduced to Level 2 of the PTO stiffness of 50000 N/m, the output power decreases by 90.75 kW. The third important influencing parameter is the damping coefficient of the PTO. The fourth important influencing parameter is the diameter of the floater. The effects of the remaining parameters on the maximum average output power are quite similar and small. But the inclination angle parameter which influences the linearity of the numerical simulation model has the least effect of 0.86 kW on the maximum average power output. When the inclination

496

X. Ji et al. / Renewable Energy 147 (2020) 487e501

Fig. 8. Wave elevation and wave excitation force for buoy and submerged body (Hs ¼ 2 m, Peak frequency 0.6283 rad/s).

angle increases from 0 to 25 , the output power is decreased by 0.86 kW. The inclination angle of the buoy is considered to affect the hydrostatic force. If the hydrostatic force is increased through the increment of the wave height, the effect of the hydrostatic force may be more prominent and play a more important role in determining the maximum average power output. As predicted, the shape of the submerged body has a significant effect on the time averaged power output of the WEC. The large viscous damping coefficient dissipates the kinetic energy and hence reduces the average power output. The added masses heavily increase the total masses and enable the WEC to resonate at a low ocean wave frequency. For the linear model, the natural resonant frequency is related to the PTO stiffness kpto and the mass of the WEC. Thereby, the increase of the WEC mass may lead to a reduced first natural frequency being less than 0.08 Hz. If the effect of the increase of the

WEC mass is less than that of the PTO stiffness kpto, it may result in the increased natural resonant frequency of the WEC being larger than 0.12 Hz. For both the situations, increase of the PTO stiffness kpto will increase the maximum average power output in the defined frequency range and consequently the large effect is shown in Table 2. The similar effects of the buoy diameter and draft of 33.62 kW and 21.35 kW respectively, on the average power output can be found. Increasing the diameter of the floater from 4 m to 6 m will increase the output power by 33.62 kW. Increasing the draft of the floater from 1 m to 2 m will decrease the output power by 21.35 kW. These two parameters can affect the hydrodynamic parameters such as the radiation damping coefficient, added mass and excitation force. In addition, the hydrostatic stiffness can also be determined by these two parameter settings. Consequently, these two parameter settings can affect the natural resonant frequency of the

X. Ji et al. / Renewable Energy 147 (2020) 487e501

497

Fig. 9. The time domain traces of the radiation damping force of buoy and instantaneous power output for Group 5 WEC in irregular wave (Hs ¼ 2 m, Peak frequency 0.6283 rad/s).

linear model through the hydrostatic stiffness. If the natural resonant frequency is in the excitation wave frequency range, the average power output will be high and if the natural resonant frequency is out of the range then the average power output will be low. This is why these two parameters are all coupled and all affect the average power output. The PTO damping coefficient bpto having the effect of 40.48 kW can also slightly affect the natural frequency. According to Equation (10), The PTO damping coefficient bpto also has an impact on energy produced. Decreasing the damping coefficient from 100,000 kg/s to 50,000 kg/s will increase the output power by 40.48 kW. The submerged depth having the effect of 23.45 kW slightly influences the average power output which has verified the statement [11] that increasing the submerged depth can also reduce the power output. This is because the submerged depth can influence the magnitude of the wave excitation force of the submerged body and thereby affect the relative motion between the buoy and the submerged body. From simulation of all the WEC groups, it can be found that Group 5 WEC device produces the maximum average power output of 216.87 kW among all the WEC groups. The optimized parameter combination of Group 5 WEC device is: 1) the inclination angle of the floating cylinder body of 25 ; 2) the floating cylinder body draft of 1 m; 3) the submerged geometry of sphere; 4) submerged depth of 20 m; 5) the PTO damping coefficient of 50000 kg/s; 6) the PTO stiffness coefficient of 100000 N/m; 7) the diameter of the cylinder floater of 6 m. Since all the parameters are coupled and affect each other, it is unlikely to

identify the other options to get a better result. Therefore, it can be claimed that in all the parameter combination settings of Table 2, the Group 5 parameter combination is the best choice. For the Group 5 WEC device in irregular waves, it can be seen from Fig. 10 (a) that the efficiency of the Group 5 (bpto ¼ 50000 kg/s, kpto ¼ 100000 N/m) WEC device does not change as the significant wave height increases for the peak periods of 8e10s but the efficiency is found to improve under the same significant wave height. The efficiency is found to increase with the significant wave height in the peak period of 11 s or 12 s. At the same significant wave height, the efficiency of the Group 5 WEC (bpto ¼ 50000 kg/s, kpto ¼ 100000 N/m) generally will increase with peak period albeit with the exceptional cases where peak periods are 10 and 11s whose efficiency is lower than that of 10 s for significant wave height smaller than 3m. Thereby, according to Fig. 10(a), the WEC device of bpto ¼ 50000 kg/s, kpto ¼ 100000 N/m, the peak period of 11s under the irregular wave excitation reaches its peak efficiency in the significant wave heights range of 2e4 m (the Australian sea state). From Fig. 10 (b), 10 (c) and 10 (d), it can be observed that the WEC devices (Groups 5 WEC device) of kpto ¼ 100000 N/m and bpto ¼ 50000 kg/s respectively (marked in plus signs) as set in Table 2, have the best efficiency among all significant wave heights and peaks periods compared with the other WEC devices where the PTO stiffness coefficient kpto are 25000 N/m and 50000N/m and PTO damping coefficient bpto is 50000 kg/s. Moreover, the efficiency variation of the WEC devices marked in the plus signs over the peak

498

X. Ji et al. / Renewable Energy 147 (2020) 487e501

Fig. 10. (a) The influence of different significant wave heights on the efficiency (bpto¼50000 kg/s, kpto¼100000 N/m); (b) The influence of different peak periods and stiffness coefficients of PTO on the efficiency (Hs ¼2 m); (c) The influence of different peak periods and stiffness coefficients of PTO on the efficiency (Hs ¼3 m); (d) The influence of different peak periods and stiffness coefficients of PTO on the efficiency (Hs ¼4 m).

period is most prominent and it can be changed from 7% (Hs ¼ 4 m, peak period ¼ 8 s) to 51% (Hs ¼ 3 m, peak period ¼ 11 s) when the peak period of irregular wave shifts from 8 to 12s. For the parameter setting of the WEC devices of bpto ¼ 50000 kg/s, kpto ¼ 25000 N/m marked in the square signs in Fig. 10(b), (c) and 10(d), it can be observed that the efficiency of the WEC devices marked in the square signs is not sensitive to the value of the peak period at Hs ¼ 2 m. The efficiency changes much largely with the peak period when the significant wave heights are increased to 3 m and 4 m. The efficiencies for the same peak period also increase with the significant wave heights. For the parameter setting of the WEC devices of bpto ¼ 50000 kg/s, kpto ¼ 50000 N/m, marked in the star signs in Fig. 10(b), (c) and 10(d), the efficiency reaches its peak value when peak period is 10 or 11 s and then it decreases with the peak period for all the significant wave heights until the peak period reaches 12s. The efficiencies for all the peak periods also increase with the significant wave heights. In addition, the PTO stiffness kpto is considered to largely affect the efficiency even in irregular waves. For a large PTO stiffness, the efficiency of the WEC

devices will be more obviously influenced by the peak period. Overall, the WEC devices with the optimal PTO parameter can achieve the efficiency from 7% to 51% in the Australian sea state, which indicates that the parameter combination for the Group 5 WEC device has the outstanding efficiency in the irregular waves as well. If the PTO stiffness kpto may be not able to reach 100000 N/m whilst the bpto is 50000 kg/s, the parameter combination of kpto ¼ 50000 N/m and bpto ¼ 50000 kg/s seems to be an alternative as they are practically low. For this scenario, the efficiency being shifted from 3%e20% can be regarded as reasonable. The classical hydrostatic force model for a buoy is a linear spring and proportional with the buoy displacement. To model this force more precisely, the non-linear model is studied and the differences of the non-linear and linear hydrostatic forces for a sphere, a cylinder and Group 5 WEC buoy are shown in Fig. 11. The deviations are increased with the displacement. The time domain traces of the displacements and velocities of the Group 5 WEC buoy coupled with a sphere in the 20 m depth for the linear and non-linear models in the irregular waves (Hs ¼ 3 m and peak period ¼ 10 s)

X. Ji et al. / Renewable Energy 147 (2020) 487e501

499

Fig. 11. Hydrostatic forces for the Group 5 WEC buoy of a sphere and a cylinder with same radius (2.5337m) at static ocean surface.

with kpto ¼ 50000 N/m and bpto ¼ 50000 kg/s are shown in Fig. 12 (a) and 12 (b). From the Archimedes’ law, the Group 5 WEC buoy shall be easier to be dragged down than to be pushed up under the influence of the wave excitation forces, due to the buoy draft volume difference in the up and down directions. Therefore, from Fig. 12 (a), it can be observed that the negative displacement amplitude of the non-linear model is higher than that of the linear model in the same wave conditions. The linear model has bigger displacement in positive direction than the nonlinear model because the resistance to go up of the linear model is not as large as the non-linear model. The deviation is small when the wave elevation is less significant but become large when the elevation amplitude is high. From Fig. 12 (b), it can be seen that the velocity for the non-linear model is generally higher than that of linear model in the positive direction, which can be explained by the impulse-momentum theorem. This is because the hydrostatic force of the nonlinear model is much larger than that of the linear model in the positive displacement direction. The negative velocities of these two models are very close to each other, because the hydrostatic forces of the linear and nonlinear models are close to each other in the negative displacement direction. Therefore, the difference of the hydrostatic forces for the non-linear and linear models has a large effect on the dynamic performance of the buoy which cannot be ignored.

10. Conclusions A two-body oscillating wave energy converter has been analysed by varying seven different system model parameters using the Taguchi method. ANSYS AQWA has been used in identifying the hydrodynamic parameters for the different parameter combination groups. Based on the frequency and time domain

simulations using the Matlab codes, the maximum average power output can be calculated for the linear and non-linear models in the frequency range from 0.08 Hz to 0.12 Hz in the Australia ocean area. It is found that the PTO stiffness coefficient and submerged body geometry have the most significant effects on the maximum average power output in the regular waves and the PTO stiffness coefficient also seriously affects the performance of WEC in the irregular waves. The best parameter combination group for the performance among the designed eight parameter combination groups is Group 5 WEC device and additionally this group also shows the maximum efficiency of 51% in the irregular waves for the significant wave heights of 2e4 m and the peak period range from 8 to 12 s. Even though the effect of the buoy inclination angle is not so prominent in this research in the regular waves, it is still necessary to conduct a further research. Because this paper mainly compares the maximum average power output within the frequency range from 0.08 to 0.12 Hz, more efforts should be made in investigating the effect of the model parameters on the energy harvesting frequency bandwidth and overall harvesting performance of the proposed WEC system in the Australia sea state. It is also necessary to reduce the viscous damping coefficient and optimize the geometry of the submerged body in order to achieve the best WEC energy harvesting performance in the desired operating frequency range. In irregular waves, the higher PTO stiffness will result in a larger sensitivity of the output power to the peak period of the Jonswap spectrum. The significant wave height has few impacts on the efficiency of the WEC with optimal PTO parameter (bpto ¼ 50000 kg/s, kpto ¼ 100000 N/m). However, the significant wave height can affect the efficiency of the WEC when kpto ¼ 50000 N/m and kpto ¼ 25000 N/m. Moreover, the Group 5 WEC device has been proved to be the best option in extracting wave energy in both the regular and irregular waves for

500

X. Ji et al. / Renewable Energy 147 (2020) 487e501

Fig. 12. (a) the displacements of Group five buoy in irregular waves (Hs ¼ 3 m, peak period ¼ 10 s and kpto ¼ 50000 N/m, bpto ¼ 50000 kg/s) for the linear and non-linear models; (b) The Velocities of Group five buoy in irregular wave (Hs ¼ 3 m, peak period ¼ 10 s and kpto ¼ 50000 N/m, bpto ¼ 50000 kg/s) for linear and non-linear model.

the Australian sea state. The non-linear hydrostatic force has essential impacts on the dynamic performance of the buoy and thus precise hydrostatic force model is necessary for the nonuniform shaped buoy in order to simulate more accurately the performance of a point absorber.

Acknowledgement Authors would like thank Australia Research Council Discovery Project grant DP170101039 for financial support.

X. Ji et al. / Renewable Energy 147 (2020) 487e501

501

Appendix. Appendix

Fig. 13. The illustration of whole procedure.

References [1] J. Falnes, A review of wave-energy extraction, Mar. Struct. 20 (4) (2007) 185e201. [2] Benjamin Drew, Andrew R. Plummer, M. Necip Sahinkaya, A Review of Wave Energy Converter Technology, 2009, pp. 887e902. ~o, F.O. Anto  nio, Joa ~o CC. Henriques, Oscillating-water-column wave en[3] Falca ergy converters and air turbines: a review, Renew. Energy 85 (2016) 1391e1424. [4] James Tedd, Jens Peter Kofoed, Measurements of overtopping flow time series on the Wave Dragon, wave energy converter, Renew. Energy 34 (3) (2009) 711e717. [5] S.J. Illesinghe, et al., Idealized design parameters of Wave Energy Converters in a range of ocean wave climates, Int. J. Mar.Energy 19 (2017) 55e69. [6] W. Dick, U.S. Patent No. 6, U.S. Patent and Trademark Office, Washington, DC, 2005, p. 857, 266. [7] Guang Li, et al., Wave energy converter control by wave prediction and dynamic programming, Renew. Energy 48 (2012) 392e403. [8] C. Liang, L. Zuo, On the dynamics and design of a two-body wave energy converter, Renew. Energy 101 (2017) 265e274. [9] A. Amiri, P. Roozbeh, R. Soheil, Parametric study of two-body floating-point wave absorber, J. Mar. Sci. Appl. 15 (1) (2016) 41e49. [10] Y. Dai, Y. Chen, L. Xie, A study on a novel two-body floating wave energy converter, Ocean. Eng. 130 (2017) 407e416. €m, et al., A resonant two body system for a point absorbing wave [11] J. Engstro energy converter with direct-driven linear generator, J. Appl. Phys. 110 (12) (2011) 124904. [12] Y. Gao, et al., A fully floating system for a wave energy converter with directdriven linear generator, Energy 95 (2016) 99e109. [13] C.J. Cargo, et al., Determination of optimal parameters for a hydraulic power take-off unit of a wave energy converter in regular waves, Proc Inst Mech Eng A J Power Energy 226 (1) (2012) 98e111. pez, P. Taveira, P. Rosa-Santos, Numerical Modelling of the CECO Wave [14] M.F. Lo Energy Converter, Renewable Energy, 2017.

[15] C. Liang, A. Junxiao, L. Zuo, Design, fabrication, simulation and testing of an ocean wave energy converter with mechanical motion rectifier, Ocean. Eng. 136 (2017) 190e200. [16] http://reneweconomy.com.au/worlds-first-grid-connected-wave-energyarray-switched-on-in-perth-77510/. [17] G. Taguchi, Taguchi on Robust Technology Development : Bringing Quality Engineering Upstream (ASME Press Series on International Advances in Design Productivity), ASME Press, New York, 1993. [18] Scott J. Beatty, et al., Experimental and numerical comparisons of self-reacting point absorber wave energy converters in regular waves, Ocean. Eng. 104 (2015) 370e386. [19] J.C.C. Henriques, et al., On the annual wave energy absorption by two-body heaving WECs with latching control, Renew. Energy 45 (2012) 31e40. [20] Elie Al Shami, Ran Zhang, Xu Wang, Point Absorber wave energy harvesters: a review of recent developments, Energies 12 (1) (2019) 47. [21] G. Tampier, L. Grueter, Hydrodynamic analysis of a heaving wave energy converter, Int. J. Mar.Energy 19 (2017) 204e318. [22] Davood Younesian, Mohammad-Reza Alam, Multi-stable mechanisms for high-efficiency and broadband ocean wave energy harvesting, Appl. Energy 197 (2017) 292e302. [23] Sighard F. Hoerner, Fluid-dynamic Drag: Practical Information on Aerodynamic Drag and Hydrodynamic Resistance, Hoerner Fluid Dynamics, 1965. [24] Z. Cheng, J. Yang, Z. Hu, L. Xiao, Frequency/time domain modeling of a direct drive point absorber wave energy converter, Sci. China Phys. Mech. Astron. 57 (2) (2014) 311e320. pez, F. Taveira-Pinto, P. Rosa-Santos, Influence of the power take-off [25] M. Lo characteristics on the performance of CECO wave energy converter, Energy 120 (2017) 686e697. ~o, P.A. Justino, Nonlinear dynamics of a tightly moored [26] P.C. Vicente, A.F. Falca point-absorber wave energy converter, Ocean Eng. 59 (2013) 20e36. [27] Johannes Falnes, On non-causal impulse response functions related to propagating water waves, Appl. Ocean Res. 17 (6) (1995) 379e389. [28] Johannes Falnes, Ocean Waves and Oscillating Systems: Linear Interactions Including Wave-Energy Extraction, Cambridge university press, 2002.