Wave Exciting Force Maximization of Truncated Cylinders in a Linear Array

: This study focuses on the determination of optimum layout conﬁgurations for a linear array of identical mutually interacting truncated cylinders. Optimum conﬁgurations correspond to those that maximize either the total heave exciting force acting on all cylinders of the array or the heave exciting force applied on pairs of cylinders within the array. For achieving this goal, we developed and applied an e ﬃ cient optimization numerical process (ONP), where a robust hydrodynamic numerical model, capable of solving the di ﬀ raction problem of the examined multi-body arrangement in the frequency domain, was appropriately coupled with a genetic algorithm solver in an integrated computational environment. Initially, the e ﬃ ciency of the ONP is demonstrated by comparing results with those of other investigations that resulted from the deployment of classical optimization methods. Then, ONP is applied for a linear array of nine cylinders for determining the optimum layout conﬁgurations under the action of the head and perpendicular to the array waves, and for di ﬀ erent maximum allowable array lengths. The resulting optimum conﬁgurations correspond to a random positioning of the cylinders within the array. Nevertheless, they are characterized by the formation of clusters of closely-positioned cylinders, which induce positive hydrodynamic interactions in terms of maximizing the exciting forces.


Introduction
Arrays of truncated cylinders can be utilized in various marine applications either as: (i) structural elements of multi-column floating structures (e.g., semi-submersibles and tension-leg platforms, ocean bridges, etc.) or (ii) as multi-body configurations for energy extraction (e.g., arrays of wave energy converters operating as point absorbers).In principle, the applications connected with items (i) and (ii) require minimization and maximization of the heave exciting forces, respectively.For ocean platforms utilized in the oil and gas industry for instance, the minimization of the environmental loading is the path to safety, while wave energy converters are designed to exploit the maximum percentage of the available energy.Clearly, the energy conversion into mechanical power is maximized in proportion with the maximization of the induced heaving forces.
In relevant situations, intense hydrodynamic interactions between the cylinders of the array are introduced as a result of multiple wave scattering, which can have a direct impact on the exciting forces applied on the cylinders and, thus, on the performance indicators connected with the required functionality of the cylinders.For a given geometry, number of cylinders within the array and specific incident wave characteristics, the aforementioned interactions and exciting forces depend strongly upon the layout configuration of the array; namely, the arrangement of the cylinders (e.g., linear arrangement with one row of cylinders, rectangular grid with multiple rows and columns of cylinders, and circular arrangement) and the center-to-center distance between the adjacent cylinders.Ultimately, an optimum layout configuration can be determined in terms of maximizing or minimizing the exciting forces applied on the cylindrical bodies within the array, in accordance with the relevant performance requirements.
For realizing such an optimum layout configuration, the multiple wave scattering problem of the mutually interacting truncated cylinders of the array should be solved in a robust and accurate way.This problem has been extensively and effectively tackled in the past by many researchers, who developed and applied relevant analytical and numerical methods aiming at: (a) the direct evaluation of the exciting forces applied on the cylinders for a specific predefined array (e.g., [1][2][3][4]), or (b) the assessment of the effect of different layout configurations on the exciting forces (e.g., [5]) or, finally, (c) the investigation of the mode trapping phenomena (e.g., [6]) and their effect on the cylinders' exciting forces (e.g., [7][8][9][10]).It should be noted that in all of the above studies the cylinders were placed at either symmetrical locations relative to the direction of the incident waves or at equally spaced positions within the examined configuration.Moreover, in [8] it was emphasized that an array of closely spaced cylinders can result in a large increase of both the free-surface elevation and the wave exciting forces, while for the case of a linear array, Chatjigeorgiou et al. [10] demonstrated that the change of the center-to-center spacing between the cylinders of the array has a significant effect on the induced surge and heave exciting forces.
Although numerous methods have been developed so far for solving the hydrodynamic problem of truncated cylinders within arrays, little attention has been given to the determination of an optimum layout configuration that maximizes or minimizes the exciting forces applied on the aforementioned bodies.In this context, the authors refer to the study by Kagemoto [11], who developed the inverse hydrodynamic interaction theory by coupling the hydrodynamic interaction theory of [2] the descent method for finding minima [12] and the sequential unconstraint minimization technique [13], in order to determine the optimum layout configuration of an array of truncated cylinders, which leads to the minimization of the total surge exciting force applied on them.Kagemoto [11] concluded that for specific wave frequencies, a non-equally spaced array configuration leads to ten times smaller total surge exciting forces, compared to the equally spaced array.It should be noted that the deployment of the aforementioned classical optimization techniques requires extensive analytical calculations to solve the optimization problem and it might impose constraints on the maximum number of independent design variables that can be considered.
For overcoming these drawbacks, alternative optimization methods can be applied, such as the evolutionary methods [14], which by utilizing an initial population of random solutions and mimicking the evolution principle of nature, result in a stochastic search.Genetic algorithms (GAs), introduced by [15], correspond to an evolutionary algorithm method inspired by the selection process of nature, where in a competition the stronger individuals survive.Selected individuals that have better scores of the objective function of the physical problem, have a better chance for reproduction and creation of offspring.Consequently, the fit individuals gradually increase the population size of the possible solutions of the optimization problem.Considering the wider area of hydrodynamics, GAs have been extensively and efficiently used by many researchers to solve numerous optimization problems, including-(a) the shape geometry optimization of a ship [16], (b) minimization of the scattered wave energy and, thus, reduction of the drift force applied on a truncated cylinder surrounded by multiple cylindrical bodies by optimizing the layout or the geometrical characteristics of the outer bodies [17,18], (c) minimization of the drift force applied on a truncated cylinder by optimizing the flexural rigidity of a concentric annular plate [19], (d) optimum hydroelastic design of flexible floating structures [20,21], and (e) determination of optimum layout configurations of wave energy converters in terms of maximizing absorbed power [22,23] or minimizing costs [24].
Energies 2020, 13, 2400 3 of 19 In the present paper, by utilizing GAs we developed an optimization numerical process (ONP) for determining the optimum layout configuration of a linear array of identical mutually interacting truncated cylinders.The term "optimum" refers to layout configurations that maximize either the total heave exciting force acting on all cylinders of the array or the heave exciting force acting on appropriately selected pairs of cylinders within the array.The developed ONP consists of two distinctive numerical components, appropriately coupled in an integrated computational environment.The first one corresponds to a hydrodynamic model (HM), which solves the diffraction problem of the examined multi-body arrangement in the frequency domain in an accurate and computationally efficient manner, by utilizing the matched eigenfunction expansion technique and the "direct" solution methodology for the numerical implementation [10,25].Accordingly, it calculates the linear heave exciting forces of the cylinders.The second component corresponds to a GAs solver (GAS), aiming at solving the aforementioned constraint optimization problems.The efficiency and the accuracy of the developed ONP is demonstrated, initially, by comparing results with those of Kagemoto [11].Then, the ONP is applied for a linear array of nine cylinders and extended results are presented focusing on the optimum layout configurations of the array and the maximum heave exciting forces for head and perpendicular to the array waves, as well as for various maximum allowable lengths of the array.

