Computational Investigation on Various Small Array Dispositions of Generic Points Absorbers

This thesis completion works it out to deepen into diverse modeling techniques for Point Absorbers. A combined frequency-time domain model is conceived, designed and developed in Matlab with FoRtran as a base., leading 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 though application of NEMOH & BEMIO solvers due to the Boundary Element Methodology. Initially, this Wave-to-motion 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. 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 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. In conclusion, this computational excursion envisions and proposes potentials fields of study, which will surely enhance new connections and link this renewable energy form.


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 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 [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 the modeling procedure for this work is schematized next: Evaluation: In Figure 1 it is noted, that the step-by-step sketch has been determined with specific focus on previous project experience in the field [9].A deeper look in Europe at some delivered 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) 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 [10].Since it is also noted, that computational capacities are growing exponentially [11], 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 [12], 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.[13] 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 [14,15].Last but not least, near-resonance conditions have to be determined for achieving a narrow operational bandwidth of the PA [16], 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 [17] 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 [18,19].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 [20][21][22].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 [19].• 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 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 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) 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 (FSI) is complex, and can be hence interpreted more abstractly, as Figure 2 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, though it also envisions future perspectives (yellow 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 2. 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, 2 reflects the hydrodynamic entry with wave amplitude ζ (in turquoise), including the most relevant parametrizable units.In the middle part, the structure presents the classical vibrational and computational approach.Reaching the top, the spreading ellipsoidal dashed curves represent characteristic force vectors, which might be tuned in advance for the FSI in whatever form research tendencies might converge into.
From the middle "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 z = ∑ F .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 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 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.On the right part of Figure ?? 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 [23][24][25], 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 angles.
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 dividef 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 [26][27][28].Another very efficient procedure turns out to be a digital filtering technique named the Prony's method [29,30], 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 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 [10], 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 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 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 DOF the form X = [ z1 , ż1 , z 1 , ..., zi , żM , 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 analytical result as method 1, but its matrix formulation can become rather cumbersome depending on the i,j 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 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 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 [13].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 [13] 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 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π [31].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, [32] 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.( 13) 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 [33].Its main characteristic is presented by the park effect [34,35]: 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 [36].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 [37].Neglecting the second subtracted term in Eq. ( 15) 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.These are presented in advance in the following Figure 3: Last figure depicts on its upper part (Subfigures a-c) 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 d-f) 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 (d) initially serves to validate the chosen model approaches of radiation and excitation in heave.Next, investigations are conducted as well on the 2-body PA from (b) and (e) for comparison with experimental results from the AWEC.Latter aparatus consists of an upper floater (body 1), which oscillates through a small gap over a larger moored spar structure (body 2).Finally, the generic vertical floating cylinder case study in (c) and (f) will lead to conclusions on the interaction of several units with variable distance l λ between them in different array layout configurations [? ], such as the ones depicted in the following Figure 4: 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 √ 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 [38,39], which refer to that triangular and squared/diamond array layouts as fairly suited for Ocean Wave Energy extraction [40].

Simulation parameters
The next endorsed Table 1 summarizes effectively the parameter configuration for the simulation runs.

Preprints
( 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 [17].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.

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 ??.Its calculation leads to the following Figure 5.
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  , 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) [17].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.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 [31], [41], [42] and [43] have been adopted during this investigation (incl.the corresponding nomenclature).The results are presented in the following: Incidentally, the last figure 7 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 CW values, presenting performance percentages of over 100% for the remaining methods in [9] .

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 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 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.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.

Single Two-body floating AWEC system
Latter results are sufficient information to simulate the second system device, which is the combined two-body PA (AWEC).This device is characterized by a small gap of 0.022 m, where strong flow interactions between the float and the spar are expected.The equation model itself is a 2 x 2 matrix system for the total mass, radiation functions C PTO and stiffnesses of float and spar in heave.The formulation of the exact equations of motion have been obtained from Bret Bosma's thesis [19], along with some the information of the parameters in Table 1.An initial approach on the model for 1 m Amplitude Waves leads to the following outcome in Figures 9 and 10: The diagrams depict on the upper left part the behavior in the displacement of both float and spar.While body 1 (float) exhibits higher oscillations in the red curve (min.:-1m,max.:1 m), this time the buoy stabilizes its movement after nearly 5 s, resulting in a regular curve that follows the free surface elevation shown on the diagram below.The spar (upper left blue curve), on the contrary, moves with less amplitude than the float due to the mooring force applied on the structure.The oscillation amplitude ranges in this case from -0.12 m to 0.12 m, and it is nearly in counter-phase with the displacement of the float.On the right Figure 10 one can see the power extraction from float and spar respectively.The spar contributes minimally to the Power absorbed, while the float peaks constantly at 3 kW, and the mean averaged Power corresponds to nearly 1.5 kW for a single AWEC unit.

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

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 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 For d 1 -d 7 one obtains from Eq. ( 13) 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   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 counterphase, exhibiting nearly same displacement amplitude ranges (-0.6 : 0.6 m).These values can be double-checked on the right graph in Figure 12, 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 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.5Hz.
Another way of interpreting the oscillation amplitude ẑ can be extracted from the following diagrams (Figs.13,15).Finding reference in Earthquake Engineering [49], 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 11 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 13 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 Figure 14: 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 15, 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 4, the two integration methods of RK4 and Nm − gen − α are compared against each other for the displacement in 1 m high waves.Next Figure 16

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

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 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 Subcase 1: beta=0

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

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.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 [33].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)λ) Summing up, it is the park effect which is referenced 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 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 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
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 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 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 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 through linear potential theory 2. Optimum PTO Damping nearly doubles the displacement of the two bodies, increasing interaction.This is verified in the literature [22,34,50].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 λ) The triangle configuration itself turns out to be the a suitable choice in terms of Wave Energy extraction for small arrays according to literature [40,50].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 recieves 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 in concordance, validating this way the research study.
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 of 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 [51].Hence, it is consistent to assert, that it is the combination of the cylindric shape with an 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 FSI.For instance, the application of Helmholtz's wave equation, combined with the existing boundary conditions and Airy theory might enhance new paths for investigation.Focusing furthermore on Einstein's energy perspective, it may be even partially possible to bypass or compensate Newton's law approach, while neglecting the related sort of uncertanties given in its relationships.Conclusively, thinking of 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.

Figure 1 .
Figure 1.Research procedure in the numerical approach

Figure 2 .
Figure 2. Fluid-Structure Interaction: Current modeling scheme and future envisioning strategy

F4Figure 4 .
Figure 4. Top View sketch for the three investigated PA array layouts 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

3. 1 . 3 .
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, Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted:

Figure 6 .Figure 7 .
Figure 6.Applied methods on capture width values for the default generic floating vertical cylinder

PreprintsFigure 8 .
Figure 8. Resulting PA motion for the regular waves case

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

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

Figure 14 .
Figure 14.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 15 .
Figure 15.Amplitude ratio ẑj ẑ0 for each j-body in relation to single body for ω e,2 , A w,1 and C Opt

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

Table 2 .
Distance varitions in the two body PA study 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 ofFigures 11 and 12: