Geometry Optimisation of a Wave Energy Converter

Abstract

The geometry optimisation of a point-absorber wave energy converter, focusing on the increase in energy absorption derived from heave forces, was performed. The proposed procedure starts by developing an initial geometry, which is later evaluated in terms of hydrodynamics and optimised through an optimisation algorithm to tune the shape parameters that influence energy absorption, intending to obtain the optimal geometry. A deployment site on the Portuguese coast was defined to obtain information on the predominant waves to assess several sea states. NEMOH and WEC-Sim (both open-source software packages) were used to evaluate the interaction between the structure and the imposed wave conditions. The results extracted and analysed from this software included forces in the six degrees of freedom. Under extreme wave conditions, the highest increase in the relative capture width between the initial and final shapes was around 0.2, corresponding to an increase from 0.36 to 0.54, while under average wave conditions, the increase only reached a value of around 0.02, corresponding to an increase from 0.22 to 0.24, as calculated through the relative capture width values.

1. Introduction

The ocean covers around 70% of the Earth’s surface [1] and is the supplier of many resources, ranging from biological to energetic assets. This energy presents itself in many forms, allowing many ways to extract and predict it. The primary energy sources in the ocean are tidal, ocean thermal, and wave energy. Waves are generated when the wind interacts with the sea’s surface, transferring the energy it has to the area of the ocean over which it acts. This energy is stored in two components, kinetic and potential, corresponding to the motion and the elevation of the water, respectively [2]. Globally, this resource is estimated to have the potential to generate 8000–80,000 TWh/yr [3]. In Portugal, a country bathed by the ocean, there is also incredible potential for this resource. In 2014, the value for this potential was estimated, concluding that the northwest coast of Portugal is more energetic than the southwest, with values for the annual average wave energy available reaching 200 MWh/m and 150 MWh/m, respectively [4].

The devices used to convert this energy—wave energy converters (WECs)—can be divided into several categories according to their working principle, of which the point absorber is the most popular, included in up to 83 projects worldwide according to the European Marine Energy Centre (EMEC) [5]. These devices are characterised by their small dimensions when compared to the wavelength of the incoming waves [6], using the relative motion of composing parts, induced by the interaction between them and the wave, to produce power [7]. Several issues commonly appear in all categories, such as the durability of the devices at sea, technological infancy, power take-off (PTO) systems which can deal with energy source variability, and also the financial and resource requirements to test devices at sea [8,9]. Device optimisation may provide answers to these challenges, and geometry optimisation has been particularly interesting from this point of view since most of the cost reductions can be related to the structure of WECs [10].

There have been several studies on the optimisation of WECs using different methods. In the case of [11], the analysis of a device consisting of a sphere-shaped buoy, a fully submerged non-buoyant sphere-shaped supplementary mass, and a PTO hydraulic unit lying on the seabed was performed. It involved applying a control technique to adjust the device’s natural period to resonate with the incident waves. The optimisation was performed by varying the diameter and the PTO’s external damping of the buoy. The results allowed us to conclude that the absorbed energy distribution was not significantly different for each configuration under the same wave conditions, which means that the deployment site is an important parameter to consider during optimisation. Shadman et al. (2018) [7] optimised a point absorber under three requirements, namely the buoy’s heave natural frequency, the resonance bandwidth, and the maximum mechanical power, using two parameters—the buoy diameter and draft. The study used the experimental design method, defining upper and lower bounds for both geometrical parameters. McCabe (2013) [12] used a genetic algorithm to optimise a WEC with fewer control points, generating bi-cubic B-spline surfaces and evaluating them through different cost functions. The results showed that considering the size and the mean power delivery in optimisation led to substantially improved performance. Another study, performed by Gilloteaux and Ringwood (2010) [13], obtained the optimisation of a point absorber considering the interaction between the geometry of the device and the control strategy adopted for it, using a deterministic algorithm to solve the problem.

An optimisation approach may address many of the challenges identified in the literature. Geometry has great potential to decrease overall costs and maintenance requirements since the ocean provides extreme conditions for WEC systems. This study aims to maximise energy absorption by the device while simultaneously increasing the durability of the buoy to decrease the need for maintenance/repair and increase the device’s lifespan. In this case, a point-absorber converter is the study’s subject, so the device’s heaving motion is considered the main driving force of energy extraction. Section 2 describes the methodology and assumptions for this study; Section 3 focuses on presenting the results from the simulations, and lastly, Section 4 includes a conclusion of the work developed.