Problem's Definition
A linear array of Q mutually interacting, identical, truncated cylinders of radius b and of draft h 2 is placed in an area of finite and constant water depth h (Figure 1).The cylinders of the array are distributed non-uniformly along a length l array defined as (Figure 1a): where l q denotes the center-to-center spacing between two adjacent cylinders q and q + 1, while l out is the distance of the external cylinders of the array from the fictitious boundaries employed for imposing an upper limit on the length of the array.
Energies 2020, 13, x FOR PEER REVIEW 3 of 19 total heave exciting force acting on all cylinders of the array or the heave exciting force acting on appropriately selected pairs of cylinders within the array.The developed ONP consists of two distinctive numerical components, appropriately coupled in an integrated computational environment.The first one corresponds to a hydrodynamic model (HM), which solves the diffraction problem of the examined multi-body arrangement in the frequency domain in an accurate and computationally efficient manner, by utilizing the matched eigenfunction expansion technique and the "direct" solution methodology for the numerical implementation [10,25].Accordingly, it calculates the linear heave exciting forces of the cylinders.The second component corresponds to a GAs solver (GAS), aiming at solving the aforementioned constraint optimization problems.The efficiency and the accuracy of the developed ONP is demonstrated, initially, by comparing results with those of Kagemoto [11].Then, the ONP is applied for a linear array of nine cylinders and extended results are presented focusing on the optimum layout configurations of the array and the maximum heave exciting forces for head and perpendicular to the array waves, as well as for various maximum allowable lengths of the array.

Problem's Definition
A linear array of  mutually interacting, identical, truncated cylinders of radius  and of draft ℎ 2 is placed in an area of finite and constant water depth ℎ (Figure 1).The cylinders of the array are distributed non-uniformly along a length   defined as (Figure 1a): where   denotes the center-to-center spacing between two adjacent cylinders  and  + 1, while   is the distance of the external cylinders of the array from the fictitious boundaries employed for imposing an upper limit on the length of the array.The array is subjected to the action of monochromatic incident waves of linear amplitude A and of circular frequency ω that propagate at angle β, relative to the horizontal O axis.All cylinders are assumed to be restrained in waves; accordingly, only the diffraction problem is considered.
Let the X coordinates of the centers of the cylinders in the global OXYZ coordinate system (Figure 1a) X q , q = 1, . . ., Q, denote the locations of the cylinders within the linear array.We sought to find the optimum values of the design variables X q , q = 1, . . ., Q, along the length l array (Equation ( 1)) or the equivalent optimum layout configuration of the cylinders for two different problems.The first problem aims at maximizing the total heave exciting force, F tot Z , acting on all cylinders of the array, while the second one aims at maximizing the heave exciting force, F subtot Z , acting on appropriately selected pairs of cylinders.The maximization of the objective functions, F tot Z and F subtot Z , of the two aforementioned constrained optimization problems is mathematically expressed as: while, additionally, for both problems the design variables, X q , q = 1, . . ., Q, are subjected to the following constraints: X q ∈ 0, l array .
Equation ( 4) mathematically expresses the avoidance of overlapping between two adjacent cylinders, while Equations ( 5) and ( 6) impose restrictions on the placement of the two external cylinders, so that Equation (1) can be satisfied.
For solving the two aforementioned optimization problems, a suitable optimization numerical process (ONP) is developed and applied.ONP consists of two distinctive numerical components appropriately coupled in an integrated computational environment.The first component corresponds to a hydrodynamic model (HM) that solves the diffraction problem of the examined multi-body arrangement, taking into account all hydrodynamic interactions between the cylinders, and, accordingly, it calculates the heave exciting forces applied on the cylinders.The second component corresponds to a GAs solver (GAS) that solves the optimization problem for the pre-defined design variables, objective function, and relevant constraints.In the following section, the two numerical components of the developed ONP as well as their coupling procedure are described.

