Investigation of new layout design concepts of an array-on-device WaveSub device

Wave Energy Converters (WECs) have not yet proven their competitiveness in the mainstream energy market. Research and development of this technology are necessary to ﬁ nd optimal solutions in terms both of energy produced and reduced cost. A WEC farm is expected to have reduced Levelized Cost of Energy (LCoE) compared to individual devices due to shared installation and grid connection costs. Studies show that energy yield of a WEC array is highly dependent on spacing and layout of the WECs. A method for selecting an optimal array layout is desirable. Here we show a comparison between 4 different design layouts of a WaveSub device with six ﬂ oats. A six ﬂ oat con ﬁ guration has been chosen because the LCoE reduces with increasing ﬂ oats per device as shown in previous research. An optimal con ﬁ guration in terms of LCoE and rated power is found for linear, rectangle, triangle and circular multi-ﬂ oat con ﬁ gurations. Parameters optimised are ﬂ oat spacing and Power Take Off (PTO) stiffness, damping and rated power. The optimisation algorithm uses a genetic algorithm combined with a Kriging surrogate model. Numerical simulations are solved in the time-domain in WEC-Sim while the hydrodynamic coef ﬁ cients are calculated in Nemoh using a linear potential ﬂ ow theory. For all geometric con ﬁ gurations, the smallest ﬂ oat spacing was the most promising because of the lower cost of the structure. In fact, the in ﬂ uence of the ﬂ oat spacing on the power produced by the device is shown to be less signi ﬁ cant than the in ﬂ uence of ﬂ oat spacing on the capital cost. Overall, the circular con ﬁ guration outperformed the other con ﬁ gurations. This study shows that layout con ﬁ gurations can be investigated with optimisation and this could be applied to other con ﬁ gurations and other WEC concepts in future. © 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Wave energy technology is a state-of-art renewable energy extraction approach that captures the energy of ocean waves to produce electricity.Wave energy is very promising because it accounts for most of the energy of the sea and it could provide a significant contribution of between 10 and 25% of the world's electricity consumption [1].The development of this technology is challenging because of the survival conditions required, the high costs and the complex optimal design process.[2] reviews various types of Wave Energy Converters (WECs) that have different working principles.The WaveSub device is an off-shore orbital absorber developed by Marine Power Systems Ltd.(MPS) characterized by a completely submerged reactor and floats [3].Its working principle can be found in Ref. [4] while a multi-float WaveSub is described in Ref. [5].Fig. 1 shows a multi-float Wave-Sub device during power production.This shows three floats in a linear configuration.Wave energy is first converted into mechanical energy thanks to the orbital motion of the floats generated by the hydrodynamic loads.Each float is connected to the reactor via four PTO lines, and the rotational motion on the reactor's pulley for each line is converted into the translational motion of the piston of a hydraulic PTO.This energy is converted into fluid power and stored in the accumulators.Finally, a hydraulic motor and generator are used to produce constant and continuous electrical power rating.A taut mooring provides a reliable and stable base to exploit the relative motion between the float and the reactor and maximise the energy generated by the floats.
This concept still needs some improvements to be competitive in the energy market.There are several barriers to its development related to the risk and reliability of the technology, LCoE, government policies and access to private investment.In this study, only the reduction of LCoE is investigated, as this is one of the most important factors for the success of a technology.A multi-float configuration has some economic advantages over a single float configuration due to a reduction in the costs of installation, grid connection and moorings.However, the energy produced by each float can be significantly affected by the presence of the other floats.Therefore, it is very important to design the layout of the floats in such a way that the interaction factor of the array is as large as possible to maximise the total amount of energy generated.To this end, different layout configurations are investigated in this work in order to optimise the performance of the floats and reduce the overall cost of the system.
A multi-float configuration can be considered as an array-ondevice, while a WEC array consists of multiple multi-float devices organized in an array.This work has some similarities with the optimisation of a WEC farm.However, the optimisation is limited to the layout of the multi-float configuration on a single device with multiple MW ratings.It is assumed that the distances between the devices are such that the shading is negligible.It is also assumed that the total nominal power of the farm is constant, so that the number of devices depends on the optimisation.
Metaheuristic algorithms are most often used for the optimisation of a WEC farm, especially when multiple optima are expected.In fact, there is a lot of literature on wave energy park optimisation using the genetic algorithm (GA).For instance, in Ref. [6], the layout and PTO of a point absorber array were optimised using MATLAB's GA optimisation tool and compared with the Parabolic Intersection (PI) method.Optimal layout of arrays of WECs with different sizes were optimised in Ref. [7] using a GA optimisation tool.The layout of a 5 point absorber WEC array was optimised using a binary GA in Ref. [8].Neural network optimisation has also been used to optimise the layout of the WECs in a size constrained farm [9] while a machine learning algorithm was applied to optimise a WEC farm of 40 devices subject to bathymetry and space constraints in Ref. [10].Metaheuristic algorithms were also compared to verify their performance, including GA, particle swarm optimisation (PSO) and Covariance Matrix Adaptation Evolution Strategy (CMA-ES).[11] compared the GA and PSO when applied to the WaveSub device in 3 different configurations.Neither algorithm in this study outperformed the other in all the configurations considered.[12,13] have shown that GA was more accurate than CMA-ES, but more computationally expensive.A new hybrid cooperative coevolution algorithm obtained larger efficiency when compared with several optimisation strategies such as CMA-ES, differential evolution (DE), particle swarm optimisation (PSO), Grey Wolf Optimiser (GWO) and Nelder-Mead simplex direct search [14].
Finally, the optimal layout of the WECs shows different conclusions between different research.[15] has shown that a square configuration of 4 heave floating-point absorbers is better in terms of power production and less sensitive to the wave direction compared with rhombus and triangular configurations.On the other hand, triangular WEC arrays with 9e25 devices whose position is defined by 2 parameters (dx,dy) are preferred in Ref. [16] as they offer the flexibility to achieve a better optimum between destructive and constructive effects.Oscillating Wave Surge Converters (OWSCs) in a farm with a staggered configuration have shown better performance than an inline configuration for a random sea analysis in Ref. [17].Further work has shown that a staggered configuration of 5 OWSCs in a "bowl" or "chevron" configuration is more efficient than a linear configuration for a model sea state [18].Unequal spacing of a linear WEC array was also investigated.It was found that the performance is higher when the devices are closer together at the front of the array, due to the positive effect of the radiated field [19].