2. Methodology

The main goals for this study were to optimise the shape of a heave point-absorber WEC, increase energy absorption, and minimise stresses on the structure. The location chosen as a possible testing location was offshore of Viana do Castelo, in northern Portugal, as represented in Figure 1. This area has been previously chosen as a testing site for other forms of WEC technology testing. The wave characteristics were taken from the Global Ocean Waves Analysis and Forecast data set from the database Copernicus.

The extracted data involved the following parameters: significant wave height (H_m0) obtained through spectral analysis, the wave peak period (T_p), and mean wave direction (θ), corresponding to the direction from which the waves are coming. The time series obtained for each parameter was analysed, and the respective statistics for the chosen period are defined in

Table 1. Statistics for the wave parameters in the analysed period.

	Hm0 (m)	Tp (s)	q (°)
Minimum (SS1)	0.54	3.87	193.19
5th Percentile (SS2)	0.78	7.20	257.60
Root Mean Square (SS3)	2.00	11.49	302.84
95th Percentile (SS4)	3.54	15.59	332.57
Maximum (SS5)	5.83	19.72	353.58

.

The extreme wave conditions chosen for this study were based on the statistics for H_m0. The 95th percentile and maximum values for this parameter were 3.54 and 5.83 m, respectively, and the value of 4.78 m was chosen as an intermediate between these two, as the 95th percentile value was found for February. The corresponding peak period was 10.82 s for θ, equal to 280.72°. Each sea state described in Table 1 is referenced with a corresponding code—SS—for the sea state and numbering for each sea state. For the average wave conditions, the RMS values for the parameters were used (corresponding to SS3). The wave spectrum for the conditions in Table 1 was calculated using the Pierson–Moskowitz spectrum through the open-source software MHKiT-MATLAB, v0.6.0 [14]. The wave spectra are presented in Figure 2.

2.1. Hydrodynamic Model

Among the different frequency-domain boundary-element method (BEM) codes available, NEMOH [15] stood out as a reliable open-source code for obtaining the required hydrodynamic coefficients as the input for WEC-Sim, as shown in published studies such as those by [16,17]. Its algorithms are run using the MATLAB (R2021a) wrapper provided by the developers. The choice of the BEM pre-processor is also based on the good predictions of the hydrodynamic parameters by NEMOH for the case of simple marine structures, as shown in [18], where the study compares the use of WAMIT, NEMOH, and HAMS. Even though HAMS, as the other open-source BEM code, showed a better performance in terms of efficiency and accuracy for more complex marine structures, in the case of the present study and the low level of complexity of the structure, NEMOH presents as a better option when also considering the integration within WEC-Sim.

The input values for NEMOH include the buoy’s mesh parameters and the wave characteristics. The initial shape analysed is represented in Figure 3. Its draft value is 10 m, and the radius is 5 m. These values remained constant throughout the optimisation process since the main goal was to evaluate different shapes.

The shape chosen to start the process was similar to a ship’s hull, as a shape that would lead to lower resistance and displacement in the water and a lower impact from viscous effects [19], and the refinement of the mesh was evaluated to assess if it had significant impacts on the results. Using a coarser mesh reduced computational effort and time. The main parameters used to define the mesh in NEMOH included an array of radial coordinates, an array of vertical coordinates, the number of points for discretisation, the number of points for angular discretisation (ntheta) (shown as the red and green circles in Figure 4), and the target for several panels (nfobj) (as seen by the panels in Figure 3). This evaluation focused on ntheta and nfobj. The characteristics of the validation simulations are described in

Table 2. Characteristics of the validation simulations.

	Simulation 1	Simulation 2	Simulation 3
ntheta	6	11	6
nfobj	500	300	300

.

In these simulations, the conditions considered were regular waves, with a 2.5 m wave height and an 8 s period. The simulation time was set to 400 s. The results are shown in

Table 3. Results from the validation process, comparing the differences between the performed simulations. The values correspond to the relative difference between each pair of simulations, as defined by Equation (1).

	Surge Forces (%)	Heave Forces (%)	Pitch Forces (%)	Energy (%)
