Study on the Resonant Behaviors of a Bottom-Hinged Oscillating Wave Surge Converter

Abstract

This paper studied the resonant behaviors of a bottom-hinged oscillating wave surge converter (OWSC) as well as the relationship of resonance with the response and capture width ratio (CWR). The time-domain dynamic equation of an OWSC in shallow water based on the boundary element method (BEM) was solved by a Python code, considering the corrected wave surface and the nonlinearities of restoring moment, drag, and friction. The unknown factors, such as wave surface corrected factor and drag coefficient, were effectively calibrated with computational fluid dynamics (CFD) method. An intermediate initial angle in free decay is appropriate for use to determine the natural period. Under regular waves, the resonance occurs near the natural period for the uniform wave amplitude, rather than the uniform wave torque amplitude, and can disappear due to the amplification of Power Take-Off (PTO) friction. Under unit-amplitude regular waves, the period of maximum CWR is relatively close to the period of maximum velocity, but far from the resonant period. Under irregular waves, no stable resonance is observed because the maximum equivalent pitch angle appears at different peak periods of wave spectra with the variation in PTO damping. When the period of a regular wave or the peak period of an irregular wave is close to the natural period, a phase hysteresis of velocity relative to wave torque always occurs.

1. Introduction

The bottom-hinged oscillating wave surge converter (OWSC), a nearshore wave energy converter, captures wave energy by oscillating under waves. It is installed in shallow water and benefits from the advantage of the amplification of horizontal water particle velocity due to shoaling effects [1]. Budal’s upper bound [2] explained that the maximum value of theoretical capture width ratio (CWR) of the flap-type absorber was double that of the axisymmetric point absorber. It was reported in [3] that the bottom-hinged OWSC captures greater mean annual power than most types of wave energy converters (WECs).

Oyster 1, a full-scale prototype of OWSC with a rated power of 315 kW, was developed by Aquamarine Power, collaborating with Queen’s University Belfast, installed at the European Marine Energy Center (EMEC) in Orkney, Scotland, in 2009, and decommissioned after its 2-year design life. The second-generation Oyster 800, with a rated power of 800 kW, was tested at EMEC in 2011, and successfully exported electricity to the national grid [4].

The hydrodynamic research of OWSCs has received attention in recent decades. Some experimental investigations on the motions and power capture of OWSCs in two-dimensional (2D) wave tank [5,6,7] and three-dimensional (3D) wave basin [8,9] were conducted by researchers. With the development of cheap and powerful computer hardware, numerical simulation is becoming more popular and there is a rapidly growing number of modelling WECs. Boundary element method (BEM) numerical simulation, based on potential flow theory, was employed in the analysis of OWSCs under waves, and validated via experiments [10]. Some hydrodynamic studies of OWSCs based on BEM can be found in [11,12,13]. BEM has the merit of efficient computation, but the nonlinearities, e.g., wave breaking [14], flow separation [15], vortex shedding [16], slamming [17], and wave overtopping, cannot be perfectly considered. A numerical wave tank (NWT), based on computational fluid dynamics (CFD) via the solving of Navier–Stokes equations, has been widely applied in the high-fidelity modeling of the interaction between waves and the motion of OWSCs [18,19,20,21,22] and other types of WECs [23]. The smoothed particle hydrodynamics (SPH) method, a purely Lagrangian meshless interpolation method [24], was used to accurately simulate the wave intertion with the rectangular OWSCs [25,26] and a cylindrical OWSC [27]. The details of these approaches (BEM, CFD, and SPH) were comparatively summarized in [28,29]. However, in the simulation of CFD and SPH, long-term motion modeling is not feasible for engineering tasks due to the high computational cost. Therefore, to ensure the accuracy, and simultaneously shorten the elapsed time, of computation, the BEM results calibrated with CFD simulations contribute to a better overall application of the hydrodynamic analysis of OWSCs.

Resonance study is one of the subsystems of hydrodynamic research. Some studies demonstrated the structural dynamic amplifications near resonance in calculating the reaction forces for offshore platforms [30], coastal bridges [31], and piles [32] subjected to the wave impacts. Taking inspiration from the above phenomenon, resonance could be prosepctive in the promotion of wave energy absorption by amplifying the dynamic response of WECs. For example, tuning a point absorber to resonance, when its undamped natural period is close to the period of a regular wave [33] or the dominant period of an irregular wave [34], can be recognized as an effective approach to derive more energy. However, a flap-type WEC, usually of relatively large size, has a different resonant mechanism to the traditional point absorber, due to its diffracting of the wave’s propagation. The resonant behaviors of an OWSC in a channel were investigated using a semi-analytical method and concluded that the resonant sloshing modes lead to increase the maximum CWR [35]. The research of an OWSC in front of a vertical wall also proved that the much higher levels of CWR were achieved when the OWSC was located at a specific distance from the vertical wall, in which the larger wave torque was produced due to an appropriate phase difference between the reflected wave and the exciting wave [36]. However, less attention has been paid to the study of resonant wave interactions with OWSCs in the open sea.

In this work, the dynamic equation of a 3D OWSC in shallow water under regular and irregular waves is formulated in Section 2, based on time-domain BEM, and the nonlinear restoring stiffness, drag, and friction are taken into account. To calibrate the BEM results in relation to two unknown factors ((1) the wave surface corrected factor altering the restoring moment and drag moment, and (2) the drag coefficient), the numerical methodology of CFD is introduced in Section 3. The study model and mesh are described in Section 4. The free decay analysis, and the relationship of resonance to the responses and CWRs under regular and irregular waves, are discussed in Section 5, Section 6 and Section 7. Finally, the conclusions are summarized in Section 8.

2. Numerical Model Based on BEM

A 3D bottom-hinged OWSC, often called a flap, of width $w$, thickness $b,$ and height $h$, captures wave power by oscillating about a horizontal hinge $H$, fixed to a base constructed on the seabed at water depth $d$, as shown in Figure 1. The distance between the axis and the still water surface is $d_a$, called the axis depth. If $d_a<h$, the flap is a surface-piercing flap, the height of which above still water is $h_f=h-d_a$, refered to as freeboard. As the water level increases to $d_a\ge h$, the flap will become fully submerged, but the top of the flap may emerge out of the water when a wave trough approaches. Points $B$ and $G$ are, respectively, buoyant center and gravity center. $\phi$ is the rotational angle. $\dot{\phi }$ and $\ddot{\phi }$ are angular velocity and acceleration, respectively.

The wave propagation direction is assumed to be unidirectional, and perpendicular to the flap (x-axis of a Cartesian coordinate system), which is an ideal scenario. Some studies on different types of coastal structures, such as low-crested structures [37], dikes [38], and bridge decks [39] demonstrated a significant reduction in wave overtopping and the wave-induced forces in relative non-normal angles between the structures and the waves (e.g., oblique waves or skewed structures). Although the effects of oblique waves is not negligible, the research of OWSCs under perpendicular waves in this work is still meaningful as the foundation for future studies, comparing the influence of oblique waves with perpendicular waves.

The incident waves are not allowed very high, based on the assumption of a linear relationship between wave amplitude and wave exciting moment amplitude. In addition, some nonlinearities, such as wave breaking, wave overtopping, flow separation, vortex shedding, and slamming, cannot be taken into account in BEM-based methodology. However, viscous drag will be simply considered (see Section 2.1.5).