Hydrodynamic Model
The hydrodynamic analysis of the array is conducted in the frequency domain and it is based on the 3D linear wave diffraction theory.The fluid is assumed to be inviscid and incompressible, and the flow to be irrotational.Accordingly, the velocity potential is introduced for describing the fluid motion.The first-order (linear) boundary value problem is formulated and solved by utilizing the matched eigenfunction expansion technique, combined with the "direct" solution methodology for the numerical implementation.This semi-analytical approach was developed and described in [10,25].For making this study self-contained, we briefly present the main aspects of the approach taken in [10,25].
Each cylinder q = 1, . . ., Q in the array defines two liquid regions-the outer region A extending to infinity and containing all other cylinders j q of the array, and the lower region B situated below each cylinder.A local polar coordinate system r q , θ q , z fixed on the center of each cylinder q is introduced Energies 2020, 13, 2400 5 of 19 with the vertical z coordinate common for all cylinders, being fixed at the bottom and pointing vertically upwards.The regions A and B are then contained within r q ≥ b, 0 ≤ θ q ≤ 2π, 0 ≤ z ≤ h and 0 ≤ r q ≤ b, 0 ≤ θ q ≤ 2π, 0 ≤ z ≤ h 1 , respectively, where h 1 denotes the vertical distance between the bottom of each cylinder and the sea bottom (Figure 1a).The complex spatial velocity potential φ with respect to the local polar coordinate system r q , θ q , z of each cylinder q should satisfy the boundary value problem consisting of the following set of equations: where ∇ 2 denotes the 3D Laplace operator and K = ω 2 /g, with g denoting the gravitational acceleration.Equations ( 8) and ( 12) are the Laplace equations, Equation ( 9) corresponds to the linearized boundary condition on the undisturbed free-surface, while Equations ( 10) and ( 13) are the bottom boundary conditions.Moreover, Equation (11) is the Neumann body condition on the wetted area of the q th cylinder and Equation ( 14) expresses the zero velocity condition on the bottom of each cylinder.As mentioned above, for each cylinder two liquid regions were defined.Denoting by φ A and φ B , the discrete potentials of the wave fields in regions A and B, respectively, it is straightforward that φ A and φ B should satisfy Equations ( 8)- (11) and Equations ( 12)- (14).Accordingly, the uniqueness of the solution requires matching the potentials φ A and φ B in the mutual surface boundary r q = b, 0 ≤ θ q ≤ 2π and 0 ≤ z ≤ h 1 .The associated continuity relations read: Equation ( 15) expresses the continuity of the velocities, while Equation ( 16) expresses the continuity of the linear pressures (or potentials).
The velocity potential φ A consists of the incident wave component φ I and the scattered components originating from all cylinders in the array.By utilizing the linear superposition principle, the total velocity potential in the outer region A is given as: where the sum term on the right-hand side describes the scattered wave field induced by all cylinders in the array, except cylinder q, and it should be regarded as an incident wave field with complicated structure that hits the cylinder q.It is noted that Equation ( 17) combines all local polar coordinate systems r j , θ j , j = 1, . . ., Q.
Assuming a local Cartesian coordinate system x q , y q fixed on the symmetrical axis of the q th cylinder and a global Cartesian coordinate system (X, Y) fixed on the sea bottom, the incident wave potential φ I with respect to x q , y q is expressed as: cosh(k 0 h) Λ q e ik 0 (x q cos β+y q sin β) where k 0 is the wavenumber, J m is the Bessel function of the first kind with integer order m, i = √ −1, and Λ q = e ik 0 (X q cos β+Y q sin β) , where X q and Y q are the X and Y coordinates of the origin of x q , y q in the global coordinate system.In the following, all velocity potentials are considered to be normalized by −igA/ω.
Continuing with the scattered components of the total velocity potential (Equation ( 17)), these quantities should satisfy the Sommerfeld far-field radiation condition for outgoing waves at infinity.Accordingly, the scattered wave component φ (q) S associated with the q th cylinder is given by: In Equation (19), , where I m and K m denote the modified Bessel functions of the first and the second kind, respectively, while the prime denotes differentiation with respect to the argument.The eigenvalues a n are obtained from the solution of the well-known transcendental equation, which reads: The eigenfunctions Z n (z) that are associated with a n and are orthogonal in z ∈ [0, h] (with the normalization constant equal to h) are given by: Finally, in Equation (19), C mn , q = 1, . . ., Q, are the unknown expansion coefficients.It could be noted that Equation (19) unifies the solution of the scattered wave component by incorporating both the wavenumber and the evanescent modes, assuming a 0 = −ik 0 .Moreover, Equation (19) satisfies Equations ( 8)- (10), as well as the Sommerfeld radiation condition, since the second independent radial solution depending on I m a n r q was omitted.

Coming back to the unknown expansion coefficients of Equation (19), C (q)
mn , q = 1, . . ., Q, these quantities are obtained by employing the Neumann condition described by Equation (11), and the continuity relations given by Equations ( 15) and ( 16).This, however, requires the expression of the scattered components for j q, in the local coordinate system r q , θ q of the randomly selected cylinder q in the array, rather than in the local system r j , θ j , j q.For this purpose, Graf's addition theorem is employed for the modified Bessel function of the second kind and for the Hankel function of the first kind, in Equation ( 19) (more details can be found in [10,25]).
Combining all of the above, the wave field outside the random cylinder q is given by the following extended form of Equation ( 17): where R jq is the radial distance between the centers of cylinders j and q and β jq is the angle formed between the distance R jq and the horizontal in the center of the j th cylinder.
Continuing with the velocity potential in the lower region B, a proper solution of φ B for the examined boundary value problem has the form: where mn , q = 1, . . ., Q, denote the unknown expansion coefficients that aer calculated together with C (q) mn , q = 1, . . ., Q, by employing Equation ( 11) (Newman condition) and Equations ( 15) and ( 16) (continuity relations).
The aforementioned expression of φ B clearly satisfies the Laplace equation (Equation ( 12)), while the satisfaction of the boundary conditions given by Equations ( 13) and ( 14) is ensured by the eigenfunctions cos(nπz/h 1 ).As also required, φ B is bounded in the origin r q = 0 (the second independent radial solution depending on K m , which is singular for zero argument was omitted), while the expression of φ B was constructed to yield finite values for both zero and very large n values.Otherwise, Equation ( 23) is always bounded, considering that in the lower liquid region r q ≤ b.
For calculating the unknown expansion coefficients F (q) mn , and C (q) mn , q = 1, . . ., Q, of Equations ( 22) and ( 23), the matched eigenfunctions technique was utilized [10].The resulting system of equations is linear and it can be solved by employing standard matrix techniques.Moreover, a finite number of harmonics m = −M, . . ., M and of vertical eigenfunctions n = 0, . . ., N for both velocity potential expressions is taken into account, depending on the accuracy of the convergence required.A small number of modes and eigenfunctions (less than five, each) is usually appropriate to obtain convergent results.This in turn ensures the realization of numerical computations in practically "no time".
Having solved the diffraction problem, the first order exciting forces are obtained by direct pressure integration on the wetted surface of the cylinders, as follows: where Y , and Z are the surge, sway, and heave exciting forces, respectively, acting on the q th cylinder of the array and ρ is the water density.Equations ( 24) and ( 25) result in robust analytical formulas.Full details of these formulas can be found in [10,25].
The amplitude of the total heave exciting force F tot Z (objective function in Equation ( 2)), applied to all cylinders of the array is calculated as follows: where Re F (q) Z and Im F (q) Z denote the real and imaginary parts of Z , respectively.On the other hand, the amplitude of the heave exciting force F subtot Z (objective function in Equation ( 3)) applied on a pair of cylinders q and j, with q j is calculated as: where Re F Z and Im F Z denote the real and imaginary parts of Z , respectively, of the j th cylinder, with j q.

