Hydrodynamic Interaction of Ocean Wave Energy Point Absorbers

This work condenses various modeling techniques for different Point Absorber configurations. A combined frequency time domain model will be developed in Matlab-FORTRAN in order to compute the displacement, velocities and the power absorbed in the heave mode. Additionally, a single buoy motion including multiple degrees of freedom will be investigated as well. Therefore, the diffraction-radiation Boundary Element Method solvers NEMOH and BEMIO will be applied in the calculation of the hydrodynamic coefficients, which will determine the solution of Newtons impulse equations of motion. Initially, the Wave to wire model will be validated through comparison with previous experimental results for a submerged cone cylinder shape (Buldra-FO3). A single, generic, vertical floating cylinder will be contemplated then, that responds to the action of the passing waves excitation. Later, two vertical floating cylinders aligned with the incident wave direction will be modeled for a variable distance between the bodies. For both unidirectional regular and irregular waves as an input in deep water, the convolutive radiation force function term will be hereby approximated through the Prony method. By changing the spatial disposition of the axisymmetric buoys, using for instance triangular or diamond shaped arrays of three and four bodies respectively, the study will focus on the interaction effects for regular waves. The results will highlight the most efficient layout for maximizing the energy production whilst providing important insights into their performance, revealing for instance displacement amplification or capture width ratios in near-resonance conditions.


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 1980's by Budal / 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.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 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 by 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 (FPSO) [6] or offshore wind energy foundations.Similarly to the way pump storage plants alternate energy production between turbomachinery and wind turbines [7], WEC's might supply Wind Energy during its production flaws, since it is a much more constant resource [8].The State of the Art of PA's is given through the following table: In Table 1 it is noted that the power specification refers to the averaged rated power per unit.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.
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 [9].Since it is also noted, that computational capacities are growing exponentially [10], 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 [11], this means, that first, hydrodynamics coefficients will be determined for the wetted geometry.Next, these parameters will be attached to the equations of motion in the model solver, which in turn will generate 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.[12] 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's 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 buoy is the only patented PA that is nearly capable of absorbing all the incident wave potential [13,14].Last but not least, near-resonance conditions have to be given for the narrow bandwidth of the PA [15], where to operate at its maximum amplitude motion.
For the current numerical study, three geometries will be under investigation.First, a single body with a cone cylinder shape employed by the SEEWEC project [16] 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 [17,18].With regard to latter PA type, different array layouts will be modeled both in the TD and in the FD, conducting the researchers to study the 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 [19][20][21].

Results
In order not to exceed the scope of this paper, only results concerning actual 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 will be discussed.Then, arrays of up to four generic PA(s) will be documented.

Cone cylinder in heave mode
The cone cylinder study is related to the FO3 project introduced in Table 1.Its calculation leads to the following Figure 1.
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 modulation of the amplitude 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).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.

Capture width in heave mode
One of the most relevant part of WEC studies is the capture width determination.Generally, there is still exists certain controversy in its calculation.Although the fundamental equation is being generally introduced in the theory (cf.section 4.2), depending on how incident wave power and absorbed power are being resolved, different variations determine the equation fraction.Herefore, different methods applied recently by relevant publishers in the field have been selected.For the specific case of the generic floating cylinder PA type, methods carried out by [22], [23], [24] and [25] have been adopted during this investigation (incl.the corresponding nomenclature).The results are presented in the following:  In Figure 2 one can distinguish either capture width values denoted with CW [in m] or capture width ratios CWR [dimensionless] for the default draught of D d,1 = 3m.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 rad s as expected after the opt.radiation damping calculation.Referring to Lewis, the PA clearly surpasses the optimum between 1 − 1.5 rad s .With relationship to the the wave frequency, the CW never exceeds a ratio greater than 0.25, 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.20%).Incidentally, the last figure 3 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 D d,2 = 5m and D d,3 = 6.8m surpass typical CW values, presenting performance percentages of over 100% for the remaining methods.

Generic Cylinder in multiple degrees of freedom
The incorporation of the BEMIO package leads for a single floating body case 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 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).It is noted, that the PTO-Damping has been adapted in advance to suit these additional degrees of freedom.The results are presented in the next graphical output: 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.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 Hz and down to min: -0.03 Hz.

Vertical floating cylinder array
For this part of the research, the generic floating cylinder remains under consideration.Used for instance in CFD simulations [26][27][28], this type of geometry has become as well predominant in the published literature on BEM calculations [29,30].From this stage on, the study will present only the results related to interacting cylinders in the configurations highlighted in Figure 14.Therefore, 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 2 shows the parametrization of the distance  11) the same eigenfrequency of ω 0 = 1.2rad / s 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 5 and 6   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 = 1rad / s and initial PTO damping C 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 6, reflecting the RAO's for each body in the situations of applying C 1 (dashed line) and C Opt (solid line) 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 is that for the initial case, the bodies shadow themselves alternately over  Finding reference in Earthquake Engineering, 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 ẑj to the amplitude of the isolated body case ẑ0 on the vertical axis.For Figure 5 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 7 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 9, 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 3 and 4 in Figure 14, the two integration methods of RK4 and Nm − gen − α are compared against each other for the displacement in 1 m high waves.Next Figure 10  In this case, the free surface elevation reacts according to Figure 5, while excitation frequency ω e,1 and initial C PTO = C 1 = C 2 = C 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 later assembled matrix inversion methods (cf. next section) have been tested herefore sequentially as well, and exactly same results are being delivered for both integration approaches.Figure 10 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 with 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 will not add significant information to the research, the study will focus next on the interaction factor q (cf. Figure 11): The diagram 11 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 β = 0 corresponding 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 (q > 1) for the incident angle comprehended between 0.6 and 1.2 Hz (light blue shape up to the red coloured one).It is noted, that symmetry appears around an angle of 90 degrees, showing another peak at 180 • (ω E = 1Hz).

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.
The capture width for the ensuing generic vertical cylinder PA ensures next, 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, concentrating the main research work.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.Further evidences can be retraced from the simulation results for two aligned cylinders, which are summarized next.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.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.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 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 the 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 Nm − gen − α) 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 4-cylinder 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.

