A Procedure to Calculate First ‐ Order Wave ‐ Structure Interaction Loads in Wave Farms and Other Multi ‐ Body Structures Subjected to Inhomogeneous Waves

: A typical assumption when performing analytical, numerical, and experimental studies in wave–structure interaction in multi ‐ body problems such as for wave farms and very large floating structures is the homogeneity of the wave field. Important interactions between the floating elements are dependent on the direction, amplitude, and phase of the waves acting on each. Then, wave homogeneity is probably unrealistic in near ‐ shore areas where these installations are to be deployed. In the present work, an existing interaction method, which allows the use of standard boundary element diffraction codes for solving the first order wave structure linear potential for each unique geometry in the problem, is shown to be able to account for inhomogeneous sea states across the domain of a multi ‐ body problem requiring only minimal modification to its implementation. A procedure to use the method to include arbitrary incoming undisturbed wave conditions at each body is presented. A verification study was done by using an artificial numerical configuration to mimic an inhomogeneous wave field in a standard diffraction code, which was used as a reference. The results obtained using the interaction ‐ method based procedure are shown to be in excellent agreement with the reference ones. Furthermore, an example of frequency inhomogeneity of the wave field in a wave farm is shown and the effects on the motion amplitudes and absorbed power are presented illustrating the applicability of the procedure.


Introduction
Computing the wave structure interaction in a multi-body problem has a wide range of applications in many offshore and coastal engineering problems configurations.It is particularly relevant for two engineering design problems: large floating structures and offshore wave energy farms.In both, an array of submerged body elements is present and wave structure multi-body interaction exists.This interaction is two-fold: (1) wave systems that are the result of scattering of incoming waves in individual submerged elements, which recursively excite all the elements of the array; and (2) wave systems that radiate away from each element as a result of the element motion, which excite the other elements of the array.Therefore, the wave excitation forces on each element differ from an element in isolation and so do the wave induced motions and derived quantities.component design.The generic floating bridge in [1] is an example of the former with its 35 floating pontoons supporting a continuous bridge girder, while the large floating hydrocarbon storage facility in [2] is an example of the latter.VLFSs also tend to include a large number of similar floating bodies such as in the floating airport design in [3,4] with tens of floating elements.Inter floater wave-structure interaction effects can lead to amplification of wave run-ups, the appearance of wave trapping modes, and excitation of natural modes.All these phenomena can significantly reduce the efficiency of the installation by impeding its operability, and even lead to catastrophic structural failures if not properly accounted for.

Array Interaction in Offshore Wave Energy Farms
Offshore wave energy exploitation is expected to be implemented in the form of wave farms composed of many wave energy converters (WECs).These farms should be compact to reduce area occupation and costs.
The layout of the farm affects the behavior of the units including the energy harvested due to the devices' hydrodynamic interactions.These can be constructive or destructive, depending on the layout of the farm and the incoming wave system.The qfactor is the ratio between the absorbed power of the farm and the power absorbed by the same number of WECs, each one in isolation.Optimized array layouts can result in a qfactor higher than 1, effectively doubling or more the power output.In ideal conditions (e.g., in narrow banded sea-states and optimally tuned devices), the q-factor can increase dramatically.Unfortunately, the q-factor is extremely sensitive to wave period and direction.Optimal configurations appear as peaks in the optimization hyperspace surrounded by vast valleys of destructive interference ones [5].Additionally, it has been proven that optimal configurations for a given wave direction are necessarily accompanied by destructive interactions for other wave headings [6].Here, it is worth noting that, to the best of the author's knowledge, wave farm studies so far do not take into account changes in the layout occurring from steady and low-frequency wave drift motions that may disturb optimal layouts.Nor do they consider variable local wave conditions throughout extensive wave farms.It is also important to bear in mind that dangerous phenomena similar to the ones described in Section 1.1 are also possible in wave farms.

