Numerical Model of Constrained Wave Energy Hyperbaric Converter under Full-Scale Sea Wave Conditions

Abstract

A 2D numerical investigation of the power absorption of a constrained wave energy hyperbaric converter (WEHC) under full-scale sea wave conditions is presented. A fully non-linear numerical model DualSPHysics, based on the coupling of a smoothed particle hydrodynamics (SPH) fluid solver with a multibody dynamics solver, is used to model the interaction between wave and WEHC sub-systems. The numerical model was first validated against experimental data for a similar device, with a good accordance between PTO position and velocity. The model is then employed to study the hydrodynamics of a constrained WEHC considering several sea states, different hydraulic power take-off (PTO) damping and breakwater geometries. It is observed that the capture width ratio (CWR) is particularly sensitive to variations in the PTO damping, although the CWR absolute maximum is less sensitive considering mild variations applied to the PTO damping. Both wave height and wave period have an important effect on the CWR. The breakwater geometry is also essential for the performance of the WEHC, with a decrease in maximum CWR of about 15% for porous breakwater. These results are necessary to understand the full-scale behaviour of WEHC.

1. Introduction

A wave energy hyperbaric converter (WEHC) is a top-mounted pitching point wave energy converter (WEC) device, which mainly consists of three different sub-systems: floater, lever arm and a complex hydraulic power take-off (PTO) system [1,2]. The floater is connected by a lever arm to a hinge in the support platform where the rotational motion of the lever arm is damped by a PTO system to absorb power in operational conditions. The PTO system consists of a hydraulic pump linked to lever arm, a hyperbaric chamber and a hydropneumatic accumulator which acts as a power storage to supply a Pelton turbine coupled to a generator, converting mechanical power into electricity. The main components of the WEHC are sketched in Figure 1.

The mechanical constraints imposed by each connection among the three sub-systems, i.e., universal joint between lever arm and floater, revolute joints between lever arm and support platform and the hydraulic pump (Figure 1), have an important influence in the power absorption of WEHC. However, most of the previous studies, e.g., [2,3,4,5,6,7], have simplified the modelling of WEHC, neglecting the description of the mechanical constraints among the various sub-systems. Additionally, the PTO system is often modelled as a linear spring and damper [8]. There are two main challenges to resolve the non-linear dynamics of WEHC and to simulate sea waves in an accurate way: (i) the numerical model should predict correctly the imposed wave motion on the floater; (ii) the non-linear feedback of floater motion on wave field under several sea states and operation modes. This paper addresses the above challenges by implementing the mechanical constraints that simulate the non-linear interaction between each sub-system and understanding its effects on the power absorption under full-scale sea wave conditions. For that purpose, the objective is providing stronger evidence on the real motion of the constrained WEHC for different hydraulic PTO dampings and breakwater geometries.

Efforts to model similar point-absorber-type WECs have been made by using both mesh-based and mesh-free computational fluid dynamic (CFD) methods. Mesh-based methods have been widely adopted to solve the governing equations by several codes such as COMFLOW [9], COMET [10] and AMAZON-SC [11], as well as in-house CFD codes [12,13,14]. Recently, [15] applied the OpenFOAM code to simulate a point-absorber-type WEC known as Wavestar, showing that the fully non-linear numerical model accurately reproduces the observed laboratory experiments of wave–floater interaction. To include PTO models, [8] developed in OpenFOAM a simulation platform for WECs by coupling a CFD-based numerical wave tank with PTO models. By comparing with other models based on the linear and partially non-linear boundary element method, their coupled model (including both the CFD and the PTO models) clearly showed that a fully non-linear numerical model is essential to resolve the wave–WEC interaction. The above listed research, among others, clearly highlighted that to accurately simulate WECs, there is a need to develop advanced models to resolve non-linear waves, physical processes of wave–WEC interaction and power absorption. However, the mesh-based method poses difficulties to include mechanical constraints. Such a method also implies solving additional equations for free-surface capturing and requires commonly expensive and complicated mesh moving algorithms to solve WECs motion. Additionally, solving mechanical constraints often requires additional complex mathematical transformations to solve non-linear complementarity problems, singularities, non-uniqueness, non-penetration and inequality constraints, which can be time-consuming when compared with solving wave–WEC interaction. Moreover, some mechanical constraints are not, in general, continuous, and such applications are not straightforward for a mesh-based method considering the ordinary differential equations that require impracticable mesh [16,17,18].

