Experimental and CFD Assessment of Harmonic Characteristics of Point-Absorber Wave-Energy Converters with Nonlinear Power Take-Off System

Abstract

The wave-energy excitation of point absorbers is highly associated with their resonant movement, and harmonic characteristics are of increasing concern in affecting resonance. However, the commonly used linearized power take-off (PTO) systems underestimate the impact of harmonics. The purpose of this study is to address the knowledge gap in assessing the contribution of hydraulic PTO systems to higher harmonic wave loads and velocities. In the present work, higher harmonics in point-absorber wave-energy converters (PA-WECs) with hydraulic power take-off (PTO) systems are investigated through both experimental and computational fluid dynamics (CFD) methods. The fast Fourier transform is used to decompose the high-order harmonics. To account for the influence of nonlinear wave–wave interaction on harmonics, the isolated PA-WEC is used as a basis for comparison with the paired PA-WECs. The influence of wave steepness is also estimated at two resonance periods. Results indicate that the additional resonance of the paired PA-WECs may be attributed to the harmonic wave loads at longer wave periods. Harmonic wave loads of paired PA-WECs typically have a more substantial impact and increase more rapidly with increasing wave steepness compared to isolated PA-WECs. Furthermore, as the wave steepness increases, there are significant enhancements in both the harmonic wave loads and heaving velocity, which strongly correlate with the instantaneous maximum hydraulic power. Consequently, our study will contribute to enhancing the maximum power output in the design of future point absorber arrays.

1. Introduction

The search for alternative energy sources to address the growing concerns over fossil fuel use and global warming has led to the exploration of ocean renewable energy. Wave energy, with a potential of 2–3 kW/m² [1], is a high-density energy source. Consequently, significant efforts have been made to commercialize its conversion. Wave-energy converters (WECs) are devices designed for this purpose, and various types have been proposed. Notably, point-absorber devices have demonstrated exceptional adaptability and effectiveness [2], rendering them the favoured selection for practical engineering applications, exemplified by W2Power [3] and Wave Star [4].

Regarding the harmonic hydrodynamics of point-absorber wave-energy converters (PA-WECs), previous results have revealed their significant role in wave elevation and motion responses. Tran et al. [5] discovered that geometric nonlinearities lead to sub-harmonic excitations that severely compromise device performance. Yang et al. [6] observed a significant contribution from the harmonic component of wave elevation and noted that this component increases performance under stationary waves around WECs. Zhang et al. [7] discovered that the harmonic component significantly influences the instantaneous power absorption of a point absorber, with higher-order harmonics being more pronounced when motion responses are predominantly influenced by pitch resonance. Haeng Sik Ko et al. [8] studied the dynamic behaviour of an asymmetric wave-energy converter and found that the third harmonic moment term is significant in both experimental and OpenFOAM results, especially at higher amplitudes. A multi-mode PA-WEC investigation by Tran et al. [9] revealed that higher harmonics in motion response can cause more frequent changes in the projected surface area, which decreases the obtained power. Therefore, harmonic components have been shown to play a crucial role in the hydrodynamic performance of wave-energy converters.

While the significant impact of harmonics on the performance of PA-WECs has been acknowledged in prior research, the current understanding lags behind that of oscillating water-column devices. In the case of the latter, harmonics have been linked to an innovative resonance mechanism that enhances power output. Zhou et al. [10] found that an accurate description of air pressure requires all five harmonic components. The wave–wave nonlinear evolution that generates higher harmonics can resonate with the chamber of an oscillating water column device, significantly affecting its efficiency [11]. The limited investigation of harmonics in PA-WECs is due to the linear simplification of the PTO system [12], which is closely associated with the hydrodynamic performance of point absorbers.

Nonlinearities stemming from the PTO system exert a significant impact on the converter’s overall performance. Windt et al. [13] performed simulations on a point-absorber device employing a linear PTO system. The results illustrate pronounced nonlinear behaviour that influences its resonance when the PTO system is optimized for energy extraction. Yu et al. [14] investigated the nonlinear wave run-up phenomenon of point-absorber devices by incrementally increasing the linear PTO load and decreasing the wave period. The results indicate that the PTO system has a notable impact on the wave run-up in front of the buoy. Cheng et al. [15] examined the nonlinear performance of the WaveRoller device and observed that the influence of the PTO system’s nonlinearity on the hydrodynamics of the WEC becomes more pronounced as wave amplitude increases. Earlier research reinforces the association between the nonlinear performance of PA-WECs and the PTO system. Nonetheless, the utilization of a linear or weakly nonlinear PTO model may result in an underestimation of their impact on hydrodynamics [16]. Consequently, the accurate analysis of harmonic wave loads necessitates the adoption of a more precise nonlinear model [17].