Simulations 1–2	1.91	0.00	0.88	0.72
Simulations 2–3	1.25	0.00	3.49	0.72
Simulations 1–3	0.64	0.00	2.64	0.00

, describing the relative differences in results between using the different ntheta and nfobj parameters, as defined by Equation (1):

$$RD=\left|{\displaystyle \frac{V_{Sim-i}-V_{Sim-j}}{V_{Sim-i}}}\right|\times 100$$

(1)

where $RD$ is the relative difference in the percentage, $V_{Sim-i}$ is the value for each variable for simulation number i, and $V_{Sim-j}$ is the value for each variable for simulation number j (different from i).

The values obtained show no significant difference when considering the total forces in heave, pitch, and surge forces, as well as in the energy absorbed by the device. Therefore, the following simulations were considered with a coarser mesh, using the parameters of simulation 3 in Table 2.

The results from the NEMOH simulation were then analysed through a BEM input/output function as part of the WEC-Sim v4.1 software package to transform them into readable inputs for WEC-Sim, building a .h5 file with them.

WEC-Sim formulates the study problem in the time domain and is a useful open-source tool for case study analysis. It is developed in MATLAB and Simulink (R2021a), with the structural model developed with the WEC-Sim Simulink library blocks. In this case, the structural model was built based on the RM3 case study provided by the developers. The PTO system was assumed to be linear.

WEC-Sim is then run using an input file with the remaining data, which includes simulation data, wave information, body data, and PTO and constraint parameters. The simulator outputs wave information (time and elevation) and body and PTO output. The results concerning body and PTO are provided in six degrees of freedom (DOF). The wave parameters considered in WEC-Sim were concordant with the NEMOH settings, using the irregular waves option and a Pierson–Moskowitz spectrum.

The use of potential flow theory in this study has limitations. Although it leads to a lower computational effort, disregarding non-linear effects, such as viscous effects, may affect the results and conclusions. Particularly, when analysing near the resonance period of the device, special consideration must be given to drag forces and viscous effects [20]. Even though the focus of this study is not optimising the resonance region for each device but rather for a particular site’s conditions, these concerns should be properly addressed. Several studies have considered this limitation and how much it affects the validity of using potential flow theory [19,20,21,22]. In a study which focuses on heaving point absorbers [19], an analysis of the viscous effects was performed with different shapes, such as cylindrical, conical, and hemispherical configurations. It was shown that the correction carried out to the damping coefficients obtained from potential flow theory to account for the viscous effects was highly dependent on the shape of the WEC. The viscous effects greatly impacted the cylindrical-shaped WEC, while the other shapes showed a lower need for correction. Particularly, the hemispherical-shaped WEC was presented as the least affected by viscous effects, followed by the conical shape with an angle of 90 degrees. Special attention is also given to cases with higher wave heights and near resonance periods, where the potential flow theory may not portray a realistic analysis of the WEC’s response, as seen in [21].

Considering the described case studies, the impact of the viscous effects and drag forces are highly dependent on the WEC’s shape and the study’s focus—whether to optimise around resonance periods or not. Since, in this study, the initial shape and its variations through the optimisation process were around a conical shape, the viscous effects were assumed to have a low impact on the overall response.

2.2. Optimisation Algorithm

The suitability of each optimisation algorithm depends on the type of study being developed, its complexity and scope, and the main goals. Considering that a wider scope was intended in this study, a genetic algorithm was chosen to optimise and find a global minimum instead of a local minimum.

In this case, the genetic algorithm from the Global Optimisation Toolbox provided in MATLAB was used [23]. This type of algorithm copies the natural selection process from biological evolution. It keeps modifying the individuals in each population and selecting them randomly to produce the next generations. The genetic algorithm uses different methods to create each population and generate each individual, considering previous generations, namely crossover, mutation, and selection. The algorithm requires the definition of a fitness function, which is then applied to the individual and attributes a fitness value to it. This algorithm provides a significant advantage in the search for global optimum values instead of focusing on local values.

The parameters chosen for the optimisation process were initially given to NEMOH, as represented in Figure 4. The green points were optimised, and the red ones were fixed.