The Direct Method
From the above, it should be evident that a proper identification of the interaction characteristics is paramount in the design spiral of VLFSs and wave farms for assuring an efficient and safe installation.Fast numerical methods are then needed and linear wavestructure diffraction codes implementing a boundary element approach (BEM) in the frequency domain are typically used, for example, WAMIT [7] and NEMOH [8] may be used when current is not considered, whereas HYDROSTAR [9] and MULDIF [10] can consider the effect of current with different degrees of accuracy.
The above codes solve the interaction problem by considering all the bodies as part of a boundary value problem for the oscillating potential flow solution-usually denominated the direct method.The computational burden can be reduced if one or two symmetry planes exist in the array of bodies, where only half or a quarter of the elements needs to be explicitly considered.However, in addition to the usually unrealistic restriction of symmetry limitation, this comes at the expense of redefining the problem into problem specific non-standard (symmetric or anti-symmetric) modes of motion.It then presents difficulties for further utilization of the results in time-domain codes (e.g., SIMO [11]) where nonlinear forces can be added under a weakly linear approach.When the number of bodies is in the order of tens or hundreds and the geometry of the bodies is not simple, the computational burden becomes prohibitive for a standard PC and inefficient for design optimization procedures when using a typical modern workstation.

Matrix Methods
In order to circumvent the aforementioned limitations, [12] developed an interaction theory for the computation of wave excitation forces for arrays of bodies that is exact within the limits of linear wave theory and accounts for both the evanescent and progressive wave systems.The method requires the knowledge of the hydrodynamic properties only for each unique body geometry as opposed to all the bodies in the direct method, giving way for significant computational savings.The interaction effects are introduced by resorting to algebraic relations that express the notion that waves scattered or radiated by one body are incident waves on the others.The procedure involves a diffraction transfer matrix, unique to each unique body geometry in the array in isolation, that transforms incident waves to the body in the corresponding scattered waves.However, both these wave systems are expressed as vectors of cylindrical partial-wave coefficients.This seriously undermines the applicability of the theory, as standard diffraction codes such as the ones mentioned in Section 1.3, which use planar waves in their formulation, cannot be used to compute the diffraction transfer matrix for bodies of arbitrary geometry.For this reason, application of this method is usually relegated in favor of the direct method, which is available in publicly available codes or commercial software.Despite this, important developments and applications of the method to arrays of bodies should be mentioned.Such are the ones by [13], [14], and [15], considering the analytic solution of arrays of cylindrical bodies; and [16] and [17], considering axisymmetric and arbitrary body geometries, respectively, to compute the diffraction matrix.
More recently, [18] presented a method to determine the necessary cylindrical partial-wave coefficients for an arbitrary geometry from information of the wave potential contained in a wet control surface surrounding the body.This method was then applied in [19] for deriving a diffraction transfer matrix from plane incident waves and where a new operator-the force transfer matrix-is formulated.The method presented in [19] effectively allows the usage of standard diffraction codes as input to the interaction methods originated by [12].It does not account for the evanescent waves in the immediate vicinity of the bodies nor is it formulated for infinitely deep water, but for practical purposes, the method holds remarkably well, as shown in [20] and [21].

Inhomogeneous Wave Environment Conditions
Another fundamental drawback of the direct method is its inability to consider prescribed arbitrary inhomogeneous wave environments as input to individual bodies in a multi-body problem, while still accounting correctly for the hydrodynamic interactions between them.On the other hand, inhomogeneous wave conditions are a realistic assumption and its effects relevant for VLFSs, as demonstrated by [22].
Accounting for inhomogeneity and diffraction interaction using standard BEM diffraction codes is stated as an open challenge by [23] on the subject of near shore VLFSs.As an example, consider the work by [24], who presented a scheme to include individual arbitrary wave conditions on individual pontoons of an end-anchored floating bridge (but neglected interaction), together with the study by [25] on pontoon interaction of the same bridge design (but who did not consider inhomogeneity).Also worth mentioning is the variable bathymetry consideration throughout a WEC array in [26]; however, the methodology was based on a non-standard dedicated BEM formulation for computing the flow.
The need for accounting for inhomogeneous wave conditions in studies of WEC farms was reported in [19].This is obviously particularly important when analyzing diffracted wave patterns in light of array layout analysis and optimization, as the coherent wave crest phases are artificial and lead to inaccurate outputs.Then, phase resolving models for wave climates become necessary as additional input if the study is purely numerical.

The Present Study
The need for considering inhomogeneous wave conditions in fast numerical simulations for some important multi-body problems has been identified.Irrespective of this, and to the best of the author's knowledge, no methodology has been published in the literature explicitly demonstrating how to achieve it, or published studies where it has been applied.
In the work herein presented, it is shown that the interaction methodology in [19] can be exploited to account for wave inhomogeneity in multi-body problems with minimal modification of the procedure.The original method characteristics are kept, but there is a need to run the matrix method for each prescribed special distribution of local undisturbed incoming wave amplitude, direction, period, and phase, throughout the domain of interest, and change the implementation to account for inhomogeneous wave headings.
The problem formulation, interaction theory, and practical procedures are presented in Section 2. The verification and an illustrative example of application is given in Section 3 for a set of five bodies, which can be considered as wave energy converters.Discussion and overall conclusions, including further suggested work, are presented in Section 4.

Methodology
The problem formulation and methodology follow the ones presented in [7,12,19], but have been adapted to the consideration of inhomogeneous environmental waves, for which some nomenclature modification is required.Derivation of equations present already in the aforementioned references are omitted here to avoid repetition.Only the final expressions are presented and modified for the specific problem of the wave inhomogeneity, where necessary.

The Direct Method
Consider the problem of N bodies submerged in a fluid that is considered inviscid and incompressible.The fluid flow in the presence of waves is described by harmonically varying velocity potentials,  , in the vicinity of each body j: where  is the angular frequency;  denotes time; and  is the complex potential at the point with cartesian coordinates  , ,  satisfying the Laplace equation: ∇  0.  only needs to be valid in the interior of the domain delimited by the still water surface, the rigid bodies' surfaces, the horizontal bottom, and an arbitrary cylindrical surface surrounding each body.The dispersion relation for linear progressive water waves is satisfied: where  is the wavelength;  is the wave number;  stands for the standard gravity acceleration; and ℎ is the water depth.The total potential is obtained from the sum of the solutions of two boundary problems targeting the solution of the radiation potential,  and the diffraction potential,  : where  is the disturbance potential, or scattering potential, due to the presence of all the bodies. is the undisturbed linear potential of a regular plane wave of amplitude  propagating through body j at an angle  with the  1,0,0 vector and is expressed by: The radiation potential is dependent on each  potential, which denotes the potential corresponding to the flow induced by a unit amplitude harmonic velocity in the  mode of the jth body.It is valid everywhere in the fluid and given by: where  is the complex amplitude in the  mode.The following body boundary conditions for the two problems are applied to the mean wet surfaces of each body: where  is the body surface normal and     • , with    denoting the unit displacement vector of the mode  .If   ∧  ,  1, … ,  , then   and   , thus forming a consistent solution valid everywhere in the fluid.The problems can then be solved using Green's second and third entities to formulate the following integral equations for the radiation and diffraction potentials: In Equations ( 7) and ( 8),   is a point in the surface  of all the bodies together and  ,   is the green function denoting the potential at  due to a point source of strength 4 located at   .Then, using the linearized Bernoulli equation, the radiation and diffraction pressures are attained.From these, motions and other response amplitude operators can be obtained (see [7] for details).Equations ( 7) and ( 8) are the integral equations used by standard boundary element method diffraction codes for the solution of the linear problem.It is evident that Equation (8) does not allow for different incoming wave amplitudes or directions to be prescribed to the different bodies.On the other hand, it is trivial to acknowledge that the problem can be solved for each body in isolation, which has been pursued by [24] and [27], but then no interaction is accounted for.Therefore, if one wishes to take advantage of these codes, an external inhomogenizing procedure must be applied while taking advantage of the computations for an isolated body: this can be achieved by the solution of an interaction problem.

The Interaction Problem
The potential of a planar wave hitting the body j can be expressed in cylindrical coordinates  ,  ,  with the origin on the center of the body, using the relation  ,    cos  , sin  , where ,  are the horizontal cartesian coordinates of the centre of body j.The representation is as follows: where    and    are vectors with an infinite number of m elements including  0.
These are defined by: where  is the mth order Bessel function of the first kind, rendering    as a vector of (incoming) basis functions, and  is the prescribed undisturbed wave phase angle at the center of body j.Note that  is only a measure of the origin of the wave when there is only one body and that it is equal to    cos   sin  for the case of homogeneous undisturbed plane waves.For the radiated and scattered potentials, one has: where    depends on the geometry of each body and    also on the mode of motion, being both infinite vectors of coefficients.   is a vector of (scattering) basis functions.
While the phase information of the undisturbed waves is prescribed in Equation ( 10), the same needs to be solved for the radiation and scattering problems.Therefore, the Hankel function is used instead of the Bessel function to include both real and imaginary parts in the definition of    : where  stands for the Hankel function of the second kind of order m.Note that the infinite vectors in Equation ( 9) through Equation ( 13) can be truncated to   , making way for an approximation with arbitrary precision, which is used herein.
From the sum of the undisturbed wave at body j with the incoming scattered waves at body j (from the other bodies), one defines the incident wave potential on body j: Introducing a transformation matrix   , so that         , Equation ( 14) becomes: It is important to note that   is the actual transformation matrix with shape 2 1 2 1 , whereas the indexing  denotes the pair of bodies to which these elements are to be applied-the elements 2 1 2 1 are specific for each pair.Using Graf's addition theorem ( [19,28]), the mnth element of the matrix comes: where: A diffraction transfer matrix   can be defined that transforms incident wave amplitudes to body j into scattered wave amplitudes from the same body.Then, the scattered wave amplitudes that are required to determine the potential in Equation ( 15) can be obtained from the solution of the following system of linear equations: where   is a vector of N elements, each element being a matrix with the shape 2 1 2 1 , corresponding to each jth body.
For the radiation problem, the formulation is similar, but replacing    with  ,   as follows: At this point, it should be clear that at no point in the interaction theory presented does the undisturbed incident wave filled at the bodies needs to be equal in amplitude, direction, or the phase equal to  .The modification relative to [19] and [12] herein implemented is to consider prescribed and possibly independent  ,  , and  for each jth body, instead of , , and  starting from Equation (9).

Excitation Forces and Hydrodynamic Coefficients
In a similar way to the diffraction transfer matrix, a force transfer matrix   can be defined as a linear operator that transforms incident linear waves to body j into linear excitation forces in body j.Then, the excitation forces due to undisturbed incident planar waves on the jth body come: where    is a  1 vector and   is a  2 1 matrix.The radiation problem solution for computing the hydrodynamic coefficients does not depend on the incoming wave system, so the formulation is kept exactly the same as in [19]: and the added mass matrix, , and the wave radiation damping matrix, , are given by:

Diffraction and Force Transfer Matrices from Standard Diffraction Codes
It is important to recall that these matrices are linear operators that transform incident cylindrical wave fields on a body into scattered wave fields, from the body, or forces, applied on the body.They are thus obviously independent of the wave field and are specified for each unique body geometry in isolation.Then, by considering a set of 2 1 incident plane wave directions on the jth body in isolation,   and   can be obtained from the systems of equations: where the body index has been dropped.In Equations ( 25) and ( 26),  , ; stands for the scattered, or incident, th cylindrical partial wave amplitude coefficient resulting from an undisturbed incoming plane wave of unit amplitude with heading angle  .In Equation (26),  , stands for the excitation force in mode  from an incoming undisturbed plane wave field with heading  .Equations ( 25) and ( 26) express a regression problem with both systems being underdetermined if  * 2 1 and overdetermined if  * 2 1. Realizing this, the actual angles to use should, in general, cover the range from 0 rad to 2 rad to avoid overfitting the coefficients within a small range.Furthermore, the larger is  * , the better the quality of the obtained diffraction and force transfer matrices' coefficients.
While  , can be determined from Equation (10), there is a need to have the information of the scattered wave amplitude coefficients to use in Equation (25).These are obtained from the following relation: can be obtained numerically from the solution of Equation ( 8) and then determining the potential at a cylindrical surface centered at the body with radius -denoted control surface from hereon.One should also bear in mind that only the progressive waves are taken into account, meaning that R should not be too close to the surface of the body for Equation ( 27) to be valid.Furthermore, the control surface of one body in the multi-body problem cannot intersect the surface of another body; on the other hand, control surfaces can intersect each other.Finally, one needs the radiated wave coefficients to use in Equation (22).These are obtained from the radiation potential distribution on the same aforementioned control surface for each mode of motion  , as follows: The derivation of Equations ( 27) and ( 28) is essentially based on application of a Fourier transform to the velocity potential and consequent exploitation of the orthogonality of the depth dependent terms obtained.Details can be found in [18] and are not covered herein.

Practical Computational Procedures
Now that the theoretical aspects have been addressed, the practical steps involved in computation of the interaction problem are explained.Furthermore, a suggested stepwise approach for exploiting the results for arbitrary prescribed inhomogeneous spectra is presented.
The following lists the procedure to define a discrete inhomogeneous field for a multi-body diffraction problem using the interaction theory described in Sections 2.2-2.4: 1. Impose an arbitrary directional wave variance spectrum at the location of each body including both amplitude and phase information.2. Optional: Define a reference location (e.g., the location of one of the bodies).3. Produce a 3-dimensional array  with shape    , where  is the number of bodies,  is the number of frequencies, and  is the number of wave headings in the entire spectral range of all bodies.
Note that this defines the inhomogeneity of the wave environment and effectively produces a locked relation of the phases, amplitudes, and directions of the incoming undisturbed waves between all the bodies.If step 2 is implemented, the complex amplitudes can be made relative to a unit amplitude wave system with origin at the reference point.
The following lists the computational steps for acquiring the necessary single body hydrodynamic data for solving the interaction problem.All steps can be carried out using standard frequency domain diffraction codes: 4. Compute the added mass, damping coefficients, and hydrostatic stiffness matrix (if motions are to be calculated) for each unique geometry in isolation.The hydrodynamic coefficients are used in Equation ( 22). 5. Compute the excitation forces on each unique geometry in isolation considering a set of at least 2M + 1 (preferably more) incident wave headings covering the range 0 to 2 and the  frequencies.The forces are used in Equation ( 26) to determine the force transfer matrix . 6. Compute the diffraction potential at a control surface surrounding each unique body geometry in isolation for a set of at least 2M + 1 (preferably more) incident wave headings covering the range 0 to 2 and the  frequencies.The potential distribution at the control surface is used in Equation ( 27) to determine  , , which then feeds Equation ( 25) to determine the diffraction transfer matrix . 7. Compute the radiation potential for each motion mode of each unique body geometry at the same control surface mentioned in the previous point.The potential distribution at the control surface is used in Equation ( 28) to determine     to be applied in Equation (22).
The following lists the steps to acquire the multi-body hydrodynamic quantities: 8. For each frequency  ,  1, … ,  : a. Solve the system of linear equations in Equation (19) to determine the diffraction    for each body and for each of the  wave headings.Here, the corresponding element of  needs to be used for each    .
b. Compute the excitation forces for each body and for each of the  wave headings using Equation ( 21).c.The excitation force on each body is the sum of the excitation forces applied on the body from all the  angles.9. Solve the system of linear equations in Equation ( 19) with the substitution in Equation (20) to determine the radiation    for each jth body and then compute the full multi-body added mass and damping coefficient matrices using Equations ( 22)-( 24).
After step 9, the full coupled hydrodynamic quantities and excitation forces have been resolved for all the frequencies.Response amplitude operators can then be computed following the usual approaches, but it is important to stress that the amplitude relations between the incoming undisturbed wave fields between all the bodies must be maintained.

Results
This section presents the verification of the methodology presented in Section 2 using a special configuration of a diffraction problem using WAMIT for the inhomogeneous wave multi-body problem.An example of the case where frequency inhomogeneity is present is also shown and analyzed.
The freely available MATLAB code in [29] was modified to include wave inhomogeneity input and computation capability and then used for obtaining the interaction problem solution's raw results in terms of hydrodynamic coefficients and excitation forces.Further pre-and post-processing and analysis were done using python scripts implemented by the author.

Verification
Consider the multi-body problem consisting of five identical bodies, each one free to oscillate only in heave (see Figure 1a).Each body is a cylinder with diameter  and draft  /2.The bottom is flat and the water depth ℎ 20 10.The bodies' positions are symmetric about the  axis.
One wishes to attain the added mass and damping matrices as well as the excitation forces on the array including the interaction effects considering different undisturbed wave conditions (amplitude, period, heading, and phase) on each body.One cannot arbitrarily introduce different incoming waves on different points in WAMIT's Equation ( 8) for the diffraction problem, but Equation ( 9) can be used by making one or more bodies act as (numerical) wave generators.In Figure 1a, body W1 is a small cylindrical body forced-oscillating with unit amplitude in heave at an angular frequency .This body's geometry is a 1/8 scaled down version of the remaining bodies.Note that Figure 1a represents a view from above and that the relative sizes of bodies 1 through 5 and W1 are in scale.Consider the set of  problems where   tanh ℎ / ,  1, . . .,  , with  / as given in Table 1.The resulting radiated wave patterns (real part of the wave elevation) for selected oscillation frequencies of W1 are shown in Figure 2. From visual observation, it is seen that the waves can be considered locally planar at a short distance from each body.Similar plots but where scattering due to presence of the other bodies is included are shown in Figure 3, where it is clear that the wave field is considerably disturbed within the wave frequency range studied, exhibiting significant array scattering phenomena.In Table 1, the prescribed incoming undisturbed wave conditions at each body center location are listed for each forced-oscillating body frequency.In each, the amplitude of the waves  and phases are relative to the motion of  .The amplitudes could have been normalized as suggested in Section 2.5 (step 3).However, these were left unchanged to deliberately have very small amplitudes of the incoming waves to verify possible numerical truncation errors within the implementation of the algebraic procedure in the interaction theory.
The values in Table 1 were obtained by running a simulation in WAMIT for the radiation problem of body  oscillating in heave.The incoming undisturbed wave amplitudes at the points corresponding to the center of each body were extracted from WAMIT's output and considered as local long-crested incident plane waves with headings as given in Figure 1a.Table 1 is the input to the inhomogeneous interaction procedure described in Section 2.5 (step 3).The mesh of each body and its panel distribution are presented in Figure 1b.The panel length is about 0.26 λ for the shortest wavelength considered.Therefore, the results for these high frequencies may not be accurate, but its relevance is negligible considering the focus of the present study.
The control surface surrounding the body, for performing steps 6 and 7 in Section 2.5, was the submerged surface of a truncated cylinder with radius 0.55 D extending from  0 to  ℎ.The surface is discretized with   panels, where  200 is the number of points in the vertical direction and  2 10 is the number of points in the angular direction.Cosine spacing in the vertical direction was applied.The body mesh, mass characteristics, and water depth were kept the same as for the cylinder case in [19].
The results obtained are shown in Figures 4-7.These correspond to the excitation forces' real and imaginary components, the added mass, and the radiation damping, respectively.Real part of the excitation forces in heave for each body."DIRECT", "IHIT", and "NO_INT" denote the solution using WAMIT with the direct method, the solution using the interaction theory with inhomogeneous wave field, and the solution for a single isolated body using the direct method and body  , respectively.F, ρ, g, and A*, stand for the excitation force, the fluid density, the standard gravity acceleration, and the heave amplitude of body W1, respectively.Figure 5. Imaginary part of the excitation forces in heave for each body."DIRECT", "IHIT", and "NO_INT" denote the solution using WAMIT with the direct method, the solution using the interaction theory with inhomogeneous wave field, and the solution for a single isolated body using the direct method and body  , respectively.F, ρ, g, and A*, stand for the excitation force, the fluid density, the standard gravity acceleration, and the heave amplitude of body W1, respectively.Figure 6.Added mass coefficients for heave motion."DIRECT", "IT", and "NO_INT" denote the solution using WAMIT with the direct method, the solution using the interaction theory, and the solution for a single isolated body using the direct method and body  , respectively.A denotes the added mass and ρ stands for the fluid density.Wave radiation damping coefficients for heave motion."DIRECT", "IT", and "NO_INT" denote the solution using WAMIT with the direct method, the solution using the interaction theory, and the solution for a single isolated body using the direct method and body  , respectively.A denotes the added mass, ρ stands for the fluid density, and ω is the angular frequency.
The verification was clearly accomplished by observing Figures 4-7, with all quantities showing excellent agreement between the interaction theory and the direct method.As a reference, the quantities relative to the bodies in isolation are also presented and demonstrate the need for properly accounting for both the diffraction under inhomogeneous wave conditions and radiation interaction within the frequency range considered.Naturally, the interactions become less important for very short waves, and wave scattering itself becomes negligible for very long waves.Both these assertions are straightforwardly verified in Figures 4-7.
An overall nominal and relative error quantification is presented in Figure 8 for the excitation forces considering the interaction theory and the bodies in isolation.The values shown are averages, but the error variation between the different bodies is marginal, meaning that the plots are good indicators of the errors.In Figure 8a, one can see that the nominal error of the interaction theory is close to zero except for the two shortest wavelengths.The corresponding relative mean errors are, respectively, 16.3% and 748% for the two highest frequencies, which compares to 39.6% and 41.67% when neglecting interactions.The reason for this error amplification could be due to numerical truncation errors.Specifically, when using the interaction theory with input from WAMIT, which has only single precision, when the excitation forces, and the quantities from which they are derived from, become very small, as is the case for these two frequencies.Another possible explanation is that the consideration of local incoming undisturbed planar waves hitting each body loses its validity for wavelengths closer to the body diameter, and the circular crests cannot be neglected.One should realize that, in both cases, the source of the error is not related to the interaction theory.
Irrespective of the origin of the error amplification, it is important to bear in mind that for wavelengths of this order (or shorter), the diffraction forces are close to zero and the radiation couplings between the bodies are likewise negligible (see Figures 4 through  7).Therefore, one can use the procedure herein described, deeming the verification excellent for practical cases.Here, it is worth to note that for the other frequencies studied, the error typically increased from about 0.1% to a maximum of 0.9% as the frequency increased, while the isolated body approximation had errors in the range of about 10-40% (see Figure 8b).It is also worth mentioning that it is always possible to simply use the nointeraction approach for wavelengths outside the range where multi-body interaction is significant.
(a) (b) . Mean error throughout the bodies for the excitation forces."IHIT ERR" and "NO_INT ERR", denote the mean error using the interaction theory with inhomogeneous wave field, and the mean error using the direct method only with an isolated body and body  , respectively.The nominal error is the magnitude of the (complex) vector connecting the solution obtained with that of the multi-body reference direct method (a); the relative error is the nominal error divided by the magnitude of the reference value (b); the values in (a) and (b) correspond to averages of these quantities.

Frequency Inhomogeneity in a WEC Farm
It is important to realize that even in the (unrealistic) case of the model used in Section 3.1 for verification, WAMIT cannot consider frequency inhomogeneity.
Consider the case of a wave farm deployed in a location where local conditions force short waves to break and dissipate in some areas, effectively acting as a low-pass filter in sections of the array.This situation is easily considered using the procedure presented herein.Using the same model of Section 3.1 for illustration purposes, one can simply remove the higher half of the incoming undisturbed wave frequency (λ/D < 6) components for bodies 1 and 4, deemed as WECs in a wave farm to impose such a condition.
In Figure 9, the motions of each body are plotted for the case similar to Section 3.1 and the one where the local frequency truncations have been applied to the incoming undisturbed waves with no Power Take Off (PTO) or viscous damping applied.Likewise, the absorbed power in both configurations is shown in Figure 10, where a PTO damping equal to each body's radiation damping at λ/D = 5.0 is applied.10 is more relevant for wave farms.In both figures, results pertaining to bodies 1 and 4 for the lower half of the frequency range are the response of the bodies solely due to the excitation forces originating from the scattering of the incoming undisturbed wave system in the other bodies, but where the presence of all bodies is considered for diffraction computation, plus the ones arising from the multi-body radiation.As expected, the results show that there are motions and absorbed power beyond the undisturbed incoming local wave system frequency range at these bodies.
It is interesting, although not surprising, to see how the absence of high frequency waves in bodies 1 and 4 led to amplification of motions and power iwo the resonant frequency for bodies 3 and 5, showing a destructive interaction that is no longer there.In fact, the difference between the sum of the product of squares of undisturbed wave amplitudes with frequency of undisturbed incoming waves between both cases shows a reduction of 43% in what can colloquially be defined as a measure of the total available wave power, whereas the power absorption is reduced by only 21%.Naturally, these assessments are but illustrative of the interaction phenomena, as the wave field and the characteristics of the wave farm used herein are not deemed to be of a real wave farm design.