2.1. Environmental Load in the Sea

2.1.1. Incident Wave Elevation and Corrected Wave Surface

Incident wave elevation at an arbitrary position, in the absence of flap, can be expressed by a sine function for regular waves, or the summarization of multi-components of sinusoidal waves for irregular waves:

$$\eta \left(t\right)=\{\begin{array}{c}{A}_i cos\left(kx-\omega t\right) \mathrm{Regular} \mathrm{waves}\\ {\displaystyle {\displaystyle \sum }_{n=1}^{N_w}}A\left({\omega }_n\right) cos\left[k_{n}x-{\omega }_{n}t+\varphi \left({\omega }_n\right)\right] \mathrm{Irregular} \mathrm{waves}\end{array}$$

(1)

where $A_i$ is the incident wave amplitude of a regular wave and $A\left({\omega }_n\right)$ is the wave amplitude of the $n$-th wave component of an irregular wave; $k$ is the wave number for a finite depth and functionally related to the dispersion relation ${\omega }^2=gk\tanh \left(kd\right)$; $k_n$ is the wave number of the $n$-th wave component; $g$ is the gravity acceleration; $\omega$ and ${\omega }_n$ are the angular frequencies; $t$ is the time; $N_w$ is the number of wave components; $\varphi \left({\omega }_n\right)\in \left[0, 2\pi \right]$ is the random phase.

For irregular waves, $A\left({\omega }_n\right)$ is determined by the wave spectrum,

$$A\left({\omega }_n\right)=\sqrt{2S_d\left({\omega }_n\right)\Delta \omega }$$

(2)

where $S_d\left({\omega }_n\right)$ is the modified wave spectrum for a finite water depth, introduced in Section 7.1; $\Delta \omega$ is the wave frequency interval of wave components.

However, the actual wave surface, with the existence of a flap, is different from the incident wave. In this study, the wave surface is corrected, as shown in Figure 2, to adjust the restoring stiffness and drag moment for more accurate results. $R_w\left(t\right)$ is the instantaneous distance between the hinge point $H$ and its projection on the corrected wave surface along the mid-surface of the flap, and this can be given by

$$R_w\left(t\right)={\alpha }_{\eta }\left[{\eta }_w\left(t\right)-d_a\right]+d_a$$

(3)

where ${\alpha }_{\eta }$ is the wave surface corrected factor; ${\eta }_w\left(t\right)$ is the instantaneous distance between point $H$ and its projection on the incident wave surface, obtained via geometric measurement according to the incident wave elevation $\eta \left(t\right)$ and the flap angle $\phi \left(t\right)$.

2.1.2. Wave Exciting Moment

The wave exciting moment $M_E\left(t\right)$, including the Froude–Krylov moment and the diffraction moment for a 3D model, can be described as

$$M_E\left(t\right)=\{\begin{array}{l}{M}_e\left(\omega \right)A_i \mathit{cos}\left[-\omega t+\psi \left(\omega \right)\right] \mathrm{Regular} \mathrm{waves}\\ {\displaystyle {\displaystyle \sum }_{n=1}^{N_w}}{M}_e\left({\omega }_n\right)A\left({\omega }_n\right) \mathit{cos}\left[-{\omega }_{n}t+\varphi \left({\omega }_n\right)+\psi \left({\omega }_n\right)\right] \mathrm{Irregular} \mathrm{waves}\end{array}$$

(4)

where $M_e\left(\omega \right)$ is the wave exciting moment amplitude for unit wave amplitude, and $\psi \left(\omega \right)$ is the phase difference between the incident wave and the wave exciting moment.

2.1.3. Radiation Damping Moment

Based on the Cummins equation [40] and the convolution integral [41], the time domain radiation damping moment $M_R\left(t\right)$, incorporating the effects of wave-making memory, is given by

$$M_R\left(t\right)=-\mu (\infty )\ddot{\phi }\left(t\right)-{{\displaystyle \int }}_0^{t}R\left(t-\tau \right)\dot{\phi }\left(\tau \right)d\tau$$

(5)

where $\mu (\infty )$ is the added inertia torque at infinite frequency; $R\left(t\right)$ is the impulse response function,

$$R\left(t\right)=\frac{2}{\pi }{{\displaystyle \int }}_0^{\infty }c\left(\omega \right) cos\left(\omega t\right)d\omega$$

(6)

where $c\left(\omega \right)$ is the radiation damping coefficient.

2.1.4. Restoring Moment

Considering the variations of the buoyant center and buoyant force due to the time-varying corrected wave surface, the nonlinear hydrostatic restoring moment $M_H\left(t\right)$ can be described as

$$M_H\left(t\right)=-\left[F_b\left(t\right)·\overline{BH}\left(t\right)-mg·\overline{GH}\right]·\mathit{sin}\phi \left(t\right)$$

(7)

where $F_b\left(t\right)$ is the buoyant force; $\overline{BH}\left(t\right)$ is the force arm of buoyancy; $m$ is the mass. Both $F_b\left(t\right)$ and $\overline{BH}\left(t\right)$ are dependent on the wet surface (below the corrected wave surface),

$$F_b\left(t\right)=\rho gw\left[h_w\left(t\right)b+\frac{\pi b^2}{8}\right]$$

(8)

$$\overline{BH}\left(t\right)=\frac{12h_w^2\left(t\right)-2b^2}{24h_w\left(t\right)+3\pi b}$$

(9)

where $h_w\left(t\right)=min\left[R_w\left(t\right),h\right]$ is the instantaneous height of the wet surface above the hinge.

2.1.5. Drag Moment

As evolved from the Morison equation [42], the nonlinear drag moment $M_D\left(t\right)$ acting on the wet surface (below the corrected wave surface) can be expressed as

$$M_D\left(t\right)=-\frac{1}{2}{C}_d\rho w{{\displaystyle \int }}_0^{h_w\left(t\right)}\left[\dot{\phi }\left(t\right)r-u_n\left(t\right)\right]·\left|\dot{\phi }\left(t\right)r-u_n\left(t\right)\right|·rdr$$

(10)

where $C_d$ is the drag coefficient; $r$ is the radial coordinate of the flap based on the hinge; $u_n\left(t\right)$ is the normal velocity of the water particles (in the absence of the flap) on the mid-surface of the flap, and is given by

$$u_n\left(t\right)=u_x\left(t\right) \mathit{cos}\phi \left(t\right)-u_z\left(t\right) \mathit{sin}\phi \left(t\right)$$

(11)

where $u_x$ and $u_z$ are, respectively, the horizontal and vertical components of the water particle velocities in the absence of the flap at a finite depth [43],

$$u_x\left(t\right)=\{\begin{array}{l}\frac{gk{A}_i}{\omega }\frac{\mathit{cosh}\left[k\left(z+d\right)\right]}{\mathit{cosh}\left(kd\right)} \mathit{cos}\left(kx-\omega t\right) \mathrm{Regular}\mathrm{waves}\\ {\displaystyle {\displaystyle \sum }_{n=1}^{N_w}}\frac{g{k}_{n}A\left({\omega }_n\right)}{{\omega }_n}\frac{\mathit{cosh}[k_n\left(z+d\right)]}{\mathit{cosh}\left(k_{n}d\right)}\mathit{cos}[k_{n}x-{\omega }_{n}t+\varphi \left({\omega }_n\right)] \mathrm{Irregular} \mathrm{waves}\end{array}$$

(12)

$$u_z\left(t\right)=\{\begin{array}{l}\frac{gk{A}_i}{\omega }\frac{\mathit{sinh}\left[k\left(z+d\right)\right]}{\mathit{cosh}\left(kd\right)} \mathit{sin}\left(kx-\omega t\right) \mathrm{Regular}\mathrm{waves}\\ {\displaystyle {\displaystyle \sum }_{n=1}^{N_w}}\frac{g{k}_{n}A\left({\omega }_n\right)}{{\omega }_n}\frac{\mathit{sinh}\left[k_n\left(z+d\right)\right]}{\mathit{cosh}\left(k_{n}d\right)}\mathit{sin}[k_{n}x-{\omega }_{n}t+\varphi \left({\omega }_n\right)] \mathrm{Irregular} \mathrm{waves}\end{array}$$

(13)

2.1.6. Power Take-Off (PTO) Moment

Ignoring the inertia and stiffness, the PTO system can be simplified as constant damping and friction, and the PTO moment $M_{PTO}\left(t\right)$ can be expressed as

$$M_{PTO}\left(t\right)=-C_{PTO}\dot{\phi }\left(t\right)-M_f\left(t\right)$$

(14)

where $C_{PTO}$ is the damping of the PTO system for wave energy harvesting; $M_f\left(t\right)$ is the frictional moment against the motion of the flap,

$$M_f\left(t\right)=\{\begin{array}{l} T_{PTO}, \dot{\phi }\left(t\right)>0\\ 0, \dot{\phi }\left(t\right)=0\\ -T_{PTO}, \dot{\phi }\left(t\right)<0\end{array}$$

(15)

where $T_{PTO}$ is a constant value of friction.

2.2. Time Domain Dynamic Equation

In the time domain analysis, the dynamic equation of the pitch motion can be written as

$$I_H\ddot{\phi }\left(t\right)=M_E\left(t\right)+M_R\left(t\right)+M_H\left(t\right)+M_D\left(t\right)+M_{PTO}\left(t\right)$$

(16)

where $I_H$ is the inertia torque around the hinge.

The workflow of solving the dynamic equation is illustrated in Figure 3. For a certain geometry of study model, the hydrodynamic coefficients $M_e\left(\omega \right)$, $\mu \left(\omega \right)$, $\mu (\infty )$, $c\left(\omega \right),$ and $\psi \left(\omega \right)$ are obtained for one time from NEMOH, an open source frequency domain BEM program. According to the inputted wave parameters (e.g., in still water, under regular waves or under irregular waves) and initial conditions of the flap (e.g., at an inclined angle for free decay or at the vertical position under waves), a Python code was developed by the authors to implement a time marching computation. In each time step, the process is as follows:

• Using the solutions of angle $\phi \left(t\right)$ and angular velocity $\dot{\phi }\left(t\right)$ at the previous time $t$, the environmental loads can be calculated by Equations (4), (5), (7), (10) and (14);
• Substitute the values of environmental loads into Equation (16) to assemble a dynamic differential equation;
• By 4th-order Runge–Kutta integration, the hydrodynamic responses of the flap at the current time ($t+dt$) can be computed, where $dt$ is the time step of integration.

2.3. Perormance Assessment

In the PTO system, the time–averaged capturing power $P_c$ is given by

$$P_c=\frac{C_{PTO}}{t_e-t_b}{{\displaystyle \int }}_{t_b}^{t_e}{\dot{\phi }}^2\left(t\right)dt$$

(17)

where $t_b$ and $t_e$ are, respectively, the beginning and end time for recording angular velocity $\dot{\phi }\left(t\right)$ in a duration.

The incident wave power of the resource per unit width $P_w$ of regular and irregular waves at a finite depth $d$ is given by

$$P_w=\{\begin{array}{c}\frac{1}{2}\rho g{A}_i{}^2{C}_g\left(\omega \right) \mathrm{Regular} \mathrm{waves}\\ \rho g{{\displaystyle \int }}_0^{\infty }{C}_g\left(\omega \right)S_d\left(\omega \right)d\omega \mathrm{Irregular} \mathrm{waves}\end{array}$$

(18)

where $C_g\left(\omega \right)$ is the wave group velocity,

$$C_g\left(\omega \right)=\frac{\omega }{2k}\left[1+\frac{2kd}{\mathit{sinh}\left(2kd\right)}\right]$$

(19)

The CWR under regular or irregular waves can be expressed as

$$CWR=\frac{P_c}{P_w·w}$$

(20)

2.4. Response of Motion

Although the actual time–history response under a regular wave is not perfectly harmonic, an equivalent sinusoidal response, based on energy conservation for a duration, can be introduced as

$$\phi \left(t\right)={\phi }_{amp} sin\left(\omega t+\theta \right)$$

(21)

where ${\phi }_{amp}$ is the equivalent amplitude of pitch; $\theta$ is the phase of the pitch. According to the kinetic energy conservation of equivalent sinusoidal motion with the actual response during the recording time, ${\phi }_{amp}$ can be calculated by

$${\phi }_{amp}=\frac{1}{\omega }\sqrt{\frac{2}{t_e-t_b}{{\displaystyle \int }}_{t_b}^{t_e}{\dot{\phi }}^2\left(t\right)dt}$$

(22)

Then, the equivalent amplitude of angular velocity can be written as

$${\dot{\phi }}_{amp}={\phi }_{amp}\omega$$

(23)

The phase angle of angular velocity ${\theta }_1$ relative to the wave exciting moment under regular waves can be estimated by measuring the averaged time difference between peak angular velocity and peak wave exciting moment in the time–history curves,

$${\theta }_1=\frac{2\pi }{N_{c}T}{\displaystyle \sum }_{n=1}^{N_c}{\left(t_{p, M_E}-t_{p, \dot{\phi }}\right)}_n$$

(24)

where $N_c$ is the number of peak values after the beginning time $t_b$; $t_{p, M_E}$ and $t_{p, \dot{\phi }}$ are, respectively, the time of peak moment and peak angular velocity. Similarly, the phase angle relative to the wave exciting moment can also be extracted with the same method, and expressed as ${\theta }_2$. Because the phase of the wave exciting moment is $-\psi \left(\omega \right)+\pi /2$, derived from Equation (4), the phase of pitch can be calculated by

$$\theta ={\theta }_2-\psi \left(\omega \right)+\pi /2$$

(25)

For a one-degree-of-freedom structure with constant inertia, stiffness, and damping, perfect resonance occurs when the period of harmonic excitation is close to its natural period. Simultaneously, the motion amplitude is maximal, and the velocity is in phase with the excitation, i.e., ${\theta }_1=0$. However, due to the strong nonlinearities of the flap under waves, the phase results are different from perfect resonance, discussed in Section 6.

Equations (21)–(25) are only available for regular waves.

Table 1. The phase angles of various items for regular waves.

Phase Angle of	Value
Incident wave elevation at $x=0$	$\pi /2$
Wave exciting moment relative to incident wave elevation	$-\psi \left(\omega \right)$
Angular velocity relative to wave exciting moment	${\theta }_1$
Pitch angle relative to wave exciting moment	${\theta }_2$
Pitch angle	$\theta$

presents the phase angles of various items for clear comparison.

To examine the response under irregular waves, according to kinetic energy conservation, the irregular response is also equivalent to the harmonic motion, the period of which is equal to the peak period, and the equivalent pitch amplitude is defined as

$${\phi }_{eq}=\frac{1}{{\omega }_p}\sqrt{\frac{2}{t_e-t_b}{{\displaystyle \int }}_{t_b}^{t_e}{\dot{\phi }}^2\left(t\right)dt}$$

(26)

where ${\omega }_p$ is the peak frequency of an irregular wave.

3. Numerical Model Based on CFD

3.1. Governing Equations

Based on the assumption of fluid incompressibility, the mass conservation equation for the NWT based on the Volume of Fluid (VOF) method [44] is given by

$$\frac{\partial \rho }{\partial t}+\nabla \left(\rho \stackrel{\rightharpoonup }{\mathit{v}}\right)=0$$

(27)

where $\stackrel{\rightharpoonup }{\mathit{v}}$ is the fluid velocity vector; $\rho =\left(1-\alpha \right){\rho }_1+\alpha {\rho }_2$ is the equivalent fluid density, calculated by air density ${\rho }_1$, water density ${\rho }_2$, and the volume fraction of water $\alpha$. The equivalent dynamic viscosity $\mu$ is also determined by $\mu =\left(1-\alpha \right){\mu }_1+\alpha {\mu }_2$, where ${\mu }_1$ and ${\mu }_2$ are the dynamic viscosity of air and water, respectively. The volume fraction can be calculated by the mass conservation equation.

The Reynolds-averaged Navier–Stokes (RANS) equation is given by

$$\frac{\partial \left(\rho \stackrel{\rightharpoonup }{\mathit{v}}\right)}{\partial t}+\nabla \left(\rho \stackrel{\rightharpoonup }{\mathit{v}}\stackrel{\rightharpoonup }{\mathit{v}}\right)=-\nabla p+\nabla \stackrel{\rightharpoonup }{\mathit{T}}+\stackrel{\rightharpoonup }{\mathit{S}}$$

(28)

where $p$ is the pressure; $\stackrel{\rightharpoonup }{\mathit{T}}$ is the stress tensor, including viscous stress and turbulent stress; $\stackrel{\rightharpoonup }{\mathit{S}}$ is the body force per unit mass (such as gravity). For incompressible Newtonian fluids and a turbulent model of eddy viscosity, $\stackrel{\rightharpoonup }{\mathit{T}}$ is defined as

$$\stackrel{\rightharpoonup }{\mathit{T}}=\left(\mu +{\mu }_t\right)\left[\nabla \stackrel{\rightharpoonup }{\mathit{v}}+{\left(\nabla \stackrel{\rightharpoonup }{\mathit{v}}\right)}^T\right]$$

(29)

where ${\mu }_t$ is the turbulent viscosity.

The dynamic equation of the flap pitch is written by

$$I_H\ddot{\phi }\left(t\right)=mg·\overline{GH}\mathit{sin}\phi \left(t\right)+M_{fluid}\left(t\right)+M_{PTO}\left(t\right)$$

(30)

where $M_{fluid}\left(t\right)$ is the total fluid moment acting on the flap surface, and is obtained by integrating the pressure and shear stress over the fluid–structure interface,

$$M_{fluid}\left(t\right)={{\displaystyle \int }}_s^{ }\left({\stackrel{\rightharpoonup }{\mathit{r}}}_s-{\stackrel{\rightharpoonup }{\mathit{r}}}_h\right)\times \left(\stackrel{\rightharpoonup }{\mathit{T}}-p\stackrel{\rightharpoonup }{\mathit{I}}\right)·\stackrel{\rightharpoonup }{\mathit{n}}ds$$

(31)

where $s$ is the fluid–structure interface; ${\stackrel{\rightharpoonup }{\mathit{r}}}_s$ is the position vector of the point on the flap surface; ${\stackrel{\rightharpoonup }{\mathit{r}}}_h$ is the position vector of the hinge; and $\stackrel{\rightharpoonup }{\mathit{n}}$ is the unit outer-pointing normal vector of the surface elements of the flap.

The motion of the flap and fluid are fully coupled by velocity and pressure on the fluid–structure interface. The velocity on the fluid–structure interface ${\stackrel{\rightharpoonup }{\mathit{v}}}_S$ is related to the motion of the flap based on the position of the hinge point,

$${\stackrel{\rightharpoonup }{\mathit{v}}}_S=\dot{\phi }\left(t\right)\times \left({\stackrel{\rightharpoonup }{\mathit{r}}}_s-{\stackrel{\rightharpoonup }{\mathit{r}}}_h\right)$$

(32)

3.2. Set-Up in OpenFOAM

The interaction of the flap with waves is simulated by the finite volume method solver overInterDyMFoam in the package OpenFOAM, an open source CFD toolbox.

In order to calibrate the BEM results, some settings in CFD should be same or similar with BEM. Therefore, the Stokes I wave model (corresponding to the linear waves in BEM) is applied for regular waves, and 100 wave components with different wave periods, wave heights, and random phases are utilized for irregular waves. Shallow water wave absorption works on both inlet and outlet.

For simulating the motion of a flap with a potential large amplitude, the overset mesh approach [21] is used. The domain is separated into two zones—a stationary main zone and a moving overset zone. The flap rotation is simulated by the motion of the whole overset zone and the cell volume weight interpolation of velocity, pressure, and volume fraction between the overlapping volume of these two zones. The flap surface is restricted by a hinged axel with a torsional spring-damper (i.e., a joint constraint releasing only Ry degree of freedom and a restraint of linear axial angular spring inputting rotational stiffness and damping) to simulate the PTO system.

A turbulence model of eddy viscosity of type standard $k$-$\epsilon$ (SKE), where $k$ and $\epsilon$ are, respectively, turbulent kinetic energy and dissipation rate, is applied to calculate the turbulent viscosity ${\mu }_t$, which will be introduced into Equation (29) for the calculation of stress tensor $\stackrel{\rightharpoonup }{\mathit{T}}$. Although the SST $K$-$omega$ model is widely applied for the simulation of wave–structure interaction [45,46], Wei et al. [18] compared the variable types of turbulence model and found that the torque on the flap due to the wave force is not sensitive to the turbulence model. Hence, the SKE model is employed in the present study, due to the relatively low computational cost [18].

In iterative processing, the pressure field is coupled to the velocity field via the PISO-SIMPLE (PIMPLE) method, a merged algorithm of the Pressure Implicit Splitting Operator (PISO) [47] and the Semi-Implicit Method for Pressure-linked Equations (SIMPLE) [48].

All other numerical solver settings were chosen based on the tutorial case of floating body in [49].

4. Model and Mesh

4.1. Study Model

The research model is a 3D surface-piercing flap, of which the $w=26 \mathrm{m}$, $b=4 \mathrm{m}$, $h=10 \mathrm{m}$, and $d_a=9 \mathrm{m}$ are referred to the Oyster 800-like model in [50]. $d=12.5 \mathrm{m}$, $I_G=9.1455\times {10}^6 \mathrm{kg}·{\mathrm{m}}^2$, and $C_{PTO}=16 \mathrm{MN}·\mathrm{m}·\mathrm{s}/\mathrm{rad}$ are given in [51]. Some parameters, $\overline{BH}=3.705 \mathrm{m}$ and $m=6\times {10}^5 \mathrm{kg}$, are estimated by the authors. $\overline{GH}=4.781 \mathrm{m}$ is calculated by the given restoring stiffness at a small rotational angle $K_0\approx F_B·\overline{BH}-mg·\overline{GH}=12.81 \mathrm{MN}·\mathrm{m}/\mathrm{rad}$ in [51].

In BEM, 40 wave cycles for the end time of simulation, $t_b=24T$ and $t_e=40T$ for the time–averaged response and CWR, and $N_c=10$ for the measurement of phases are imposed for regular waves. We used 100 peak periods as the simulation time, and $t_b=20T_p$, $t_e=100T_p$ ($T_p$ is the peak period of wave spectrum) and $N_w=100$ for wave components from 0.1 to 4.6 rad/s in irregular waves. To eliminate the effects of random phase, the solutions of equivalent pitch amplitude and CWR for five arrays of random phases have been averaged.

To reduce the computational cost in CFD, one half of the wave basin is built as the computational domain of NWT, as shown in Figure 4a, since the geometry, restriction of the flap, and boundary conditions are symmetric. The lengths of the upstream and downstream in the main zone are set as 5 times the flap width, and the width of the main zone (one half of the full model) is 3 times the flap width. The height of the initial air zone is 11 m. The distance between external cubic surface of the overset zone and the flap surface, $d_{overset}$ as shown in Figure 4b, should contain multi-layer grids. We define the number of grid layer as $N_{GL}$. A small value of $N_{GL}$ could result in the divergence of computation due to inaccurate interpolation in the overlapping region. However, for more layers (i.e., a large value of $N_{GL}$), the implement of dynamic overset grid algorithm [52] in each iteration is time consuming. The influence of the value of $N_{GL}$ will be discussed in Section 4.2.2.

Due to the expensive computational cost of long-duration modeling, the duration of simulation under regular waves is 5–13 times the legth of the wave period, according to the stability of the response. Usually, more cycles in a short wave are necessary. The method of extracting pitch amplitude results in CFD is similar to that in the time domain BEM when the response is steady. The recording period for the computation of pitch amplitude is the last 1–3 periods in regular waves.

4.2. Mesh

4.2.1. BEM Grids

Only the grids on the wet surface should be generated in BEM. The quadrilateral grids are applied, and the element sizes are almost uniform. Three different element sizes (coarse, medium, and fine meshes), as shown in

Table 2. The mesh statistics and computational cost (on 1 processor) for three configurations of grids in BEM.

Mesh	Mean Cell Size (m)	# of Wet Surface Cells above the Hinge	# of Divisions on the Bottom Round	# of Cells	Computational Cost
Coarse	≈0.45	16 × 52 × 23	12	3.9K	$\approx 5 \mathrm{h}$
Medium	≈0.35	20 × 65 × 30	16	6.4K	$\approx 19 \mathrm{h}$
Fine	≈0.3	24 × 75 × 35	20	8.8K	$\approx 49 \mathrm{h}$

, are tested for mesh independence by comparing the results of the hydrodynamic coefficients, which are calculated by Nemoh program on a desktop of 3.4 GHz CPU with 16 processors and 32GB RAM. Unfortunately, Nemoh does not provide the parallel computation with multi-processors. The normalized root mean square error (NRMSE) is utilized to estimate the agreement of the results of a certain mesh configuration (coarse mesh or medium mesh) with the results of fine mesh. In Figure 5, each hydrodynamic coefficient seems to be insensitive to the mesh size. Finally, the medium grids, as shown in Figure 6, are selected for the case studies. The calculated value of added inertia torque at infinite frequency $\mu (\infty )=28.94\times {10}^6 \mathrm{kg}·{\mathrm{m}}^2$.

4.2.2. CFD Grids

The mixed grids consisting of cut-cell elements [53] for the main zone and hexahedral elements for the overset zone are employed in CFD, as shown in Figure 7. To understand the generation of cut-cell grids, the meshes in the main zone are created as follows:

• The overall Cartesian meshes with almost uniform size are constructed;
• The transition region and refined region (near the water surface and overset zone) are cut off once, and their grid size bisects;
• Cut off the meshes in the refined region again, so the grid size of refined region is a quarter of outermost Cartesian region.

Thus, the grids near the overset zone and the water surface are fine for the accurate simulation of flap motion and wave propagation. The length of refined region in x-direction is set to more than $3 \mathrm{h}$ to guarantee a larger coverage area than the range of flap motion. The cell size in the overset zone should be approximately same as the grid size in refined region of the main zone, as this will diminish the interpolation error between these two zones [23].

In the overset zone, the boundary layer (very fine mesh near the structure surface), which affects the calculation of tangential fluid shear stress on the structure surface, is neglected in this work, because the shear stress is not essential relatively to the normal stress for the prototype OWSC [18]. The simulations in [54] indicate that the shear forces on the fixed flap calculated are less than 1% of the total surge force.

Three different element sizes (the number of grid layer between external cubic surface of overset zone and the flap surface is uniformly set to $N_{GL}=10$), as shown in

Table 3. The mesh statistics and computation cost (on 16 processors) for three configurations of grids in CFD.

Mesh	Main Zone			Overset Zone			Total # of Cells	Computational Cost
	# of Water Surface	# of Top & Bottom	# of Cells	# of Flap Surface above Hinge	# of Divisions on Flap Rounds	$\# \mathbf{of} \mathbf{Cells} ({\mathit{N}}_{\mathit{G}\mathit{L}} = 10)$
Coarse	468 × 140	117 × 35	0.94M	12 × 24 × 19	16	34.1K	0.98M	$\approx 5 \mathrm{h}$
Medium	636 × 192	159 × 48	1.96M	16 × 33 × 26	20	58.7K	2.02M	$\approx 28 \mathrm{h}$
Fine	868 × 264	217 × 66	3.83M	20 × 44 × 34	24	88.5K	3.92M	$\approx 76 \mathrm{h}$

, are tested for mesh independence by comparing the wave elevation at 20 m in front of the flap and the time–history total fluid moment $M_{fluid}\left(t\right)$ (for 55 s simulation time) on a fixed flap under a regular wave ($A_i=1 \mathrm{m}$ and $T=17.5 \mathrm{s}$) in CFD. The mean cell sizes of the refined regions for coarse, medium and fine mesh configurations are 0.56 m, 0.41 m and 0.30 m, respectively. Because of the limitation in BEM, the wave height is not allowed very high (see Section 2). With the purpose of calibrating for BEM, the unit-amplitude wave is choosen. The selected wave period is a representative period close to the resonant period of the flap (see Section 5) with greatest interests in this study.

The NRMSE is applied to estimate the agreement of the results of current mesh with the results of fine mesh. A more than 5% NRMSE related to the results are found in Figure 8 from a comparison between coarse and fine grids. However, the result for the medium grid is very close to that of the fine grid, in which the total number of elements is almost double of the medium grid. In addition, the computation for medium mesh is much faster than fine mesh (see Table 3). Accordingly, the medium grids with acceptable accuracy and efficient computation is selected for the case studies.

To assess the influence of the grid layers between the external cubic surface of overset zone and flap surface, four different numbers of grid layer ($N_{GL}=5$, $N_{GL}=10$, $N_{GL}=15$ and $N_{GL}=20$, as shown in Figure 9) are tested to assess the effects of overset zone dimension. The solver with the current settings cannot produce the converged result of flap motion in the case of 5 layers ($N_{GL}=5$), possibly caused by the amplified interpolation error (mentioned in Section 4.1).

The comparison of time–history responses of a flap without a PTO system ($C_{PTO}=0$ and $T_{PTO}=0$) under a regular wave ($A_i=1 \mathrm{m}$ and $T=17.5 \mathrm{s}$) for 10 layers, 15 layers and 20 layers exhibits a negligible difference, as shown in Figure 10. Similarly, a slight influence of overset zone size can also be found in [21].

The computational costs of three different $N_{GL}$ on a desktop (see Section 4.2.1) using 16 processors are approximately 90 h ($N_{GL}=10$), 95 h ($N_{GL}=15$) and 105 h ($N_{GL}=20$), respectively, for a 55 s simulation time. Because of a shorter time consumption of computation in the smaller number of grid layer and insignificant influence on the motion of flap, $N_{GL}=10$ are adopted. Figure 11 illustrates the selected mesh model in CFD for the case studies.

5. Free Decay Analysis

In the free decay analysis of a flap with an initial inclination under still water, $M_E\left(t\right)=0$, $\eta \left(t\right)=0$, $T_{PTO}=0$, $\dot{\phi }\left(0\right)=0$, $\ddot{\phi }\left(0\right)=0$ and a nonzero $\phi \left(0\right)$ are set.

The response of free decay in still water is studied with different initial angles. According to the time–history curve of response extracted from OpenFOAM, as shown in Figure 12, the natural period can be estimated by measuring the time distance between the first and second valley values of the response. Without friction ($T_{PTO}=0$), the averaged periods of response in cases with different initial angles and either an undamped flap ($C_{PTO}=0$) or a damped flap ($C_{PTO}=16 \mathrm{MN}·\mathrm{m}·\mathrm{s}/\mathrm{rad}$) are, respectively, 17.57 s or 19.15 s. The measured period of response lengthens when configuring PTO damping. In actuality, the natural period without additional declaration usually defaults as the natural period of an undamped flap without friction for the judgment of resonance, i.e., 17.57 s.

In BEM, the actual water surface, affected by the moving flap’s radiating wave, is difficult to determine. We assume that the water surface is always horizontal in free decay analysis, i.e., ${\alpha }_{\eta }=1$ in Equation (3). The free decay responses for different drag coefficients in BEM are calibrated with CFD, as shown in Figure 13. It is noticed by observation from the periods of these signals that the case of a 30 deg. intial angle matches relatively well with OpenFOAM, as shown in Figure 13b. However, a longer period with a small initial angle (15 deg.) and a shorter period with a large initial angle (45 deg.) are observed in Figure 13a,c, respectively. The authors believe that these inaccurate results in BEM are caused by the assumption of a horizontal water surface.

To understand the distinctions between BEM and CFD in free decay analysis, Figure 14 shows some typical frames in the CFD of a flap without a PTO system at three different initial angles. The water surface near the flap is elevated in several frames, and the elevation seems quantitatively unpredictable. Theoretically, the values of the radiation damping moment, the restoring moment, and the drag moment are disturbed by the change in the realistic wet surface, but cannot be accurately considered in BEM.

We scan the values of $C_d$ from 2 to 6 to find a proper drag coefficient with a minimum NRMSE of BEM against CFD in the case of a 30 deg. intial angle. As shown in Figure 15a, $C_d=3.4$ leads to the lowest NRMSE. The comparison of time–history free decay response between BEM ($C_d=3.4$) and CFD is illustrated in Figure 15b. For engineering applications of the free decay analysis based on BEM, an intermediate initial angle (e.g., 30 deg.) with a selection of $C_d=3.4$ for the current model is more appropriate.

For a 40 s simulation time of a 30 deg. initial angle free decay response, the computational cost on the same desktop (see Section 4.2.1) is about 30 s for the Python code based on BEM (on 1 processor), except the elapsed time of extracting hydrodynamic coefficients for one time from Nemoh. However, it takes about 40 h for the CFD simulation with 16 processors in OpenFOAM.

6. Under Regular Waves

In this section, the responses of the flap under regular waves are studied. The range of wave period is from 5 s to 23 s. The unit-amplitude waves and the varying wave heights corresponding to the uniform wave exciting moment are considered, respectively. The range of PTO damping is from 0 to 64 $\mathrm{MN}·\mathrm{m}·\mathrm{s}/\mathrm{rad}$, and the range of PTO friction is from 0 to 3 $\mathrm{MN}·\mathrm{m}$. However, the expanded range of PTO parameters can be included in some study cases with additional illustrations.

6.1. Calibration of BEM with CFD

There are two undecided factors in the computation of the response under waves in BEM: wave surface corrected factor ${\alpha }_{\eta }$ and drag coefficient $C_d$. We calculate the NRMSE of the pitch amplitude in BEM with respect to the CFD solutions (the results of equivalent pitch amplitude in CFD are extracted according to Equation (22)) of both the undamped ($C_{PTO}=0$) and damped ($C_{PTO}=16 \mathrm{MN}·\mathrm{m}·\mathrm{s}/\mathrm{rad}$) flaps, without friction and under unit-amplitude regular wave conditions ($T=5 ~ 22.5 \mathrm{s}$ with interval 2.5 s).

Figure 16 indicates that ${\alpha }_{\eta }=0.16$ and $C_d=5.4$ are well calibrated for the time domain BEM results in the study’s model. It should be noted that the corrected wave surface is not a realistic water surface, and it only numerically influences the computations of restoring moment and drag moment for the sake of better accuracy.

With the calibrated factors, a good agreement of BEM with CFD in the whole range of wave periods from 5 s to 23 s is shown in Figure 17. However, BEM slightly overestimated the response in short periods due to the constant drag coefficient, which was Keulegan–Carpenter ($KC$) and Reynolds ($Re$) number-dependent [55]. Using a larger value of the non-constant drag coefficient would seem to be necessary for a high $Re$ number under short waves in future work. Given its acceptable accuracy and efficient computation, the calibrated BEM will be adopted in the next subsections.

Taking an example of a simulation under a unit-amplitude regular wave, the computational cost to run the Python code based on BEM (on 1 processor) on the same desktop (see Section 4.2.1) is about 4 min for 40 wave cycles. However, it takes about 190 h for 6 wave cycles in OpenFOAM with 16 processors.

6.2. When Does the Maximum Pitch Occur?

The responses of the flap under regular waves of uniform amplitude ($A_i=1 \mathrm{m}$) and uniform wave exciting moment (e.g., $M_e\left(\omega \right)A_i=10 \mathrm{MN}·\mathrm{m}$), assessed by varying the wave amplitude for each period, are studied to determine the maximum pitch amplitude.

Figure 18a shows the wave amplitudes for different wave periods with respect to $M_e\left(\omega \right)A_i=10 \mathrm{MN}·\mathrm{m}$. When the wave period is longer than 15 s, the wave amplitude is greater than 1 m (the freeboard height $h_f=1 \mathrm{m}$). Based on the assumption of a linear relationship between wave amplitude and wave exciting moment amplitude, the wave torque of a high wave in BEM may be imprecise, due to wave overtopping. Figure 18b shows a comparison of the pitch amplitude of a flap without a PTO system ($C_{PTO}=0$ and $T_{PTO}=0$) between BEM and CFD, confirming the reasonable precision of BEM related to the slight influence of wave overtopping under such intermediate waves (the wave amplitude slightly exceeds the freeboard height).

Figure 19 shows the pitch amplitude results of a no-friction flap with various values of PTO damping $C_{PTO}$ under regular waves of uniform amplitude and exciting moment amplitude. The maximum pitch amplitudes (solid point on each curve, between 16 s and 19 s in Figure 19a) with different $C_{PTO}$ values occur within a small range of wave periods from 17.6 s to 18.3 s, close to the natural period (17.57 s is obtained from the free decay analysis in CFD). This means that the flap resonates near the natural period. However, there are no obvious maximum values of pitch amplitude near the natural period with uniform wave exciting moment amplitude, as shown in Figure 19b. Therefore, the resonant conditions of regular waves, in which an incident wave near the natural period produces the maximum response over wave periods, should be restricted by employing the same amplitude of incident wave elevation, rather than the same amplitude of wave torque.

Figure 20a shows the effects of friction (with an expanded range) on the period of maximum pitch amplitude of a damped flap ($C_{PTO}=16 \mathrm{MN}·\mathrm{m}·\mathrm{s}/\mathrm{rad}$). The increasing value of friction breaks the resonance rule (the maximum response occurs near the natural period) by moving the period of maximum pitch amplitude towards the period of peak wave exciting moment, which is near 7 s for the current model, as shown in Figure 20b.

6.3. Phase Characteristics Near and Far from Resonance

We have studied the phase angle of angular velocity relative to the wave exciting moment (${\theta }_1$ in Equation (24)) of a flap with different PTO damping and friction values under unit-amplitude regular waves. A negative value of ${\theta }_1$ means that the angular velocity lags behind the wave exciting moment, and a positive value denotes the advance of angular velocity. As mentioned in Section 2.4, the velocity is in phase with the wave torque, i.e., ${\theta }_1=0$ represents the occurrence of resonance in a one-degree-of-freedom structure with constant inertia, stiffness, and damping. It is evident from Figure 21 that the values of ${\theta }_1$ are not zero near the natural period (orange circle). Therefore, the standard of ${\theta }_1=0$ for the assessment of resonance in the hydrodynamics of a flap is unreliable, due to the strong nonlinearities of the restoring and drag moments.

A higher value of $C_{PTO}$ or $T_{PTO}$ generally produces a trend of diminishing hysterical angles of angular velocity (absolute value of a negative ${\theta }_1$) within a shorter period than the natural period, and it causes increases in hysteretic angle or reductions in advanced angle (value of a positive ${\theta }_1$) over a longer period. However, the intersections of these phase vs. period curves always occur near the natural period.

To understand the phase characteristics near to or far from resonance, Figure 22 shows the time–history curves of the ratios of the wave exciting moment and angular velocities for a no-friction flap with various PTO damping during one wave cycle, over periods of 12 s, 17.5 s (near the natural period), and 23 s. The ratios are defined as below with a steady response,

$$R_E\left(t\right)=\frac{M_E\left(t\right)}{\max \left(\left|M_E\left(t\right)\right|\right)}$$

(33)

$$R_v\left(t\right)=\frac{\dot{\phi }\left(t\right)}{\max \left(\left|\dot{\phi }\left(t\right)\right|\right)}$$

(34)

where $R_E\left(t\right)$ and $R_v\left(t\right)$ are, respectively, the ratio of the wave exciting moment and the ratio of the angular velocity.

With various values of PTO damping, there are obviously different phases of angular velocity relative to the wave exciting moment in a wave period of 12 s, but the distinctions are slight over a wave period of 17.5 s, close to the natural period. Under a long wave of period 23 s, the angular velocity of an undamped flap ($C_{PTO}=0$) precedes the wave exciting moment slightly, shown as the red dotted line in Figure 22c, but the angular velocities are almost in phase with the wave torque for various values of PTO damping from 16 to 64 $\mathrm{MN}·\mathrm{m}·\mathrm{s}/\mathrm{rad}$.

6.4. Relationship of Maximum CWR with Resonance

To understand the relationship between the maximum CWR and the resonance, Figure 23a shows the CWRs of a no-friction flap with various degrees of PTO damping under unit-amplitude regular waves. The period of maximum CWR is apparently shorter than the natural period 17.57 s, and the difference increases almostly with the increase in PTO damping. Figure 23b,c exhibits the curves of time–averaged capturing power vs. period and angular velocity amplitude vs. period, respectively. For each value of PTO damping, the period of maximum capturing power is identical with the period of maximum angular velocity amplitude. However, the period of maximum CWR is closer to the period of maximum angular velocity amplitude than to the period of maximum angle amplitude (near the natural period), as listed in

Table 4. The periods of peak values of various items under unit-amplitude regular waves.

$\mathbf{PTO} \mathbf{Damping} {\mathit{C}}_{\mathit{P}\mathit{T}\mathit{O}}$ $(\mathbf{M}\mathbf{N}\mathbf{·}\mathbf{m}\mathbf{·}\mathbf{s}/\mathbf{r}\mathbf{a}\mathbf{d})$	The Periods of Max.				Difference of the Peak Periods
	$\mathbf{Angle} \mathbf{Amplitude} {\mathit{\phi }}_{\mathit{a}\mathit{m}\mathit{p}}$	CWR	$\mathbf{Capturing} \mathbf{Power} {\mathit{P}}_{\mathit{c}}$	$\mathbf{Angular} \mathbf{Velocity} \mathbf{Amplitude} {\dot{\mathit{\phi }}}_{\mathit{a}\mathit{m}\mathit{p}}$	$\mathbf{From} \mathbf{Max}. \mathit{C}\mathit{W}\mathit{R}$$\mathbf{to} \mathbf{Max}. {\mathit{\phi }}_{\mathit{a}\mathit{m}\mathit{p}}$	$\mathbf{From} \mathbf{Max}. \mathit{C}\mathit{W}\mathit{R}$$\mathbf{to} \mathbf{Max}. {\dot{\mathit{\phi }}}_{\mathit{a}\mathit{m}\mathit{p}}$
16	18.1 s	14 s	14 s	14 s	22.65%	0
32	18.3 s	10 s	13 s	13 s	45.36%	23.08%
48	18.2 s	8 s	10 s	10 s	56.04%	20.00%
64	18 s	8 s	9 s	9 s	55.56%	11.11%

. This implies that the highest CWR does not occur simultaneously with resonance.

7. Under Irregular Waves

7.1. Wave Spectrum and Modification

Goda [56] introduced the JONSWAP spectrum $S\left(\omega \right)$, formulated from the significant wave height and the peak period,

$$S\left(\omega \right)=16{\beta }_J{\pi }^4\frac{H_s^2}{T_P^4}\frac{1}{{\omega }^5}exp\left(-\frac{20{\pi }^4}{T_p^4{\omega }^4}\right)·{\gamma }^{exp\left[-{\left(T_p\omega /2\pi -1\right)}^2/2{\sigma }^2\right]}$$

(35)

where $H_s$ is the significant wave height and $T_p=2\pi /{\omega }_p$ is the peak period, and

$${\beta }_J\approx \frac{0.06238\left(1.094-0.01915ln \gamma \right)}{0.230+0.0336\gamma -0.185{\left(1.9+\gamma \right)}^{-1}}$$

(36)

$$\sigma =\{\begin{array}{l} 0.07, \omega \le {\omega }_p\\ 0.09, \omega >{\omega }_p\end{array}$$

(37)

where $\gamma$ is the peak enhancement factor, the range of which is 1–7, with a mean value of 3.3 [57].

When a severe wave in an offshore region propagates into the nearshore area, the energy will be attenuated by wave-breaking and seabed friction. In shallow water, the original spectral values should be modified to describe the reduction in wave energy according to water depth. The modified wave spectrum $S_d\left(\omega \right)$ [58] is given by

$$S_d\left(\omega \right)=S\left(\omega \right)\cdot \frac{k^{-3}\cdot \partial k/\partial \omega }{k_{\infty }{}^{-3}\cdot \partial k_{\infty }/\partial \omega }$$

(38)

where $k_{\infty }$ is the wave number at infinite depth.

Figure 24a shows a comparison between the original and modified JONSWAP spectra for an irregular wave ($H_s=6 \mathrm{m}$ and $T_p=17.5 \mathrm{s}$) in 12.5 m-deep water, and Figure 24b shows the time–history elevation (generated by Equation (1)) of a modified irregular wave with the specific array of random phase $\varphi \left({\omega }_n\right)$ used for the validation of BEM. Due to the assumption of linear relationship between wave exciting moment and wave amplitude in BEM, the significant wave height is also not allowed very high. After correcting the spectral values, the highest wave elevation is a little more than 1 m (the freeboard heigh of the flap). It implicates that a higher significant wave height might cause the overestimated results because of ignoring the wave overtopping in BEM.

7.2. Validation of BEM against CFD

The wave surface corrected factor ${\alpha }_{\eta }=0.16$ and drag coefficient $C_d=5.4$ described in Section 6.1 are continuously applied in irregular waves. Figure 25 shows that the time–history response in BEM matches well with that in CFD, using the same wave components and the identical array of random phases within the time–history wave elevation of Figure 24b. Due to the relatively small wave elevation in the modified wave spectrum, the linear assumption of wave exciting moment vs. wave amplitude is still applicable for such an intermediate irregular wave. Considering the linear assumption of wave exciting moment vs. wave amplitude in BEM, a relatively small significant wave height $H_s=2 \mathrm{m}$ with lower nonlinerities will be applied in the next case studies.

For the 120 s simulation time under an irregular wave, running the Python code based on BEM (on one processor) spends about 4 min on the same desktop (see Section 4.2.1), and the elapsed time for CFD with 16 processors is about 192 h.

7.3. When Does the Maximum Equivalent Pitch Occur?

Under modified irregular waves with the same significant wave height $H_s=2 \mathrm{m}$, the equivalent pitch amplitudes ${\phi }_{eq}$ of each peak period for various values of PTO damping $C_{PTO}$ are calculated. The maximum equivalent pitch amplitude of a flap without a PTO system ($C_{PTO}=0$ and $T_{PTO}=0$) appears in a peak period longer than the natural period, shown via the black dotted lines in Figure 26a, but the period of maximum ${\phi }_{eq}$ decreases significantly with the increase in $C_{PTO}$. This reveals that the maximum equivalent pitch amplitude cannot consistently occur near the natural period, which differs from the results under regular waves. This can be explained via the response under an irregular wave, achieved from the linear combination of the results derived using multi-components of regular waves.

In a specifc JONSWAP spectrum, the proportion of longer waves relative to the peak period is less than that of the shorter waves due to the asymmetry of spectral values, and the transfer function (i.e., the response for unit-amplitude regular waves in Figure 19a) under shorter waves thus dominates the synthetic solution. Figure 19a shows that the difference between the transfer function in a short period and its maximum value near the natural period reduces with the increase in $C_{PTO}$. In other words, the short-period transfer function for a large $C_{PTO}$ contributes relatively more to the equivalent pitch amplitude than a small $C_{PTO}$. That is why the maximum equivalent pitch amplitude for a larger $C_{PTO}$ appears within a shorter peak period.

As shown in Figure 26b, the impact of the friction on the peak period of maximum equivalent pitch amplitude is similar to that under regular waves, as discussed in Section 6.2.

7.4. Phase Characteristics near the Natural Period

Figure 27 shows the time–history wave exciting moment and angular velocities for different values of PTO damping and friction when the peak period is close to the natural period, $T_p=17.5 \mathrm{s}$. It is found that the angular velocity lags behind the wave exciting moment (orange arrows), and the phase difference reduces with the increase in either PTO damping $C_{PTO}$ or friction $T_{PTO}$, as opposed to the near-constant phase difference (hardly affected by $C_{PTO}$ and $T_{PTO}$) near the natural period under regular waves, shown in Figure 21.

7.5. CWR vs. Peak Period

Figure 28 shows the curves of CWR vs. peak period for various values of PTO damping. Similar to the results under regular waves, the maximum CWR here appears in a peak period lower than the peak period of the maximum equivalent pitch amplitude, as quantitatively compared with Figure 26a, and it declines with the increase in PTO damping.

8. Conclusions

Based on BEM, the time domain dynamic equation of a 3D bottom-hinged OWSC in shallow water, considering the nonlinear restoring moment, drag moment, and friction, is solved by a Python code importing the hydrodynamic coefficients obtained from NEMOH. The wave surface corrected factor ${\alpha }_{\eta }$ and drag coefficient $C_d$ are calibrated with the CFD results in OpenFOAM. For an Oyster 800-like flap, the factors ${\alpha }_{\eta }=1$ and $C_d=3.4$ are applicable for free-decay analysis in still water, and ${\alpha }_{\eta }=0.16$ and $C_d=5.4$ are suitable in regular and irregular waves.

In the analysis of free decay, an intermediate angle—30 deg. in this study— is recommended as the initial inclination in BEM for the estimation of the natural period, because a small initial angle overestimates the natural period and a large initial angle underestimates it.

Under regular waves, one criterion when judging resonance is the period of maximum pitch amplitude with uniform wave amplitude, rather than the uniform amplitude of the wave exciting moment. A large value of friction destroys the resonant condition and drives the period of maximum pitch closer to the period of peak wave torque. Perfect resonance cannot be achieved, because a hysteresis of angular velocity relative to the wave exciting moment near the natural period is always present, and is almost independent of PTO damping and friction. Maximum CWR does not coexist with resonance but relatively near the largest angular velocity, the period of which is lower than that of maximum pitch (resonance).

Under irregular waves of the same significant height, the peak period of maximum equivalent pitch amplitude reduces with the increase in PTO damping. This implies that resonance cannot exist stably at a peak period near the natural period. The angular velocity also lags behind the wave torque within the current regime of PTO damping and friction when the peak period is close to the natural period. Here, one discrepancy from the results under regular waves is that the hysterical angle of angular velocity relative to wave torque relies on either PTO damping or friction. In addition, the curves of CWR vs. peak period for different PTO damping exhibit a similar trend to the results under regular waves.