A promising alternative to overcome these mesh-based limitations is the mesh-free method. This method is favoured for modelling non-linear sea waves, since a mesh is not required when solving the Lagrangian form of the governing equations. Among the more up-to-date mesh-free methods, the smoothed particle hydrodynamics (SPH) has been widely adopted to study complex wave–WEC interaction. This method has been properly validated to study point-absorber-type WECs (e.g., [19,20,21]). SPH is able to simulate unsteady and non-linear flows, extreme deformations and complex topological evolutions, since the free surface is captured naturally [22,23]. However, given the dynamic nature of the SPH method, mechanical constraints are also difficult to implement in SPH form. Regarding this limitation, refs. [24,25,26] have successfully developed a coupled numerical model to simulate wave–WEC interactions. This coupling approach consists of an augmented SPH model of DualSPHysics [27] with the multibody solver of Project Chrono [28]. This is specifically designed to support the simulation of complex fluid–structure interaction considering its mechanical constraints. Hence, the SPH is used to model wave–WEC interaction, and multibody solver is used to model the flap–mechanical constraints interaction. The top-mounted pitching point absorber has also been modelled using DualSPHysics [29]. This device has a working principle similar to WEHC, however, in [29], a mechanical constraint between lever arm and floater and the effects of breakwater geometry were not considered.

Following the success of the SPH coupling approach for WEC modelling, this paper presents an analysis of the power absorption of a constrained WEHC under full-scale sea wave conditions. The details of the numerical model (coupling procedure between DualSPHysics and Project Chrono) are presented in Section 2. The numerical model validation is performed in Section 3. In Section 4, the characteristics of the capture width ratio (CWR) and PTO system are studied in both optimal and sub-optimal conditions, and the effects of wave height and wave period are analysed. Finally, Section 5 summarises the main findings.

2. Numerical Model

The coupling of an SPH fluid solver of DualSPHysics [30] with a multibody dynamics solver of Project Chrono [28] is used. This allows for efficient and accurate modelling of complex fluid–multibody interactions with any mechanical constraints, as described in [25] and [26]. The fluid–multibody interactions (e.g., non-linear waves’ interaction with floaters and the interaction between floater, lever arm and PTO system) are mainly solved in three steps [26] that are described below.

In the first step, DualSPHysics solves the following continuity and Navier–Stokes equations in SPH form [31] with the density diffusion term [32]:

$$\frac{d{\rho }_i}{dt}=\sum _j{m}_j{\overrightarrow{v}}_{ij}·{\nabla }_i{W}_{ij}+2\delta h{c}_0\sum _j\left({\rho }_{ji}^T-{\rho }_{ij}^H\right)\frac{{\overrightarrow{r}}_{ij}·{\nabla }_i{W}_{ij}}{{\overrightarrow{r}}_{ij}^2}\frac{m_j}{{\rho }_j}$$

(1)

$$\frac{d{\overrightarrow{v}}_i}{dt}=-\sum _j{m}_j\left(\frac{P_i+P_j}{{\rho }_i{\rho }_j}+{\Pi }_{ij}\right){\nabla }_i{W}_{ij}+\overrightarrow{g}$$

(2)

where $\rho$ is the density; t is the time; m is the mass; ${\overrightarrow{v}}_{ij}={\overrightarrow{v}}_i-{\overrightarrow{v}}_j$, with ${\overrightarrow{v}}_i$ and ${\overrightarrow{v}}_j$ being, respectively, the velocity vector of the target particle i and neighbouring particle j; ∇ is the gradient operator; P is the pressure; $c_0$ is the reference speed of sound; g is the gravitational acceleration; ${\rho }_{ji}^T$ and ${\rho }_{ij}^H$ are, respectively, the total and hydrostatic part of the density according to [32]; and $W_{ij}=W\left({\overrightarrow{r}}_{ij},h\right)$ is the kernel function, which depends on the distance ${\overrightarrow{r}}_{ij}={\overrightarrow{r}}_i-{\overrightarrow{r}}_j$ between i and j particles and characteristic smoothing length h. In this work, the Quintic Wendland kernel suggested by [33] is employed. $\delta$ is the delta-SPH coefficient that controls the intensity of the diffusive term. This coefficient is recommended for applications with open and moving solid boundaries, where a non-complete kernel is inevitably present, in order to reduce numerical diffusion [30]. In a strong dynamic case, this coefficient contributes with a force that may be several orders of magnitude smaller than the pressure and viscous terms, not contributing to a significant degradation of the solution. In Equation (2), $\Pi$ is the viscous term described using the artificial viscosity proposed by [31].

