MODELING LINK FOR OCEAN WAVE ENERGY OF POINT ABSORBERS : A MATHEMATICAL FRAMEWORK

This compendium presents new mathematical techniques for modeling Point Absorbers. A combined frequency-time domain framework is developed. It is used to simulate the energy generated by the wave farms. With Matlab and Fortran as a base, this leads to obtain physical variables of primary importance, namely position, velocity and power to energy net balance relationships of absorption. Integration of different degrees of freedom with heave as main executable leads in turn to a single buoy motion focus. Acquisition of the needed hydrodynamic coefficients is provided through application of potential field solvers with Boundary Element Methodology background. Initially, this Wave-to-motion model is validated by comparison with previous experimental results for a floating cone cylinder shape (Buldra-FO3). A single, generic, vertical floating cylinder is contemplated then, that responds to the action of the passing regular waves excitation. Later, two equally sized vertical floating cylinders aligned with the incident wave direction are modeled for a variable distance between the bodies. For both unidirectional regular and irregular waves as an input in deep water, we approximate the convolutive radiation force function term through the Prony method. By changing the spatial disposition of the axisymmetric buoys, using for instance triangular or rectangular shaped arrays of three and four bodies respectively, the study delves into motion characteristics for regular waves. The results highlight efficient layouts for maximizing the energy production whilst providing important insights into their performance, revealing displacement amplificationand capture width-ratios, while deriving in possible interpretations of scenarios related to the known park effect. These terms are encompassed by the novelty of a new conceptual Post-Processing methodology in the field, which leads to obtain an optimal distance for the separated bodies with effective energy absorption in a regular wave regime. The main objective is to generate a tendency within the hydrodynamic field of study, which is the Wave to motion perspective. More generally, this computational excursion envisions and depicts potential fields of study, which will surely enhance new connections and link this renewable energy form. Therefore, this research delves first into the historical and technical background on Ocean Wave Energy. Next, it is in the section regarding Materials and Methods, where boundaries and related equations are introduced step by step, together with latter mentioned case scenarios, and their corresponding configuration parameters. A separate section frames then the scope of results, while finally, there is an ensuing discussion and conclusions for evaluation assessment.


Introduction
Among the large list of patented Ocean Wave Energy Converter (WEC) types, the Point Absorber (PA) technology has become predominant in the field according to [1].Studied intensively since the oil crisis during the late 1970's by Budal and Falnes [2] in Scandinavia and by Evans [3] in parallel in the United Kingdom, focus was set upon the Mathematical & Physical Modeling of such devices among others.They are best described by the fact that their characteristic length is much smaller than the surrounding wavelength from the passing waves.Nowadays, several countries all over the world are investigating the performance of this technology [4].PA's can be floating or bottom-mounted structures, with a Power-Take-Off (PTO) mechanism.They vary depending on the oscillation mode of the geometry.Most of the devices use mainly the vertical displacement (heave) of the structure induced by the incident waves to activate the electricity generation either through hydraulics or directly using magneto-inductive components.By positioning and connecting several units close to each other, the system is addressed as a multipoint absorber [5].Another possibility is to deploy lines or "matrices" of single PA's in order to cover a larger area, while absorbing more energy.These so-called Arrays might become the future of Ocean Wave Energy, since they not only produce more energy in theory, but they can also be adapted to existing floating structures, such as Floating Production Storage and Offloading vessels [7] or offshore wind energy foundations.Similarly to the way pump storage plants alternate energy production between turbomachinery and wind turbines [8], WEC's might supply Wind Energy during its production flaws, since it is a much more constant resource [9].The State of the Art of the modeling procedure for this work is schematized next: In Table 1 it is noted, that the step-by-step table has been determined with specific focus on previous project experience in the field [10].Procedure steps A-E involve a pre-/ solver / post-processing structure, that ranges from 1 − 4 point absorbers.It also includes certain optimization perspectives at the end stage E as a feedback loop, which might in turn lead to restart at A again.
A deeper look in Europe at some implemented projects (incl.real sea testing) reflects the following list Wedge Global buoy (Gran Canaria-Spain) Pipo Systems buoy (Gran Canaria-Spain) OPT buoy array (Spain) Searev (France) Tecnalia-BIMEP point absorber (Basque country-Spain) Corpower buoy (Sweden) Seabased buoy (Sweden) ISWEC (Italy) AWEC (USA) At multipoint absorber level, two main projects are highlighted, namely: Wave Star Energy (Denmark-Europe) Albatern WaveNET (Scotland-Europe) Emphasis is set upon the insular base, which surrounds the Canary Islands, due to • location-based proximity: either european or outer borders • functionality: analogy to present study conditions