Recently, the arrangement of PA-WECs has become a topic of interest, as the interaction between buoys is significant for performance analysis. Agamloh et al. [18] used CFD technology to simulate a PA-WEC array. Results indicate that buoy #1 experiences a smaller displacement than an isolated floater, while buoy #2 has an increased amplitude. Sun et al. [19] conducted hydrodynamic tests on a floating PA-WEC and found that intense buoy interaction can significantly decrease the optimal wave period for the wave-energy converter. According to Li et al.’s [20] framework, which combines CFD with a multi-body dynamic structural solver, there is a persistent transient disruption in the natural frequency of the floater’s pitch due to complex interactions. The intensive wave–structure interaction and wave–wave nonlinear evolution have proven to be significant sources of harmonic excitations that affect the performance of PA-WECs. Therefore, a more thorough understanding of the harmonic hydrodynamic characteristics of PA-WECs, especially when considering buoy interaction, is crucial for performance analysis.

The main objective of the present work is to examine the impact of harmonic components on the hydrodynamics of PA-WECs. Previous linearization of the PTO system has resulted in the underestimation of harmonic hydrodynamic characteristics, hindering further understanding of the resonance mechanism. Therefore, in this work, we include a hydraulic PTO system to consider the nonlinearity arising from its effect on harmonic characteristics. Second, we use paired PA-WECs to preliminarily consider the effect of buoy interaction on harmonic wave loads. Both experimental and CFD-based numerical investigations of the harmonic characteristics of paired PA-WECs were carried out to provide a more convincible estimation. The contribution of this work is twofold: it provides a hydraulic PTO system model that can consider nonlinearity in hydrodynamics for both experimental and numerical models, and it offers a thorough investigation into the harmonic characteristics of PA-WECs with buoy interaction taken into consideration.

The paper is organized into five sections. Section 2 introduces the concept of a pair of PA-WECs, and discusses the physical model and the CFD-based model used in the study. In Section 3, a detailed comparison between the numerical model and physical test results is presented to validate the accuracy of the numerical model. Section 4 analyses the harmonic components of the model and investigates the effects of wave steepness and buoy interaction. Finally, the paper concludes by summarizing the findings in Section 5.

2. Mathematical Theories and Paired PA-WECs Model

2.1. Problem Description of Paired PA-WECs

The isolated PA-WEC, as shown in the left portion of Figure 1a, consists of an absorber and a hydraulic PTO system. The origin of the global coordinate system was established at the gravity centre of the point-absorber device. For paired PA-WECs, the origin was established midway between the two point absorbers, with the x-axis directed towards the wave elimination bank and the y-axis determined by the right-hand rule.

The hydraulic PTO system comprises a working cylinder, an accumulator, and a PVC hose. The working cylinder is located near the point absorber, while the accumulator is situated on the trailer. The two cylinders are connected by a PVC hose. The operation can be divided into two stages. During the upward heave of the point absorber, as indicated by the blue line in Figure 1, the liquid in the working cylinder is compressed and flows into the accumulator. This conversion transforms wave energy into the potential energy of the load. During the downward heave of the point absorber, seawater is drawn into the working cylinder while preventing the conversion of wave energy into potential energy as the load remains stationary, preventing the conversion of wave energy into potential energy. Check valves are used to prevent water backflow. At the end of a test condition, the gate valve is opened to clear out the water in the accumulator.

The paired PA-WECs are a combination of two isolated PA-WECs, as illustrated in Figure 1b. Both isolated and paired PA-WECs are constrained to heave motion only. The usage of the paired PA-WEC is to consider the effect of hydrodynamic interaction between point absorbers on harmonic characteristics.

Each absorber in the paired PA-WECs is subjected to the hydrodynamic force (F_hydro), the PTO system force (F_PTO) and its gravity (G). The motion equation of the floater can be expressed as Equation (1):

$$M_i{\ddot{z}}_i=F_{hydro}^i+F_{PTO}^i+G_i$$

(1)

where M_i is the mass of the ith floater in paired PA-WECs, and $\ddot{z}$ is the heaving acceleration.

The heave-response-amplitude operator (RAO_hea) of each PA-WEC is illustrated in Equation (2):

$$RA{O}_{hea}=\frac{z_{amp}}{H}$$

(2)

where H represents the incident wave height and z_amp represents the amplitude of heave motion.

The power generated by the hydraulic PTO system is referred to as the hydraulic power of the PA-WEC. The instantaneous hydraulic power, defined as the product of PTO system force and heaving velocity, has an average value expressed in Equation (3):

$$P_1=\frac{1}{\tau }{\displaystyle \underset{0}{\overset{\tau }{\int }}-F_{PTO}(t)\dot{z}(t)dt}$$

(3)

The overall power of the PA-WEC is defined as Equation (4):

$$P_2=\frac{1}{\tau }{\displaystyle \underset{0}{\overset{\tau }{\int }}mgV(t)dt}$$

(4)

where m represents the PTO system load, τ is the selected time interval, and V is the load lifting velocity.

The incident wave power is defined as Equation (5):

$$P=\frac{\rho g}{2}{(\frac{H}{2})}^2\frac{\omega }{2k}(1+\frac{2kh}{\mathrm{sh}2kh})L_c$$

(5)

where sh is the hyperbolic sine function, L_c is the characteristic length, ω is the incident wave circular frequency, k is the wave number, and h is the tank’s water depth.

The efficiency of the paired PA-WECs can be defined as Equation (6):

$${\eta }_i=\frac{{\displaystyle \sum _{j=1}^2{P}_i^j}}{P},i=1,2$$

(6)

where i = 1 refers to the hydraulic efficiency, η₁, and i = 2 refers to the overall efficiency, η₂; j is the number of floaters.

The first three harmonic components of the variable of interest of the current paired PA-WEC model, under monochromatic wave excitation, can be described as Equation (7):

$$\begin{array}{l}x(t)=2[B_1(j\omega ,A)\cos (\omega t+{\theta }_1)+B_2(j\omega ,A)\cos (2\omega t+{\theta }_2)\\ \begin{array}{cc}& \end{array}+B_3(j\omega ,A)\cos (3\omega t+{\theta }_3)+\cdots ]\end{array}$$

(7)

where B_i is the amplitudes of harmonic components defined by Rugh’s general formula, and θ_i is the phase.

2.2. Experimental Setup

2.2.1. Layout of the Experiment

The model test for both isolated and paired PA-WECs was conducted in the towing tank of the National Ocean Technology Center, as depicted in Figure 2. The origin of the experimental coordinate is the same as the numerical setup. The tank measured 18 m in width, and 130 m in length, and had a water depth of 4.5 m. The PA-WEC featured a cone-bottom floater shape optimized by a previous investigation [21] using a viscous-corrected potential flow theory model. Figure 3 shows the physical model of the point absorber, while

Table 1. Main dimensions of the wave-energy converter.

Item	Value	Unit
Diameter	0.4	m
Draft	0.25	m
Absorber spacing	0.5	m
Displacement	14.65	kg

lists the main dimensions of the isolated and paired PA-WECs. The dimensions of the absorber are sourced from a floating wave-energy converter, as documented in [19]. The distance between the paired absorbers is defined by the distance from one’s centre to the other’s.

Constrained by a prismatic joint, the PA-WEC can only move in the vertical translation direction. A guyed displacement sensor, YM02-P, was installed between the truss and the point absorber. The heaving velocity was obtained by differentiating the displacement measurements from YM02-P. An ultra-small strain gauge, LCM201, was installed between the PTO system and the point absorber to measure the PTO force exerted on it. Two wave probes were installed before and after the wave-energy converter. The wave probe has a sampling frequency of 20 Hz, while the displacement sensor and strain gauge have a sampling frequency of 50 Hz.

2.2.2. Physical Test Conditions

Table 2. Overview of the physical test schedule.

Number	Item	Model	PTO Load	Wave Period
1	Resonance period	Isolated PA-WEC	0	1.2–2.2 s
2	Optimal PTO load	Isolated/Paired PA-WEC	3.48–14.76 kg	1.2–2.2 s
3	Power performance		12.12 kg	1.2–2.2 s

shows the schedule for the physical test. Before conducting a comparison, we estimated the fundamental characteristics of the PA-WEC device, such as its heave resonance period and optimal PTO load, through tests No. 1 and No. 2. During test No. 1, the PTO system was connected to the point absorber without external loading, enabling us to determine the resonance period. Figure 4 presents the results: the heaving resonance period of the isolated PA-WEC without PTO load is 1.8 s, at which point the absorber achieves maximum heave amplitude.

Previous research has demonstrated that PA-WEC performance is significantly influenced by PTO system damping [22]. Damping of the present system is closely related to the PTO system load. Consequently, determining the optimal load, achieved through test No. 2, is essential.

During test No. 2, we varied the load of the PTO system while keeping the incident wave period constant. We evaluated wave periods ranging from 1.2 s to 2.2 s to determine the optimal load for each period, using maximum overall efficiency (η₂) as the selection criterion. As shown in Figure 5, the optimal PTO system load for both the isolated and paired PA-WECs is 12.12 kg, while the suboptimal load for the pair of PA-WECs is 9.24 kg. Since our aim is not to conduct an extensive analysis of the optimal load in this study, we present only the results for an incident wave period of 1.8 s in the figure. More detailed results can be found in [19].

After establishing the fundamental characteristics of the PA-WEC, we conducted test No. 3 to examine its power performance. The results of this test will be validated using the CFD model in Section 3.

2.3. Numerical Hydraulic PTO System Simulation

As aforementioned, the current PTO system is single-acting and cannot effectively produce work when the point-absorber device moves downward. Consequently, the force generated by the PTO system is defined as a piecewise function, with $\dot{z}=0$ serving as the threshold for segmentation, as shown in Equation (8):

