HYDRODYNAMIC ASSESMENT ON MOTION PERFORMANCE CHARACTERISTICS OF GENERIC WAVE ENERGY POINT ABSORBERS

This work condenses various modeling techniques for different Point Absorber configurations. An alternating frequency time domain model is implemented in MatlabFORTRAN in order to compute the displacement, velocities and the power absorbed in the heave mode for both single and multiple body configurations. Coupling of different degrees of freedom are merely contemplated with regard to a single buoy motion. NEMOH and BEMIO solvers are applied in the solution of Newtons second law according to the Boundary Element Methodology. Initially, this Wave-to-Wire model is validated through 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. In deep water, we approximate the convolutive radiation force function term through the Prony method. Using for instance triangular or diamond shaped arrays of three and four bodies respectively, the study delves into the interaction effects for regular waves. The results highlight the most efficient layout for maximizing the energy production whilst providing important insights into their performance, revealing displacement amplification-, capture width-ratios and the commonly known park effect.


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 . Studied intensively since the oil crisis during the late 1970's by Budal and Falnes [Budal and Falnes(1975)] in Scandinavia and by Evans [Evans(1976)] 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 [Pecher and Kofoed(2017)]. PA's can be floating or bottom-mounted structures, with a Power-Take-Off (PTO) mechanism that varies 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 [Marquis et al(2010) Marquis, Kramer, and Frigaard]. 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 (FPSO) [De Rouck(2006)] or offshore wind energy foundations. Similarly to the way pump storage plants alternate energy production between turbomachinery and wind turbines [Merino and Veganzones(2012)],, WEC's might supply Wind Energy during its production flaws, since it is a much more constant resource [Stoutenburg et al(2010)Stoutenburg, Jenkins, and Jacobson]. The State of the Art of some relevant PA's is given through the following table:  Table 1 it is noted, that the power specification refers to the averaged rated power per unit. The choice of the devices has been determined 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 five projects correspond to single Point Absorbers, while the three remaining products represent Multipoint Absorbers. From the total number of introduced projects, only three are being tested right now, namely the AWEC, Wave Star and WaveNET. For a broader perspective on latter mechanisms, further global documentation can be found in [Magagna and Uihlein(2015)].
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 [Li and Yu(2012)]. Since it is also noted, that computational capacities are growing exponentially [Koomey(2011)] , 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-FD-TD analysis will take place. For a so called wave to wire (W2W) model [Garcia-Rosa(2014)] , 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 FD or TD.
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. [Falnes(2002)] 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 [Cruz and Salter(2006), Mei(1989) Bosma(2013)] . With regard to latter PA type, different array layouts will be modeled both in the TD and in the FD, 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 [Evans and Porter(1997), McIver(2002), Garnaud and Mei(2010)]. It is noted however, that the added value of the present work relies on ensuring a valid and novel approach in the Post-Processing assesment 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 Regarding the floating system, further limitations are introduced next:  Mooring, viscous and damping forces have been determined empirically, using previous thesis results from Bosma [Bosma(2013)]..  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 either unidirectional regular or irregular waves are used as an input to the model. Linear potential or so called Airy theory applies to the previous mentioned sea states, mainly based on small amplitude wave heights with ratios such that ℎ < 1corresponding to shallow water and < 1referring to deep water. Here, ℎ represents the water depth, and 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 FD. These are added mass ( ), radiation damping ( ℎ ) and excitation force ( ) as sum of diffraction ( ) and Froude-Krylov forces ( ). Additionally, one obtains in the TD the values for at → ∞ and Impulse Response Functions (IRF) for the radiation force function ( )(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

HYDRODYNAMIC EQUATIONS OF MOTION
The following Figure 1 presents the solving procedure for a generic hydrodynamic problem of a single floating body in waves: . 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 ∑ = + + ℎ + + , 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 Second term is the radiation force, which is described generally by the equation   The third term inserts the hydrostatic restoring force into the sum of forces, yielding with ℎ as the hydrodynamic stiffness. Fourth term considers the application of a mooring force on the floating system, such that where stands for the mooring coefficient. Last term includes the PTO-force, which is simply reflected as a damping mechanism in the form whereas PTO represents the PTO damping coefficient. On the right part of Figure [fig:1] there is the spectral distribution ( ), 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 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 [Perez and Fossen(2007), Fossen(1994), Faltinsen(1990)],namely represents the wave amplitude, while | ( )| is the computed magnitude of the diffractive force term for the excitation frequency and both and ∠ correspond to wave phase angles.
For multiple bodies, the equations of motion are introduced explicitly for the sake of clarity as follows Once the different coefficients have been determined and loaded into the program, the program finds solutions for displacement, velocities and accelerations both in TD and FD for a given initial value problem of position and velocity. For the 2 order partial differential equation in the adapted MSD form, one calculates iteratively the integrals of acceleration and velocity through a   problem order reduction in the TD. 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 [Li and Yu(2012)], both a classical 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 (ℎ 5 ), while the second one is characterized by being unconditionally stable for > 1 4 and = 1 2 with second order accuracy. 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 ([eom]). Method 1 inverts the total mass = + ,∞ term and leads to ̈= −1 (∑ ). Then, one only needs to integrate the acceleration twice to get the remaining variables. In the second method, the mathematically considered coefficient matrix becomes much more dense, since it does not only include the total mass, but it also incorporates PTO damping and hydrodynamic stiffness.
Remark is set upon the fact, that analytically, the research has made it possible to determine as well the RAO for two interacting bodies. This is the reason, why the approach is addressed as an alternating FD-TD model, where Having introduced the main hydrodynamic variables and numerical schemes, it is worth to present next the most relevant Ocean Wave Energy equations. They form part of the Post-Processing of results, which follow 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 by the PTO mechanism, and the energy related to it through the integral = 1 ∫ 0 . The power absorbed follows directly from the relationship where is the exciting power and the radiated power after [Falnes(2002)] . 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 per meter crest length, otherwise called Wave power level for real sea waves: Here, deep water is assumed and the variables correspond to the group velocity = 4 , significant wave height and the energy period . Note, that represents the mean wave height of the highest third of the waves.
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 . It is quoted as capture width and for an axisymmmetric WEC oscillating merely in heave one can determine it from the maximum absorbed power , = 2 , such that = 2 [Babarit (2015)]. . While regarding the given sea-state spectrum, one notices that the peak frequency corresponds to the most energetic sea-state. Ideally, it is possible to set the natural frequency of the PAsystem equal to . 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, [Falcao(2010)] 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 = 0and only regarding the real quantities from the FD formulation of the equations of motion. Replacing in latter equation both the complex velocity term = ( ) and the relationship $U = F_e / 2 C_{h}$ leads to the second condition in Eq.([eq:11]) 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 [Child(2011)]. Its main characteristic is presented by the park effect [Babarit(2013), Cruz(2008) Based on the FD, the q-factor presents the ratio of the total power in the array arr to the single power generated by one device multiplied by the number of bodies . If < 1it 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 [   there is a another direct approach, such that $q = P_{\text{arr,max}}(\omega) / N_b P_{\text{iso,max}}(\omega)$, 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, ( )is the full radiation matrix, ( )the excitation force component in the FD and ( ) the complex velocity vector according to [Bellew et al(2013) Stallard]. Neglecting the second subtracted term in Eq. ([eq:13]) for the case where the velocity is in phase with the excitation force (optimum condition) arr becomes arr,max and latter relationships may be applied.

3.1.CASE STUDIES
Following the past introduced theory, the next section will introduce different case scenarios. These are presented in advance in the following Figure [    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 xzsymmetry plane needs to be modeled. 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) leads to conclusions on the interaction of several units with variable distance between them in different array layout configurations [Babarit et al(2008) Babarit, Borgarino, Ferrant, and Clement],, such as the ones depicted in the following Figure [fig:3]: Cases 1-3 in the sketch represent an aligned (1), a triangular (2) and a 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 √11. 25 , 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 [Borgarino et al(2012)] , which conclude that triangular and diamond array layouts are optimally suited for Ocean Wave Energy extraction.

SIMULATION PARAMETERS
In the next endorsed Table [tab:2] the parameter configuration for the simulation runs is detailed. . This becomes part of the first analyzed task as described in the next section.

RESULTS
In order not to exceed the scope of this paper, only results concerning current research tendencies are outlined here. For outcome related to basic displacement and power absorbed values in the TD, 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 the FO3 project introduced in Table [tab :1]. Its calculation leads to the following Figure [fig: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 approx. 16.67max. 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) . 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 parametrization of geometrical and hydrodynamic variables, as it will be seen next.  [dimensionless] for the default draught of ,1 = 3 . It is noted, that the black curves have been colored differently, since they present the rather expected curve shape acc. to the representative RAO in heave mode. The green curve on the two right diagrams presents the optimal incident efficiency for opt. radiation damping acc. to Evans and Lewis, respectively. It can be appreciated, that after Evans, there is clearly one single intersection close to 1 as expected after the opt. radiation damping calculation. Referring to Lewis, the PA clearly surpasses the optimum between 1 − 1.5 . With relationship to the the wave frequency, the CW never exceeds a ratio greater than 0.27, while there is a significant agreement in the methods presented by Babarit and Price on the one hand, and deviation on methods by Evans and Lewis on the other(max. 40%).

CAPTURE WIDTH IN HEAVE MODE
Incidentally, the last figure Figure6 in this section shows the most significant methods for determining the capture width and its ratio after Babarit and Evans with variable draught in the generic floating cylinder case. De facto, the greater the submerged part, the greater the resulting efficiency after Evans method. However, draughts ,2 = 5 and ,3 = 6.8 surpass typical CW values, presenting performance percentages of over 100% for the remaining methods.

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. 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 = 0.5 as an input, the motion is restricted to the resulting degrees of freedom surge indexes: (11), heave (33) and pitch(55). 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, 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.17 , min.: −0.17 , 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 research, the generic floating cylinder remains under consideration. Used for instance in CFD simulations [Agamloh et al(2008)], this type of geometry has become as well predominant in the published literature on BEM calculations [Eriksson(2005), Falnes and Hals(2012)]. From this stage on, the study presents only the results related to interacting cylinders in the configurations highlighted in Figure [fig:3].

TWO GENERIC CYLINDERS (ALIGNED)
Therefore, two separated, aligned bodies with the wave direction are exposed initially to regular waves of 0.5 m amplitude. The distance is varied according to the wavelength , which is calculated for deep water waves through the relationship $\omega_E=1 \hspace{0.1cm} \text{Hz} \Rightarrow T=6.28 \hspace{0.1cm} \text{s}$, hence $\lambda \cong 1.56 T^{2} = 61.59 \hspace{0.1cm} \text{m} $. The following Table [     For sake of conciseness, merely one distance calculation is outlined. It refers to distance 3 , which is half the wave length between the two cylinders. The results are observed in the following set of Figures [fig:8] and [fig:9]:The diagrams on the left upper part present the floaters response to the incoming waves on the lower part for an excitation frequency ,1 = 1 and initial PTO damping 1 . As it can be seen, both devices move in counterphase, exhibiting nearly same displacement amplitude ranges (-0.6 : 0.6 m). These values can be double-checked on the right graph in Figure [fig:rao3], reflecting the RAO's for each body in the situations of applying 1 (dashed line) and Opt (solid line) respectively. For ,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 is that for the initial case, the bodies shadow themselves alternately over the frequency range. In near-resonance conditions on the other hand, it is body 2 which shadows body 1 predominantly over the frequency range of 1.25 − 1.5 .
Another way of interpreting the oscillation amplitude ̂can be extracted from the following diagrams (Figs. 10 -12).
Amplitude ratio 0 for each j-body in relation to single body for ,1 , ,1 and 1 Figure 8:   [Hernández et. al., Vol.4 (Iss.4)   Finding reference in Earthquake Engineering [Hong-Wen Maa et al (2006)], 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 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 to the amplitude of the isolated body case 0 on the vertical axis. For Figure [fig:10] this means, that having 0 = 0.64 , the amplitude amplification ratio becomes 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 Figure [   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 = 0.15, (with 1 0 = 1.4and 2 0 = 1.1) and close to = 0.2, (with 1 0 = 1.35and 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 [fig:12], with ten operative working points, where a significant amplification ratio at = 0.75 can be appreciated. Here, 1 0 ≈ 2 0 ≈ 2.5, achieving maximum amplification of nearly five times the incident wave amplitude ,1 .
Amplitude ratio 0 for each j-body in relation to single body for ,2 , ,1 and Opt

THREE GENERIC CYLINDERS (TRIANGULAR DISPOSITION)
For the interaction of cases 3 and 4 in Figure [fig:sketch2], the two integration methods of 4and − − are compared against each other for the displacement in 1 m high waves. Next Figure [fig:figampamp3tri] reflects for the triangular configuration the agreement of results:

FOUR GENERIC CYLINDERS (DIAMOND DISPOSITION)
In this case, the free surface elevation reacts according to Figure [fig:8], while excitation frequency ,1 and initial PTO = 1 = 2 = 3 values remain constant. This time, the diagrams are split into parts, emphasizing the correlation between the integration methods in all three of them. The only significant change to the previous 2-body calculation is the time step, which becomes 0.5 s instead of the initially 0.1 s used in all cases before. It is additionally noted, that the two assembled matrix inversion methods have been tested here for sequentially as well, and exactly same results are being delivered for both integration approaches. Figure  [fig:figampamp3tri] also reports important information about the displacement. Body 1 starts to oscillate and reaches a higher amplitude than the incoming wave (ranges from -0.9 m : 0.9 m). Comparatively, bodies 2 and 3 move synchronously over time, but delayed with respect to body 1. Their amplitude reveals a displacement smaller than the surface elevation amplitude, with an operative range of -0.45 m to 0.45 m. For case 4, a similar analysis has taken place. However, since results for the displacement of the four floaters do not contribute with significant information to the research, the study focuses next on the interaction factor q (cf. Figure  [fig:14]):   The diagram [ fig:14] reveals the q-factor as a function of the wave frequency and the incident wave angle for the given case study of a diamond configuration including 4 floating cylinders. The wave angle shifts its direction counter-clockwise, with initial = 0corresponding to a wave train from the right plane. It becomes clear, that the curve peaks at nearly 1 Hz for the excitation frequency. Also, there is beneficial interaction ( > 1 ) for the incident angle comprehended between 0.6 and 1.2 Hz (light blue shape up to the red colored one). It is noted, that symmetry appears around an angle of 90 degrees, showing another peak at 180 ∘ ( = 1 ). However, there is an unexpected tendency of approaching zero for considerably small frequencies. The expectation should be hereby, though, that the factor would rather tend to zero instead [Wolgamot, H.(2012)].
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 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. For a more detailed view of the operational bandwidth covered, the reader is addressed to the thesis currently being carried out by the author, where e.g. irregular waves are treated. Despite of this, several aspects can be highlighted from the commented simulation runs. For the three introduced problems of varying frequency and PTO damping ( ,1 , 0 and 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. In the distance   range = 0: 0.5interactions 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 [Mei(2012),Newman (2001)]. 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 = 0 . For this case, the displacement nearly quadruplicates the wave amplitude, as the RAO shows. As expected for intermediate distances, the bodies move in counterphase (0.5 , 1.5 ... (n-0.5) ), while they are in phase for odd values of ( , 3 ... (2n-1) ) The triangle configuration itself turns out to be an optimal choice in terms of Wave Energy extraction for 3 bodies according to literature ][][. This research concludes that F1 facing the incoming wave gets 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 that need to be studied more carefully. Moreover, the two applied integration methods (RK4 and − − ) are in concordance, validating this way the research study.
Conclusively, it is the park effect which is pointed out for showing the interactions in the 4cylinder diamond configuration. Although the results present a different perspective 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 axis (from the top coming incident angle of 90 ∘ ) confirms the expectations for 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 perspective for future research is encouraging. 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 from this study remains the subject of presenting results for irregular waves, where it becomes really difficult to understand phenomena such as shadowing effects.

CONCLUSIONS & RECOMMENDATIONS
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 W2W model has been provided.   With respect to the Single Cylinder analysis it can be said that both TD and FD 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 understand better how interactions work. Additionally, it is the use of packages such as BEMIO or OPENWEC, which open new perspectives for the investigation of multiple DOF.
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 by linear potential theory 2) Optimum PTO Damping nearly doubles the displacement of the two bodies, increasing interaction. This is verified in the literature ][][. 3) As it is expected, for intermediate wave lengths the bodies move in counterphase (0.5 , 1.5 .. (n-0.5) ), while they are in phase for odd values of ( , 3 .. n ) For the three and four bodies, not only dispositions are being proposed for absorption in specific distance cases. Also, different integration methods are applied for validation purposes, as well as a new diagram form for the q-factor in a regular waves regime.
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 of last but not least, the improvement of capture width ratios by means of the first three preceding terms.