• proven concept: based on prototype technical reports
The first seven projects correspond to single Point Absorbers, while the two remaining products represent Multipoint Absorbers.From the total number of introduced projects, only two are being tested right now, namely the Wave Star and WaveNET.For a broader perspective on latter mechanisms, further global documentation can be found in [1].
However, certain drawbacks have to be considered before deploying such converters in the sea.The first issue corresponds to the high corrosive, harsh environment, which is the ocean.Then, there is the need for large infrastructures, such as mooring lines, both to anchor the floater to the bottom and to transmit the electricity to shore if the PTO is in the PA itself.Also, there is always a risk of encountering extreme wave conditions, which may affect the device during operation.All of these arguments determine the importance of minimizing costs during design and test phases, while considering possible worst case scenarios.Computational methods employed during the course of this research are a consistent approach to achieve the mentioned purposes, specially with respect to hydrodynamics.
The number of different applied numerical methods in Ocean Wave Energy is increasing steadily [11].Since it is also noted, that computational capacities are growing exponentially [12], it is a matter of time before the scientists will be able to describe real fluid phenomena with certain accuracy.For the given cases, a combined BEM-frequency-time domain analysis will take place.For a so called wave to wire (W2W) model [13], this means, that first, hydrodynamics coefficients are determined for the wetted geometry.Next, these parameters are attached to the equations of motion in the model solver, which in turn generates amplitudes of oscillation, velocities and power (energy) absorbed by the PA's in either frequency domain or time domain.In short, given an input wave incidence, the codes compute pressure and velocities on the panels (FEM), and determine characteristic hydrodynamic parameters.Latter ones can be employed to solve the governing equations of motion, which in turn, depending on the chosen domain (either time or frequency), will provide insights on the motion displacement patterns.
As in any other related field, optimization is the keyword to refine the overall performance of the technology.From the hydrodynamic point of view, there are mainly three arguments that have to be taken into account in order to obtain maximum energy.The first one refers to the oscillation mode of the device.[14] determined that combining heave and surge, the device will absorb more than 80% of the available Wave Energy.Hence, one has to find suitable geometries and PTO systems that can superpose the different degrees of freedom.The geometry itself turns out to be determinant, and thus becomes the second argument.The Salter's duck is the only patented device mechanism, that is theoretically nearly capable of absorbing all the incident wave potential [15,16].Last but not least, near-resonance conditions have to be determined for achieving a narrow operational bandwidth of the PA [17], from which to obtain its maximum amplitude motion.
For the current numerical study, two geometries will be under investigation.First, a single body with a cone cylinder shape employed by the SEEWEC project [18] will validate the present model approach against experimental results.Then, the study will base on the generic vertical floating cylinder case, such as the one proposed by Seabased [19,20].With regard to latter PA type, different array layouts will be modeled both in the time and in the frequency domain, conducting the researchers to deal with basic hydrodynamic interactions of such heaving devices for various distances and dispositions.Thereby, a restriction on regular incident waves will give indications on the extent the devices might operate with greater power absorption capability.If not, then causes for the losses will be accounted for through existing phenomena such as shadowing or near-trapping effects [21][22][23].It is noted however, that the added value of the present work relies on ensuring a valid and novel approach in the Post-Processing assessment of the PA motion characteristics.

Materials and Methods
Throughout the course of this research, certain conditions have been met for the incoming waves and more specifically for the floating systems.This involves choosing the correct physical and mathematical background for the underlying theory behind.Both subjects are treated in the following section.

Assumptions
According to the wave mechanics theory, following characteristics are being considered: • Linear Potential Airy theory • Various Depths (ranges from shallow water to deep water) • Unidirectional regular and/or irregular waves (specified on own Phd thesis) Regarding the floating system, further limitations are introduced next: • Mooring, viscous and damping forces have been determined empirically, using previous thesis results from Bosma [20].• The PTO force is estimated randomly for the generic heaving cylinder, since no such an existing physical device has been found with published damping characteristics.
It is noted, that basically, unidirectional regular waves are used as an input to the model, as stated in 1 .Irregular superpositions are only mentioned due to the own PhD thesis carried out by the author.Linear potential or so called Airy theory applies to the previous mentioned sea states, mainly based on small amplitude wave heights H with ratios such that H D h << 1 corresponding to shallow water and H L << 1 referring to deep water.Here, D h represents the water depth, and L is the wave length.For the boundaries, classical periodicity and Sommerfeld radiation appear at inlet and outlet respectively.
Thereby, the open-source Diffraction-Radiation solver NEMOH (BEM) is used for the acquisition of the hydrodynamic coefficients in the frequency domain.These are added mass (A r ), radiation damping (C h ) and excitation force (F e ) as sum of diffraction (F d ) and Froude-Krylov forces (F f k ).Additionally, one obtains in the time domain the values for A r at ω → ∞ and Impulse Response Functions for the radiation force function k r (t)(retardation kernel).All these terms become part in the solution of Newton's equations of motion, whereas specific terms have to be first integrated in convolutive form, for instance the radiation force F r .