Materials and Methods
Throughout the course of this research, certain conditions have been met for the incoming waves and specifically for the floating systems.They are being summarized in the following section.

Assumptions
According to the Wave Mechanics, 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.
• 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 thereby 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 H with ratios such that H

Hydrodynamic Equations of Motion
The following Figure 12 presents the solving procedure for a generic hydrodynamic problem of a single floating body in waves: The Fluid-Structure Interaction (FSI) is split into parts.On the left branch tree of the workflow, one can see first 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 z = ∑ F .This procedure assumes that the surrounding water medium can be replaced by springs and stiffness's, 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 Second term is the radiation force, which is described by the equation The third term inserts the hydrostatic restoring force into the sum of forces, yielding Fourth term considers the application of a mooring force on the floating system, such that Last term includes the PTO-force, which is simply reflected as a damping mechanism in the form On the right part of Figure 12 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 chosen for this research according to [31][32][33], 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 .
One has to solve M-equations of motion for M given bodies in N degrees of freedom according to the following matrices • Total mass: Sum of structural mass plus added mass Hereby, the indices i and j denote respectively the M and N boundaries for the given given problem, while δ ij corresponds to the Kronecker delta.On the right hand side (RHS) of the equation, the resulting vectors are introduced by the forces respectively in the form F = [F ,1 ...F ,N ] T depending on the chosen domain.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 [34][35][36].Another very efficient procedure turns out to be a digital filtering technique named the Prony method [37,38], 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 TD and FD 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 of the problem 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 [9], both a classical Runge-Kutta method of 4th order (RK4) and an implicit Newmark Generalized-α method (Nm-gen-α) have been implemented for 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 β >= 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 in Newton's Equation of Motion (4.2).Method 1 inverts the total mass M t 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 is also incorporates PTO damping and hydrodynamic stiffness's.This way, the equation leads to the generalized form X = A −1 b, whereas the state vector becomes the form X = [ z1 , ż1 , z 1 , ..., zM , żM , z M ] T and b = [F e,1 − F r,1 , ż1 , z 1 , ..., F e,M − F r,M , żM , z M ] T .On the RHS, initial values are being inserted for z 1 , ż1 , ..., z M , żM .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 [M x N] 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 analytical result as method 1, but its matrix formulation can become rather cumbersome depending on the M,N combination.
For the FD, 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 It is noted that analytically, the research has made it possible to determine as well the RAO's for two interacting bodies.
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 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 straightforward through the relationship where P e is the exciting power and P r the radiated power after [12].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 [12] for real sea waves: Here, deep water is 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.
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π [22].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 maximum energy can be obtained while the system is tuned to the wave.In theory, [39] 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 FD 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.( 11) through reformulation of the RAO motion.
Arrays of WEC's represent the future 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 [40].Its main characteristic is presented by the park effect [41,42]: Based on the FD, 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 [43].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 FD and U(ω) the complex velocity vector according to [44].Neglecting the second subtracted term in Eq. ( 13) 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 order in the variables appearance, the next section will introduce results for different case scenarios.These are presented in the following Last figure depicts on its upper part (Subfigures a-b) the sectional sketches of the given cases.One can notice that only the wetted 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 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 l λ between them in different array layout configurations [45], such as the ones depicted in the following ( 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.25r, 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 [46,47], which conclude that triangular and diamond array layouts are optimally suited for Ocean Wave Energy extraction.

Simulation parameters
In the next endorsed Table 3 the parameter configuration for the simulation runs is detailed.According to the table all three 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 thesis [16].This becomes part of the first calculated task.

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

Figure 4 .
Figure 4. Resulting PA motion for the regular waves case :

Figure 6 .
Figure 6.RAO's for C 1 and C Opt for a variable frequency range

Preprints
(www.preprints.org)| NOT PEER-REVIEWED | Posted: 19 November 2016 doi:10.20944/preprints201611.0102.v1 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 rad/s.Another way of interpreting the oscillation amplitude ẑ can be extracted from the following diagrams (Figs.7, 9).

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

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

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

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

Figure 11 .
Figure 11. Park effect for variable frequency in the diamond configuration of 4 floaters 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 Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 19 November 2016 doi:10.20944/preprints201611.0102.v1length.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 (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 TD the values for A r at ω → ∞ and Impulse Response Functions (IRF's) for the radiation force function k r (t)(retardation kernel).All these terms become part in the solution of Newtons impulse equations of motion, whereas specific terms have to be first integrated in convolutive form, for instance the radiation force F r .

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

Table 1 .
Point Absorber Projects

Table 2 .
Test For d 1 -d 7 one obtains from Eq. (