Genetic Algorithms Solver
The GAS applied in the present paper corresponds to a search evolutionary numerical method that is based on the process of natural evolution [14], where the fittest individuals of an initial population representing possible random solutions are selected for reproduction, in order to generate offspring of the next population.An individual of a population corresponds to a specific set of values of the design variables, X q , q = 1, . . ., Q, and it is considered to be a candidate solution of the optimization problem.The size of the population (equal to the number of individuals) is defined at the beginning of the GAS operators and remains constant during GAS execution.
Iteratively and at each step of the optimization process, the GAS scores each individual of the current population by computing its fitness value, with the use of the relevant objective function (Equation (26) or Equation ( 27)).Moreover, the GAS considers the individuals of the current population to create the individuals of the next one by employing specific GAS operators.Constraints and bounds for the design variables are also imposed (Equations ( 4)-( 7)) and needed to be satisfied when creating individuals.Part of the individuals that have better scores are selected as the "elite" and are maintained within the next population without any change, while some individuals are selected as "parents" for the generation of the new individuals (selection operator of GAS); in this way the fit solutions are kept in the next population.The "parents" are used to generate new individuals by-(a) making random changes to a single "parent" in order to maintain diversity within the population and prevent premature convergence (mutation operator of GAS) or (b) combining the design variables entries of a pair of "parents" to create fitter offspring (crossover operator of GAS).The same process continues iteratively until one of the predefined stopping criteria is met.Accordingly, the GAS is terminated and the individual with the best fitness value among the individuals of the last generated population presents the optimum solution.In the present paper, GAS is numerically realized by using the Optimization Toolbox™ R2019b [26] of MATLAB (R2019b, Natick, MA, USA) [27].

Optimization Numerical Process (ONP)
The HM and GAS numerical components described in the previous sections are coupled within MATLAB R2019b [27] in order to solve the optimization problems that are mathematically expressed by Equations ( 2)-( 7).This coupling is implemented by developing and applying a suitable ONP in the present paper, as shown in Figure 2. mutation, and crossover).The individuals of the new population are used as input to the HM in order to calculate new values of   () ,  = 1, … , , and accordingly new values of    or    .The individuals of the new population are scored and the stopping criteria are checked again.The aforementioned process continues iteratively until the stopping criteria are satisfied.When satisfaction of these criteria is achieved, the individual of the last generated population with the best score is assigned as the optimum solution of the constrained maximization problem examined.

ONP Validation
In order to validate the accuracy and the efficient development of the ONP, this process is initially applied in order to compare results with those of Kagemoto [11].The latter author utilized the inverse hydrodynamic interaction theory and developed analytical expressions for the objective function, to determine the optimum arrangement of cylindrical bodies in terms of minimizing the total surge or the heave exciting force applied on arrays of cylinders.
Herein, we compare the optimum solutions that minimize the total surge exciting force, for two linear arrays of truncated cylinders with, respectively  = 3 and  = 4, identical cylinders of radius  = 0.25 m and of draft ℎ 2 = 0.5 m.The length of the array,   , is taken to be equal to 2.0 m and 3.0 m, respectively, with   = 0.0 m (Figure 1a).Both arrays are placed in an area of water, depth ℎ = 1.0 m, and are subjected to the action of head monochromatic waves ( = 0 o in Figure 1a) of  = 5.54 rad/s.In the work of Kagemoto [11], the external cylinders of the array were considered to be fixed at the two ends of   .Accordingly, for each array, the optimum solutions correspond to the values of the  coordinates of the internal cylinders in the array.Based on the above, the optimization problem for the array with  = 3 is formed as follows: subjected to: The objective functions 26) and ( 27) respectively), the design variables X q , q = 1, . . ., Q, the constraints (Equations ( 4)-( 7)), the size of the populations, and the definition of GAS operators are provided as an input to the ONP.A random initial population is then created by the GAS for the design variables.The individuals of the initial population are used as input to the HM, along with the values of the remaining required physical quantities, i.e., radius, b, draft, h 2 , water depth, h, wave frequency, ω, and incident wave direction, β (Figure 1).HM is executed for each individual, F (q) Z (Equation ( 25)) for each CYLq, q = 1, . . ., Q, is obtained and the value of the objective function, F tot Z or F subtot Z , is calculated by utilizing Equation (26) or Equation (27), respectively.Then, the fitness of each individual to the objective function is evaluated, each individual is accordingly scored and the stopping criteria related to those scores are examined.If the stopping criteria are not satisfied, a new population of individuals is generated by the GAS with the use of the relevant GAS operators mentioned in the previous section (i.e., selection, mutation, and crossover).The individuals of the new population are used as input to the HM in order to calculate new values of F (q) Z , q = 1, . . ., Q, and accordingly new values of F tot Z or F subtot Z .The individuals of the new population are scored and the stopping criteria are checked again.The aforementioned process continues iteratively until the stopping criteria are satisfied.When satisfaction of these criteria is achieved, the individual of the last generated population with the best score is assigned as the optimum solution of the constrained maximization problem examined.

ONP Validation
In order to validate the accuracy and the efficient development of the ONP, this process is initially applied in order to compare results with those of Kagemoto [11].The latter author utilized the inverse hydrodynamic interaction theory and developed analytical expressions for the objective function, to determine the optimum arrangement of cylindrical bodies in terms of minimizing the total surge or the heave exciting force applied on arrays of cylinders.
Herein, we compare the optimum solutions that minimize the total surge exciting force, for two linear arrays of truncated cylinders with, respectively Q = 3 and Q = 4, identical cylinders of radius b = 0.25 m and of draft h 2 = 0.5 m.The length of the array, l array , is taken to be equal to 2.0 m and 3.0 m, respectively, with l out = 0.0 m (Figure 1a).Both arrays are placed in an area of water, depth h = 1.0 m, and are subjected to the action of head monochromatic waves (β = 0 o in Figure 1a) of ω = 5.54 rad/s.In the work of Kagemoto [11], the external cylinders of the array were considered to be fixed at the two ends of l array .Accordingly, for each array, the optimum solutions correspond to the values of the X coordinates of the internal cylinders in the array.Based on the above, the optimization problem for the array with Q = 3 is formed as follows: minimize F tot X (X 2 ), subjected to: X q − X q−1 ≥ 0.5, q = 2 and 3, On the other hand, the optimization problem for the array with Q = 4 is formed as follows: minimize subjected to: In Equations ( 28) and (30), F tot X corresponds to the amplitude of the total surge exciting force applied on all cylinders of the array, which is calculated as follows: where Re F (q) X and Im F (q) X denote the real and imaginary parts of F (q) X (Equation ( 24)), respectively.Table 1 shows the comparison of the results obtained using the present ONP with the corresponding ones of [11].It is evident that both the calculated optimum solutions and the minima of F tot X (normalized in terms of ρgπbh 2 AQ as in [11]) agree favorably with those in [11], which proves the efficiency and the accuracy of the developed ONP in the present paper.The observed very small differences are attributed to the fact that ONP includes a search evolutionary method for solving the optimization problem, while in [11] the minimum of the objective function was derived analytically (existence of an explicit mathematical formula).
Table 1.Comparison of the optimum solutions between the present ONP and Kagemoto [11].

Examined Cases
The ONP developed in this study is applied for solving the optimization problems described with Equations ( 2)-( 7) along with Equations ( 26) and ( 27), for a linear array of Q = 9 identical truncated cylinders (Figure 1a with Q = 9) of b = 1.0 m and non-dimensional draft h 2 /b = 1.0, placed in a liquid area of non-dimensional constant depth h/b = 10.0 (Figure 1b).In total, 18 optimization cases are considered and solved (Table 2).Cases EC1.1-EC1.8correspond to the optimum layout configurations of the array that maximize F tot Z (i.e., the total heave exciting forces acting on all the cylinders of the array) for different combinations of l array /b and l out /b values, under the action of either head or perpendicular to the array waves (β = 0 o and 90 o , respectively, in Figure 1a).On the other hand, EC2.1-EC2.10correspond to optimum layout configurations of the array that maximize F subtot Z , namely, the heave exciting forces applied on appropriately selected pairs of cylinders within the array, as well as in the middle cylinder of the array (in the latter case F subtot Z ≡ F (5) Z ).For each different objective function F subtot Z , two different combinations of l array /b and l out /b are considered, along with the action of the perpendicular to the array waves.For all cases examined, the array is subjected to monochromatic waves of ω = 2.525 rad/s, which coincides with the heave natural frequency of a cylinder within the array at h/b = 10.0 [28], while during the implementation of the GAS the non-dimensional design variables, X q /b, q = 1, . . ., 9, are defined to have values up to their first decimal.Moreover, for the cases with l out /b = 0, the application of Equations ( 5) and ( 6) leads to constant values of X 1 /b and X 9 /b equal to 0 and l array /b, respectively; accordingly, the external cylinders of the array, CYL1 and CYL9, are assumed to be fixed at the two ends of l array .
For applying the ONP, the following options were assigned to the GAS-(a) the population size was defined equal to 1000, while the double vector was selected as the population type, (b) a rank fitness scaling was applied, (c) a crossover fraction was taken to be equal to 0.8 with the intermediate function, and (d) an adaptive feasible mutation was selected.In order to ensure that the GAS predicts the global maximum of the physical problem correctly, multiple different runs for each case of Table 2 were made.It is also noted that for each examined case, a different number of GAS operator iterations was required in order to achieve convergence of the optimization algorithm.In general, the cases with β = 90 o required a smaller number of GAS operator iterations (∼ 150), compared to the cases with β = 0 o , where ∼ 430 iterations were required to obtain the optimum solution.The average computational time for one iteration (1000 runs of the HM) was approximately equal to 3 min, using a standard PC (Personal Computer) with 128 GB RAM (Random Access Memory) and Intel ® Xeon ® Silver 4110 CPU@ 2.1 GHz 2.1 GHz (2 processors).

Optimum Solutions for EC1.1-EC1.8
Table 3 shows the results (optimum values of the design variables X q , q = 1, . . ., 9, normalized in terms of b, and maximum values of F tot Z , normalized in terms of ρgb 2 A) for the cases EC1.1-EC1.8 of Table 2.For these cases the optimum layout configuration of the array in the X − Y plane is shown schematically in Figure 3, while, Figures 4 and 5 show the effect of l array /b on the maximum values of F tot Z and on the optimum X q /b values, respectively, for a given β direction.
The results in Table 3 and Figure 3 clearly demonstrate that irrespective of the examined l array /b and β combinations, the optimum solutions correspond to a random positioning of the cylinders within the array, characterized by the existence of unequal center-to-center distances between the adjacent cylinders and of non-symmetrical features.The only distinctive characteristic identified is the formation of clusters (sub-groups) of one, two, or more (up to six) closely-positioned cylinders, depending on l array /b and β.These clusters are distributed at various distances between each other, along the length of the array.Under the action of perpendicular to the array waves (EC1.3-EC1.4,EC1.7-EC1.8), a smaller number of clusters is realized, compared to the corresponding cases with β = 0 o , while, for all cases, except that of EC1.3, the clusters tend to be situated close to the ends of the examined l array .The latter feature is not observed for β = 0 o , where the clusters are distributed along the whole available length of the array.From a physical point of view, all of the above illustrate that the closely-positioned cylinders within the clusters induce intense hydrodynamic interaction effects, which are positive in terms of maximizing F tot Z .
Table 3. Optimum values of X q /b, q = 1, . . ., 9 and corresponding maximum F tot Z / ρgb 2 A values for EC1.1-EC1.8. the formation of clusters (sub-groups) of one, two, or more (up to six) closely-positioned cylinders, depending on   / and .These clusters are distributed at various distances between each other, along the length of the array.Under the action of perpendicular to the array waves (EC1.3-EC1.4,EC1.7-EC1.8), a smaller number of clusters is realized, compared to the corresponding cases with  = 0 o , while, for all cases, except that of EC1.3, the clusters tend to be situated close to the ends of the examined   .The latter feature is not observed for  = 0 o , where the clusters are distributed along the whole available length of the array.From a physical point of view, all of the above illustrate that the closely-positioned cylinders within the clusters induce intense hydrodynamic interaction effects, which are positive in terms of maximizing    .Regarding the maximum values of    (Table 3, Figure 4), it is evident that for a given   / value, the maximum total heave exciting force is significantly smaller for the cases EC1.1-EC1.2 and EC1.5-EC1.6, where  = 0 o (Figure 4a) compared to the cases EC1.3-EC1.4,EC1.7-EC1.8(Figure 4b), respectively, where  = 90 o .This could be attributed to the fact that, for head waves, the incident wave energy is successively reduced towards the leeward section of the array.Moreover, for  = 0 o , the successive increase of   / results in a smooth increase of    maxima, since larger values of   / enable the reduction of shadow effects between the clusters of the cylinders.This is not observed for  = 90 o , where the change of   / does not have any significant impact on the Regarding the maximum values of F tot Z (Table 3, Figure 4), it is evident that for a given l array /b value, the maximum total heave exciting force is significantly smaller for the cases EC1.1-EC1.2 and EC1.5-EC1.6,where β = 0 o (Figure 4a) compared to the cases EC1.3-EC1.4,EC1.7-EC1.8(Figure 4b), respectively, where β = 90 o .This could be attributed to the fact that, for head waves, the incident wave energy is successively reduced towards the leeward section of the array.Moreover, for β = 0 o , the successive increase of l array /b results in a smooth increase of F tot Z maxima, since larger values of l array /b enable the reduction of shadow effects between the clusters of the cylinders.This is not observed for β = 90 o , where the change of l array /b does not have any significant impact on the maximum values of F tot Z .