In order to complete the set of governing equations, it is necessary to define an equation for pressure. A trivial implementation is to use an equation of state to relate the pressure to a slight variation in density, so-called Weakly Compressible SPH (WCSPH). This method requires a very small time step in order to solve the sound speed and depends on the delta-SPH term to stabilise the density field. The following Tait equation is used in DualSPHysics to determine fluid pressure based on particle density [22]:

$$P_i=\frac{c_0^2{\rho }_0}{\gamma }\left[{\left(\frac{{\rho }_i}{{\rho }_0}\right)}^{\gamma }-1\right]$$

(3)

where ${\rho }_0=1000$ kg/m${}^3$ is the reference density, $\gamma =7$ is the polytropic constant and $c_0$ is the sound speed at ${\rho }_0$, given as $c_0=c\left({\rho }_0\right)=\sqrt{\partial P/\partial \rho }{|}_{{\rho }_0}$ [30].

In DualSPHysics, a rigid-body dynamic is calculated based on the basic Newton equations, and the discretisation consists of summing the contributions from each SPH particle, as

$$M\frac{d\overrightarrow{V}}{dt}=\sum _{k\in b}{m}_i\frac{d{\overrightarrow{v}}_k}{dt}$$

(4)

$$I\frac{d\overrightarrow{\Omega }}{dt}=\sum _{k\in b}{m}_k\left({\overrightarrow{r}}_k-{\overrightarrow{R}}_0\right)\times \frac{d{\overrightarrow{v}}_k}{dt}$$

(5)

where M is the mass and I is the inertial tensor of the rigid body b, ${\overrightarrow{R}}_0$ is the position of the centre of mass of a rigid body, $\overrightarrow{V}$ is the linear velocity and $\overrightarrow{\Omega }$ is the angular velocity. These vectorial quantities are computed at every time step by summing the contributions of each boundary particle k, belonging to body b.

In the second step, the linear $d\overrightarrow{V}/dt$ and angular $d\overrightarrow{\Omega }/dt$ acceleration to be applied in the centre of mass of a rigid body b are transferred to the Project Chrono library. During this time step, this library updates the motion, considering the given mechanical constraints and using the multibody dynamic model. The configuration of a rigid multibody system is described by using generalised coordinates $\overrightarrow{q}={[{\overrightarrow{R}}^{\mathrm{T}}+{\overrightarrow{\Theta }}^{\mathrm{T}}]}^{\mathrm{T}}$, where $\overrightarrow{R}$ is the translational and $\overrightarrow{\Theta }$ is the rotational coordinates that define each body in the system frame [16].The dynamics of rigid bodies are characterised by a system of two differential algebraic equations that relate the time derivative of generalised coordinates and velocity through a linear transformation and the equation of the force balance that ties the inertial forces to the applied and constraint forces in the following form [18]:

$$\frac{d\overrightarrow{q}}{dt}=\overrightarrow{L}\left(\overrightarrow{q}\right)\overrightarrow{V}$$

(6)

$$M\left(\overrightarrow{q}\right)\frac{d\overrightarrow{V}}{dt}={\overrightarrow{F}}_a(\overrightarrow{V},\overrightarrow{q},t)+\sum _{j\in \Phi }\left({\gamma }_{j,n}{\overrightarrow{D}}_{j,n}+{\gamma }_{j,u}{\overrightarrow{D}}_{j,u}+{\gamma }_{j,w}{\overrightarrow{D}}_{j,w}\right)$$

(7)

where $\overrightarrow{V}$ is the velocity, M is the inertia matrix, ${\overrightarrow{F}}_a$ is the applied forces, $\overrightarrow{L}$ is the Jacobian matrix, and ${\gamma }_{j,n}{\overrightarrow{D}}_{j,n}$, ${\gamma }_{j,u}{\overrightarrow{D}}_{j,u}$ and ${\gamma }_{j,w}{\overrightarrow{D}}_{j,w}$ are the constraint forces at contact point j in normal n and tangential u and w direction [34]. Equations (6) and (7) are solved using a variable time step in Project Chrono [34,35,36,37].

In the third step, position, linear and angular velocity of the rigid body are transferred back to the DualSPHysics to update the particles that form the rigid body with the information transferred from the Project Chrono library (i.e., $\overrightarrow{V}$, $\overrightarrow{\Omega }$ and ${\overrightarrow{R}}_0$). The velocity of each particle of a rigid body is given by:

$${\overrightarrow{v}}_i=\overrightarrow{V}+\overrightarrow{\Omega }\times \left({\overrightarrow{r}}_i-{\overrightarrow{R}}_0\right)$$

(8)

The rigid bodies are also defined with a set of discrete particles [38]. The new modified Dynamic Boundary Condition (mDBC) approach is applied to describe the boundary particles [39]. This mDBC approach overcomes some of the drawbacks of the default Dynamic Boundary Condition (DBC) approach in DualSPHysics [40], resulting in less dissipation in flows near boundaries and reducing unphysical pressure and density values for the boundary and adjacent fluid particles. The motion of a floating body is described by net force on each boundary particle as the sum of the contributions of all surrounding fluid particles. The motion of a floating body is given by integrating Equation (8) in time, using the Symplectic scheme [41]. A variable time step established by the Courant–Friedrich–Lewy condition is considered [42]. This is based on the force per unit mass, viscous diffusion term and Courant–Friedrich–Lewy (CFL) condition. Lastly, DualSPHysics updates the position, velocity, density and pressure of all particles of the computational domain at the new time step.

3. Numerical Model Validation

Since no detailed data are available for WEHC, the numerical model was first validated against experimental data of [43] for a Wavestar point-absorber-type device. This selected device has a working principle similar to WEHC, based on the interaction between the three different sub-systems: floater, lever arm and hydraulic PTO system. The floater is connected by a fixed arm to a hinge in the support platform where the rotational motion is damped by a hydraulic pump to absorb power in operational conditions. In [43], the experimental tests were carried out using a 1:5-scale Wavestar model. The schematic of the computational domain with the Wavestar model is shown in Figure 2. The channel is 24 m-long and has depth in the flat-bottom section of $d=3$ m. A passive dissipative beach located at the end of the channel with a slope of ${10}^{\circ }$ is used. The wave was generated in DualSPHysics using a piston-type wavemaker with active wave absorption system (AWAS) to absorb any reflected wave energy coming from the impact with the beach and Wavestar [44]. The Wavestar model was located 14 m from the wavemaker with an initial submergence of 0.4 m. The geometry and dimensions of the Wavestar model are given in Figure 2b.

In this study, four monochromatic wave conditions were considered.

Table 1. Wave conditions considered for the validation of Wavestar.

	Test 1	Test 2	Test 3	Test 4
H (m)	0.1	0.15	0.25	0.25
T (s)	1.4	1.4	1.4	2.8
$kA$ (-)	0.103	0.154	0.257	0.07

presents the wave height H, wave period T and wave steepness $kA$, with k being the local wavenumber and A the wave amplitude of each test.

The hydraulic PTO system was modelled as an external force applied to the lever arm. This force was imposed 1 m from the revolute joint between lever arm and support platform at ${52}^{\circ }$, as shown in Figure 2b. As defined in [43], the applied PTO force is modelled as a linear spring–damper system:

$$F_{PTO}\left(t\right)=D\dot{z}\left(t\right)+Kz\left(t\right)$$

(9)

where the first term, $D\dot{z}\left(t\right)$, is the linear damping force and second term, $Kz\left(t\right)$, is the restitution force (linear spring) applied on the lever arm by the hydraulic pump, D and K being the damping and stiffness coefficients of the linear spring–damper PTO system, respectively. Similar to [43], $K=0$, and four different damping coefficients were used $D=0$, 50, 100 and 200 Nms.

Simulation accuracy is estimated by quantifying the discrepancy between the numerical and experimental Wavestar motion, including heaving position z and velocity $\dot{z}$. Before performing these comparisons, a convergence analysis was carried out. The numerical resolution in DualSPHysics is defined using the initial interparticle distance $dp$. Four different $dp$ were considered: $dp=H/5$, $H/10$, $H/15$ and $H/20$. Similar to several applications of SPH models to simulate waves, a computational domain with $dp=H/10$ shows a good agreement between numerical results and experimental data [26]. This discretisation allows a good balance between computational cost and simulation runtime. The $dp$ chosen in this work is $dp=0.01$ m, which corresponds to the resolution of $dp=H/10$ for the smallest wave height tested (Test 1). With this resolution, about 624,500 particles are generated in DualSPHysics. The average time consumed was about 1 h to solve 50 s of physical time, running on an Nvidia RTX 3060.