The fitness function was developed with the following distribution: 15% for the decrease in the total force in the surge DOF; 15% for the decrease in the pitch DOF; 30% for the increase in the power absorbed by the PTO system; and 40% for the increase in the total force in the heave DOF. This distribution was defined according to the study’s goals, as the heave forces the focus of the process and maximises the energy absorption. To reduce energy losses in other DOFs, the fitness function also considers decreasing the total forces in the surge and pitch. The fitness value attributed to each solution defines the optimisation course. The variables optimised are calculated by Equations (2) and (3) [24]:

$$P_{pto}=K_{pto}{X}_{rel}{\dot{X}}_{rel}+C_{pto}{{\dot{X}}_{rel}}^2$$

(2)

where $P_{pto}$ is the instantaneous power absorbed by the PTO system, $K_{pto}$ is the PTO stiffness, $C_{pto}$ is the PTO damping, $X_{rel}$ is the relative motion between two bodies and ${\dot{X}}_{rel}$ is the relative velocity between two bodies.

$$m\ddot{X}=F_{exc}\left(t\right)+F_{rad}\left(t\right)+F_{pto}\left(t\right)+F_v\left(t\right)+F_{me}\left(t\right)+F_B\left(t\right)+F_m\left(t\right)$$

(3)

In the above equation, $m$ is the body mass matrix, $\ddot{X}$ is the acceleration vector of the device, $F_{exc}\left(t\right)$ is the wave excitation force/torque vector, $F_{rad}\left(t\right)$ is the wave radiation force/torque vector, $F_{pto}\left(t\right)$ is the PTO force/torque vector, $F_v\left(t\right)$ is the damping force/torque vector, $F_{me}\left(t\right)$ is the Morison element force/torque vector, $F_B\left(t\right)$ is the net buoyancy restoring force/torque vector, and $F_m\left(t\right)$ is the mooring connection force/torque vector. A synthesis of the workflow is presented in Figure 5.

3. Results and Discussion

After the mesh validation, the optimisation simulations were run. Three conditions were tested: (1) optimisation under irregular waves (extreme conditions), (2) optimisation with a cylinder as the initial shape, and (3) results from testing under average wave conditions.

3.1. Irregular Waves—Extreme Conditions

Choosing a more extreme scenario helps determine the maximum energy extraction values while showing how the device performs in a harsher environment. The settings for the algorithm were a maximum of 100 generations and a population size of 10. The maximum of stall generations was also defined and given a value of five. The genetic algorithm stops if there is no improvement in the fitness function results after the number of stall generations. The initial shape’s mesh is represented in Figure 3. The solution which gave the best result of the fitness function was achieved at the 957th iteration, with a shape as described in Figure 6, and the results from this improvement are shown in

Table 4. The results from the optimisation process under extreme conditions. The variation is considered the relative difference between the initial and final shapes’ results, as defined by Equation (5).

	Surge Forces (N)	Heave Forces (N)	Pitch Forces (N)	PTO Power (W)	Energy (kWh)
Initial shape	4.00 × 104	5.25 × 105	8.46 × 104	3.71 × 105	41.26
Final shape	4.40 × 104	6.09 × 105	8.95 × 104	5.64 × 105	62.72
Variation	10.00%	16.00%	5.79%	52.02%	52.02%

.

The fitness values are shown in Figure 7. This figure shows the evolution of the fitness values along the number of analysed generations. The curve defines the mean fitness values, and the dots define the best fitness values of each generation.

The main comparison between the two shapes was made considering the root mean square of each parameter. There was an increase in the total force in the three DOFs analysed (surge, heave, and pitch) and in the PTO-absorbed power (Figure 8) and energy. The relative difference shown in Figure 7 is calculated similarly to Equation (1), with the following correction:

$$RD=V_{FinalP}-V_{InitialP}$$

(4)

where $V_{FinalP}$ corresponds to the values of the final population and $V_{InitialP}$ corresponds to the values of the initial population.

Although one of the goals was to decrease the forces in surge and pitch, given that the increase in heave was higher and around ten times the magnitude of the former, the results are considered positive.