Case
Energies 2020, 13, x FOR PEER REVIEW 13 of 19 Continuing with the effect of   / on   ,  = 1, … , , (Figure 5), it is evident that for a given P th cylinder, the successive increase of   / in the case of  = 0 o leads to larger values of   /,  = 1, … ,9, especially for  ≥ 5.For  = 90 o , this trend is only observed for the last three cylinders of the array, namely, for  ≥ 7.Moreover, under the action of head waves a linear correlation between  and   /,  = 1, … , 9, for each examined   / value is observed.This is also shown in Figure 6, where   /  ,  = 1, … , 9, is plotted as a function of .By implementing linear fitting to the aforementioned data, a linear trend line along with its equation is obtained and it is additionally included in Figure 6.It is evident that the formula shown in Figure 6 does not correspond to a direct of the solution of the examined optimization problem (Equations ( 2) and ( 4)-( 7)), but to an approximation of the relevant optimum solutions obtained for  = 0 o , by utilizing the present ONP.In the case of  = 90 o (Figure 5b), a linear dependence of   / with  exists only for EC1.3 and EC1.7 (  /=40 and 32, respectively) restricting the use of any possible fitting scheme of the data related to  = 90 o .Continuing with the effect of   / on   ,  = 1, … , , (Figure 5), it is evident that for a given P th cylinder, the successive increase of   / in the case of  = 0 o leads to larger values of   /,  = 1, … ,9, especially for  ≥ 5.For  = 90 o , this trend is only observed for the last three cylinders of the array, namely, for  ≥ 7.Moreover, under the action of head waves a linear correlation between  and   /,  = 1, … , 9, for each examined   / value is observed.This is also shown in Figure 6, where   /  ,  = 1, … , 9, is plotted as a function of .By implementing linear fitting to the aforementioned data, a linear trend line along with its equation is obtained and it is additionally included in Figure 6.It is evident that the formula shown in Figure 6 does not correspond to a direct outcome of the solution of the examined optimization problem (Equations ( 2) and ( 4)-( 7)), but to an approximation of the relevant optimum solutions obtained for  = 0 o , by utilizing the present ONP.In the case of  = 90 o (Figure 5b), a linear dependence of   / with  exists only for EC1.3 and EC1.7 (  /=40 and 32, respectively) restricting the use of any possible fitting scheme of the data related to  = 90 o .Continuing with the effect of l array /b on X q , q = 1, . . ., Q, (Figure 5), it is evident that for a given q th cylinder, the successive increase of l array /b in the case of β = 0 o leads to larger values of X q /b, q = 1, . . ., 9, especially for q ≥ 5.For β = 90 o , this trend is only observed for the last three cylinders of the array, namely, for q ≥ 7.Moreover, under the action of head waves a linear correlation between q and X q /b, q = 1, . . ., 9, for each examined l array /b value is observed.This is also shown in Figure 6, where X q /l array , q = 1, . . ., 9, is plotted as a function of q.By implementing linear fitting to the aforementioned data, a linear trend line along with its equation is obtained and it is additionally included in Figure 6.It is evident that the formula shown in Figure 6 does not correspond to a direct outcome of the solution of the examined optimization problem (Equations ( 2), ( 4)-( 7)), but to an approximation of the relevant optimum solutions obtained for β = 0 o , by utilizing the present ONP.In the case of β = 90 o (Figure 5b), a linear dependence of X q /b with q exists only for EC1.3 and EC1.7 (l array /b = 40 and 32, respectively) restricting the use of any possible fitting scheme of the data related to β = 90 o .Figure 6.  /  as a function of  for EC1.1-EC1.2 and EC1.5-EC1.6.

Optimum Solutions for EC2.1-EC2.10
Table 4 shows the results (optimum values of   /,  = 1, … , 9 , and maximum values of    /( 2 )) for the cases EC2.1-EC2.10 of Table 2, where only the action of waves with  = 90 o is considered.Moreover, Figure 7 schematically shows the corresponding optimum layout configuration of the array in the  −  plane.For all examined   / values, the optimum solutions correspond again to random placements of the cylinders within the array, with unequal center-to-center distances between adjacent bodies.The formation of clusters of closely-positioned cylinders, distributed along the length of the array is also observed.In general, the maximization of    for a given pair of cylinders P th and P th is realized by placing one or more adjacent bodies at small distances from the P th and/or the P th cylinder.It is also worth noting that for both examined   / values, the largest maxima of    occur for EC2.9 and EC2.10 (Table 4), where maximization of the total heave exciting force applied on CYL4 and CYL6 of the array is implemented.For the rest of the examined pairs of cylinders, the    maxima successively decrease and the smallest values are observed for cases EC2.3 and EC2.4 (maximization of total heave exciting force applied on the outer cylinders, CYL1 and CYL9, of the array).X q /l array q =0.121374 -0.09308 Figure 6.X q /l array as a function of q for EC1.1-EC1.2 and EC1.5-EC1.6.

Optimum Solutions for EC2.1-EC2.10
Table 4 shows the results (optimum values of X q /b, q = 1, . . ., 9, and maximum values of F subtot Z / ρgb 2 A ) for the cases EC2.1-EC2.10 of Table 2, where only the action of waves with β = 90 o is considered.Moreover, Figure 7 schematically shows the corresponding optimum layout configuration of the array in the X − Y plane.For all examined l array /b values, the optimum solutions correspond again to random placements of the cylinders within the array, with unequal center-to-center distances between adjacent bodies.The formation of clusters of closely-positioned cylinders, distributed along the length of the array is also observed.In general, the maximization of F subtot Z for a given pair of cylinders q th and j th is realized by placing one or more adjacent bodies at small distances from the q th and/or the j th cylinder.It is also worth noting that for both examined l array /b values, the largest maxima of F subtot Z occur for EC2.9 and EC2.10 (Table 4), where maximization of the total heave exciting force applied on CYL4 and CYL6 of the array is implemented.For the rest of the examined pairs of cylinders, the F subtot Z maxima successively decrease and the smallest values are observed for cases EC2.3 and EC2.4 (maximization of total heave exciting force applied on the outer cylinders, CYL1 and CYL9, of the array).Table 4. Optimum values of X q /b, q = 1, . . ., 9 and the corresponding maximum  Figure 8 shows the variation of the optimum values of   /,  = 1, … , 9, for   / = 40 (Figure 8a) and   / = 60 (Figure 8b), as a function of .By comparing the same symbols of the two subfigures, it can be concluded that for the middle cylinder CYL5 or for a given pair of P th and P th cylinders, the consideration of a larger   / for maximizing the corresponding    , leads to a more linear relationship between   / and q.For a better insight on the heave exciting forces,   () ,  = 1, … , 9, and    for the optimum solutions of cases EC2.1-EC2.10 are presented in Table 5.For each case, the values of   () written in bold are related to the cylinders that have been considered to form the relevant    .Accordingly, they correspond to the largest values of the heave exciting forces among all cylinders, within the array.These values are also plotted in Figure 9a.By comparing the pair of cases, where the same objective function was considered (e.g., cases EC2.1 and EC2.2, cases EC2.3 and EC2.4,etc.), it is evident that the change of   / has a minor effect on the largest   () ,  = 1, … , 9, values (an increase or decrease of up to 2%), except for cases EC2.7 and EC2.8, where the increase of   / leads to a 7% increase of   (3) .With regards to    (Table 5 and Figure 9b), as resulted from the optimum layout configurations of cases EC2.1-EC2.10, the largest total heave exciting force for   / =40 is observed for EC2.5, where maximization of    applied on CYL2 and CYL8 of the array is performed.On the other hand, for   /=60 the case of EC2.8 (maximization of    applied on CYL3 and CYL7 of the array) leads to the largest value of    .Nevertheless, the aforementioned    values are smaller, compared to those obtained for the relevant cases EC1.3 and EC1.4 (Table 3, Figure 4b).Figure 8 shows the variation of the optimum values of X q /b, q = 1, . . ., 9, for l array /b = 40 (Figure 8a) and l array /b = 60 (Figure 8b), as a function of q.By comparing the same symbols of the two subfigures, it can be concluded that for the middle cylinder CYL5 or for a given pair of q th and j th cylinders, the consideration of a larger l array /b for maximizing the corresponding F subtot Z , leads to a more linear relationship between X q /b and q.
Energies 2020, 13, x FOR PEER REVIEW 15 of 19 Figure 8 shows the variation of the optimum values of   /,  = 1, … , 9, for   / = 40 (Figure 8a) and   / = 60 (Figure 8b), as a function of .By comparing the same symbols of the two subfigures, it can be concluded that for the middle cylinder CYL5 or for a given pair of P th and P th cylinders, the consideration of a larger   / for maximizing the corresponding    , leads to a more linear relationship between   / and q.For a better insight on the heave exciting forces,   () ,  = 1, … , 9, and    for the optimum solutions of cases EC2.1-EC2.10 are presented in Table 5.For each case, the values of   () written in bold are related to the cylinders that have been considered to form the relevant    .Accordingly, they correspond to the largest values of the heave exciting forces among all cylinders, within the array.These values are also plotted in Figure 9a.By comparing the pair of cases, where the same objective function was considered (e.g., cases EC2.1 and EC2.2, cases EC2.3 and EC2.4,etc.), it is evident that the change of   / has a minor effect on the largest   () ,  = 1, … , 9, values (an increase or decrease of up to 2%), except for cases EC2.7 and EC2.8, where the increase of   / leads to a 7% increase of   (3) .With regards to    (Table 5 and Figure 9b), as resulted from the optimum layout configurations of cases EC2.1-EC2.10, the largest total heave exciting force for   / =40 is observed for EC2.5, where maximization of    applied on CYL2 and CYL8 of the array is performed.On the other hand, for   /=60 the case of EC2.8 (maximization of    applied on CYL3 and CYL7 of the array) leads to the largest value of    .Nevertheless, the aforementioned    values are smaller, compared to those obtained for the relevant cases EC1.3 and EC1.4 (Table 3, Figure 4b).For a better insight on the heave exciting forces, F Z , q = 1, . . ., 9, and F tot Z for the optimum solutions of cases EC2.1-EC2.10 are presented in Table 5.For each case, the values of F (q) Z written in bold are related to the cylinders that have been considered to form the relevant F subtot Z .Accordingly, they correspond to the largest values of the heave exciting forces among all cylinders, within the array.These values are also plotted in Figure 9a.By comparing the pair of cases, where the same objective function was considered (e.g., cases EC2.1 and EC2.2, cases EC2.3 and EC2.4,etc.), it is evident that the change of l array /b has a minor effect on the largest F (q) Z , q = 1, . . ., 9, values (an increase or decrease of up to 2%), except for cases EC2.7 and EC2.8, where the increase of l array /b leads to a 7% increase of Z .With regards to F tot Z (Table 5 and Figure 9b), as resulted from the optimum layout configurations of cases EC2.1-EC2.10, the largest total heave exciting force for l array /b = 40 is observed for EC2.5, where maximization of F subtot Z applied on CYL2 and CYL8 of the array is performed.On the other hand, for l array /b = 60 the case of EC2.8 (maximization of F subtot Z applied on CYL3 and CYL7 of the array) leads to the largest value of F tot Z .Nevertheless, the aforementioned F tot Z values are smaller, compared to those obtained for the relevant cases EC1.3 and EC1.4 (Table 3, Figure 4b).

Table 5. F (q)
Z and F tot Z for the optimum solutions of EC2.1-EC2.10(values given are non-dimensional in terms of ρgb 2 A).