The numerical and experimental time series of PTO position and velocity for Test 1 with different D are compared in Figure 3.

As can be observed, a good agreement was obtained between numerical results and the experimental data. However, a slight difference in both position and velocity can be found, especially for lower damping values (more amplitude of heaving motion). It is also observed that as the damping increases, the position and velocity decrease, as expected. The same behaviour was observed for other wave conditions and damping values; therefore, such results are not presented here. Figure 4 shows the time series of the PTO position and velocity for Test 2 with $D=0$ and 200 Nms. All numerical results are consistent with the experimental data.

The time series for Test 3 and 4 with $D=200$ Nms are compared in Figure 5. Similarly to Test 1 and 2, good agreement was found between the experimental data and the numerical results. Test 3 achieves slightly smaller differences with the experimental results than Test 4, even for the same wave height. This is due to a small lag in the wave period. However, the velocity appears to be slightly lagging behind when compared with the experimental velocity, with almost the same amplitude.

The time series of Wavestar motions obtained in the numerical model with different PTO values show a satisfactory agreement with the experimental data in the physical model tests. For a more quantitative comparison, the accordance between the numerical and experimental data was quantified using the average normalised root mean error (NRMSE) in amplitude.

Table 2. Average NRMSE for the different wave states with $D=0$.

	Test 1	Test 2	Test 3	Test 4
z (m)	10%	8%	7%	9%
$\dot{z}$ (m/s)	11%	6%	5%	13%

shows NRMSE for the four tests with $D=0$. Overall, the NRMSE is less than 10%, although in some conditions it is slightly higher, caused by the small difference in phase between the numerical results and experimental data. For all wave conditions, the NRMSE decreases with increasing D.

4. Results and Discussion

4.1. Simulation Set-Up

The power absorption of the WEHC device is investigated under full-scale sea wave conditions in the south-western coast of Portugal. The location of the western breakwater of the Port of Sines was considered for the present study. In order to correctly predict the power extracted from the waves, the first step is to correctly predict the incoming waves. For this reason, the waves chosen to simulate the device were based on the work of [45], who developed a detailed library of the most common waves and sea states along the Portuguese coast. Both regular and irregular sea states, and different operation modes and design parameters, were considered, including hydraulic PTO damping and breakwater geometries. Two different breakwater geometries were modelled: one as a straight wall with no porosity and another made by a double layer of cubic blocks with a porosity of around 40%, laying on an impermeable slope bed.

A schematic of the computational domain is shown in Figure 6. The channel is 240 m-long and has a depth in the flat-bottom region of $d=30$ m. The WEHC device is located at the right-hand end of the channel with an initial submergence of 0.67 m. The geometry and dimensions of the device are presented in Figure 6b.

The computational domain was discretised with a $dp=0.1$ m, i.e., $dp=H/10$ for the smallest wave height tested. The simulations were also performed on an Nvidia RTX 3060. The average computational time was about 5 h for a resolution of 564,374 particles and 300 s of physical time.

The hydraulic PTO system is modelled as a single-piston effect with the external force given as

$$F_{PTO}\left(t\right)=\left\{\begin{array}{ccc}{C}_{PTO}& \mathrm{if}& z\left(t\right)<0\\ 0& \mathrm{if}& z\left(t\right)\ge 0\end{array}\right.$$

(10)

where $C_{PTO}$ is the constant force applied 2 m below of the revolute joint between lever arm and support platform. The friction torque on the universal joint between lever arm and floater and the revolute joint between lever arm and support platform was modelled using the Coulomb friction coefficient of 0.16 [26].

4.2. Power Absorption under Regular Sea State

The hydrodynamic response of a WEHC device under a regular sea state is analysed considering the capture width ratio (CWR), defined as the ratio between the average extracted power $W_{out}$ and the average incoming wave power $W_{inc}$, which is

$$\mathrm{CWR}=\frac{W_{out}}{W_{inc}}$$

(11)

The $W_{out}$ is calculated as

$$W_{out}=\frac{1}{t_f-t_0}{\int }_{t_0}^{t_f}{F}_{PTO}\left(t\right)\dot{x}\left(t\right)dt$$

(12)

where $\dot{x}$ is the velocity of the point where the external force is applied (Figure 6b), and $t_0$ and $t_f$ are the initial and final time, respectively.

According to linear wave theory, $W_{inc}$ is calculated as

$$W_{inc}=\frac{1}{8}\rho g{H}^2{c}_g$$

(13)

where $c_g$ is the wave group velocity, and for this study, it is calculated as

$$c_g=\frac{\omega }{2k}\left(1+\frac{2kd}{sinh\left(2kd\right)}\right)$$

(14)

where $\omega =2\pi f$ is the angular wave frequency, d is the water depth, f is the wave frequency and k is calculated by ${\omega }^2=gktanh\left(kd\right)$.

Figure 7 compares the variation of CWR with $C_{PTO}$ for six different sea states presented in

Table 3. Regular wave conditions considered in the numerical tests.

	State 1	State 2	State 3	State 4	State 5	State 6
H (m)	2	2	2	2.5	2.5	2.5
T (s)	7	9	11	7	9	11
$kA$ (-)	0.083	0.054	0.040	0.104	0.067	0.050

. These sea states’ characteristics are frequent on the south coast of Portugal [45].

A large variation of CWR with PTO damping can be observed in Figure 7. These constraints, imposed by the hydraulic PTO system, have a non-negligible influence on the power absorption of WEHC. As the $C_{PTO}$ increases, the CWR reaches its maximum. The curve shows almost similar trends for the different sea states tested, however, the maximum position is envisaged. For example, for $H=2$ m and $T=7$ s, the maximum CWR occurs with $C_{PTO}=175$ kN, and for $H=2$ m and $T=9$ s, it occurs with $C_{PTO}=225$ kN. A large variation of CWR is observed when the $C_{PTO}$ increases due to the rise of the latching effect of the hydraulic PTO system, where the floater remains stationary for as long as the mass forces are unable to overcome the resisting force introduced by the PTO system. In general, the CWR decreases when the wave period increases. For $H=2$ m and $T=7$, the maximum CWR ≈ 32%, and for $H=2$ m and $T=11$, the maximum CWR ≈ 19%. A similar trend is also observed when increasing the wave height, for example, for $H=2.5$ m and $T=7$, the maximum CWR ≈ 23%. With increasing wave height, the maximum CWR occurs for a large $C_{PTO}$. With increasing wave period, the CWR curve becomes flatter, and the maximum occurs at a higher value of $C_{PTO}$ for a longer period. This behaviour is caused by the non-linear waves interacting with the sea bottom that produces larger excitation forces, enlarging the bandwidth of the maximum CWR. A direct consequence of this behaviour is that it allows a simple control strategy for the pressure of the hydraulic PTO system based on wave period measurements and predictions from swell forecast models.

Figure 8 provides instantaneous snapshots of the velocity field in the vicinity of the WEHC device for $H=2$ m and $T=7$ s with straight breakwater. As expected, the flow field is characterised by a strong dynamical velocity field, with large acceleration near the floater. This causes differences in the turbulent field, characterised by the increase in the turbulent kinetic energy near the free surface and close to the floater.

As seen in Figure 8, the breakwater is a straight wall with no porosity, however, the breakwaters tend to be built of shaped rocks or cement blocks. This method of construction allows for more dissipation of wave energy on the breakwater, ensuring a better protection of the onshore structures. Given the proven capability of the DualSPHysics–Project Chrono coupling for simulating fluid–structure interaction [46], one additional breakwater geometry was tested in this study in order to obtain results on the influence of more realistic breakwaters. The variation of CWR with $C_{PTO}$ for six different sea states with $H=2$ and 2.5 m, $T=7$, 9 and 11 s presented in Table 3 are shown in Figure 9. Similar trends of CWR are found for different breakwater geometries, however, the maximum position is slightly envisaged. For example, for $H=2$ m and $T=7$ s with straight breakwater, the maximum CWR occurs with $C_{PTO}=175$ kN, and for porous breakwater, it occurs with $C_{PTO}=150$ kN. In general, the CWR reduces by about 5% from straight breakwater. The reflected waves are reduced due to the porosity of the breakwater. However, its effect on the excitation forces is much smaller than the incident counterpart, so it has less effect on CWR. The breakwater porosity has more effects when H increases; indeed, complex interactions arise when H increases, mainly for higher T, and the excitation forces are much affected.

In order to illustrate the influence of the breakwater porosity and geometry on the flow field, the detailed instantaneous velocity field in the vicinity of the WEHC device for $H=2$ m and $T=7$ s with porous breakwater is given in Figure 10.

The geometry shows the dissipation of the energy from incoming waves. Moreover, energy dissipation also causes a reduction in wave reflection, which in turn reduces the WEHC position motion, leading to a lower CWR. Even though this effect has a negative impact on energy extraction, this is closer to how breakwaters are built.

4.3. Power Absorption under Irregular Sea State

The ability of the model to reproduce more realistic wave conditions was tested by simulating the behaviour of the WEHC device for irregular wave conditions. The following JONSWAP (Joint North Sea Wave Observation Project) spectral shape was used [47]:

$$S_{\eta }\left(f\right)=\alpha g^2{\left(2\pi \right)}^{-4}{f}^{-5}exp\left[-\frac{5}{4}{\left(\frac{f_p}{f}\right)}^4\right]{\gamma }^{exp\left[{\displaystyle -\frac{{(f-f_p)}^2}{2{\sigma }^2{f}_p^2}}\right]}$$

(15)

where $\alpha =0.291$, $f_p$ is the peak frequency, $\gamma =3.3$ is peak enhancement factor, $\sigma =0.07$ for $f\le f_p$ and $\sigma =0.09$ for $f>f_p$.

On the basis of the position time series, $x\left(t\right)$, and free surface elevation, $\eta$, the power density spectra of the WEHC motions $S_x$ and free-surface $S_{\eta }$ were calculated for each simulation. These power spectra were computed using Welch’s technique, in which each power spectral density is estimated by dividing the data into overlapping segments, computing a modified periodogram for each segment and averaging the periodograms. The Hann window with segment length of 1024 points and overlap of 50% was used.

As can be observed in Figure 11, the spectra also properly reflect the spectral peaks, i.e., where most of the energy is concentrated and, thus, the frequency range of the greatest interest. Three sea states were considered with the significant wave height $H_s=2.5$ m and peak period $T_p=$ 7, 9 and 11 s.

The time series of each sea state for 100 s $<t<$ 300 s is shown in Figure 12. Similarly to the regular sea state, an increase in CWR with the peak period is observed. For $H_s=2.5$ m and $T_p=$ 7 s, a CWR of about 25% was found. The CWR for irregular waves is evaluated by Equation (11), where the average incident wave power was calculated as

$$W_{inc}=\rho g{\int }_0^{\infty }{c}_g\left(f\right)S_{\eta }\left(f\right)df$$

(16)

Figure 13 shows the power spectrum for three sea states considering significant wave height $H_s=2$ m and peak period $T_p=$ 7, 9 and 11 s. The spectra also properly reflect the spectral peaks.

Figure 14 shows the time series of each sea state for 100 s $<t<$ 300 s. As for a regular sea state, an increase in CWR with the peak period is observed. For $H_s=2$ m and $T_p=$ 7 s, a CWR of about 30% was found.

5. Conclusions

This paper presents a numerical investigation of the power absorption of a constrained wave energy hyperbaric converter (WEHC) under full-scale sea wave conditions. A fully non-linear numerical model DualSPHysics was used to model the interaction between waves and each WEHC sub-system. The numerical model was first validated against the available experimental data of a similar device, a 1:5-scale Wavestar model. A good agreement was found between the experimental data and the numerical results, with an average normalised root mean error (NRMSE) of about 10%.

The numerical model was then employed to perform a study of the non-linear hydrodynamics of a constrained WEHC under regular and irregular sea states, and for different operation modes and design parameters, including hydraulic power take-off (PTO) damping and breakwater geometry. For regular waves, it was observed that the capture width ratio (CWR) is particularly sensitive to the variations in the PTO damping, although the absolute maximum CWR is less sensitive to the slightly lower variations applied to the PTO damping. It was also found that the influence of the wave height and wave period have an important effect on the CWR. For irregular waves, a similar trend on the CWR was observed, with CWR increasing with decreasing peak period and an increase in significant wave height. In general, the efficiency of the system is very dependent on the mechanical constraints (PTO system) and on the breakwater characteristics. These results show the potential and effectiveness of this type of numerical modelling approach for understanding the non-linear dynamics of WEHC in a full-scale application.