The body’s position is also shown in Figure 9. It is important to note that there are two y-axes in the body position graph, the left corresponding to the wave elevation and the right to the body’s position. While the wave elevation ranges from −5 to 5 m, the body’s position, with the graph at its centre of gravity, only ranges from −4.5 to 0.5 m, about half of the range of the wave elevation.

Table 4 also shows an increase of 52.02% in energy absorption, which is a more important factor than the impacting forces. The relative difference in Table 4 is calculated similarly to Equation (1) but with the corrections shown in Equation (5):

$$RD=\left|{\displaystyle \frac{V_{Sim-i}-V_{Sim-f}}{V_{Sim-i}}}\right|\times 100$$

(5)

where $V_{Sim-i}$ is the value for the initial population (initial shape), and $V_{Sim-f}$ is the value for the final population (final shape).

After analysing these results, the impact of the first shape fed to the algorithm was also investigated to see if the final shape was similar or even the same as the one in this simulation.

3.2. Cylinder Case Study

The influence of the first shape chosen for the algorithm run was studied to evaluate whether this initial configuration led to significantly different results under the same conditions. A cylinder shape was chosen, as shown in Figure 10a. It is important to consider, as mentioned in Section 2.1, that cylindrical shapes are shown to be impacted by viscous effects, which are not included in this study and may affect the results. The wave conditions remained the same as in the previous simulations, as did the fitness function (extreme conditions). The settings for this algorithm, when run, were also the same as the previous algorithm, with a maximum of 100 generations, five stall generations, and a population size of 10. The best results were obtained after the 665th iteration, with the shape shown in Figure 10b.

The final shape is similar to the one resulting from the previous run. The results were slightly lower than the previous ones regarding power absorption (Figure 11) but showed an improvement in impact forces. Although the cylinder shape showed a poorer performance when compared to the initial shape in the previous section, the genetic algorithm still reached similar results.

Table 5. Results from the optimisation process for the cylinder case study. The variation is considered the relative difference between the cylinder and final shapes’ results, as defined by Equation (5).

Shape	Surge Forces (N)	Heave Forces (N)	Pitch Forces (N)	PTO Power (W)	Energy (kWh)
Cylinder	1.08 × 105	3.46 × 105	2.16 × 105	1.53 × 105	17.01
Final	4.30 × 104	5.97 × 105	9.33 × 104	5.45 × 105	60.60
Variation	60.19%	72.54%	56.81%	256%	256%

represents the results of this simulation.

The variable values in this algorithm, when run, showed a significant improvement from the cylinder shape, with decreased surge and pitch forces and increased heave forces. In terms of energy, there was a remarkable increase of 256% due to the energy values of the cylinder being significantly lower than the shape introduced in the first algorithm used. As shown in Table 4, the initial shape showed an energy value for that run of around 41 kWh, whereas the cylinder shape had a value of around 17 kWh.

3.3. Average Wave Conditions

After testing the impact of the chosen first shape, the results from the first optimisation run were tested under average wave conditions. The values for the wave conditions are described in the Methodology section. Evaluating the device’s performance in average wave conditions provides a more realistic assessment since it is the scenario it will face most of the time. The results from the simulations are presented in

Table 6. Results from the optimisation processes’ final shape under average wave conditions. The variation is the relative difference between the initial and final shapes’ results, as defined by Equation (5).

Shape	Surge Forces (N)	Heave Forces (N)	Pitch Forces (N)	PTO Power (W)	Energy (kWh)
Initial	4.26 × 104	1.69 × 105	8.03 × 104	4.29 × 104	4.77
Final	4.88 × 104	1.70 × 105	7.89 × 104	4.75 × 104	5.28
Variation	14.55%	0.59%	1.74%	10.72%	10.72%

.

Comparing these results with the extreme wave conditions, it is possible to verify a significantly lower difference (Figure 12) since the variation in energy values was five times lower. There is also an increased impact of the pitch and surge total forces, and the difference in the heave total force was slight. Compared to the previous simulations, the values for pitch and surge total forces appeared higher when compared to heave forces, which meant that their expression in the movement of the device increased. This may lead to a lower level of extracted energy since the considered PTO system extracts energy solely from heave forces.

3.4. Resource Evaluation

Evaluating the resulting shapes in terms of efficiency requires comparing the PTO-absorbed power, P_PTO, and wave power. It is possible to calculate the power density of the wave, P_dens, using Equation (6) (Guillou, 2020) [25]:

$$P_{dens}={\displaystyle \frac{\rho g^2}{64\pi }}{H_s}^2\alpha T_p [{\displaystyle \frac{\mathrm{k}\mathrm{W}}{\mathrm{m}}}]$$

(6)

and then use the value of P_wave:

$$P_{wave}=P_{dens}W \left[\mathrm{k}\mathrm{W}\right]$$

(7)

where W corresponds to the device’s width—10 m—to calculate the relative capture width, RCW, according to Equation (8).

$$RCW={\displaystyle \frac{P_{PTO}}{P_{wave}}}$$

(8)

Calculating P_dens requires values for H_m0 and T_p. For each of the sea states analysed, the wave power density value per unit length of the wavefront is presented in

Table 7. Parameter definition and results for wave power density.

	Hm0 (m)	Tp (m)	Pdens (kW/m)	Pwave (kW)
Extreme Wave Conditions	4.78	10.82	104.09	1040.90
Average Wave Conditions	2.00	11.49	19.38	193.80

. The value for α was assumed to be 0.86 since the chosen type of wave spectrum was Pierson–Moskowitz [26].

The RCW values were calculated using the values for the absorbed PTO power for each of the simulations provided in the tables from previous sections. The results are shown in

Table 8. Relative capture width values for the optimisation results.

	PPTO (kW)	Pwave (kW)	RCW (%)
Initial shape—Extreme conditions	371.00	1040.90	36
Final shape—Extreme conditions	564.00	1040.90	54
Initial shape—Average conditions	42.90	193.80	22
Final shape—Average conditions	47.50	193.80	24
Cylinder shape	153.00	1040.90	15
Final shape—Cylinder case study	545.00	1040.90	52

. There is a significant increase in RCW from initial to final shapes under extreme wave conditions, both with a more conical initial shape and a cylinder, with the latter suffering a higher increase. Regarding the results under average wave conditions, the values are both low and slightly increased between the initial and final shapes. This may be due to these simulations’ high pitch and surge force values, as mentioned in the previous section. The difference between the cylinder case study’s initial and final shapes was significantly high because the initial shape results were poor compared to the first case study, which did not use a cylinder as the initial shape. When the final shapes from these two scenarios were compared, it was possible to observe a similar energy value.

Although the increase in RCW is extremely positive, the values may differ when tested in actual conditions since these simulations do not consider the non-linear aspects of wave–structure interactions. Since the optimisation process is considered to have a high wave height, important effects may be disregarded by the potential flow theory, as discussed in Section 2.1. Comparing the obtained values with other studies, it is possible to verify that the simulations with average wave conditions are the ones that approach the average values for efficiency for heaving wave energy converters, even though some may reach values as high as 51% [27].

4. Conclusions

The optimisation process produced satisfactory results in terms of absorbed power. For the same operating conditions, the energy absorbed increased by around 52% in the final geometry. This increase was smaller for the average wave conditions—around 11%—which may be because the final shape was adjusted to the most severe wave conditions, and, therefore, an optimisation process dedicated exclusively to these conditions should be considered, as the wave conditions had an impact on the results. However, it was impossible to reduce the effect of surge and pitch forces, as intended, to reduce possible stress in the device and improve its robustness. The existence of different sea states in a particular site can be incorporated into the optimisation process using recorded data and statistical analysis to differentiate the probability of occurrence for each sea state. Upon defining these values, the different sea states can be evaluated in each of the optimisation algorithm’s iterations, and using the probability of occurrence for each sea state, their impact on the results can be attributed.

This attribution can be incorporated into the fitness function to consider not only one sea state but several, depending on the likelihood of their occurrence. Further work should also include the motion of other degrees of freedom. The focus on the PTO system can also improve the device’s performance in terms of power production and should, therefore, be considered as the next stage in the optimisation process. A process focused on the co-design of the WEC could also enhance power optimisation and should be considered in future work.

The variation in the initial shape did not result in relevant changes since the final geometry resulting from the optimisation process reached a similar absorbed power starting from a cylinder and a conical shape. The relative capture width values obtained may be higher than expected since non-linear effects were not considered. The optimisation process was automated and successfully implemented, allowing its use in other studies.