Hydrodynamic Equations of Motion
The Fluid-Structure Interaction is complex, and can be abstracted even more, as Figure 1 shows: Taken the PA system as a whole software application itself, one does not only provide fundamental simulation theory and the main focus of this work (navy blue rounded shape (0)) related to the present investigation, but it also envisions future perspectives (yellow dashed ellipsoids), while foreseeing new research fields (grey nodal shapes).It is merely an abstraction of a black-box system, which includes systematically further optimization strategies (confer to violet variable NACA profiles on the bottom corners of (0) in Figure 1.Assuming (0) as the starting unit confronting the waves, the scheme is extended hereby into a small triangular array (objects I and II), which is expected to perform efficiently.Moreover, from the base to the top of the sketch, 1 reflects the hydrodynamic entry with wave amplitude ζ (in turquoise), including the most relevant parametric units.In the middle part, the structure presents the classical vibration and computational approach.Reaching the top, the spreading ellipsoidal dashed curves represent characteristic force vectors, which might be tuned in advance for whatever focus research tendencies might converge into.
From the main "EXE" workflow, one can first deduce, how the hydrostatic and hydrodynamic parts are introduced by the BEM calculation.In the middle part, a modified Mass-Spring-Damper (MSD) equation is solved, which is a second order ordinary differential equation in the form M t z = ∑ F , whereas M t involves the total mass of the floating body.Additionally, simple (double) dotted variables denote first (second) order time derivatives of position, velocity and position respectively.In this case, the z formulation refers to the heave component in vertical cartesian coordinates.Moreover, this procedure assumes that the surrounding water medium can be replaced by multiple spring and stiffness coefficients, which are then inserted into the discretized equation.Furthermore, the sum of forces becomes ∑ F = F e + F r + F h + F m + F pto , from which one can assume the single components information.First item on the RHS is the excitation force, which can be formulated through the following convolution integral where τ represents the causality time span, such that some previous stage is taken out of memory in advance.Second term is the radiative term, which is described generally by the equation with ż as the vertically displacing velocity component of the radiating body in motion.
The third term inserts the hydrostatic restoring vector into the sum of forces, yielding with K h as the hydrodynamic stiffness.Fourth term considers the application of a mooring force on the floating system, such that where K m stands for the mooring coefficient.Last term includes the PTO-force, which is simply reflected as a damping mechanism in the form whereas C PTO represents the PTO damping coefficient.There is the spectral distribution S(ω), which serves as the last input to the model.One can calculate out the free surface elevation η, which is then convoluted with the diffraction force F d to determine the heave excitation force.It is noted at this stage, that a simple equation form for the excitation force in regular waves has been considered for this research according to [24][25][26], namely F e (t) = A w |F e (ω E )| cos(ω E t + φ + ∠F e ).Here, A w represents the wave amplitude, while |F e (ω E )| is the computed magnitude of the diffractive force term for the excitation frequency ω E and both φ and ∠F e correspond to wave phase shift of initial free surface elevation state and exciting angle respectively.
For multiple bodies, the equations of motion are introduced explicitly for the sake of clarity with respect to the two-body system as follows Hereby, each variable is divided by ς = η, whereas the denominator term With regard to the radiation force determination, different alternatives have been used up to date in order to approximate the integral.The common approach is to use a rational function fit, such as done by [27][28][29].Another very efficient procedure turns out to be a digital filtering technique named the Prony's method [30,31], which is implemented in this case for an exponential function approximation of the convolution integrals of radiation.
Once the different coefficients have been determined and loaded into the program, the program finds solutions for displacement, velocities and accelerations both in time domain and frequency domain for a given initial value problem of position and velocity.For the 2 nd order partial differential equation in the adapted MSD form, one calculates iteratively the integrals of acceleration and velocity through an order reduction in the time domain.Integration methods need to be applied in order to be able to solve the derivatives.Following the state-of-the-art in numerical modeling of WEC's [11], both a classical explicit Runge-Kutta method of 4th order (RK4) and an implicit Newmark Generalized-α method (Nm-gen-α) have been implemented during the conducted research.First one presents an accuracy through a local truncation error on the order of O(h 5 ), while the second one is characterized by being unconditionally stable for β G >= 1  4 and γ G = 1 2 with second order accuracy.The variables β G and γ G are numerical parameters that control both the stability of the method and the amount of numerical damping introduced into the system by the method [52].For γ G = 1 2 there is no numerical damping; for γ G >= 1 2 numerical damping is introduced.For this study in particular, when solving the equations of motion for more than one floating body, two additional matrix inversion methods have been used in Matlab.They are needed when solving the acceleration terms for accessing the Equation of Motion (2.2).Method 1 inverts the total mass M t = M + A r,∞ term and leads to z = M −1 t (∑ F).Then, one only needs to integrate the acceleration twice to get the remaining variables.In the second method, the mathematically considered coefficient matrix A becomes much more dense, since it does not only include the total mass, but it also incorporates PTO damping and hydrodynamic stiffness.This way, the equation leads to the generalized form X = A −1 b, whereas e.g. the state vector becomes in a single degree of freedom the form X = [ z1 , ż1 , z 1 , ..., zi , żi , z i ] T and b = [F e,1 − F r,1 , ż1 , z 1 , ..., F e,i − F r,i , żi , z i ] T .On the RHS, initial values are being inserted for z 1 , ż1 , ..., z i , żi .Then the model calculates the force absorbed and so determines the b vector.Now, it is straightforward to make the inversion of the A matrix [k x k] on the left hand side for obtaining z, ż and z the same way it worked for the previous Mass-Matrix Inversion method.It leads to exactly the same analytic result as method 1, but its matrix formulation can become rather cumbersome depending on the i,j combination.
For the frequency domain, the Response Amplitude Operator (RAO) represents the magnitude of the displacement Ẑ divided by the incoming wave amplitude ς as a function of frequency.For one body with a linear PTO damper as an example, the equation yields Remark is set upon the fact, that analytically, the research has made it possible to determine as well the RAO for two interacting bodies.For this reason, the approach is addressed as an alternating frequency-time domain model, where Having introduced the main hydrodynamic variables and numerical schemes, it is worth to present then the most relevant Ocean Wave Energy equations.They form part of the Post-Processing of results, which follows consequently from the previous steps of the program workflow.
For the Wave Energy absorption of floating bodies, one has to rely on the behavior of the power absorbed P a by the PTO mechanism, and the energy related to it through the integral E = 1 t t 0 P a dt.The power absorbed follows directly from the relationship where P e is the exciting power and P r the radiated power after [14].More generally, the average of the power absorbed yields In terms of available Ocean Wave Energy it is more appropriate to calculate the Wave Energy flux E f per meter crest length, otherwise called Wave power level [14], which is an approximation for real sea waves [6], such that Here, intermediate or shallow waters are assumed and the variables correspond to the group velocity c g = gT 4π , significant wave height H s and the energy period T e .Note, that H s represents the mean wave height of the highest third of the waves and S represents the generalized spectral distribution of the present sea-state.
Another relevant characteristic of a WEC measures the power absorption capability of a Wave Energy device through the ratio of time averaged absorbed power to Energy flux P a E f .It is quoted as capture width c w and for an axisymmmetric WEC oscillating merely in heave one can determine it from the maximum absorbed power P a,max = E f L 2π , such that c w = L 2π [32].Specifically, latter author reveals the so-called capture-width ratio (CWR), stated through the following relationship where D c represents the characteristic dimension in the device, the diameter in the present case.While regarding the given sea-state spectrum, one notices that the peak frequency ω p corresponds to the most energetic sea-state.Ideally, it is possible to set the natural frequency ω n of the PA-system equal to ω p .In other words, shift the motion RAO to the incident sea-state, so that maximum energy can be obtained, while the system is tuned to the wave.In theory, [33] determines the following two resonant conditions for the j-body: The first condition is determined from the undamped free motion of the device by setting F e = 0 and only regarding the real quantities from the frequency domain formulation of the equations of motion.Replacing in latter equation both the complex velocity term U = (iω) Ẑ and the relationship U = F e / 2C h leads to the second condition in Eq.( 14) through reformulation of the RAO motion.
Arrays of WEC's represent the perspective of Ocean Wave Energy, since they enable the scaling of a single floating unit into larger power levels.A number of new effects and variables appear in this field, since little is known practically on the interaction of this kind of layout, where several bodies are positioned close to each other in diverse configurations [34].Its main characteristic is presented by the following expression for the park effect [35,36] q(ω) = P arr (ω) Based on the frequency domain, the q-factor presents the ratio of the total power in the array P arr to the single power generated by one device P iso multiplied by the number of bodies N b .If q < 1 it results in destructive interference, while a ratio > 1 corresponds to a constructive, positive effect on the array.It is noted, that others prefer to call this effect interaction factor, such as [37].They conclude for axisymmetric devices, that there is a another direct approach, such that q = P arr,max (ω) / N b P iso,max (ω), which is applied as well during this study.
Whenever one wants to find the power absorption for an array, the following equation based on linear wave theory is required: where * defines the conjugate transpose, B(ω) is the full radiation matrix, F e (ω) the excitation force component in the frequency domain and U(ω) the complex velocity vector according to [38].Neglecting the second subtracted term in Eq. ( 16) for the case where the velocity is in phase with the excitation force (optimum condition) P arr becomes P arr,max and latter relationships may be applied.

Case Studies
Following the past introduced theory, the next section will introduce different case scenarios (Figure 2): Last figure depicts on its upper part (Subfigures a-b) the sectional sketches of the given cases.One can notice that only the wet part of the bodies is shown, for the reason that it is the only relevant part in the BEM calculation.On the lower part (Subfigures c-d) merely one half of the given bodies is exposed.This is due to the fact, that all bodies are axisymmetric, and for the initial configuration only one half around the xz -symmetry plane needs to be modeled in order to obtain the hydrodynamic coefficients.It is noted that in all three cases, the incoming wave direction is the positive x-axis.For a single body, the cone cylinder case in (a) and (c) initially serves to validate the chosen model approaches of radiation and excitation in heave.Finally, the generic vertical floating cylinder case study in (b) and (d) will lead to conclusions on the interaction of several units with variable distance l λ between them in different array layout configurations [32], such as the ones depicted in the following Figure 3: Cases 1-3 in the sketch represent an aligned (1), a triangular (2) and a squared (or diamond) shaped (3) farm of two, three and four bodies respectively (floaters F1-F4).In all three cases the incident angle is equal to 0 degrees, since the incident wave corresponds to the positive x-axis direction.Furthermore, the distances have been chosen according to the radius of a single floater.This way, f becomes  keeping the bodies more stretched in the x-direction than in the y one.The reason for choosing latter cases relies on the results determined by [39,40], which refer to that triangular and squared/diamond array layouts as fairly suited for Ocean Wave Energy extraction [41].

Simulation parameters
The next endorsed Table 2 summarizes effectively the parameter configuration for the simulation runs.
According to the table both cases are calculated in Scale 1:1.Only the first one is then sized down for Post-Processing through Froude scaling laws to 1:15.9, which is the scale of the experimental setup in de Backer's (GdBck) thesis [18].This becomes part of the first analyzed task as described in the next section.From this step on, the researcher will use the results delivered by the BEM code NEMOH according to the geometries chosen next.It is noted though, that for multiple bodies, a sensitivity analysis of the time step takes place, and values are placed linearly in between 0.1s and 0.5s on a 0.05s step base, resulting in 0.5s as best case scenario.Reference here for has been previous work done in thesis delivered by [18,29].

Results
In order not to exceed the scope of this paper, only results concerning current research tendencies for regular waves are outlined here.For outcome related to basic displacement and power absorbed values in the time domain, the author recommends further reading on his own written thesis, currently being assembled.The following analysis is being subdivided into a function of the number of bodies under study.First, the single cases of the cone cylinder and the generic floating cylinder are discussed.Then, the behavior of arrays of up to four generic PA(s) are documented.

Cone cylinder in heave
The cone cylinder study is related to latter introduced FO3-Buldra project.Its calculation leads to the following Figure 4.
For a range of 15-30 s the graph reflects on the one hand the heave oscillation amplitude of the cone cylinder buoy undergoing regular wave excitation in laboratory measurements (cf.dashed blue curve).The amplitude goes from -0.1 to 0.1 m following a clear sinusoidal pattern.On the other hand, the solid red curve shows the results from the current NEMOH-Matlab model for this time range.It becomes clear, that the current model has not achieved final convergence at this stage, since there is a  certain amplitude modulation over the interval.All over the peaks and troughs, the numerical results deliver higher amplitudes than the experimental results.In general, there is ≈ 16.67% max.deviation from the experiments displacement.However, the phase is nearly matched by the numerical results, which in the original thesis delivered as well higher amplitudes than the test setup.Hence, one can compare the numerical solution with the one obtained by de Backer (cf.dashed magenta curve) [18].It can be noticed, that its numerical solution also peaks over the experimental results, exhibiting a sort of amplitude modulation, which behaves in a very similar fashion to the presently obtained solution.

Generic cylinder in heave
In addition to the previous, the following generic cylinder case serves as a base for an immersion into the parametric of geometrical and hydrodynamic variables, as it will be seen next.

Capture width in heave mode
One of the most relevant part of WEC studies is the capture width determination.The argumentation for the incorporation of the capture width at this stage relies on several reasons.Mainly, it replaces the usual displacement representation as primary output.The fundamental equation is introduced in section 2.2.Depending on how incident wave power and absorbed power are being resolved, different variations determine this equation fraction.For the specific case of the generic floating cylinder PA type, methods carried out by [32], [42], [43] and [44] have been adopted during this investigation (incl.the corresponding nomenclature).The results are presented in the following diagram:  Incidentally, the last figure 6 in this section shows the most significant method for determining the capture width and its ratio after Evans with variable draught in the generic floating cylinder case.De facto, the greater the submerged part, the greater the resulting efficiency gets.However, it is noted, that draughts D d,2 = 5m and D d,3 = 6.8m surpass typical C w values, presenting performance percentages of over 100% for the studied methodologies in [10] .

Multiple degrees of freedom
The incorporation of the BEMIO package for a single floating body case leads to a full calculation of its hydrodynamic coefficients, together with some transformations in the state-space modeling part through the use of the NEMOH tool.While the first package is a Python oriented program for the array calculation in multiple degrees of freedom, the second application contains a FORtran based code for the hydrodynamic parameter acquisition, among other capabilities.Further details are outlined in the following table 3: Table 3 reflects a qualitative assessment ranging from one plus to three plus units depending on its performance capabilities.The criteria behind is simply the flexibility in order to be able to adapt the code to specific problems, as well as its performance in terms of quality and time consumption.OpenWEC, e.g. is the newest code, and reveals itself pretty handy in terms of user interface parametrization, but presents lacks of adaptability to modify the available package.On the opposite, NEMOH has proven to be able to cope with WAMIT performance [54], while providing multiple ways to redefine procedures and therefore obtains the highest evaluation mark.
This allows to modify the existing equation solver for the fundamental cylinder point absorber problem addressed throughout the course of the research.Considering mainly regular waves with an amplitude of A w = 0.5m as an input, the motion is restricted to the resulting degrees of freedom surge (indexes 11), heave (33) and pitch (55).The double numbered notation confers to the previously introduced i,j index allocation within the corresponding matrix elements in the equation of motion.Hence, the subindexes point out values contained in their diagonals.It is noted, that the PTO-Damping has been adapted in advance to suit these additional degrees of freedom.The results are displayed in the next graphical interface:  As it can be seen on Figure 7, the results remain consistent for the heave mode, since nothing with regard to the integration method (time step or sim.time) has been modified.Additionally, one obtains a surge component of max.0.17m, min.−0.17m, which explains why the buoy shifts horizontally during excitation and resulting oscillation.Also, there follows a consequent pitching moment along the rotation axis of the buoy, revealing values up to max.0.03 rad and down to min.-0.03 rad.

Multiple floating bodies in heave
For this part of the paper, the generic floating cylinder remains under consideration.Used for instance in CFD simulations [45][46][47], this type of geometry has become as well predominant in the published literature on BEM calculations [48,49].From this stage on, the study presents only the results related to interacting cylinders in the configurations highlighted in Figure 3.

Two generic cylinders (aligned)
At this stage, two separated, aligned bodies with the wave direction are exposed initially to regular waves of 0.5 m amplitude.The distance d is varied according to the wavelength λ, which is calculated for deep water waves through the relationship ω E = 1 Hz ⇒ T = 6.28 s, hence λ ∼ = 1.56T 2 = 61.59m.The following Table 4 shows the parametrization of the distance    14) the same eigenfrequency of ω 0 = 1.2Hz for each one of the bodies in every distance case.Again, z, ż, z, RAO's and P a are computed for both bodies.For sake of conciseness, merely one distance calculation is outlined.It refers to distance d 3 , which is half the wave length between the two cylinders.The results are observed in the following set of Figures 8a and 8b: The diagrams on the left upper part present the floaters response to the incoming waves on the lower part for an excitation frequency ω e,1 = 1Hz and initial PTO damping C 1 .As it can be seen, both devices move in counter-phase, exhibiting nearly same displacement amplitude ranges (-0.6 -> 0.6 m).These values can be double-checked on the right graph in Figure 8b, reflecting the RAO's for each body in the situations of applying C 1 (initial damping-dashed line-) and C Opt (optimal radiation damping-solid line-acc.to 14) respectively.For ω e,1 , the RAO curves intersect for both damping cases, which means both bodies are oscillating with same amplitude.However, the main difference in both RAO cases is that for optimum damping, the amplitude in the displacement nearly doubles the values of the initial damping.Another interesting phenomena of the RAO graphs happens during the initial case, where the bodies shadow themselves alternately over the frequency range.In near-resonance  Finding reference in Earthquake Engineering [50], a method has been determined for variable distances, which delivers amplitude amplification ratios for each body with respect to the isolated body case.These ratios are not be confused with RAO's though their curve shapes are particularly similar.Limitation to these novel ratios are given generally for irregular waves, since the amplitude in the displacement varies continuously due to the superposition of waves.The graphs reflect for a fixed frequency and variable distances on the horizontal axis the ratio of the maximum amplitude for each body ẑj to the amplitude of the isolated body case ẑ0 on the vertical axis.For Figure 8a this means, that having ẑ0 = 0.64m, the amplitude amplification ratio becomes ẑ1 ẑ0 = 0.938 (body 1) and ẑ2 ẑ0 = 0.984 (body 2), which are the points on the vertical dashed line for d λ = 0.5.On Figure 8 one can also observe, that there is amplification in the movement of the devices, if both bodies oscillate above the horizontal dashed limit of ẑj ẑ0 = 1.For instance, the distances d 2 and d 4 show amplification ratios above one for both bodies.For the first mentioned distance value of d λ = 0.25, ẑ1 ẑ0 = 2.35 and ẑ2 ẑ0 = 1.05, meaning that although the first body receives clearly more energy, the second still gets some motion amplification.This also happens to be the case for d λ = 1, where ẑ1 ẑ0 = ẑ2 ẑ0 = 1.1.Moreover, the interaction tends to decrease with the distance, and the shadowing effect alternates between the bodies according to the diffracted and radiated waves.So, the question arises for what will happen if the PTO damping is optimally fitted to the energy period.The results are given in the next Latter graph reflects only ten operative points from the twelve studied distance cases.This is due to diverging simulations for the initial time step of 0.1 s.However, one can deduce a similar behavior to latter curves, having amplified displacements close to d λ = 0.15, (with ẑ1 ẑ0 = 1.4 and ẑ2 ẑ0 = 1.1) and close to d λ = 0.2, (with ẑ1 ẑ0 = 1.35 and ẑ2 ẑ0 = 1).Again, the curve tends to decay for increasing distances.Further studies need to be done for an optimum excitation frequency in near resonance conditions with the same optimum damping in order to validate this part of the study.Its result is shown in Figure 10, with ten operative working points, where a significant amplification ratio at d λ = 0.75 can be appreciated.Here, ẑ1 ẑ0 ≈ ẑ2 ẑ0 ≈ 2.5, achieving maximum amplification of nearly five times the incident wave amplitude A w,1 .For the interaction of cases 2 and 3 in Figure 3, the two integration methods of RK4 and Nm − gen − α are compared against each other for the displacement in 1 m high waves.