$$F_{PTO}=\left\{\begin{array}{cc}-p_1{S}_1-F_c{\dot{z}}^3-f_v\dot{z}& \dot{z}\ge 0\\ k_{dam}\dot{z}& \dot{z}<0\end{array}\right.$$

(8)

where S₁ represents the piston working area, F_c represents the Coulomb friction force [23], f_v represents the viscous friction coefficient, k_dam is the damping factor and $\dot{z}$ denotes the heaving velocity.

Assuming that the flow in the tube is stable, we can use the Bernoulli equation to describe the pressure. The pressure exerted on the piston in the working cylinder, p₁, is associated with the PTO system load as Equation (9):

$$p_1=k_{cor}(\frac{mg}{S_2}+\rho g(h_2-h_1))$$

(9)

where k_cor is the correction factor to consider deviations caused by dynamic pressure, m is the PTO system load, and h₁ and h₂ are the height of the piston and tailer as in Figure 1.

The coefficients k_dam and k_cor can be obtained from experimental data, and the fitted values shown in Figure 6 are used in the numerical model setup. This reveals that the damping term in Equation (8) is associated with wave periods and a minimum damping factor can be obtained at 1.8 s. Additionally, the correction factor increases as the PTO system load decreases and also achieves its minimum at 1.8 s. Therefore, the PTO system shows a nonlinear relationship with either heaving velocity or system load and wave periods.

2.4. CFD Approach

2.4.1. Flow Field Model

To explore the complex hydrodynamic performance of the PA-WEC’s wave–structure interaction, a CFD model is created based on the STARCCM+ [24]. The finite volume method (FVM) [25] is commonly used in constructing CFD models. This involves discretizing the computational domain into control volumes, with fluid properties stored at the volume centroid. Temporal domain discretization is then performed to investigate time-domain quantities. The governing equations for the model are the continuity equation and Navier–Stokes equations, as shown in Equation (10):

$$\begin{array}{l}\nabla \cdot u=0\\ (\frac{\partial }{\partial t}+u\cdot \nabla )u=-\nabla (\frac{p}{\rho }+gz)+\nu {\nabla }^{2}u\end{array}$$

(10)

where $u$ represents velocity, $p$ represents pressure, and $\nu$ represents the kinematic viscosity coefficient.

As wave–structure interaction is a flow problem involving multiple phases, the volume of fluid (VOF) method [26] is used to model the free surface. The volume of fluid (VOF) method has a distinct advantage over the marker-and-cell (MAC) method in simulating nonlinear phenomena like wave breaking and run up. This is because the VOF model tracks the water volume fraction, which is shown in Equation (11):

$$\frac{\partial \alpha }{\partial t}+\nabla \cdot (u\alpha )+\nabla \cdot ((u-u_a)\alpha (1-\alpha ))=0$$

(11)

where $\alpha$ represents the water volume fraction, and $u_a$ represents the air phase velocity.

There are discrepancies in the choice of turbulence models, with k-ε, k-ω, k-ω SST, and Reynolds stress models being options. However, based on successful applications to solve wave-breaking [27] and run-up [28] problems, the k-ω SST model is recommended for this study. The k-ω SST model [29] uses the blending function to compute the turbulence flows by coupling the k-ω model in the near walls and the k-ε model in the free shear flows. It can be written as Equations (12)–(14).

$${\mu }_t=\rho a_{1}k\frac{1}{\max (a_1\omega ,b_{1}S{F}_2)}$$

(12)

$$\frac{\partial k}{\partial t}+U_j\frac{\partial k}{\partial x_j}=\min (P_k,10\beta *k\omega )-\beta *k\omega +\frac{\partial }{\partial x_j}[(\mu +{\mu }_t{\sigma }_k)\frac{\partial k}{\partial x_j}]$$

(13)

$$\begin{array}{l}\frac{\partial \omega }{\partial t}+U_j\frac{\partial \omega }{\partial x_j}=\rho \alpha \frac{\omega }{k}\min [G,\frac{c_1\beta *\omega }{a_1}\max (a_1\omega ,b_1{F}_{23}S)]-\beta {\omega }^2\\ +\frac{\partial }{\partial x_j}[(\mu +{\mu }_t{\sigma }_{\omega })\frac{\partial \omega }{\partial x_j}]+2(1-F_1)\frac{{\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}\end{array}$$

(14)

where ${\mu }_t$ is the eddy viscosity, S denotes the magnitude of strain rate tensor and P_k is the production kinetic energy, in which the Kato–Launder correction is included (P_k = G = 2μSΩ). F₁ and F₂ are the blending functions and the value of model constants are a₁ = 0.31, b₁ = 1.0, β∗ = 0.09, α₁ = 5/9, α₂ = 0.44, σ_k1 = 0.85, σ_k2 = 1.0, σ_ω1 = 0.5, σ_ω2 = 0.856, β₁ = 0.075, β₂ = 0.0828, σ_k = F₁σ_k1 + (1 − F₁)σ_k2, σ_ω = F₁σ_ω1 + (1 − F₁)σ_ω2, β = F₁β₁ + (1 − F₁)β₂, α = F₁α₁ + (1 − F₁)α₁.

2.4.2. Computational Domain and Boundary Conditions

The computational domain depicted in Figure 7 has an identical width and water depth as that of the physical tank. To enhance computational efficiency, the numerical tank’s length is equal to five times the wavelength of the incident wave [28]. The figure clearly illustrates the boundary conditions: the inlet is a velocity inlet, the outlet is a pressure outlet, and the other boundaries are non-slip walls. Additionally, one forcing zone and one wave-elimination zone are located at one and two wavelengths from the inlet and outlet, respectively.

2.4.3. Numerical Wave Generation

The wave-making system in the CFD model serves as a momentum source. Waves are generated by the wave elevation and velocity at the velocity inlet. To account for nonlinear wave generation, the fifth-order Stokes theory, as described in Fenton’s work [30], is adopted. The resulting free surface profile can be expressed as Equation (15):

$$\begin{array}{l}\eta =\frac{1}{k}{\displaystyle \sum _{m=1}^q{\displaystyle \sum _{n=1}^m{b}_{mn}}}{(\frac{kH}{2})}^m\cos (n(kx-\omega t))\\ u=\sqrt{\frac{g}{k}\tanh (kd)}{\displaystyle \sum _{m=1}^5{(\frac{kH}{2})}^m{\displaystyle \sum _{n=1}^{m}n{a}_{mn}}}\cosh (nk(z+d)\sin (n(kx-\omega t))\\ v=0\\ w=\sqrt{\frac{g}{k}\tanh (kd)}{\displaystyle \sum _{m=1}^5{(\frac{kH}{2})}^m{\displaystyle \sum _{n=1}^{m}n{a}_{mn}}}\cosh (nk(z+d)\sin (n(kx-\omega t))\end{array}$$

(15)

where $\eta$ is the wave elevation; $u,v,w$ are velocity components; and $\omega ,k,d,H$ are circular frequency, wave number, water depth and wave height, respectively.

A wave-elimination zone is placed near the outlet boundary to eliminate reflected waves as in Figure 7. This is achieved by adding a damping term to the vertical velocity, as shown in Equation (16):

$$S=\rho (f_1+f_2\left|w\right|)\frac{e^{\kappa }-1}{e^l-1}w,\kappa ={(\frac{x-x_s}{x_e-x_s})}^q$$

(16)

where $f_1,f_2,q$ are damping model parameters, and $x_s,x_e$ are the start and end points of the elimination region, respectively.

To avoid wave reflection or viscous effect dominated wave-elevation reduction, a forcing zone is implemented at the inlet. This is achieved by adding a source term [31] to the momentum equation, as shown in Equation (17):

$$\begin{array}{l}{q}_{\varphi }=-\gamma \rho (\varphi -\varphi *)\\ \gamma =-{\gamma }_0{\cos }^2(\frac{\pi x}{2})\end{array}$$

(17)

where $q_{\varphi }$ is the additional source term, $\gamma$ is the force coefficient, $\varphi$ is the present potential, and ${\varphi }^{*}$ is the terminal potential.

2.4.4. 6 DoF-Solver

The motion responses of each point absorber are obtained using a dynamic fluid–body interaction (DFBI) solver that is included in the software. The buoy motion is represented by the near-body grid velocity, as shown in Equation (18):

$$v_{buoy}=w_{grid}\times r+v_{grid,t}$$

(18)

where $v_{grid}$ is the grid velocity, $w_{grid}$ is the rotational velocity, r gives the node position, and $v_{grid,t}$ is the translation velocity.

The calculation was performed by a supercomputer, whose CPU is AMD EPYC 7H12 @ 2.6 GHz, single core has 128 threads, memory is 512 G, and operating system is CentOS 7.

3. Convergence Study

3.1. Wave Generation Test

In this section, we ensure the numerical model’s independence to guarantee the convergence of results before comparing them with the experimental outcomes.

The mesh model, as shown in Figure 8, is primarily divided into three regions: the background region, the near-body region, and the off-body region. To achieve a more precise simulation of the point-absorber device’s motion, the mesh in the off-body region is refined. This region has dimensions of eight times the floater’s diameter by twice its diameter by nine times the wave height. Additionally, refinement of the mesh around the free surface is conducted to ensure accurate simulation of the waves.

First, the refining of the free surface is discussed and

Table 3. Mesh parameters for free surface refining.

Mesh Type	Nodes in a Wavelength	Aspect Ratio	Grids Number	Time Step
Corse	25	1/2	18,528	0.025733
Medium	50	1/4	115,358	0.012867
Fine	100	1/8	1,163,566	0.006672
Extra Fine	135	1/16	8,477,988	0.003336

outlines four types of mesh for refinement. As indicated in Figure 9, the primary distinction in wave elevation is concentrated on the crests. The extra fine mesh type, which has a total of 8 million grids, requires more time for numerical simulation. However, since the fine mesh type yields comparable results to the extra fine mesh type and aligns with the results obtained from the Stokes fifth-order regular wave theory (as shown in Figure 9), it is chosen for constructing the present CFD model.

3.2. Mesh Convergence Tests

The parameters of the near-body region are also addressed. When the isolated point absorber is excited under regular waves without additional damping, mesh independence can be observed through the gradual stabilization of the heaving amplitude.

Table 4. Mesh independence test parameters.

Mesh Type	Time Step/s	Iterations in a Step	Boundary Layer		Grid Length		Grid Numbers
			Layers	Height/m	Off-Body/m	Near-Body/m
Mesh1	0.005	10	5	0.004	0.04	0.02	1,537,152
Mesh2	0.005	10	10	0.004	0.02	0.01	2,126,418
Mesh3	0.005	10	20	0.004	0.01	0.005	6,134,235
Mesh4	0.005	10	40	0.004	0.005	0.0025	31,982,971

provides detailed parameters for various mesh models. Iterations within a step represent the inner iterations of a time step, and residuals can be decreased by increasing the number of iterations.

As shown in Figure 10, the floater’s heave motion initially experiences rapid growth before eventually converging to a specific amplitude. The relative error is less than 0.5% between the last two types of mesh and the time step. Consequently, there are only minimal discrepancies between Mesh 3 and 4, but the latter will require more time for numerical simulation. Hence, Mesh 3 is selected for the current simulation.

3.3. Time Step Convergence Tests

Similarly, the convergence test is conducted to determine the proper time step.

Table 5. Mesh independence and time step independence test parameters.

Mesh Type	Time Step/s	Iterations in a Step	Boundary Layer		Grid Length		Grid Numbers
			Layers	Height/m	Off-Body/m	Near-Body/m
Mesh3	0.001	20	20	0.004	0.01	0.005	6,134,235
Mesh3	0.003	15	20	0.004	0.01	0.005	6,134,235
Mesh3	0.005	10	20	0.004	0.01	0.005	6,134,235
Mesh3	0.01	5	20	0.004	0.01	0.005	6,134,235

provides detailed parameters for time steps of 0.01 s, 0.005 s, 0.003 s, and 0.001 s, respectively.

As shown in Figure 11, the floater’s heave motion initially experiences rapid growth before eventually converging to a specific amplitude. The relative error is less than 0.5% between the last two types of time steps. Step 4 will require more time for numerical simulation. Hence, Step 3 is selected for the current simulation.

3.4. CFD Model Validation

The performance estimated by the CFD model of the isolated and paired PA-WECs is presented in Figure 12, and the results are compared with the physical test. The heaving amplitude (RAO_hea) and hydraulic efficiency (η₁) of the isolated point absorber reached the maximum level once, whereas the paired device reached it twice. In terms of the first resonance period, the RAO_hea and η₁ of the isolated point absorber were higher than those of the paired WEC. However, as the wave period exceeded the isolated point absorber’s resonant period, its efficiency and heave response were significantly reduced to a low level. On the other hand, the paired PA-WECs experienced another peak at a period of 2.2 s.

Phase analysis was conducted to further investigate the reason behind the second peak in the performance of the paired PA-WECs, as shown in Figure 13. In this figure, ξ refers to the elevation of both the wave and point absorber. When considering the isolated PA-WEC, it behaved as a “wave follower”, heaving upward and downward in phase with the incident wave. This behaviour is typical when away from the resonance period. However, an apparent phase lag was observed in the elevation of the point absorber in the array, revealing that the second peak was caused by another resonance primarily due to the harmonic components.

Good agreement is achieved between the physical tests and the CFD model, and detailed validation of the heaving velocity in the time domain is presented in Figure 14. The maximum experimental velocity is generally lower than the maximum velocity predicted by the numerical model, mainly due to the frictional works not being fully included in the CFD model. Additionally, some irregular variations were observed in the experimental results, which may be associated with the prismatic joint or pulling force from the displacement sensor.

The comparison of the PTO system force of CFD and experimental results in the time domain is shown in Figure 15. If F_PTO is greater than zero, it corresponds to the linear damper actor, while values below zero demonstrate the hydraulic characteristics of the PTO system.

4. Results and Discussion

4.1. Higher Harmonic Components in Hydrodynamics

A more comprehensive understanding of the velocity of the PA-WEC device is presented as follows. The fast Fourier transform (FFT) technique is used to analyse the harmonic components. Figure 16 depicts the heaving velocity spectrum where only the velocity of floater #1 in the pair is shown because the velocity of floater #2 is the same as floater #1. It reveals that both isolated and paired PA-WECs have three main peaks corresponding to the wave frequency, double wave frequency, and triple wave frequency, respectively. The low-level residuals, defined as $\mu =\left|{\dot{z}}_{sim}-{\dot{z}}_{exp}\right|$ [10], indicate that the CFD simulation can capture the low-level residual values associated with the harmonic components.

Figure 17 displays the decomposition of the heaving velocity at 1.8 s, while Figure 18 illustrates the decomposition at 2.2 s. The first harmonic component plays a major role in the heaving velocity due to wave excitation. Essentially, the second and third components mainly contribute to the maximum velocity. As high heave velocity is usually accompanied by a large PTO force, the harmonic components are crucial for optimal power performance. For both isolated and paired PA-WECs, the second harmonic heaving velocity takes a greater part as the period increases, while the proportion of the third harmonic velocity always exceeds that of the second. Nevertheless, a notable disparity between the isolated and paired PA-WECs becomes evident at T = 2.2 s: the second harmonic velocity of the isolated PA-WEC resembles that of the third component, whereas the third harmonic velocity of the paired PA-WECs continues to surpass the second-order component. This may be due to intensive interaction, which is a high-frequency phenomenon that excites the heave motion of paired PA-WECs. Additionally, at a period of 1.8 s, the third component in heaving velocity achieves maximum values of 0.022 m/s and 0.017 m/s for the isolated and paired PA-WECs, respectively. These maximums reach values of 0.023 m/s and 0.024 m/s for both the isolated and paired PA-WECs, respectively, at a period of 2.2 s. Compared to the isolated PA-WEC, the paired PA-WECs have a greater amplitude at a period of 2.2 s, while the isolated PA-WEC has a greater amplitude at 1.8 s.

The amplitudes of the first three harmonic wave loads ($F_{hydro}$) on both isolated and paired PA-WECs are illustrated in Figure 19. A fair agreement is achieved between the CFD model and experimental data. Both isolated and paired PA-WECs experience an increasing trend in the first harmonic component ($F_{hydro}^1$), with the paired PA-WECs experiencing a smaller load, possibly due to more severe viscous effects. The isolated PA-WEC exhibits a relatively higher peak at 1.8 s, while the paired PA-WECs have peaks at both 1.8 s and 2.0 s. The first peak is related to the resonance peak, while the second peak is not. This may be primarily due to the interference caused by a substantial increase in double and triple frequency wave loads. Additionally, the $F_{hydro}^3$ of the pair are generally greater than those of the isolated PA-WEC, possibly due to intense buoy interaction. Therefore, it is necessary to consider harmonic hydrodynamics when investigating performance, as it may concern resonance motion in the regime of low-frequency waves.

4.2. The Influence of Increasing Wave Steepness on Higher Harmonics

The previous section discusses the hydrodynamics associated with both the wave period and the number of floaters. This section examines the impact of wave steepness (kA) on hydrodynamics.

Figure 20 illustrates the heave RAO of both isolated PA-WEC and paired PA-WEC. Only floater #1 of the pair is released, as floater #2 possesses an identical RAO. As kA increases, RAO_hea initially rises rapidly before gradually increasing. The turning point differs slightly between isolated and paired PA-WECs but is significantly influenced by the wave period. At 1.8 s (kh = 5.59), the turning point is kA = 0.15, while at 2.2 s (kh = 3.75), it is 0.124. At 1.8 s, the heave motion of an isolated PA-WEC exceeds that of paired PA-WECs, but at 2.2 s, it becomes lower. This can be attributed to another resonance experienced by paired PA-WECs at 2.2 s. At T = 1.8 s, RAO_hea improves from 0.56 to 0.65, an increase of 15.3% for the isolated PA-WEC, and an increase of 22.7% for the paired PA-WECs. For T = 2.2 s, the RAO_hea increases by 21.6% and 9.1%, respectively, for the isolated and paired PA-WECs. The results show that RAO_hea increases more rapidly at non-resonance periods.

The harmonic wave loads on both isolated and paired PA-WECs are illustrated in Figure 21. The wave loads on the isolated PA-WEC are larger than those on the paired PA-WECs, because the performance of the isolated version is better than the paired at the first resonance. Specifically, the second harmonic wave loads ($F_{hydro}^2$) account for a maximum of 9.8% of total wave loads, and the third harmonic ($F_{hydro}^3$) is 4.9% for the isolated PA-WEC. Regarding the paired PA-WECs, the $F_{hydro}^2$ accounts for 10.3% at a kA of 0.25, while the $F_{hydro}^3$ is evaluated to be 6.8% under the same condition. Although the paired PA-WECs experience less wave excitation, the intensive wave–structure interaction increases the proportion of higher harmonic wave loads.

As the kA increases from 0.1 to 0.25, the $F_{hydro}^2$ experiences a 3.14-fold increase, and the $F_{hydro}^3$ increases by 61.5% for the isolated PA-WEC. For the paired PA-WECs, the $F_{hydro}^2$ experiences a 3.18-fold increase, and the $F_{hydro}^3$ experiences a 1.05-fold increase.

With the wave period increasing to 2.2 s, the paired PA-WECs are found to outperform the isolated PA-WEC concerning both $F_{hydro}^1$, $F_{hydro}^2$, and $F_{hydro}^3$ as demonstrated in Figure 22. The $F_{hydro}^1$ increases linearly with increasing wave kA. The $F_{hydro}^2$ account for a maximum of 8.5% of total wave loads, and the $F_{hydro}^3$ is 3.9% for the isolated PA-WEC. Regarding the paired PA-WECs, $F_{hydro}^2$ accounts for 11.4% at a kA of 0.25, while $F_{hydro}^3$ is evaluated to be 5.2% under the same condition. As kA increases from 0.07 to 0.17, the second harmonic wave loads experience a 2.78-fold increase, and the third harmonic component experiences a 1.67-fold increase for the isolated PA-WEC.

In summary, the higher harmonic components can contribute to as much as 17.1% of the total wave loads when kA is 0.25. When compared to isolated PA-WECs, the paired PA-WECs consistently encounter a higher proportion of harmonic wave loads, potentially attributable to the intensive interaction between buoys leading to wave–wave evolution.

The influence of kA on harmonic heaving velocity is discussed based on the decomposition of the paired PA-WECs’ heaving velocity, as shown in Figure 23. It can be observed that increasing kA benefits the second component. The maximum second harmonic velocity increases from 0.129 m/s to 0.026 m/s as kA increases from 0.12 to 0.25. The second component has a sharper peak and a flatter trough, resulting in a higher maximum and a smaller minimum for the total heaving velocity.

4.3. Power Performance of the PA-WECs

Since the effects of wave steepness and buoy interaction on the harmonic velocity and wave loads of PA-WECs have been thoroughly discussed, it is important to investigate their impact on power performance to gain a deeper understanding of harmonic hydrodynamics.

The influence of the buoy’s interaction on power performance can be detected using the q-factor criterion defined as Equation (19):

$$q-factor=\frac{{\displaystyle \sum _{i=1}^2{P}_{paired}}}{2\times P_{isolated}}$$

(19)

where P refers to the averaged and maximum hydraulic power.

The q-factor of the WEC under two different periods is depicted in Figure 24. The averaged hydraulic power of paired PA-WECs consistently outperforms that of isolated PA-WECs at T = 2.2 s, but only exceeds it under higher wave steepness at T = 1.8 s. This coincides with the previously mentioned resonance periods, where isolated PA-WECs experience resonance at 1.8 s while paired PA-WECs experience resonance at 2.2 s. The q-factor for maximum hydraulic power also follows this pattern, though its variation is significant. The maximum q-factor falls below the average q-factor at 1.8 s but exceeds it at 2.2 s, indicating that differences in hydraulic power performance are primarily due to discrepancies in maximum power output.

The fast Fourier transform (FFT) method is used to decompose the power output of a wave-energy converter, as shown in Figure 25. Like wave loads and heaving velocity, power also consists of three components. When comparing the harmonic components of isolated and paired PA-WECs, we find that the first and second components of isolated PA-WEC exceed those of paired PA-WECs at 1.8 s while the reverse is true at 2.2 s. This observation aligns with the motion responses discussed earlier.

For isolated PA-WEC, the maximum second harmonic hydraulic power can reach 4.67 W with kA = 0.25 and T = 1.8 s, while the third component can achieve 1.26 W under the same condition. The current contribution of the second and third harmonic components to the total hydraulic power is 50.8% and 13.7%, respectively. For paired PA-WECs, floaters perform better at T = 2.2 s with kA = 0.17, where the second harmonic hydraulic power can reach 5.07 W and the third component can achieve 1.63 W. Their contribution can be up to 57.4% and 18.5%, respectively. When compared with a small wave steepness condition, the second component’s contribution increases by 12.2% and 7.1% for the isolated and paired PA-WECs, respectively. The contribution of the third component increases by 4.3% and 4.8% for isolated and paired PA-WECs, respectively. Notably, both isolated and paired PA-WECs have their second harmonic power surpassing their first-order components. Figure 26 shows that the instantaneous hydraulic power of PA-WECs experiences two peaks within a cycle, which is commonly observed in point-absorber systems during nonlinear investigations [28]. This explains this phenomenon.

Additionally, the more significant role of higher harmonics in the hydrodynamics of paired PA-WECs can be attributed to the nonlinear wave evolution between two buoys. Figure 27 shows a noticeable difference in wave elevation detected by wave probes at the left and middle of the paired PA-WECs, particularly at harmonic frequencies. The higher amplitude at harmonic frequencies for the middle wave probe suggests that the interaction between buoys excites the production of harmonic wave evolution, which resonates with point absorbers and significantly affects performance.

5. Conclusions

Higher harmonics hydrodynamics have received less attention due to the linearization of the PTO system. However, the nonlinear wave–structure interaction and wave–wave evolution arising from paired PA-WECs have strengthened the harmonic phenomenon, necessitating a more thorough analysis. In this study, a hydraulic PTO system was incorporated into the experimental and CFD model of paired PA-WECs. A comprehensive analysis of higher harmonic wave loads, heaving velocity, and power performance was conducted to examine the relationship between harmonic characteristics and hydrodynamics.

Research indicates that higher harmonic wave loads can contribute to as much as 17.1% of the total wave load, specifically related to the additional resonance observed in the paired PA-WECs. With an increase in wave steepness, the harmonic wave loads experienced by paired PA-WECs exhibit a more rapid ascent compared to those of an isolated PA-WEC, attributed to intensive interaction among buoys. Furthermore, owing to the pronounced peaks and shallower troughs at higher harmonic velocities, the heaving velocity of the resulting absorber accelerates more rapidly and decelerates more gradually, thereby exerting a substantial impact on its power performance. The second-order and third-order hydraulic power can account for up to 57.4% and 18.5% of total power output, respectively, due to the PTO system’s nonlinear characteristics and the influence of harmonics. Harmonic components in the wave elevation between two floaters are more evident, implying significant harmonic characteristics of paired PA-WECs.

In conclusion, assessing the performance of paired PA-WECs, particularly in long wave periods with higher wave steepness, highlights the importance of harmonic characteristics.