Case
F Z F Z F Z F Z F Z F Z F Z F Z F   value occurs when    is maximized for the pair of cylinders that includes the specific P th body.For the outer cylinders of the array ( = 1 and 9), the second largest value occurs for EC2.5 and EC2.6, where    is maximized for the pair of the two adjacent cylinders (i.e.,  = 2 and 8).Same conclusions can be derived for CYL2 and CYL8, where the second largest value occurs for EC2.3 and EC2.4 (maximization of    for the pair of the two outer cylinders, i.e.,  = 1 and 9).For Z , q = 1, . . ., 9 , for each case; and (b) F tot Z as resulted from each optimum layout configuration.
Finally, Figure 10 shows the values of F Z , q = 1, . . ., 9, applied on each cylinder of the array for each of the optimum layout configuration of cases EC2.1-EC2.10.For a q th cylinder, the largest Z value occurs when F subtot Z is maximized for the pair of cylinders that includes the specific q th body.For the outer cylinders of the array (q = 1 and 9), the second largest value occurs for EC2.5 and EC2.6, where F subtot Z is maximized for the pair of the two adjacent cylinders (i.e., q = 2 and 8).Same conclusions can be derived for CYL2 and CYL8, where the second largest value occurs for EC2.3 and EC2.4 (maximization of F subtot Z for the pair of the two outer cylinders, i.e., q = 1 and 9).For the rest of the cylinders, a clear trend can not be observed.However, it is evident that the consideration of different pairs of cylinders for maximizing the corresponding F subtot Z has a direct impact on the values of F (q) Z for each q th cylinder of the array.the rest of the cylinders, a clear trend can not be observed.However, it is evident that the consideration of different pairs of cylinders for maximizing the corresponding    has a direct impact on the values of   () for each P th cylinder of the array.

Conclusions
In this paper, we developed and applied an ONP for determining the optimum layout configuration of a linear array of identical mutually interacting truncated cylinders.The optimum configurations maximize either the total heave exciting force,    , acting on all cylinders of the array or the heave exciting force,    , acting on appropriately selected pairs of cylinders within the array.ONP consists of two coupled distinctive numerical components capable of solving the diffraction problem of the examined multi-body arrangement in the frequency domain and the examined optimization problem, by utilizing GAs.The efficiency and the accuracy of the proposed ONP was demonstrated by the very good agreement of results with those of Kagemoto [11].The developed ONP was applied for the case of a linear array consisting of nine identical truncated cylinders.In total, 18 optimization cases were considered and solved under the action of either head or perpendicular to the array waves, and for different maximum allowable lengths of the array.
The results illustrated that irrespectively of the objective function definition, the examined length of the array and the incident wave direction, optimum solutions corresponded to a quite random positioning of the cylinders within the array with unequal center-to-center distances between adjacent bodies.However, a common feature observed for all optimum layout configurations was the formation of clusters (sub-groups) of closely-positioned adjacent cylinders.Under the action of head waves, the clusters were distributed along the whole available length of the array, while for perpendicular to the array waves, a smaller number of clusters was realized, which tended to be situated close to the ends of the examined length of the array.
For the optimization cases, where the maximization of    was sought, the maximum values of    under head waves (  = 0 o ) were significantly smaller, compared to those obtained for perpendicular to the array waves.Moreover, for  = 0 o , the increase of the array's length resulted in a smooth increase of the total heave exciting force maxima, while, for all examined lengths of the array, a linear correlation was observed between the cylinders' number and their optimum locations within the array.Regarding the optimization cases, where maximization of    was sought, the consideration of the pair of cylinders situated on the two sides of the middle cylinder in the objective function led to the largest maximum value of    .The optimum layout configurations of the aforementioned cases also resulted in    values that were smaller, compared to those obtained by the direct maximization of the total heave exciting force.Finally, the consideration of different pairs of cylinders for maximizing the corresponding    had a direct impact on the values of   () for each P th cylinder of the array.

Conclusions
In this paper, we developed and applied an ONP for determining the optimum layout configuration of a linear array of identical mutually interacting truncated cylinders.The optimum configurations maximize either the total heave exciting force, F tot Z , acting on all cylinders of the array or the heave exciting force, F subtot Z , acting on appropriately selected pairs of cylinders within the array.ONP consists of two coupled distinctive numerical components capable of solving the diffraction problem of the examined multi-body arrangement in the frequency domain and the examined optimization problem, by utilizing GAs.The efficiency and the accuracy of the proposed ONP was demonstrated by the very good agreement of results with those of Kagemoto [11].The developed ONP was applied for the case of a linear array consisting of nine identical truncated cylinders.In total, 18 optimization cases were considered and solved under the action of either head or perpendicular to the array waves, and for different maximum allowable lengths of the array.
The results illustrated that irrespectively of the objective function definition, the examined length of the array and the incident wave direction, optimum solutions corresponded to a quite random positioning of the cylinders within the array with unequal center-to-center distances between adjacent bodies.However, a common feature observed for all optimum layout configurations was the formation of clusters (sub-groups) of closely-positioned adjacent cylinders.Under the action of head waves, the clusters were distributed along the whole available length of the array, while for perpendicular to the array waves, a smaller number of clusters was realized, which tended to be situated close to the ends of the examined length of the array.
For the optimization cases, where the maximization of F tot Z was sought, the maximum values of F tot Z under head waves (β = 0 o ) were significantly smaller, compared to those obtained for perpendicular to the array waves.Moreover, for β = 0 o , the increase of the array's length resulted in a smooth increase of the total heave exciting force maxima, while, for all examined lengths of the array, a linear correlation was observed between the cylinders' number and their optimum locations within the array.Regarding the optimization cases, where maximization of F subtot Z was sought, the consideration of the pair of cylinders situated on the two sides of the middle cylinder in the objective function led to the

Figure 1 .
Figure 1.Geometry of the examined linear arrangement and definition of basic quantities: (a)  −  plane and (b)  −  plane.

h 2 Figure 1 .
Figure 1.Geometry of the examined linear arrangement and definition of basic quantities: (a) X − Y plane and (b) Y − Z plane.

Figure 2 .
Figure 2. Flow chart of the developed optimization numerical process (ONP) in the present paper for calculating the optimum arrangement of cylinders in a linear array, in terms of maximizing the heave exciting force.

Optimization Numerical Process (ONP)
Figure 2. Flow chart of the developed optimization numerical process (ONP) in the present paper for calculating the optimum arrangement of cylinders in a linear array, in terms of maximizing the heave exciting force.

Table 2 .
Details of the examined optimization cases.