Three generic cylinders (triangular disposition)
Next Figure 11  It can retraced from 11, there is a matching of results for the two integration methods applied with regard to body 1. Bodies 2 and 3 experience a minimum delay of second decimals in the Newmark approach, which can be related to the difference in time step employed, or simply due to the matrix inversion calculation.However, amplitudes seem to be approximated properly, so that the methods are also being applied on four separated bodies next.

Four generic cylinders
Two main separate subcases are analyzed in the following subsections, namely the squared and diamond disposition according to the layout from case 3.For each subcase, the heave displacement (acc.to previous applied integrational methods), is presented selectively [case 3].

Squared disposition
Being characterized by the fact, that the perpendicular distance between each body is equal to f (acc.to the rombus diagonal), the layout is taken from the Case studies, such that Floaters 1-4 appear counter-clockwise, starting from the east-left position.This way, 9 different incident angles are tested (β = 0, 25, 45, 60, 90, 120, 135, 155, 180), and the bold highlighted cases are detailed in the following.Reason therefore is, that there is some expectation in advance, of at least what should be the amplitude tendency within the bodies response [B,1] to [B,4] Subcase 1: beta=0 [Figures 12 and 13] Motion amplitude assessment: ( ẑB,2 ∼ = ẑB,4 ) << ( ẑB,3 > ẑB,1 ) Subcase 2: beta=45 [Figures 14 and 15   Generally speaking, the motion characteristics of previous assessment is understandable according to the incidence angle variation and involved radiation and diffraction effects.The comparison of integrative methods applied works it out to conceive, which might be a suitable result for the displacement amplitude in heave.However, it is noted with regard to the integral methods compared during latter case scenario, that there exists a certain deviation of nearly 10-17% in body's 3 displacement, which is not present for the rest of the bodies 1-4

Diamond disposition
For case 3 in diamond, a similar analysis has taken place.However, the dimensions are here assumed strictly from 3. Despite of this, since results for the displacement of the four floaters do not contribute with newer significant information to the research, previous work reflects additional information for such a case, where deeper knowledge on the interaction factor is demanded [10].

Discussion
Throughout the course of the past section, the credibility of the present modeling approach has been stated.For the initial cone-cylinder PA case analyzed it is found out, that it is very well suited to compare both numerical and experimental solutions to a specific given problem.Experience tells, that for increasing simulation time intervals (e.g.40-100 s) the displacement will tend to converge into the expected regular wave pattern given by Griet de Backer's experiments.Next, it is the capture width for the generic vertical cylinder PA, which ensures, that the case is analog to what can be found in the literature.The research concludes, this study can be easily parametrized into configurations of two, three or four interacting floaters with various frequencies and/or distances.It is being noted again, that the Post-Processing of results has been reduced to specific distance and frequency cases, in order not to exceed the scope of this publication.Despite of this, several aspects can be highlighted from the commented simulation runs.For the three introduced problems of varying frequency and PTO damping (ω e,1 ,ω 0 and C Opt ), the computation diverges at specific distances for the chosen initial time step of 0.1 s.However, a variation of the time step, choosing for instance 0.5 s achieves convergence.This generates a conflict in the sense, that the explicit scheme used initially, presents only preliminary accuracy for 1-2 buoys.However, for three or more floating bodies, the implicit scheme shows at least more suitable results for 0.5 s time step, while the explicit scheme with any time step < 0.5 s does not converge at all.The before mentioned conflict relies on the experience, that generally speaking, explicit integration methods require smaller time steps than implicit ones, leading to less stability.In the distance range d λ = 0 − 0.5 interactions seem to be more intense for all three cases.These might be due to near-trapping effects.Afterward, the curves tend to decay with distance.This behavior has been documented in the literature [34].There is an existing proportionality between the increase of wave amplitude and the displacement of the bodies, which is explained through linear potential theory.With optimum damping, the amplitude in the displacement for a single body increases significantly.This increase is maximal, when ω e = ω 0 .For this case, the displacement nearly quadruplicates the wave amplitude, as the RAO shows.As expected for intermediate distances, the bodies move in counter-phase (0.5λ, 1.5λ ... (n-0.5)λ), while they are in phase for odd values of λ (λ, 3λ ... (2n-1)λ) Summing up, the perspective for future research is enlightening.For a single cylinder, optimum fitting should be regarded primarily with respect to draft adjustment and PTO damping.Not only could there be application of CFD studies, but also experimental testing to any of the given cases.Moreover, capture width ratios should be determined, in order to be able to calculate the system(s) efficiency.For the problems of 3-4 cylinders, one can foresee how greater array configurations may be designed according to the presented principle.Finally, it is the q-factor study which could be optimized for Arrays, by carrying out simulations for various distances and a fixed optimum frequency.Main issue during this study remains the subject of presenting specific results for irregular waves due to the lack of field data for this PA types, where it becomes really difficult to understand phenomena such as shadowing effects, etc. Besides, the displacement discrepancy of floaters 2-4 in Figures 13 -17 requests deeper knowledge on the applied numerical packages, and are therefore out of this research scope.

Conclusions
During the course of this research, a number of milestones have been achieved.First of all, for both unidirectional regular or irregular waves, a solid, flexible, transparent and computationally efficient W2M model has been provided.
With respect to the Single Cylinder analysis it can be said that both time and frequency domain solutions have been verified and correlate for different wave amplitudes, excitation frequencies and PTO Damping values.Although the case has been elucidated for both regular and irregular waves as an input, for the multiple bodies scenario the restriction to simple regular waves as an input has been kept in order to better understand how interactions work.Additionally, it is the use of packages such as BEMIO or OPENWEC, which open new perspectives for the investigation of multiple degrees of freedom.
For the Two-Body System, further conclusions can be retraced from the simulation results, which are summarized in the following: 1. Correlation between the increase of wave amplitude and the displacement of the bodies, which is explained through linear potential theory 2. Optimum PTO Damping nearly doubles the displacement of the two bodies, increasing interaction.This is verified in the literature [23,35,51].3.As it is expected, for intermediate wave lengths the bodies move in counter-phase (0.5λ, 1.5λ .. (n-0.5) λ), while they are in phase for odd values of λ (λ, 3λ .. n λ) The triangle configuration itself turns out to be the a suitable choice in terms of Wave Energy extraction for small arrays according to literature [41,51].Out of this recent references it is also worth to investigate results for greater arrays.In the present research, it can be appreciated, that F1 facing the incoming wave receives an amplitude amplification in displacement which nearly doubles the wave amplitude.Besides, the two remaining cylinders in the background keep moving synchronously, with a displacement minimally lower than the incident wave.The only possible reason points towards floaters 2-3 feeding back at F1 through radiation processes, which need to be studied more carefully.Albeit, the two applied integration methods (RK4 and Nm − gen − α) are relatively in concordance, validating this way the research study.Despite this fact, it is recommendable to rather take advantage of implicit schemes only while modeling multiple bodies for the reasons exposed during the discussion section.Elsewhere, the application of explicit schemes might reveal itself more appropriate in terms of accuracy.
Last but not least, it is the park effect which is pointed out for showing the interactions in the 4-cylinders There is, despite some discrepancy in the amplitude results, quite an agreement within the symmetrical boundaries for this configuration.For beta 0 and 180, there is an opposed behavior within the floater numbers, which is expected though.Also for beta=45, the side floaters to the wave front get more input than the ones fronting the wave.The 90 degrees incidence case merely validates the previous cases, since it is again floaters 1 and 3 located laterally to the wave train, the ones that are significantly amplified.

• Diamond
Although the results present a different shape to what has been done up to now in the field, there appears to be constructive interference for a wide range of combined frequencies and incident angles.Moreover, symmetry around the perpendicular from the top coming incident angle of 90 • confirms the expected behavior of the four floaters.It presents nearly the same interaction for the cases, where the wave train originates either from the right (0 degrees) or the left (180 degrees) of the diamond configuration.Furthermore, there is a decaying tendency for an increasing frequency, both facts that have been dealt with before in the given references.
The complexity of the studied case leaves a broad margin for future investigation and optimization strategies.Hydrodynamic aspects, which are worth to set focus on in the near future, are e.g.geometry modification strategies, application of mooring configurations, behavioral study of shadowing / near-trapping effects, consequences of extreme waves events and last but not least, the improvement of capture width ratios by means of the first three preceding characteristics.With regard to this, it is noted, that recent research reveals and confirms the right choice of the generic PA geometry during the present investigation [53].Hence, it is consistent to assert, that it is the combination of the cylindrical shape with a functional bottom shape, which will surely provide future insights about the performance optimization within the point absorber technology field.
Moreover, with respect to the fundamental theory, it could be of advantage to assume and combine new approaches.From the physical point of view, quantum mechanics, string or even fractal theory might be able to delve into detail within the Fluid-Structure Interaction.These are the grey knuckles in 1 involving e.g.Fractal and quantum theory [55] .They are of primary importance, since studying fluids up to date has revealed itself a task of looking deep into the details, as CFD or Smooth Particle Hydordynamic simulations show.Besides, the application of existing formulations, such as the Helmholtz's wave equation, combined with the existing boundary conditions and Airy theory might enhance new paths for investigation.Focusing furthermore on the energy perspective, it may be even partially possible to bypass or compensate the balance of forces.While neglecting the related sort of uncertainties given in its relationships and assuming the whole as a computational approach, neural networks might set up an innovation breakline both for the estimation of statistics and probability of the PA's overall performance.In other words, it might be determining to consider the system wave to motion in computer science terms a Big Data unit [56] with regard to real sea conditions.Latter one might enable to draw knowledge from a waves database history, resulting in being able to link different affected mathematical and physical scenarios, such as irregular waves handling, wave state prediction or extreme events execution.

Figure 3 .
Figure 3. Top View sketch for the three investigated PA array layouts

Figure 5 .
Figure 5. Comparative of applied methods on capture width determination for the default generic floating vertical cylinder

Figure 6 .
Figure 6.Capture width as a function of variable draught

Figure 9 .
Figure 9. Amplitude ratio ẑj ẑ0 for each j-body in relation to single body for ω e,1 , A w,1 and C 1

Figure 9 :Figure 10 .
Figure 10.Amplitude ratio z j z 0 for each j-body in relation to single body for ω e,1 , A w,1 and C Opt

Figure 11 .
Figure 11.Amplitude ratio ẑj ẑ0 for each j-body in relation to single body for ω e,2 , A w,1 and C Opt

Figure 12 .
Figure 12.Displacement for the triangle configuration of bodies 1-3 over two different integration methods

Table 3 .
Hydrodynamic packages assessment

Table 4 .
Distance variations in the two body PA study