Aims of this work
This paper aims to optimise the layout of the array-on-device float configuration of the WaveSub device in terms of Levelised Cost of Energy (LCOE).Four main configurations were considered for this work: linear, rectangular, triangular and circular.For each layout, three main float spacings have been investigated from the smallest practicable value to the largest (described as a multiple of the radius of the absorber).The number of floats was set to 6 because the LCOE reduces with an increased number of floats and because multi-float configurations larger than 6 become less feasible due to the size (E [20].A genetic algorithm supported by the Kriging surrogate model was used to optimise the PTO parameters for each layout and float spacing.The Kriging model increases the convergence speed to the optimal solution as shown in (E [20] where it was used to optimise a linear WaveSub configuration.Parallel computation has been used to reduce the computational time.

Methodology
The numerical methodology is based on different steps to numerically simulate the wave device.First a CAD/mesh software is used to obtain the mesh of the different hydrodynamic bodies describing the wave device, namely the floats and the reactor.Salome-Meca [21] has been used because it is open-source and a conversion of the mesh format to Nemoh [22] can be easily implemented.Nemoh is used to obtain the hydrodynamic coefficients based on linear potential flow theory and a Boundary Element Method (BEM).Linear potential flow theory assumes an irrotational and incompressible fluid, small wave amplitudes compared to the wavelength and a small motion of the device compared to its characteristic dimension to have negligible viscous forces.Therefore, the Navier-Stokes equation is reduced to the Bernoulli equation: where 4 is the fluid potential, p is the pressure, r is the density of the fluid, g is the acceleration of gravity and z is the elevation.In particular, the linear theory neglects the quadratic term of this equation.
The fluid potential ( 4) is defined as where V is the flow velocity.The fluid potential can be expressed also as the sum of the incident (4 i ), scattered (4 s ), and radiation potential (4 r ).
The BEM is then applied using Green's function to solve the radiation and diffraction problem and considering the boundary conditions on free surface, bottom and body.
Nemoh is able to take into account the hydrodynamic interaction.The interaction between the hydrodynamic bodies considers the six degrees of freedom (DoF) of each hydrodynamic body for both the diffraction and radiation wavefields.The excitation (F ex ðf Þ) and radiation forces (F rad ðf Þ) acting on each hydrodynamic body can be expressed as where n is the unit vector normal to the body surface.Hydrodynamic interaction of the radiation force is accounted for using the added mass coefficient (A ðf Þ2R NBNDoF ÂNBNDoF ÂNF ) and the radiation damping coefficient (B ðf Þ2R NBNDoF ÂNBNDoF ÂNF ) where N B , N DoF and N F are the number of hydrodynamic bodies, degrees of freedom and frequencies respectively.The hydrodynamic coefficients resulting from this calculation are subsequently used for the dynamic system simulation in WEC-Sim [23].The equation of motion is solved in the time-domain for the 6 of freedom of each body based on the Cummins equation [24]: where m is the mass matrix, A ∞ is the infinite-frequency added mass matrix, X is the translational and rotational displacement vector of the body, K is the radiation impulse response function, F ex is the wave excitation force and torque vector, F v is the quadratic viscous drag force and torque vector, F res is the net buoyancy restoring force and torque vector, F PTO is the PTO force and torque vector and F m is the mooring force and torque vector.The total PTO force on each float is the sum of the forces of the four PTO lines, while the total mooring force on the reactor is the sum of the forces of each mooring line.The PTO force for each PTO line is expressed as where k PTO and c PTO are the PTO linear stiffness and damping coefficients, x is the PTO line linear position and v is the PTO line linear velocity.The mooring force is expressed only as a stiffness term.PTO and mooring lines also include universal joints that provide three rotational degrees of freedom at the lower and upper ends of each line.Finally, the power produced is determined by post-processing the simulation data.The total power produced (P tot ) is calculated as T sim (7) where n f is the number of floats, N PTO is the number of PTO lines, T sim is the simulation time, Dt and N t are the timestep and their number respectively.In this calculation, the ramp time used to achieve a gradual increase to the full flow condition is subtracted from the total simulation time and only the full operating condition is considered.
The optimisation methodology of (E [20] was modified for this research.The main difference is that the hydrodynamic computation is calculated before the optimisation rather than during the process.There are two different optimisers nested within each other, as shown in Fig. 2. The outer process computes the power matrix using a timedomain simulation in WEC-Sim.The second inner optimiser is not affected by the time-domain simulation and is used to optimise the rated power in a post-processing step.In this second step, two calculations are performed: The energy matrix is determined as the product of the power matrix and the occurrence matrix, and a commercially relevant WaveSub array cost calculation is performed.The rated power limits the maximum possible power from the sea states (increasing LCOE) while simultaneously limiting the PTO cost (reducing LCOE).The overall rated power of the farm is fixed, so the number of devices changes as a consequence.The optimiser determines the lowest LCOE in this step and returns it to the outer time-domain optimiser.The LCOE calculation is described in more detail in Section 2.2.4.
The design variables considered in the outer optimisation are float spacing, PTO spring and PTO damping, while in the inner optimisation the rated power is considered.The float spacing is discretised between 3 possible values, while the PTO parameters are optimised within a certain range.They are defined as follows: The optimisation of the design parameters has been simplified by associating a different optimisation run to each float spacing value, so that the optimisation algorithm of the first optimiser only considers two variables.The range of PTO parameters is based on previous experience of MPS but a wide range is specified to maximise the search space.The optimisation algorithm used for the first optimiser is the genetic algorithm combined with the Kriging surrogate model as described in (E [20].The population is considered as 5 times the number of parameters to be optimised (10 individuals), while the number of generations is limited to 5.More precisely, each individual of the population is described by the PTO parameters.The Kriging surrogate model has been combined with the genetic algorithm to reduce the number of simulations required to reach the optimum because it can be used to estimate the simulation result where not simulated.In particular, the next population of the genetic algorithm is determined by an elitism parameter, where new individuals are selected from the best ones of the surrogate model.The elitism factor is chosen very high (0.9) because the optimisation function has been found to be convex and with a single optimum.For this purpose, the Oodace toolbox was used [25].The efficiency of the surrogate model can be obtained using the Normalized Root Mean Square Error (NRMSE) and the maximum error (MAX), as also previously described in (E [20].They are evaluated for a certain number of testing points (n p ) and are expressed as where y i and b y i are the real and predicted value of the response.This algorithm has been chosen because it is derivative-free and a global minimum search algorithm.Furthermore, it can be easily used in a parallel computation because the individuals of each population optimisation step can be run in parallel.The code has been run in MATLAB/Simulink on a PC based cluster.Computational time is in the order of a few days.The second optimiser uses a simple genetic algorithm from Ref. [26].The population is 100 and the number of generations is 30.The computational time required for the second optimiser is on the order of a few seconds, which is almost negligible compared to around an hour for the WEC-Sim simulations.

Computational model set-up 2.2.1. Device design numerical inputs
The design specification of the device considers three hydrodynamic entities: a spherical float with a mass of 80% of the displaced mass; a reactor consisting of a tubular steel frame connected to the float's PTO lines; several air-filled buoyancy tanks connected to the reactor.Two mooring cables at each corner of the reactor (one vertical and the other angled) give stability to the reactor.The numerical optimisation process allows the design team to investigate different multi-float configurations at very low cost.However, to increase confidence in the numerical results, validation through tank testing is required.The most important design parameters are listed in Table 1.
The design of the reactor mass and the net buoyancy and the dimensions of the buoyancy tanks is dependent on the float spacing and number and based on some preliminary relationships provided by MPS.Some details can be found in Appendix A. The schematic of Fig. 2. Optimisation scheme for each multi-float configuration and each float spacing.

Table 1
The design parameters for the float and the reactor.

Float
Number  The mooring layout consists of a number of oblique and vertical lines to increase the reliability of the mooring system.Indeed, the redundancy of the mooring lines helps to reduce the consequences of the failure of a mooring component.In the linear, rectangular and triangular configurations, each corner of the reactor is connected to both an oblique and a vertical mooring line while in the circular configuration there are four oblique and four vertical lines.More specifically, the mooring attachments on the reactor for the circular configuration are located in an imaginary inscribed square to the external circle.The oblique and vertical lines share the same vertical preload tension.The mooring layout for each configuration in the xy plane is shown in Fig. 7.For more details on the configuration of the mooring line, see Appendix A.

Simulation parameters
The parameters for the numerical simulations using Nemoh and WEC-sim are shown in Table 2. Frequency range and number are based on previous experience of frequency-domain simulations using Nemoh, while the simulation time, ramp time, solver and timestep are reasonable choices for time-domain simulations.The simulation time without including the ramp time is 30 min, the typical recording length suggested for irregular waves [27].The ramp time is used to gradually increase the flow until the full-wave conditions are reached [23].

Sea states simulation
Belmullet, the Atlantic Marine Energy Test Site (AMETS), Ireland [28], was selected as the location for this optimisation.This area has a water depth of 100 m, which corresponds to the design depth of the WaveSub device.A wave data analysis for the years 2013e2017 from the Berth A wave buoy in Belmullet is carried out to obtain information on wave spectra and occurrence matrix, as shown in Fig. 8.A similar methodology is also described in (E [20].The full scatter matrix considers significant heights between 0.25m and 15.5m and a spacing of 0.25m, while energy periods between 4s and 22s and a spacing of 0.5s are considered.However, in order to achieve an acceptable runtime for the optimisation, three main wave spectra in three main ranges of energy periods were chosen.These correspond approximately to energy periods below 9s, between 9s and 11s and above 11s.The central range between 9s and 11s represents the energy periods that are the most energetic.Therefore, the optimised device will be mostly tuned to produce the maximum amount of energy in this region.It is important to note that the most energetic periods have been chosen for the design and not the most frequently occurring ones, as this paper is primarily concerned with maximising the energy produced by the device.The spectra simulated in the optimisation are shown in Table 3.The spectra used in the simulations are the average raw spectra while an approximated theoretical spectrum is given as a comparison in Table 3.The significant heights of the spectra have been normalized to 1 m for the simulations in WEC-Sim, while the significant heights in the post-processing step are the original ones of the full occurrence matrix.Previous simulation experience with the WaveSub device has shown that all forces considered in the numerical model of the dynamic system are linear, except for minor drag terms, and the power for moderate waves depends almost linearly on the square of the significant wave height.Therefore, a full numerical simulation was performed only for the different energy periods and for a single significant wave height.This is followed by a post-processing step to compute the power and energy matrices.The power variation with Hs has been obtained in the post-processing as where P and H are the power and significant heights, respectively, while the indices s and p represent the simulation and postprocessing step, respectively.The power is also capped at the rated power, while for wave heights larger than 10 m a survival condition without generation is assumed.The energy matrix is then the product of the power and the occurrence matrices.Finally, regarding the wave direction, two different types have been used.These are a single wave direction and a uniform normal wave direction distribution with the numerical settings already presented in Ref. [4].The distribution is based on the average standard deviation of the wave directional distribution (28.95 ) in Belmullet.

LCOE calculation
LCOE is an important indicator of technology performance.The accuracy of the calculation of energy production and the wave farm costs is of primary importance to the industrial partner in this work.The calculation of the total cost used in the optimisation is based on confidential information provided by MPS.The cost calculation takes into account a WEC farm with a fixed size of 500 MW.For this calculation, the experience from the wind farm sector is used [29] as well as a detailed cost calculation for some device components.A simplified cost parameterisation was employed to allow optimisation without requiring a detailed design for each variant.Only the most significant cost components are considered, which leads to uncertainties in the accuracy of the cost estimate, but provides a robust and efficient metric for comparing the different configurations.
The cost calculations in this paper refer to a WEC array as a realistic choice for industrial applications.The CAPEX is considered as the sum of the cost of the structure, PTO, mooring, Balance of System (BoS) and installation costs.All these costs are proportional to the number of devices.The specific cost per device can be described as follows: The structure cost (C st ) is calculated as where M R and M F are the mass of the reactor and the float, c R and c F are the specific material costs for the reactor and the float.
The calculation of the hydraulic PTO cost considers specific costs of the components of the primary transmission.Then the sum of the costs for the secondary transmission and energy storage is proportional to the rated power (P rated ) using a constant specific value (c PTO ).The hydraulic PTO cost (C PTO ) can be expressed as where C PTOp is the cost of the primary transmission system of the PTO.
The mooring cost (C mo ) is based on the cost of the mooring lines, anchors (C anchors ) and winches (C winch ).It is expressed as where N lines is the number of mooring lines, L line and c l are the length and specific cost per metre of the mooring lines.
The Balance of System cost (C BoS ) is proportional to the rated power of the device: where c BoS is the specific cost related to the BoS.
The installation cost per device is fixed.
Finally, the OPEX cost is a fraction of the CAPEX while the decommissioning cost (DECOM) per device is a constant value.More details of the cost calculation are confidential.
The LCOE is expressed as where PVF OP , PVF D and PFV E are the present value factors related to the OPEX, DECOM and annual energy production from the farm (EP year ) for a project lifetime of 25 years.EP year is calculated considering an availability factor, conversion and transmission efficiencies.

Mesh convergence
A mesh convergence study is undertaken and this is described in Appendix B.

Drag term
The hydrodynamic coefficients from the BEM calculation in the frequency-domain in Nemoh do not take viscous drag into account.However, the time-domain simulation in WEC-Sim can include quadratic drag if a drag coefficient (C d Þ is defined by the user.It is expressed as where A d is the characteristic area for the drag calculation.Drag gives a more realistic representation of fluid viscosity, but the accuracy of this effect depends on a good choice of the drag coefficient.A drag coefficient of 0.4 has been used for the floats considering the literature.[30] shows very low drag coefficients for a sphere for a low KC number, while in Ref. [31] a value of 0.4 is suggested for a laminar steady flow.This last value has been considered reasonable compared to Ref. [30] and was therefore chosen.The cylindrical buoyancy tanks that are part of the reactor have a low velocity and consequently a low KC number.[32] shows an increase in the drag coefficient for a lower KC number for smooth and rough cylinders and so a value of 2 has been chosen.Rotational drag is included for the tanks of the reactor while it is assumed to be negligible for the floats (spheres).The characteristic areas of the drag rotational moments is calculated considering a  simplified rectangular cross-sectional area of the tanks and calculating the integral as follows for each tank: where l min and l max are the minimum and maximum distances between the rectangle cross-sectional area of each tank and the rotational axis of the reactor, L TankÀts is the transversal size of the cross-sectional area of the tanks, v and u are the linear and angular velocities about the axis of rotation.

The linear 6 float design
The results of the optimisation for the linear configuration with the smallest float spacing is shown in Fig. 9.The red data points represent the WEC-Sim simulations of the model using the genetic algorithm used to create the surrogate model shown as a surface plot.Each subplot represents a different parameter and each generation of the genetic algorithm converges the simulations towards the optimal minimum LCOE.The results shown are the LCOE, the average total power, the average power for each of the three simulated sea states, the rated power and the capacity factor.To maintain confidentiality, each parameter has been normalized to the maximum value obtained during optimisation (R MAX ).The relationship between the simulation results (R s ) and the normalized results (R n ) is expressed as The average total power is calculated based on the rated power, the occurrence and the power matrix.The average power for each sea state is calculated from each wave spectra described in Section 2.2.3.Results show that there is generally a single optimum.The LCOE plot (Fig. 9a) shows a very good convergence to the minimum found.The optimisation of the stiffness and damping of the PTO has led to a reduction of the LCOE by about 50% compared to the maximum value found.Fig. 9b has shown that the maximum value of the average power has a very similar optimum to the LCOE, as the PTO stiffness and damping have a negligible influence on the cost of the device.The optimum LCOE is influenced by the rated power (Fig. 9f) and the capacity factor (Fig. 9g) and requires high values for both.A high rated power reduces the structural cost due to a smaller number of devices in the WEC farm, while a high capacity factor increases the average power produced for that specific rated power.The power produced by the sea states shows that the optimum moves towards low stiffness for larger wave periods.The optimal design parameters are similar to the values of the second sea state.
All optimal results are summarized in Table 4.Each column in Table 4 refers to single direction (SD) or multidirectional (MD) waves in combination with three float spacings (FS1, FS2, FS3).The optimal values of rated power, energy production, capacity factor, total CAPEX and LCOE have been normalized to the maximum value obtained from the three float spacings and the two types of directional cases.Optimal damping tends to decrease with an

Table 3
The spectra used in the simulations and the sea states used in the post-processing step for the estimation of the energy matrix.increase in float spacing for both the single direction and multidirectional cases.The optimal stiffness decreases slightly with float spacing in the single direction case, while it is less predictable in the multidirectional case.The float spacing does not have a clear and significant influence on the energy produced by the device.In particular, the single wave direction and the multidirectional waves have an opposite effect on the total energy.The capacity factor generally increases slightly when the float spacing is increased.CAPEX tends to increase with a larger float spacing for both the wave cases, and since CAPEX structural cost is a major cost driver, consequently LCOE also increases with the float spacing.OPEX and decommissioning costs are negligible compared to CAPEX.They each account for about 3% of CAPEX.In fact, the OPEX costs for the WaveSub device are estimated to be relatively low.The NRMSE and the MAX error are also shown in Table 4 and are calculated considering the predictability of the surrogate model built on the first 4 generations for the last generation.The NRMSE of the surrogate model shows reasonable predictive capabilities for most of the cases as defined in Ref. [33].

The rectangular 6 float design
Table 5 shows the main results for the rectangular configuration.The optimal stiffness and damping generally decrease with the float spacing.However, the damping coefficient increases for the largest float spacing and the single direction case.The energy produced for the single direction case decreases with the float spacing while it is fluctuating for the multidirectional case.In the single direction case, the capacity factor increases with float spacing, while an opposite trend is observed in the multi-directional case.CAPEX and LCOE generally increase with float spacing, but the second float spacing for the multidirectional case shows a slight reduction in both values compared to the lowest float spacing case.Structure cost is also the main cost component for the rectangular configuration.The surrogate model shows good predictive performance based on the NRMSE and MAX error, demonstrating good optimisation convergence.

The triangular 6 float design
Table 6 shows the main results for the triangular configuration.The optimal stiffness and damping generally decrease with the float spacing.However, the damping coefficient increases for the largest float spacing and the multidirectional case.The energy produced increases with float spacing, while the capacity factor decreases for both the wave configurations.CAPEX and LCOE both increase with float spacing.The main difference in the cost components compared to the other configurations is the mooring cost.Since this configuration has three corners instead of four the mooring cost is reduced.The optimisations also converge quite well for this configuration, as shown by the prediction performance of the surrogate model using the NRMSE and the MAX error.

The circular 6 float design
Table 7 shows the main results for the circular configuration.The optimal stiffness and damping show a decrease from float spacing 1 to float spacing 2. Instead, there is generally a slight increase in the design parameters between float spacing 2 and 3.
The energy produced decreases with float spacing, while the capacity factor increases in the single direction case, while it oscillates in the multidirectional case.CAPEX and LCOE both increase with float spacing.The cost components ratio is quite similar to the linear and rectangular configurations.The NRMSE and the MAX error also show good optimisation convergence for this case.

Visualization and power check of the optimal case
In this section we show a visualization and power check for each configuration and lowest float spacing of the optimal PTO springdamping for a selected representative regular wave.The optimal PTO spring-damping are described in the previous Tables 4e7, while the sea state is a regular wave with a wave height of 3 m and a wave period of 10 s and a single wave direction.Figs.10e13 and Tables 8e12 show the total wave amplitude generated by the perturbed wave field and the power produced by each float, respectively.Table 8 shows that the power generated by float1 to float6 generally decreases for the linear configuration due to the shadowing effects of the previous floats facing the incoming wave.In fact, if interaction effects between the floats are not considered in a first approximation, each float generates a parabolic wave amplitude field [6], giving the following floats a smaller total wave amplitude from which to extract energy.A larger power produced by a float is also associated with a larger reduction in the total amplitude behind the float, as it is possible to notice for example when comparing float 1 and float 5 in Fig. 10.
Figs. 10e13 shows that the first three linear, rectangular and triangular configurations have a main parabolic shape of the wave field centred on the front float, while the circular configuration has several local shapes for each float.The incident wave amplitude and the values below 1.2 m have been marked in the figures to make the parabolic shapes more visible.We can see that the shadowing effect is more pronounced for the linear, rectangular and triangular configurations, while for the circular configuration it is more limited to the rearmost float.This means that the front float in the first three configurations is likely to have more power than the rear ones.In addition, some of the floats in the rectangular, triangular and circular configurations have some symmetry with respect to the x-axis.As a result, the overall amplitude pattern and the power generated depend only on the x-position.Some of the floats of these configurations have the same x-coordinate, so the power produced is the same.
Table 9 shows similar power produced between floats 1 and 6, floats 2 and 5, and floats 3 and 4 in the rectangular configuration.Floats 1 and 6 are the first floats to face the incoming wave and they produce the largest power.The hydrodynamic interaction between  the floats leads to an increase in power in the back row (floats 3 and 4).Table 10 shows that similar power is generated between floats 2 and 6 and between floats 3 and 5 in the triangular configuration.Float 1 produces about twice the power compared to the other floats and is the main influence of the wave field on the total amplitude.The floats in the back of float 1 receive a lower total wave amplitude due to the energy extraction by float 1 and therefore produce less power (shadowing effect).
Finally, Table 11 shows similar performance between floats 2 and 6 and between floats 3 and 5 for the circular configuration.Floats 2 and 6 produce the largest power, because of the positive interactions with the nearer floats.Floats 3, 4 and 5 are mostly shadowed by floats 1, 2 and 6 and so their power is reduced.

Comparing configurations
The overall results for all geometries can now be compared.Fig. 14 shows the LCOE against the rated power for the different multi-float configurations obtained by the single wave direction and multi-directional wave cases.The circular and rectangular configurations have the highest rated power compared to the other ones.The lowest LCOE is generally given by the lowest float spacing for each configuration.The linear configuration is outperformed by the other configurations.The lowest float spacing of the circular configuration has the minimum LCOE (about 12% LCOE reduction compared to the worst configuration), followed by the lowest float spacing of the triangular configuration.The MAX error is generally very low.The float spacing 3 of the rectangular configuration for the single wave direction case has the largest error, but is still adequate for comparison.Fig. 14b shows that the inclusion of multi-direction waves does not change the overall conclusion.

Discussion
The genetic algorithm, supported by the Kriging surrogate model with a high elitism factor, has found rapid convergence to the optimum.In fact, the optimisation function has been shown to be convex with a single optimum.The surrogate model has played a very important role in reducing the number of simulations to reach the optimum.In this optimisation, a high elitism factor was used, while a lower factor is required when the optimisation function has multiple optima.
Optimal PTO stiffness and damping were determined for each configuration in terms of LCOE.However, the PTO parameters have no direct influence on the cost calculation, so the second stage optimisation could also maximise average power.This choice has been considered insignificant as the simulation time required for the postprocessing block is negligible compared to the simulation in WEC-Sim.
LCOE has been shown to depend on both rated power and capacity factor.High values for both factors are required for optimal LCOE.A fixed total farm rated power has been used, so a high rated power reduces the structure cost due to the lower number of devices required.On the other hand, the optimal capacity factor cannot be optimised directly.Instead it is determined by the optimal average power and rated power.
LCOE is also strongly influenced by CAPEX and especially by structure cost.A configuration with a larger float spacing is penalised by higher structural costs, while the total energy from the WEC farm is similar for all the float spacing cases.The influence of the float spacing on the total energy is not very clear, but in general relatively small.
Optimal PTO stiffness and damping generally decreases with increasing float spacing, especially for the linear, rectangular and triangular configurations.This optimal value could be related to both the float spacing and the sea states considered in the optimisations.Finally, it is worth noting that the stiffness and damping were optimised to a single value for all 3 sea states based on the total annual energy of the energy matrix.However, a more advanced control would probably have an optimal spring-damping for each sea state with a possible further reduction of the LCOE for each device.
This work has some limitations, such as the use of a time-domain model instead of a frequency-domain model and the linearisation of the PTO force.Usually, the optimisation of a WEC farm is performed with a frequency-domain solver, which significantly reduces the required computational time (see for example [6,7,13,18,34]).The hydrodynamic quadratic drag can also be linearised and included in a frequency-domain model [35e37].Excellent agreement between the quadratic drag and the statistical linearised version is shown in Ref. [37], while an example of application to a self-reacting two body device is given in Ref. [36].However, when non-linear effects such as PTO and mooring are considered, a time-domain model is usually used because it is more accurate than a frequency model.For example, the differences between a linear and a nonlinear PTO were compared in Ref. [38], while the mooring effect was investigated by comparing time and frequency domain approaches in Ref. [39].The time domain has also been preferred for optimising the geometry of a WEC in several works such as in Refs.[40e42].In this work, a timedomain simulation was preferred as it allows the inclusion of nonlinear effects, such as the nonlinear quadratic drag and the impulse response convolution term of the hydrodynamic radiation force.In addition, the number of floats is limited to 6, which allows the use of a time-domain approach, whereas a frequency approach is more likely to be required for a larger number of floats.Moreover, most geometries chosen for the frequency domain optimisation of a WEC array consider a single body or two bodies as point absorbers, while a multi-body system with more than two bodies and 6 of freedom has not been considered yet [43] and would probably require the verification of the numerical model.On the other hand, WEC-Sim has been verified and validated in a large number of works [44,45] and is a powerful tool for the study of multibody systems due to the use of the Simscape library [46].WEC-Sim is also coupled with Mechanics Explorer [47] that enables real time debugging of the simulation to identify errors in the numerical set-up and numerical instabilities.Therefore, the methodology presented in this study was considered the most appropriate to increase confidence in the results, which are very important for the industrial partner of this work.
The implementation of other non-linear effects, including nonlinear hydrostatics and Froude-Krylov, a dynamic mooring model and a realistic hydraulic PTO are already implemented in the WEC-Sim code, making such a model particularly useful for upgrades compared to a frequency-domain model, and can be used in a more detailed design phase.The hydraulic PTO has also been simplified as a linear PTO and this is probably the most important assumption of the model.However, we must bear in mind that this model refers to a concept design phase where high accuracy is not yet so important.As the machine develops further, the hydraulic PTO needs to be designed in more detail, e.g. the number of cylinders, which could greatly influence the results.A confidential exchange with MPS has also shown that a control system that can use a sufficient number of hydraulic cylinders could replicate the linear stiffness-damping behaviour fairly accurately.However, these nonlinearities are not considered in this work, as it affects the computational time significantly (more than 5e10 times).

Conclusions
Different layouts of multi-float configurations of the WaveSub device have been investigated in this work.The baseline linear configuration is compared with the rectangular, triangular and circular configurations.For each configuration, a medium fidelity time-domain model has been created considering hydrodynamic coefficients from linear potential flow theory and quadratic drag.The model is based on open-source software (Salome-Meca, Nemoh and WEC-Sim) that allows complex multibody systems to be efficiently created and hydrodynamic interactions to be taken into account.Salome-Meca and Nemoh are used to create the mesh and obtain the hydrodynamic coefficients needed for timedomain simulation in WEC-Sim.Nemoh is also used to check the convergence of the hydrodynamic coefficients for different sizes of the mesh panels.In particular, this work investigates multibody systems consisting of six floats connected to a submerged reactor via a series of PTO lines.The collaboration with MPS has allowed the conceptual design of the different layouts to be developed that take into account their industrial experience, such as the use of buoyancy tanks and the estimation of the various costs associated with the device.
A linear representation of the dynamic system, apart from a few terms such as the quadratic drag and the memory convolution term of the radiation force, was preferred in order to reduce the computational time required for the optimisation.The optimisation was carried out considering an efficient optimisation algorithm combining a genetic algorithm and a Kriging surrogate model, which is particularly suitable for computational expensive simulations.For the time-domain simulations, three sea states representative of the Belmullet site were selected, while in a postprocessing step an extended occurrence matrix is used to estimate the energy produced under the assumption of linearity.The optimisation was performed both for the case of a single wave direction and for a more realistic case with multi-directional waves.In particular, the comparison between the configurations leads to a similar conclusion between the two different directional cases.The results show that float layout can be optimised to improve the power generated by a multi-float configuration.However, a possible design of a multi-float configuration must also take into account the associated costs.Therefore, the results here are optimised for LCOE rather than energy produced.Generally, the lowest float spacing of each configuration has demonstrated the lowest LCOE.It is clear that the most important result of this study is the demonstration that the circular configuration is most optimal.The better performance of the circular configuration is not completely unexpected.In this configuration, the rear floats are generally less shadowed while the two floats in the second row might have some energy benefit from the parabolic effect of the front float.Furthermore, other studies have shown that the linear configuration can be outperformed by other configuration types [18].
The numerical model could be improved in future studies in terms of drag, the PTO and mooring.The estimated drag coefficients with constant value are sufficient for small waves but could be improved for large float orbits.In fact [20], shows that the drag coefficient depending on the Keulegan-Carpenter number is generally similar for a small wave.However, the nonlinear effects are expected to increase for a larger wave, so a constant drag coefficient is likely to be more uncertain.A linear spring-damper is used for the PTO, which approximates the PTO control force.However, more complex effects such as the end-stop force are not taken into account.The mooring is also simplified as a spring with high stiffness, as it has no particular influence on the dynamics of the device.The accuracy of this modelling should be adequate for the taut mooring expected for this type of device.In a previous study, the stiffness of the mooring was investigated and found to be fairly linear [20].The main issue of this simplified modelling of the mooring is that does not take into account some possible but quite rare slackline events.
Then the optimisation can be improved to account for the tuning of the PTO parameters for each sea state.This will further reduce the LCOE of the configurations.For a more detailed investigation of the Belmullet site, an extended occurrence matrix needs to be included.Non-linear effects that depend on Hs could also be considered in a further study.However, it is expected that a CFD model would be used to simulate especially extreme events.
Finally, the cost model is limited to industrial experience gained from prototype deployments and will become more robust as more devices are built.However, the comparison between configurations takes into account all significant cost effects, including structure and PTO, these are modular and there is reasonable confidence in the estimates.Therefore, the comparison between configurations are valid and can be used to inform design and commercial decisions.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The mooring lines of the linear, rectangular and circular configurations generate an angle of 45 in the orthogonal views in the xy, xz and yz planes.The triangular concept with multiple floats has three mooring lines that create an angle of 45 in the planes described by the vertex of each triangle and the bisector of the corner angle.The reactor and seabed mooring attachments (R and S index of the next equation, respectively) for the two corners of the backside of the device can be defined as follows: The attachment points relative to the corner at the front of the device can be defined as follows: x M À x S ¼ z M À z S y M ¼ y S ¼ 0 (24)

Float
spacing depends on a multiple of the float radius (R) and three discrete values are used: FS1 ¼ 7R, FS2 ¼ 8R, FS3 ¼ 10R.PTO linear spring (kPTO): [5e1000] kN/m PTO linear damping (cPTO): [50e1000] kNs/m Rated power: [0 100] MW the four different multi-float configurations is shown in Figs.3e6.The linear configuration has six floats in a line perpendicular to the wavefront.The rectangle has a 2x3 layout with two floats facing the waves.The three rows of the triangle consist of 1, 2 and 3 floats respectively, with the front vertex facing the waves.Finally, the circular configuration consists of six floats evenly distributed around a diameter.All of these configurations have the same configuration of PTO lines where the attachment points to the reactor are inscribed a square, as shown in Figs.3e6.The PTO lines of the linear, rectangular and triangular configurations generate an angle of 45 in the orthogonal views in the xy, xz, yz planes.In the circular multi-float concept, the attachment points of the PTO lines to the reactor are located along two circles, as shown in Fig.6.

Fig. 3 .
Fig. 3. Schematic of the linear configuration for the 3 float spacing values.

Fig. 4 .
Fig. 4. Schematic of the rectangular configuration for the 3 float spacing values.

Fig. 5 .
Fig. 5. Schematic of the triangular configuration for the 3 float spacing values.

Fig. 6 .
Fig. 6.Schematic of the circular configuration for the 3 float spacing values.

Fig. 8 .
Fig. 8. Diagram showing the process used to obtain the power matrix of the device.

Fig. 9 .
Fig. 9. Kriging surrogate model of the main results of the linear multi-float configuration for the lowest float spacing considered and a single wave direction.Red data points represent the simulation of the numerical model.

Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .
Fig. 10.Total wave amplitude generated by the lowest float spacing of the linear configuration for a regular wave with H ¼ 3m and T ¼ 10s.

Fig. 14 .
Fig. 14.Comparison of the LCOE/power results between the different multi-float configurations and between the single direction (SD) and the multi-directional (MD) cases.The error bars refer to the MAX error (lower and upper bounds) of the surrogate model on the last generation.Results are normalized to the maximum LCOE and rated power of the single direction case.

Figure B. 1
Figure B.1 Mesh hydrodynamic convergence relative to the linear 6 float configuration.The float refers to the float1 of the schematic of the linear multi-float design shown in Fig. 3.

Figure B. 2
Figure B.2 Mesh hydrodynamic convergence relative to the rectangular 6 float configuration.The float refers to the float1 of the schematic of the rectangular multi-float design shown in Fig. 4.

Figure B. 3
Figure B.3 Mesh hydrodynamic convergence relative to the triangular 6 float configuration.The float refers to the float1 of the schematic of the triangular multi-float design shown in Fig. 5.

Figure B. 4
Figure B.4 Mesh hydrodynamic convergence relative to the circular 6 float configuration.The float refers to the float1 of the schematic of the circular multi-float design shown in Fig. 6.

Table 2
Main simulation parameters for Nemoh and WEC-Sim.

Table 4
The optimal PTO parameters and main results of a linear multi-float configuration.

Table 5
The optimal PTO parameters and main results of a rectangular multi-float configuration.

Table 6
The optimal PTO parameters and main results of a triangular multi-float configuration.

Table 7
The optimal PTO parameters and main results of a circular multi-float configuration.

Table 8
Normalized power produced by the lowest float spacing of the linear configuration for a regular wave with H ¼ 3m and T ¼ 10s.

Table 9
Normalized power produced by the lowest float spacing of the rectangular configuration for a regular wave with H ¼ 3m and T ¼ 10s.

Table 10
Normalized power produced by the lowest float spacing of the triangular configuration for a regular wave with H ¼ 3m and T ¼ 10s.Floats are numbered clockwise from Float1 facing the incoming waves.