Discussion
In the present study, it was shown that the interaction theory originally presented in [12], and modified in [19] to allow the use of standard diffraction codes for solving the linear potential, can also be used when accounting for inhomogeneous sea states across the domain of a multi-body problem without modifying its fundamental supporting algebraic relations.A procedure to include the inhomogeneous sea state in the implementation of the theory was presented, which avoids the need to modify the inner core functions of the code available in [29] considerably.
The new procedure should allow researchers to accurately include wave field inhomogeneity in numerical models of wave farms and accurately assess the sensitivity of the results to the perturbations on otherwise perfectly coherent long-crested wave components.These perturbations are realistic and will significantly influence the results of many published studies concerning optimal wave farm layout and wave farm power output in general, which, to the best of the author's knowledge, have not considered such perturbations.In this respect, one should bear in mind that coastal or nearshore structures are more prone to be hit by inhomogeneous wave fields arising from shore reflections and bathymetric variations.Furthermore, near-shore wave climates tend to exhibit shorter peak periods, therefore exhibiting a relative higher energy for waves within the wavelength range where interaction effects are non-negligible.
Additionally, dangerous modes of operation including wave-trapping modes are dependent on the direction, frequency, and amplitude distribution at each body composing the structure or array of bodies.For the specific case of VLFSs such as floating bridges or wave farms composed of WECs connected through a structure, there is a need to properly account for inhomogeneity and interaction among floating elements, which can possibly excite natural flexural modes that are otherwise neglected if full diffraction interaction is not accounted for.
The focus of the present paper is to demonstrate the capability of the interaction theory to cope with inhomogeneous wave conditions and the results are shown only for regular waves.However, the generalization to irregular waves is straightforward and can be implemented following the steps listed in Section 2.5.
It is also worth noticing that the absence of the consideration of evanescent wave modes in the present theory did not show any effect on the results for the frequency range where multi-body interactions are non-negligible.
Regarding future work, perhaps the most interesting immediate step would be to include the output regarding mean second order drift forces.Work on delineating a possible strategy has been initiated by the author, with the considerable increase in complexity being the biggest challenge.

Figure 1 .
Figure 1.Layout and geometry: (a) Positions of the bodies (distance values denote multiples of D); (b) WAMIT mesh for each body.

Figure 2 .
Figure 2. Radiated wave patterns from forced oscillations of body W1.

Figure 3 .
Figure 3. Wave patterns from forced oscillations of body W1 including scattering from the other bodies.

Figure 4 .
Figure 4. Real part of the excitation forces in heave for each body."DIRECT", "IHIT", and "NO_INT" denote the solution using WAMIT with the direct method, the solution using the interaction theory with inhomogeneous wave field, and the solution for a single isolated body using the direct method and body  , respectively.F, ρ, g, and A*, stand for the excitation force, the fluid density, the standard gravity acceleration, and the heave amplitude of body W1, respectively.

Figure 7 .
Figure 7.Wave radiation damping coefficients for heave motion."DIRECT", "IT", and "NO_INT" denote the solution using WAMIT with the direct method, the solution using the interaction theory, and the solution for a single isolated body using the direct method and body  , respectively.A denotes the added mass, ρ stands for the fluid density, and ω is the angular frequency.

Figure 9 .
Figure 9. Heave motions with only wave radiation damping."1-5" and "1*-5*" denote full wave field and local frequency truncated wave field at bodies 1 to 5, respectively.ξ denotes the heave motion amplitude and A stands for the local incoming undisturbed (full wave field) wave amplitude at each body.

Figure 10 .
Figure10.Power absorbed with PTO damping applied."1-5" and "1*-5*" denote full wave field and local frequency truncated wave field at bodies 1 to 5, respectively.P, ρ, ω, g, and A, stand for the power absorbed, the fluid density, the angular frequency, the standard gravity acceleration, and the local incoming undisturbed (full wave field) wave amplitude at each body, respectively.

Figure 9
Figure 9 is more interesting for the general case of multi-body interaction, while Figure10is more relevant for wave farms.In both figures, results pertaining to bodies 1 and 4 for the lower half of the frequency range are the response of the bodies solely due to the excitation forces originating from the scattering of the incoming undisturbed wave system in the other bodies, but where the presence of all bodies is considered for diffraction computation, plus the ones arising from the multi-body radiation.As expected, the results

Funding:
This research was funded by the Research Council of Norway under the Ekstraordinaere GM 2020 funding attributed to SINTEF Ocean.Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.Data Availability Statement: Not applicable.