Surrogate modeling of hydrodynamic forces between multiple floating bodies through a hierarchical interaction decomposition

Journal of Computational Physics 408 (2020) 109298

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Surrogate modeling of hydrodynamic forces between multiple floating bodies through a hierarchical interaction decomposition Jize Zhang a , Alexandros A. Taflanidis a,b,∗ , Jeffrey T. Scruggs c a b c

Department of Civil and Environmental Engineering and Earth Sciences, University of Notre Dame, United States of America Department of Aerospace and Mechanical Engineering, University of Notre Dame, United States of America Department of Civil Engineering, Department of Electrical Engineering, University of Michigan at Ann Arbor, United States of America

a r t i c l e

i n f o

Article history: Received 22 June 2019 Received in revised form 1 December 2019 Accepted 27 January 2020 Available online xxxx Keywords: Hydrodynamic coefficients Multi-body interaction Wave-energy converters Surrogate modeling Many-body expansion Hierarchical interaction decomposition

a b s t r a c t Efficient estimation of wave interactions between multiple floating bodies (diffraction forces and radiation coefficients) is of importance to many applications, with a key one being the layout optimization of arrays of wave energy converter. The computational complexity for this estimation dramatically increases for configurations with large number of bodies. To address this challenge, a data-driven surrogate modeling implementation is discussed in this paper. The foundation of the approach is an innovative application of the many-body expansion principle that overcomes the curse of dimensionality for the surrogate model development. Instead of using a single surrogate model to predict the hydrodynamic characteristics of the multi-body configuration, multiple surrogate models corresponding to clusters with fewer bodies are employed. These lower-order surrogate models can be developed at a substantial smaller computational cost, especially for the first terms of the many-body expansion that contribute dominantly to the total hydrodynamic characteristics. Additional enhancements for the surrogate modeling implementation pertain to the characterization of the surrogate model input and output, exploiting symmetry and invariance principles for the hydrodynamic problem. Using Kriging as surrogate model, the approach is demonstrated for predicting the hydrodynamic characteristics for an array of vertical axisymmetric bodies, but it is extendable to similar multi-body interaction problems or other surrogate modeling techniques. Several numerical case studies are finally presented, showcasing the computational efficiency and prediction accuracy of the proposed method, as well as its scalability to arrays of large size, and the advantages it can offer within the context of layout optimization for wave energy converters. © 2020 Elsevier Inc. All rights reserved.

1. Introduction The estimation of the interaction characteristics among a large ensemble of bodies plays an important role in many physics applications, such as gravitational astrophysics [1], flow simulation [2], quantum mechanics [3], nuclear physics [4] and molecular dynamics [5]. Direct numerical simulation for such problems can be computationally demanding, stemming

*

Corresponding author. E-mail address: a.tafl[email protected] (A.A. Taflanidis).

https://doi.org/10.1016/j.jcp.2020.109298 0021-9991/© 2020 Elsevier Inc. All rights reserved.

2

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

from the need to account for the interactions between all bodies. This has motivated researchers to seek alternative approximate approaches that can reduce the overall computational burden [6–8]. Although the specifics of these approaches depend on the characteristics of the governing equations, their computational complexity in general increases with the number of bodies. This paper specifically discusses the wave interaction between multiple obstacles, a classical many-body interaction problem [9] with practical importance to various fields including acoustics [10], electromagnetism [11], seismology [12] and hydrodynamics [13]. In particular, it focuses on hydrodynamic interaction (diffraction forces and radiation coefficients) between identical floating bodies, though the proposed methodology is extendable to similar problems. The calculation of such type of hydrodynamic interactions is of interest to many offshore and coastal applications [14,15], and its relevance has dramatically increased in the past two decades due to the rising interest in the deployment of arrays of wave energy converters (WECs) [16,17]. For wave energy extraction applications, arrays of oscillating floating resonant bodies are widely acknowledged as one of the most promising implementations [18–20], since their deployment can reduce cost of energy production/transfer and can improve efficiency by increasing the effective capture width of wave crest, compared to the collective widths of the individual devices in the array. To fulfill this efficiency improvement, the array layout needs to be carefully selected [20] to positively exploit the inter-WEC dynamic coupling which arises as a consequence of fluid-structure interaction. This in turn has created a strong interest in the estimation of the hydrodynamic interactions among WECs. The typical modeling assumption for such applications is linear wave theory, whereas hydrodynamic interaction forces are estimated for harmonic motion [21], with the objective being the calculation of the characteristics of the frequency response function for the WEC array to support its design [22,23]. Computational solution techniques for the multi-body hydrodynamic interactions in this setting fall into two main categories: (i) approximate methods, such as the Point Absorber (PA) method [23] and the Plane Wave method [24,25], that invoke a simplification of the multi-body interaction and can provide a computationally efficient estimation of hydrodynamic characteristics albeit at a limited accuracy; and (ii) exact methods, that can in principle achieve any desired accuracy but at a substantial computational burden. Exact methods can be further categorized into semi-analytical methodologies based on multi-scattering (MS) approximation, which may involve either iterative [6] or non-iterative implementation [14, 26] (also referenced as direct-matrix method), and numerical methods such the Finite Element Method (FEM) [27] or the Boundary Element Method (BEM) [28]. While numerical approaches have arguably broader scope, the MS approach is often preferred because of the relatively higher computational efficiency as well as the valuable physical insights it can offer [21]. Additional discussions on relevant numerical techniques for WEC hydrodynamic interaction estimation can be found in [21]. All aforementioned exact solution techniques face computational challenges especially with a larger number of bodies and/or complex boundary domains. This cost bottleneck is particularly relevant for WEC array formation design, where repeated estimations are needed for different array layouts in order to identify an optimal configuration [26,29]. This realization has incentivized the adoption of advanced numerical schemes to improve the computational efficiency, such as the use of the fast multipole method [8,30] and the Block Toeplitz structure in matrix inversion [31], or of reduced-order modeling techniques that project the governing equations to a lower dimensional subspace [32,33]. An alternative approach for circumventing this challenge is to exploit some hierarchical approximation scheme, such as grouping of the larger array to smaller clusters [34,35], or ignoring certain interactions in the near [36] or far field [37], which essentially trade the solution accuracy for computational efficiency. This paper investigates a different path to reduce the computational cost. It adopts a surrogate modeling implementation that replaces the computationally intensive simulator of the hydrodynamic interaction characteristics with a fast-to-evaluate emulator, developed by regression or interpolation over a set of training data capturing the input/output relationship of the original, high-fidelity simulator [38]. After the surrogate model is constructed, it can facilitate computationally efficient predictions of the hydrodynamic interactions (outputs) for all desired inputs, supporting optimization or other relevant tasks, such as sensitivity analysis or uncertainty quantification, that would require multiple calls to the expensive simulator. To date, a limited number of papers have examined the use of surrogate modeling techniques in the context of multibody hydrodynamic interaction for WEC applications [29,39]. These existing approaches have adopted a specialized formulation, with surrogate models approximating the WEC power absorption performance instead of individual hydrodynamic interaction coefficients, even though approximating these coefficients and subsequently estimating the absorbed energy based on them is known to offer substantial benefits in terms of prediction accuracy [40]. More importantly, these efforts did not offer comprehensive solutions for controlling the surrogate model input dimensionality, a feature that creates substantial computational challenges when considering arrays with large number of bodies. The study in [39] was limited to very simple transformations considering symmetry of the interaction problem, offering ultimately only minor dimensional reduction, whereas a simplified heuristic dimensional reduction was adopted in [29], assuming that only the two nearest neighbors interact with each body, which evidently cannot accurately capture the interaction effects among the whole array. This work proposes a novel surrogate model implementation framework for predicting the hydrodynamic interaction characteristics for arrays of identical floating bodies. To accomplish a systematic dimensional reduction of the surrogate model input, a hierarchical interaction decomposition scheme is developed based on the many body expansion (MBE) principle [5]. The key idea is that the hydrodynamic coefficients considering the total interaction can be expressed as a series summation of contributions from the interactions of clusters consisting of an increasing number of bodies, starting with all single bodies, then examining all combinations of two bodies and up to the array consisting of all bodies. If this expansion converges rapidly, the exact solution should be accurately reconstructed even when truncating the summation after only a few terms. For the surrogate modeling framework, MBE can be used as a dimensional reduction technique, by developing

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

3

Fig. 1. Geometric characteristics in Cartesian coordinates (x- y-z system) for N body array hydrodynamic interactions.

surrogate models to predict the additive interaction effects of clusters with fewer bodies, instead of working directly with the original many-body system. Similar concept of combining surrogate modeling with MBE principles has been shown promising in the particular field of molecular dynamics simulations [41–43]. For modeling hydrodynamic characteristics of floating bodies, we further demonstrate that the additive interaction effects of a cluster, stemming explicitly from the addition of a new body compared to its sub-cluster(s), can be efficiently modeled through a proper characterization of the surrogate model input and output, exploiting symmetry and invariance principles [43–45]. The interaction order is truncated when the magnitude of the additive interaction effects becomes sufficiently small. Kriging (or Gaussian process regression) [46] is adopted as the surrogate model in this work, although the overall formulation extends to many other type of surrogate models. In the numerical study, we demonstrate that the MBE-based surrogate models can lead to orders of magnitude speed up with small error in predicting hydrodynamic interactions of floating bodies compared to conventional solvers. The surrogate assisted layout optimization of WEC arrays for optimal power absorption is also investigated. The remainder of this paper is organized as follows. In the next section, the underlying physical problem is introduced along with discussion of relevant solution techniques. In Section 3, Kriging is reviewed and the hierarchical interaction decomposition surrogate modeling implementation is presented in detail. Section 4 discusses relevant numerical examples whereas Section 5 demonstrates the implementation for WEC layout optimization. 2. Problem formulation: hydrodynamic interaction for multiple bodies Consider a group of N rigid, identical bodies floating in still water of depth d as shown in Fig. 1. In the surrogate model formulation discussed in the illustrative applications later, the assumption of vertical axisymmetric bodies will be further utilized, but the overall framework can be applied regardless of the body geometry. A global Cartesian coordinate system (x, y , z) is defined with the z-axis pointing upwards with z = 0 corresponding to the seabed. The array layout is fully characterized by the 2-by-N dimensional matrix W = [w1 w2 · · · w N ], whose pth column w p = [x p y p ] T is a 2-dimensional vector representing the center coordinate of the pth body in the xy-plane Cartesian coordinate system. All bodies are subject to progressive incident wave propagating along direction with angle β relative to the x-axis, with angular frequency ω and wave number k satisfying the dispersion equation ω2 = gktanh(kd), where g is the gravitational constant. For vertical axisymmetric bodies, the coordinate system can be rotated so that β = 0o , and therefore the incident wave angle does not need to be treated as a separate parameter in the formulation. Assuming small amplitude displacements and inviscid, incompressible and irrotational flow, the fluid motion can be described through the linearized wave theory with the time-dependent total velocity potential given by [9]:

  (x, y , z, t ) = Re φ(x, y , z)e −iωt

(1)

√

where i = −1, Re{.} corresponds to the real part of the argument inside brackets and the spatial (time-independent) term φ(x, y , z) is governed by the Laplace’s equation (∇ 2 φ = 0) within the fluid, subject to Neumann boundary conditions on the body surface, the seabed at z = 0 and the free surface at z = d, and the Sommerfeld radiation condition. For the full descriptions of boundary conditions, one can refer to [6]. It is convenient to decompose the total potential φ as a summation:

φ = φ D + φ R = (φ I + φ S ) + φ R

(2)

where φ R represents the radiation wave potential due to the motion of the bodies absent of incident waves and φ D is the solution to the wave diffraction problem for the still bodies, with φ I being the incident wave potential and φ S the scattering wave potential due to the presence of the bodies. The term φ R characterizes the interaction due to the independent

4

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

oscillation of rigid bodies in the different directions/modes of motion (1st-6th mode being, respectively, surge, sway, heave, roll, pitch and yaw), whereas φ D describes the interaction due to the incident wave diffracted on the bodies acting as obstacles (no movement). Knowing the potential, hydrodynamic pressures on the bodies can be obtained by the linearized Bernoulli’s equation and hydrodynamic forces can be calculated by integrating these pressures over the surface of each body [6,47]. For practical applications, like the WEC implementation discussed in the introduction, we are not interested in the potential field per se, but in hydrodynamic characteristics representing these resultant forces that are exerted on the bodies [23]. This is separately performed for the diffraction and radiation problem yielding, respectively: (i) the wave excitation force on the p pq pth body along the ith mode of motion F i , a complex quantity; (ii) the hydrodynamic reaction forces f i j acting on the q

pth body in the ith mode of motion due to oscillation of the qth body in the jth mode of motion with amplitude ξ j . The latter forces can be expressed as [48]: pq

pq

f i j = ω2 (ai j +

i

ω

pq

q

b i j )ξ j

(3) pq

pq

which leads to the well-known definition of the added mass ai j and added damping b i j coefficients, both real quantities dependent on the frequency ω . The added mass and damping coefficient are symmetric with respect to the modes of motion and the array bodies, and invariant under any rotation or translation of the coordinate system. Diffraction forces are invariant under translation vertical to the wave propagation direction, symmetric with respect to the wave propagation direction, whereas translation of the coordinate system along the wave propagation direction by D simply creates a phase shift of the diffraction forces by exp(i · k · D ). Ultimately, the necessary hydrodynamic characteristics to describe the problem, p pq pq corresponding to our quantities of interest QoI, are force F i and the added mass ai j and damping b i j coefficients, forming for each mode of motion the diffraction force vector and the added mass and damping matrices, respectively. For WEC applications we are commonly interested in oscillating bodies in the heave mode [26] (properly constrained through tethers to allow significant movement only along this direction), restricting the relevant motion to i = j = 3. For notational simplicity we will drop the dependence on i and j in the next section, with the understanding that the estimated characteristics can correspond to any directions of motion. The resultant output vectors, collecting all relevant unique hydrodynamic terms for all bodies (i.e., excluding symmetric, thus redundant, elements for added mass/damping), will be denoted as F, a and b, respectively. Estimation of these characteristics given the potential fields is discussed in detail in [6,47] and is not further elaborated here. Ultimately, for a given frequency ω all aforementioned QoIs can be expressed as some functions explicitly dependent on the array layout W [21] and for the diffraction problem additionally on the wave incident angle β . Appendix A briefly reviews the multiple scattering (MS) solution method [6,10,13] for the hydrodynamic interaction solution, and discusses its relevant physical concepts as well as computational challenges. As discussed in the introduction, MS and other relevant exact solvers have an associated high computational burden. To alleviate this burden, one can always employ some low-fidelity numerical version of the solver, for example, MS with fewer scattering order or fewer retained eigenfunctions (details in Appendix A), but the compromise on solution quality might be substantial. Instead, we consider here a data-driven surrogate modeling approach to address this computational challenge. 3. Surrogate modeling for hydrodynamic characteristics using a hierarchical interaction decomposition Surrogate models provide a simplified mathematical approximation to the input/output mapping of computationally expensive functions [49]. They are constructed in a data-driven manner, by interpolating or regressing over a limited number of input/output data from the expensive function, formally denoted as training (or support) points or simulation experiments. Generating training data and calibrating the surrogate model can take considerable computational effort, but after construction, the surrogate model can provide predictions with negligible numerical complexity compared to the original expensive function. In this setting, it becomes computationally affordable to adopt high-fidelity numerical solvers for obtaining the training data. Kriging (or Gaussian process regression) [46] is adopted as the surrogate model choice in this work, for its highly versatility in approximating complex functions and the ability to provide a local approximation of the surrogate model prediction error [50], but the overall formulation is applicable to other commonly used surrogate models. Before discussing further the surrogate model implementation, a quick overview of Kriging is first offered. 3.1. Review on Kriging Let u ∈ Rnu and v (u) ∈ R denote the input and output of the function to be approximated. For illustration clarity, we assume a scalar output v for now, with the extension to multivariate outputs discussed at the end of this section. The fundamental assumption of Kriging is that the output function v(u) is one realization of a Gaussian process (GP) coupled with a global regression model [46,51]. Formulation of Kriging requires first selecting the regression basis functions h(u) and the correlation function R (u, u |s) for the GP, the latter dependent on the hyper-parameters s to be calibrated. Let U = [u1 . . . u Nk ] T denote the training input and V = [ v (u1 ) . . . v (u Nk )] T the corresponding output, Kriging approximates the response v at u as a Gaussian distribution conditioned on such observation [46]:

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

 v (u)|U, V ∼ N v˜ (u), σ 2 (u)

5



(4)

where N stands for Gaussian distribution with mean a and variance b and ∼ denotes “distributed according to”. The predictive mean v˜ (u) is the Best Linear Unbiased Prediction (BLUP) for the complex function [46] whereas the predictive variance σ 2 (u) offers a local estimation of the variability of the surrogate model predictions (i.e., approximation of the Kriging error). Expressions for v˜ (u) and σ 2 (u) may be found in [46]. With respect to the surrogate model characteristics, in this study, we adopt constant basis function, h(.) = 1, the generalized exponential as correlation function [51], and perform the calibration of the hyper-parameter vector s through Maximum Likelihood Estimation [52]. Note that in the illustrative examples considered later, the surrogate model performance was found to be insensitive to these choices. For the multivariate output case (v(u) ∈ Rn v >1 ), since we do not anticipate strong correlation between output components in the context of multi-body wave interactions, the recommendation in [53] is followed to formulate a separate Kriging model for each output component. Common regression basis functions and correlation function are selected to improve the computation efficiency, but the hyper-parameters s need to be repeatedly calibrated/calculated for each component. Approximation to the multivariate output vector v can be finally obtained as a direct extension to the scalar case, augmenting the predictive mean and variance for all components. 3.2. Direct surrogate modeling implementation One surrogate modeling implementation option for the hydrodynamic characteristics estimation is to adopt a direct formulation, regarding the array layout coordinate u(W) = [x1 , x2 , . . . , x N , y 1 , y 2 , . . . , y N ] ∈ R2N as the surrogate input, and all components of F, a and b as the surrogate output. This approach however faces a number of challenges: (a) The surrogate model suffers from the curse of dimensionality with large-scale arrays (N  1). The amount of training data (N k ) needed to construct an accurate surrogate model increases exponentially with the input dimensionality [54], creating significant computational challenges. Assuming each output component is independently approximated by a separate surrogate model, the output dimensionality also creates challenges, albeit less critical, due to the need to tune multiple surrogate models. (b) For a new examined array size N, a separate surrogate model needs to be developed. Even if available, surrogate models corresponding to different array sizes provide no information for this creation. (c) This approach fails to incorporate knowledge about physical characteristics of the input/output relationship, for example, some symmetry or invariance properties of the solutions [43–45], including the fact that all QoIs are invariant to the ordering of the bodies, or that the potential is radially symmetric for vertical axisymmetric bodies [47]. In what follows, a different surrogate model development approach is examined, with the foundation being a hierarchical interaction decomposition based on the many-body expansion (MBE) [5] principle. Though the emphasis is placed on addressing challenges (a) and (b), this development also efficiently addresses challenge (c). Before discussion the surrogate model development, a brief discussion on MBE is first offered. 3.3. Review of the many-body expansion principle Let ψ(W) represent some hydrodynamic interaction characteristic of interest, which can be any component of F, a or b, that considers the entire-array interaction (thus ψ is a function of the entire array geometry W). MBE decomposes ψ(W) to contributions from constituent cluster of n bodies, with n ranging from 1 to N:

ψ(W) =

N  p

ψ(w p ) +

N −1  N  p

ψ(w p , wq ) +

q> p

N −2 N −1  N   p

ψ(w p , wq , wr ) + · · ·

(5)

q> p r >q

where ψ is the additive interaction effect for ψ considering n bodies. This term can be obtained by computing the total interaction of the n-body cluster, and subtracting all single-body contributions and additive effects from its distinct subsets. An alternative MBE expression may be established by replacing the multi-body subscript notation in Eq. (5) with the use of notation {W}ln to represent the lth distinct n-body cluster for the N-body array W, which helps rewriting the summation in Eq. (5) with regard to all unique n-body combinations:

ψ(W) =

N  l =1

N !/[2!( N −2)!]

ψ({W}l1 ) +

 l =1

N !/[3!( N −3)!]

ψ({W}l2 ) +



ψ({W}l3 ) + ...

l =1

Under this representation, the additive interaction effect for the lth n-body cluster ψ({W}ln ) can be expressed as:

(6)

6

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

⎡ ⎤ n!/[2!(n−2)!] n n    ψ({W}n ) = ψ({W}n ) − ⎣ ψ({{W}n }r1 ) + ψ({{W}n }r2 ) + · · · + ψ({{W}n }nr −1 )⎦ l

l

l

l

r =1

=

ψ({W}ln ) −

n 

⎡

r =1

n!/[s!(n−s)!]

⎣

s =1



l

r =1

⎤

(7)

ψ({{W}ln }nr −s+1 )⎦

r =1

which can also use the total interaction effects of clusters (ψ ) rather than additive ones ( ψ ):

ψ({W}ln ) = ψ({W}ln ) − =

n −1 

⎡ ⎣(−1)s

n 

ψ({{W}ln }nr −1 ) + · · · + (−1)n−1

r =1 n!/[s!(n−s)!]



n 

ψ({{W}ln }r1 )

r =1

⎤

(8)

ψ({{W}ln }nr −s )⎦

r =1

s =0

For the two and three body clusters, that will be used later in this work, these expressions simplify, respectively, to

ψ(w p , wq ) = ψ(w p , wq ) − ψ(w p ) − ψ(wq )

(9)

ψ(w p , wq , wr ) = ψ(w p , wq , wr ) − ψ(w p , wq ) − ψ(w p , wr ) − ψ(wq , wr ) + ψ(w p ) + ψ(wq ) + ψ(wr ) = ψ(w p , wq , wr ) − ψ(w p , wq ) − (w p , wr ) − (wq , wr ) − ψ(w p ) − ψ(wq ) − ψ(wr )

(10)

Depending on the interested hydrodynamic characteristics, not all clusters need to be examined in the summation of Eq. (5) or Eq. (6). For quantities involving only the pth body (i.e. ψ = F p , a pp or b pp ), only clusters including the pth body should be considered (evidently other clusters provide no contribution). Similarly, for quantities involving two bodies (i.e. ψ = a pq or b pq ), only clusters that include both bodies will need to be examined. Let ψ p (W) denote the quantities involving only the pth body and let ψ pq (W) denote the quantities involving both the pth and qth bodies. Eq. (6) is modified to

ψ p (W) =

N   l =1

pq



I w p ∈ {W}l1 ψ({W}l1 ) +

ψ (W) = 0 +

N !/[2!( N −2)!]

 l =1

I [w p wq ] ∈

N !/[2!( N −2)!]



l =1

{W}l2





I w p ∈ {W}l2 ψ({W}l2 ) + ...

ψ({W}l2 ) +

N !/[3!( N −3)!]



I [w p wq ] ∈

{W}l3



(11)

ψ({W}l3 ) + ...

l =1

where I [.] corresponds to the indicator function being one if the condition insides the brackets holds and zero else. Ultimately ψ p (W) has contributions from the pth body alone, as well as (N − 1) two-body subsets including body p, ( N − 1)( N − 2)/2 three body-clusters including body p and so forth. For ψ pq (W), its contributions come only from the two-body cluster including both bodies, (N − 2) three-body clusters including both bodies and so forth. Same modification as in Eq. (11) applies to the additive effects in Eq. (7) and Eq. (8). As demonstrated above, MBE is solely based on combinatorial techniques such as the principle of inclusion/exclusion. Its idea of representing the total interaction as summations of contributions from increasing number of bodies is not unique, and can be actually viewed as a special case of the high dimensional model representation (HDMR) technique [55,56] (or the analysis of variance (ANOVA) in statistics [57]), which treats the multiple input variables analogously to how MBE deals with physical objects, e.g., atoms in molecular simulations and WECs in our work. HDMR has a broader application scope than MBE: it can deal with any high-dimensional function without requiring that the system can be decomposed to sub-systems and the outputs from such sub-systems be available. On the other hand, this generalizability of HDMR also comes with a cost: the decomposition terms in HDMR are constructed completely under mathematical principles such as the orthogonality condition, and thus lack the clear physical interpretation of their MBE counterparts. As a result, for systems where MBE can rapidly converge for particular physical reasons, HDMR might not be able to perform similarly well. In an analogous manner to our work, there has also been approaches to combine HDMR with surrogate modeling [58–60]. We refer them to readers for solving problems with systems that are not decomposable via MBE. Although MBE is exact upon the order of N, the convergence behavior of MBE with respect to the cluster size might be problematic for some applications [43,61], but this shall not be an issue for the examined hydrodynamic interaction problem, with the following justification coming from considering the exact solution obtained through infinite scatteringorder summation in MS discussed in Appendix A. For MBE truncated at the n-body level, because of its additive nature, the contributions from all orders (up to infinity) of scattering waves whose path is confined within n-body clusters has been included, while the truncation error is only due to not accounting for scattering waves that visit more than n bodies (i.e., having a MS scattering order of at least n). An alternative way to think of this is that the newly created additive potential due to a new body in the cluster (compared to all lower order clusters) needs to visit all bodies in the n-body cluster and so has by default a scattering order of at least n. On the other hand, for MS truncated at the order of m, its error directly comes from the inability to capture any waves with the scattering order higher than m. Comparing both sources of errors,

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

7

we can see that at the same truncation level (n = m), MBE include richer information than MS from its knowledge of higher scattering-order contribution among the n-body clusters and their subsets. Thus, the convergence rate of MBE is even faster than MS, whose favorable convergence characteristics have been already demonstrated [6,47]. It is therefore expected that MBE can offer an accurate reconstruction to the exact solution even when truncated at some value N b that is lower than the typically adopted max scattering order for MS. It should be pointed out that MBE can be applied with any desired numerical approach for estimating the hydrodynamic characteristics, not necessarily constrained to MS. The provided MS-related arguments simply justify the convergence properties of MBE for the wave interaction problem. 3.4. Surrogate modeling using hierarchical interaction decomposition 3.4.1. Motivation and overview of the formulation The MBE formulation of Eq. (5) [or Eq. (11)] motivates a different approach for the surrogate model implementation. Separate surrogate models can be established to approximate the additive effects on the right side of Eq. (5) from n = 2 up to n = N b , the intended MBE truncation order. The total hydrodynamic characteristics are ultimately obtained by summing up contributions from all n-body approximations along with single body contributions ψ(w p ), whose exact value can be very efficiently calculated as further discussed in Appendix B. This alternative surrogate model implementation addresses key challenges discussed in Section 3.2. Most importantly, both input and output dimensionality are significantly reduced when truncating the MBE expansion after first few terms (corresponding to few-body clusters), as these clusters are anticipated to have bigger contributions to the total response. It leads to substantial computational benefits; the more impactful surrogate models (corresponding to ψ for small n) can be constructed with high accuracy at a reduced computational burden, since their input and output will be low dimensional. This is also beneficial, as it allows to obtain higher-fidelity surrogate training data, for example with higher scattering order N m in MS as discussed in Appendix A, thanks to the significantly higher computational efficiency of MS for fewer-body arrays. Finally, each surrogate model development only depends on the considered cluster size in MBE expansion (n), and therefore is independent of the entire array size (N). This enables the use of surrogate model approximation even for large arrays (N  1), and also ensures that new surrogate models are only needed if the MBE truncation order N b is increased. Beyond the MBE implementation, the surrogate model development needs to carefully incorporate knowledge about the physical characteristics of the input/output relationship. Although this is independent of the array size or whether MBE is utilized, it can be applied more seamlessly for fewer-body arrays. The fundamental principles that need to be considered are the following. (a) Invariance and symmetry principles need to be incorporated, as discussed in Section 2 for the diffraction forces and the radiation interaction coefficients. Specifically, this information can be used to reduce the surrogate model input domain. Taking the two-body diffraction case as an example, the surrogate model simply needs to consider the following class of two-bodies: the second body can be restricted into the first quadrant or the upper half-plane with the first body at the origin, as arbitrary two-body configurations can be mapped to this class. Note that this input domain needs to be carefully defined to avoid extrapolations for predicting over all possible configuration geometries. (b) The response is permutation invariant with respect to the body orderings. Combining this with the above invariance and symmetry principles, a single surrogate model can be used to predict multiple hydrodynamic characteristics. For example, the surrogate model approximating F 1 ([w1 w2 ]) can also provide prediction to F 2 ([w1 w2 ]), through proper permutations so that w2 becomes the first body. This also affects the proper definition of surrogate input domain to avoid extrapolations. For the two-body diffraction case discussed previously, the possible placement domain for the second body [see discussion in point (a) above] is the entire half-plane if we only develop a surrogate model to approximate F 1 , while it would had been the first quadrant if both forces were separately approximated. Ultimately, the minimal independent hydrodynamic characteristics that need to be approximated are one diffraction force, and for radiation coefficients (added mass and added damping matrices) one diagonal (single body) and one off-diagonal element (pair of bodies). This drastically reduces the output dimension as it is now independent of the cluster size (n). (c) The best input characterization might be output dependent. Instead of relying on the Cartesian coordinate system, alternative array geometry descriptions can be adopted by considering the pairwise distances and angles between bodies. This choice has no impact on the exact solver but can drastically affect the surrogate model performance, as certain characterizations might lead to simpler (e.g., smoother or more linear) input-output relationships and facilitate more accurate surrogate model prediction [62]. The most appropriate input characterization can be selected through physical understanding of the problem or by adopting an entirely data-driven approach: examine different options and adopt the one providing the highest accuracy according to some validation criteria. The integration of these principles with the MBE decomposition leads finally to the surrogate model framework discussed next.

8

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

3.4.2. Hierarchical interaction surrogate modeling Two key components need to be described for the hierarchical interaction surrogate model implementation: the development of the surrogate models for the additive effects of the n-body cluster, and the use of these surrogate models to provide the hydrodynamic interaction predictions for arbitrary arrays. These components are discussed next for a general formulation. Appendix B contains detailed information on the adopted development approach for the numerical examples presented later. Let W define the considered geometric domain for the multi-body configuration. The first step is to obtain the singlebody solutions. For the considered hydrodynamic characteristics, it is sufficient to numerically solve them once for a body at the origin, since any other placement can be analytically mapped to that. For the n-body cluster, assuming that all lowerorder surrogate models (for clusters with fewer bodies) have been developed, the surrogate model for the additive effects ψn can be established through the following steps: Step D.1: Characterize the output and input for the surrogate model. At this stage, multiple characterization choices can be considered, with final selection made in Step D.4. For our problem of interest, the input/output will be separately defined for the diffraction and radiation problems. The notation vnc and unc , with superscript c being dif (diffraction) or rad (radiation), is used herein to distinguish the surrogate output and input. A transformation gnc that maps the Cartesian coordinates of n bodies to the surrogate input unc will also be defined. It is expected that the higher-order surrogate input unc can be chosen by extending the one for the precedent order (unc −1 ) with additional inputs to describe the placement of the added body. Step D.2: Define the proper surrogate input domain U nc to eliminate extrapolated predictions for any n-body configuration within W , exploiting symmetries and invariances. Define also the output transformation Tnc that maps the output for the n-body configuration in W to the corresponding configurations in U nc . Step D.3: Generate N k experiments for the surrogate input unc within the defined domain U nc via a space-filling approach (e.g., Latin hypercube sampling). For each experiment, obtain its surrogate output vnc (unc ) through Eq. (7) [or Eq. (8)] considering only the relevant clusters (contributing to the specific output). For the lower-order additive effects ( ψr
 Nk  i =1

CC = 

 Nk  i =1

v (ui ) −

v (ui ) −

 Nk

i =1

 Nk

v (ui )/ N k

v (ui )/ N k i =1



2 

v˜ −i (ui ) − Nk i =1



 Nk

i =1

v˜ −i (ui ) −

v˜ (ui )/ N k



 Nk

v˜ (ui )/ N k i =1

2

(12)

where v˜ −i (ui ) represents the predictive mean of the training input (ui ) after removing it from the available observations (training points). Closed form expression for v˜ −i (ui ) may be found in [52]. If multiple input characterization choices have been examined in Step D.1, one can decide on the most appropriate one based on LOOCV accuracy. Step D.5: If the accuracy is insufficient, obtain additional experiments and go back to Step D.3. In this case advanced design of experiment approaches can be used, exploiting the already established surrogate model to guide the selection of new experiments [46]. If the accuracy is deemed sufficient, the surrogate model development can be completed. This development process is successively performed for different values of n ranging from 2 to N b . The value of N b can be adaptively determined by monitoring the MBE convergence properties and surrogate model approximation quality. As n increases, greater computational effort is required to obtain a specified prediction accuracy level due to the curse of dimensionality, while the amplitude of additive effects is expected to decrease due to the MBE convergence. This naturally suggest that the MBE expansion can be stopped at the order where a sufficient prediction accuracy cannot be obtained within a reasonable computational effort (i.e., when the training data size N k reaches a predefined limit). Considering the commonly used max scattering order N m = 4 for large dimensional WEC arrays, the truncation value is set to the three-body level (N b = 3) in the illustrative implementation and in Appendix B. The logic behind this choice has been explained in Section 3.3: since MBE converges faster than MS, it should be able to use a smaller truncation order. This selection will be further shown to be sufficient to achieve high prediction accuracy for the WEC applications of interest. Finally, the surrogate approximated MBE expansion truncated at the order of N b can be established as:

˜ W) = ψ(

N  l =1

N !/[2!( N −2)!]

ψ({W}l1 ) +

 l =1

˜ W}2 ) + · · · + ψ({ l

N !/[ N b !( N − N b )!]

 l =1

˜ W} b ) ψ({ l N

(13)

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

9

˜ W}n ) where only the clusters contributing to the specific output need to be considered [modification of Eq. (11)] and ψ({ l is the surrogate approximated additive effects for the lth n-body cluster {W}ln . The detailed definitions and evaluation procedures for those quantities can be referred to Appendix B. 4. Case study for prediction of hydrodynamic characteristics In this section, the surrogate modeling framework for the multi-body hydrodynamic interaction estimation is validated. The development and verification of two-body and three-body surrogate models is discussed first and subsequently the established surrogate models is used to predict the hydrodynamic characteristics of different multi-body array configurations. The consider bodies are identical cylinders oscillating in heave direction only, with radius r p = 3 m and mass m p = 2.54e5 kg, corresponding to a draft of D r = 9 m and period of resonance 6.9 s whereas water depth is d = 50 m. The ‘exact’ hydrodynamic interaction solver is MS coupled with matched eigenfunction expansion. Numerical details for the approach are included in [6,47]. The eigenfunction series are truncated at twenty terms for the fluid region adjacent to the cylinder and fifty terms for the fluid region beneath the cylinder, and unless otherwise specified, the max scattering order is set to N m = 5 to ensure MS convergence [47]. The array deployment domain (W ) is set as a rectangular domain with a length of 300 m along the x-direction and 600 m in the y-direction, similar to what has been used in previous WEC layout optimization literature [20,26]. 4.1. Two-body and three-body surrogate models We first establish additive-effect surrogate models from n = 2 and n = 3 bodies based on the approach discussed in Appendix B. To see the effect of varying wave frequencies, three frequencies ω = 0.87, 1.14 and 1.62 rad/s are chosen, corresponding to dimensionless wave number parameters of kr p = 0.2, 0.4 and 0.8, respectively. To create consistent comparison between further examined options, no specialized design of experiment approach is adopted for the surrogate model development; Latin hypercube space filling sampling is used with a pre-determined training data size of N k = 250 × 2nu to achieve satisfactory accuracy. This dimensional-dependent selection of training data size follows the general guideline in [65]. We use test-sample validations with N t = 1000 samples to verify the surrogate model accuracy. The validation metrics include the previously used correlation coefficient (CC), which is identical to the LOOCV version in Eq. (12) but replaces N k with N t and the LOOCV prediction v˜ −i (ui ) with the ordinary prediction v˜ (ui ). The normalized root mean square error (NRMSE) is also used:

 Nt

NRMSE(%) = 100 ×

i =1

[ v (ui ) − v˜ (ui )]2 / N t

(14)

v max − v min

where v max and v min is the maximum and minimum output among all test samples. Such normalization is introduced to facilitate comparison across different cases, and lower NRMSE percentage indicates higher accuracy. The same surrogate model development procedure is repeated under all three wave numbers. In Tables 1–2, the prediction accuracy is reported. Table 1 Test-sample validation results for the two-body surrogate model (N k = 500 for diffraction and 250 for radiation). kr p

Metric

Re( F 1 )

Im( F 1 )

a11

a12

b11

b12

0.2

CC NRMSE

1.00 0.07%

1.00 0.04%

1.00 0%

1.00 0%

1.00 0%

1.00 0%

0.4

CC NRMSE

1.00 0.04%

1.00 0.09%

1.00 0.01%

1.00 0.001%

1.00 0%

1.00 0%

0.8

CC NRMSE

1.00 0.13%

0.994 0.25%

1.00 0.07%

1.00 0.008%

1.00 0.008%

1.00 0.001%

Table 2 Test-sample validation results for the three-body surrogate model (N k = 2000 for diffraction and 1000 for radiation). kr p

Metric

Re( F 1 )

Im( F 1 )

a11

a23

b11

b23

0.2

CC NRMSE

0.872 1.29%

0.867 0.79%

0.975 0.53%

0.997 0.29%

0.989 0.41%

0.997 0.22%

0.4

CC NRMSE

0.825 1.35%

0.833 2.35%

0.963 0.66%

0.965 0.61%

0.981 0.59%

0.963 0.68%

0.8

CC NRMSE

0.553 2.96%

0.472 4.35%

0.901 1.06%

0.967 0.53%

0.926 1.03%

0.946 0.99%

10

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

Table 3 Test-sample validation results using alternative inputs for diffraction problem. kr p

Metric

Re( F 1 )

Im( F 1 )

0.2

CC NRMSE

0.588 1.92%

0.451 3.26%

0.4

CC NRMSE

0.524 2.57%

0.346 4.49%

0.8

CC NRMSE

0.184 3.45%

0.109 5.05%

Table 4 Test-sample validation results using alternative inputs for radiation problem. kr p

Metric

a11

a23

b11

b23

0.2

CC NRMSE

0.851 1.46%

0.925 0.86%

0.927 1.66%

0.956 1.29%

0.4

CC NRMSE

0.867 0.88%

0.931 1.09%

0.966 1.01%

0.960 0.96%

0.8

CC NRMSE

0.877 1.32%

0.906 1.04%

0.971 0.55%

0.923 1.16%

Across almost all cases, we can observe excellent surrogate model prediction accuracy from high (close-to-one) CC and low NRMSE (near 0%) values. The only relatively challenging task is predicting the three-body additive contribution to waveexcitation forces as shown at the third and fourth column in Table 2. Even in that case, the NRMSE values (never exceeding 5%) indicate good overall predictive capabilities (as a reference, 10% is a NRMSE threshold recommended for adequate accuracy in [66]). The relative difficulty encountered for the three-body diffraction problem should be attributed to the higher input dimensionality and stronger nonlinearity. However, as we will demonstrate later in Section 4.2, such accuracy level is still sufficient to yield near identical prediction of wave excitation forces, partly because those 3-body additive terms have smaller amplitude and are less important according to the convergence trait of MBE. Further accuracy improvement can be accommodated if an adaptive DoE is adopted [46], which, though, falls outside the scope of this study as discussed earlier. One can also observe a slightly decreased prediction accuracy when the wave number (k) increases. This should be attributed to the fact that, though the spatial design domain (W ) is fixed, the wave number impacts the non-dimensional distance between bodies (this non-dimensionalization can be enforced by multiplying with the wave number [14]) and has the effect of an equivalent expansion of the non-dimensionalized domain, making the surrogate model training more challenging. The use of alternative input/output definitions for the three-body surrogate model is also investigated. First the focus is on the input choices. Adopting convenient definitions, for (i) the diffraction problem, the Cartesian coordinates of the second and third bodies relative to the first one are used as the input, whereas for (ii) the radiation problem the pairwise distance between all bodies are adopted. Following terminology of Appendix B, the new input vectors are di f



u3 = x21



urad 3 = l12

y 21 x31 y 31

l13

l23

T

;

di f

T

;



di f

T y 21 x31 y 31 T   x231 + y 231 x223 + y 223

g3 ([w1 w2 w3 ]) = x21

g3 ([w1 w2 w3 ]) =



x221 + y 221

(15) (16)

For the radiation problem this choice ultimately corresponds to the replacement of the angle information θ213 (the bearing angle between vectors connecting the first body to the other two) with distance l23 . Results under alternative input definitions are reported in Tables 3–4. We observe considerable accuracy losses when comparing to results in Table 2. This performance discrepancy demonstrates the advantage of using polar coordinates to define surrogate inputs, as well as the inclusion of the angle information rather than the pairwise distance for the radiation input. Such advantage should be attributed to the fact that the chosen input variables represent more direct characterizations of the underlying scattering and radiation processes [6]. For alternative outputs, the selection of a different off-diagonal radiation output is examined, using [.]12 element of the radiation matrix instead of the initially selected [.]23 element. The other off-diagonal element [.]13 shows practically identical behavior as [.]12 under our chosen input due to permutation invariance. The surrogate model prediction accuracy is reported in Table 5. Comparing to Table 2, the results show that our adopted output [.]23 features a significantly higher accuracy than the alternative one [.]12 . The advantage of such input/output selection can be actually explained by examining the potential field that contributes dominantly to such output in the MS process. The potential originates from the zero-order radiation on the third body, then gets scattered by the first body and finally excites the second body. One can see that the three

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

11

Table 5 Test-sample validation results for alternative radiation outputs. kr p

Metric

a12

b12

0.2

CC NRMSE

0.981 0.97%

0.972 1.23%

0.4

CC NRMSE

0.931 1.49%

0.934 1.77%

0.8

CC NRMSE

0.710 1.08%

0.755 2.07%

chosen input variables (two pairwise distances and one bearing angle from the other two body to the first one) provide a natural description of this process. The discussion in the previous two paragraphs demonstrate the effectiveness offered by the data-driven selection of the proper input and output characterizations. Different potential candidates can be examined and the most appropriate one is finally chosen. The choice can be further intuitively rationalized leveraging the physical description of the problem; but even the data-driven validation alone can be leveraged to guide the selections. A final remark is warranted about the computational efficiency gain in obtaining the training data (Step D.3 of Section 3.4.2), offered by using the surrogate model predictions when subtracting lower-order cluster contributions. When evaluating the three-body additive effects for the training data, thanks to the high accuracy of the two-body surrogate models (Table 1), rarely the exact MS solver had to be used for the two-body combinations. In all cases, it contributed to a two-fold (for diffraction) and three-fold (for radiation) computational saving per three-body experiment. This strategy would have been even more beneficial when dealing with higher-order cluster contributions, or more expensive numerical solvers. 4.2. Predictions for multi-body configurations The established surrogate models are now used to predict hydrodynamic characteristics for different multi-body configurations. A fixed wave number setting kr p = 0.4 is considered, following the typical setup in [20,26]. The proposed approach (further referred as Hierarchical) is compared to an alternative surrogate model formulation and a lower-fidelity MS solver: 1. Direct: the direct surrogate modeling implementation discussed in Section 3.2. To facilitate a fair comparison to the Hierarchical implementation, symmetry and invariance principles with respect to the input are also adopted to reduce the dimensionality. The surrogate model input for the diffraction problem is defined by the locations of all bodies relative to the first body considering the incident wave angle β : udi f = [l12 , . . . , l1N , θ12 − β, . . . , θ1N − β] ∈ R2N −2 , following the same notation as Appendix B. For radiation problem, leveraging the rotational invariance the angular information can be further simplified, leading to urad = [l12 , . . . , l1N , θ13 − θ12 , . . . , θ1N − θ12 ] ∈ R2N −3 . As for output, all components of F, and the upper triangle elements in a and b (lower parts are not needed according to symmetry) are separately approximated. The number of training data is always fixed at the maximum size used in Hierarchical, i.e., N k = 2000. The surrogate input domain is chosen with the goal of not extrapolating for any configuration inside the whole deployment domain W . 2. LF (Low-Fidelity) MS: a faster, yet less accurate, MS solver. It reduces the max scattering order to 3 to match the max MBE expansion order and the eigenfunction series to 10 for the fluid region adjacent to the cylinder and 20 for the fluid region beneath the cylinder, contrasting the choice of 5, 20 and 40 in the exact MS. After constructing the surrogate models, our first task is to provide a simple assessment of how good they perform. Therefore, the initial focus will be on visualizing the approximation quality of those methods for two classical topological configurations taken from [36] and [67], as shown in Fig. 2. These configurations correspond to different array size (N = 4 or 5). In addition, their wave incident angle and characteristic spacing between devices are also adjustable to allow examining a wider range of different configurations. In Fig. 3–6, the interaction effects for selective non-dimensional hydrodynamic characteristics from different methods are presented for varying wave incident angle (not needed for radiation due to its symmetry) and characteristic length. Results for other hydrodynamic characteristics demonstrate nearly identical trends. The non-dimensionalization is established following the recommendations in [6,47]:

F¯ = F/ρ g π r 2p a¯ = a/ρπ r 3p b¯ = b/ωρπ r 3p

(17)

where ρ corresponds to the fluid density. In certain scenarios, the results for direct surrogate modeling are not shown, because it yields extremely poor predictions that would prohibit a good resolution for the result comparison among other approaches if plotted in the same figure.

12

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

Fig. 2. Characteristic 4 and 5 body array topological configurations.

Fig. 3. For the 4-body array, the change of chosen diffraction-related hydrodynamic coefficients (the second body wave excitation force F¯ 2 ) versus varying characteristic lengths L (with β = 0), including (a) its real part and (b) its imaginary part, as well as its change versus varying incident angles β (with L = 5r p ), including (c) its real part and (d) its imaginary part.

Results in Fig. 3–6 show that the Hierarchical approach facilitates good agreement with the exact solution over the whole range of varying characteristics lengths/angles. LF-MS sometimes over or under-estimates radiation, but provides adequate accuracy for forces, though at a significant higher computational cost than Hierarchical. The direct surrogate modeling approach fails to offer even marginal predictive capabilities despite of the large number of experiments it used, affirming the substantial benefits from the hierarchical surrogate modeling implementation. To further examine how the established approach performs on larger-scale arrays with a wider variety of configuration geometries (e.g., asymmetric and/or unevenly spaced layouts), two array cases of N = 10 and 20 bodies are examined. This should cover the common size ranges for WEC arrays, which is usually chosen below 10 in existing studies [26,68] and hardly goes beyond 20 [36]. For each case N t = 100 different random array configurations are considered. Each configuration is obtained by randomly placing bodies inside W following uniform distribution. The approximated hydrodynamic charac˜ from both Hierarchical and LF MS are then compared to the actual values F, a and b from the exact MS. teristics (F˜ , a˜ , b) Direct surrogate modeling implementation is not included due to its inability to achieve sufficient accuracy. Approximation errors are quantified using the Frobenius norm (||.|| F ) of the error vectors/matrices:

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

13

Fig. 4. For the 4-body array, the change of chosen radiation-related hydrodynamic coefficients versus varying characteristic lengths L (with β = 0), including (a) diagonal added mass coefficient a¯ 11 , (b) off-diagonal added mass coefficient a¯ 14 , (c) diagonal added damping coefficient b¯ 11 , and (d) off-diagonal added damping coefficient b¯ 14 .

Fig. 5. For the 5-body array, the change of chosen diffraction-related hydrodynamic coefficients (the second body wave excitation force F¯ 2 ) versus varying characteristic lengths L (with β = 0), including (a) its real part and (b) its imaginary part, as well as its change versus varying incident angles β (with L = 5r p ), including (c) its real part and (d) its imaginary part.

14

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

Fig. 6. For the 5-body array, the change of chosen radiation-related hydrodynamic coefficients versus varying characteristic lengths L (with β = 0), including (a) diagonal added mass coefficient a¯ 22 , (b) off-diagonal added mass coefficient a¯ 24 , (c) diagonal added damping coefficient b¯ 22 , and (d) off-diagonal added damping coefficient b¯ 24 .

Table 6 Average approximation errors for 10 and 20 body random arrays.



Array size

Approach

err(F)

err(a)

err(b)

N = 10

Hierarchical LF-MS

1.03% 0.82%

0.19% 4.05%

0.31% 1.78%

N = 20

Hierarchical LF MS

2.27% 0.96%

0.46% 4.05%

0.66% 2.34%

√

err(F) = || (F − F˜ )(F − F˜ )∗ || F /|| FF∗ || F err(a) = ||a − a˜ || F /||a| F err(b) = ||b − b˜ || F /||b|| F

(18)

where ∗ stand for the complex conjugate. This norm takes the square root of the summed element-wise square errors. For the force, the focus is on the amplitude of the (complex) error, which is a standard approach to assess accuracy for complex quantities [69]. Table 6 reports the average (over the 100 examined configurations) relative errors as percentages, with smaller values suggesting better agreement. From these results, we can see that both approaches achieve very good agreement with high accuracy even for large-scale arrays as the relative error never exceeds 5%. Hierarchical obtains better (for radiation) and comparable (for diffraction) approximation accuracy against LF-MS. The bigger error for diffraction stems from the lower accuracy of its three-body cluster surrogate models, as discussed in the previous section (see results in Table 3). As mentioned earlier, a potential remedy for this challenge is to employ some advanced DoE strategy to further improve accuracy. We also notice that the error for Hierarchical seems to be more array size-dependent than that for LF-MS. This is an expected trend under the combinatorial nature of MBE: even if the errors from additive contributions are small, they can accumulate to non-negligible total errors especially for large-scale systems. To facilitate a more direct understanding on this characteristic, let us simply assume the prediction error in each n-body additive contribution [ ψ({W}n ) in Eq. (6)] to be the same, denoted as εn (the superscript .n denotes the dependency of that error on the number of bodies). If all combination of bodies contribute to ψ , ˜ W) can be written as: the total error between ψ(W) and its prediction ψ(

˜ W) = N ε 1 + [ N ( N − 1)/2]ε 2 + [ N ( N − 1)( N − 2)/6]ε 3 + · · · ψ(W) − ψ(

(19)

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

15

Table 7 Average approximation errors for 10 body random arrays with different MS scattering order. MS scattering order

1

2

3

4

err(F) err(a) err(b)

9.63% 2.06% 33.9%

1.19% 0.26% 4.62%

0.15% 0.03% 0.36%

0.02% 0.01% 0.08%

Fig. 7. Computational running time versus array size for different numerical methods.

where the multiplication factor of εn rapidly increases with the number of bodies N, and thus even a small error can be amplified significantly when N is large. Similar trend also holds even if ψ is some specific hydrodynamic characteristic that not all combinations contribute to (e.g., only combinations that include body 1 contribute to F 1 ), which is also evident in Eq. (B.20) in Appendix B. The same-error assumption is of course not realistic, but the above discussion stresses nevertheless the importance of reducing errors εn when using the proposed hierarchical surrogate model implementation for large arrays, something that can be achieved with an adaptive DoE scheme. We also provide evidence on how fast MBE converges for our problem utilizing this random array setting. The combinatorial nature of MBE prohibits its use on the exact solver on this scale (for example, a single 10-body array requires at least solving 45 two-body arrays and 120 three-body arrays). Instead, we monitor the convergence of MS with respect to the max scattering order. As we have already discussed in Section 3, MBE components converge faster than MS, allowing MS convergence to be used as approximate indicator (effectively a lower bound) for the MBE behavior. We record the average error rate for 10-body arrays for the estimation using the MS (exact) solver with different max scattering orders (from 1 to 4) in Table 7. From the table, we see that the error rate reaches near-zero and the improvement are rarely significant for terms beyond third order; this means that our use of MBE up to the three-body cluster should be able to achieve good convergence. In fact, a major reason behind the proposed method’s good performance is due to the rapid convergence of MBE: the additional contributions from clusters with more than three bodies are negligible for the various configurations and wave number settings examined in this paper. For cases where it will be necessary to account for such higher-order contribution, our proposed method is still theoretically applicable, but its practical implementation can be more challenging: constructing the corresponding surrogate model for increasing cluster size n will inevitably require more training data, and will face more difficulty to find the best input/output characterization. Thus, we recommend such MS based error-rate convergence test as a preliminary check: the proposed surrogate modeling scheme is anticipated to work well, as long as the error rate corresponding to the maximum MS scattering order of three is indeed tolerable for the intended application. Finally, we report the typical running times of Hierarchical, LF MS and the exact MS solver for different array sizes in Fig. 7. All measurements were carried out on a PC with 2.0 GHz Xeon processor and 28GB RAM. Hierarchical always offers several orders-of-magnitude speedup compared to MS solvers. Even for the 20-body arrays, it only requires 0.38 s per run, which is 100 times faster than LF-MS (40 s) and nearly 2000 times faster than the exact MS (650 s). Considering the already demonstrated high accuracy, this speedup shows that the proposed hierarchical surrogate modeling scheme presents a promising direction to remove the computational limitations of using exact hydrodynamic interaction solvers in practical applications, such as the WEC layout optimization problem discussed next in Section 5.

16

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

Fig. 8. The q-factor versus the second body’s location for 2-WEC arrays estimated by different numerical schemes.

5. Application to WEC layout formation This section investigates the surrogate model implementation to support the WEC layout formation optimization. We adopt the same setting from a similar study [67] that has examined the hydrodynamic interaction impact on the WEC performance: the impedance matching approach is utilized for controlling the WEC array response (phase angles and amplitude of vibration), leading to the optimal power absorption given by [23]:

P max (W) =

1 ∗ −1 F B F 8

(20)

where F is the column vector representing the complex amplitudes of the exciting forces, defined earlier in Section 2, and B is the matrix of radiation-damping coefficients (composed of elements of b vector). This formulation does not explicitly use the added mass coefficients, which would have been needed in alternative optimal control formulations [70]. The interaction potential between WECs is measured by the q-factor, defined as the ratio of the power absorption from the N-WEC array to the one from N independently operated WECs:

q(W) = p

P max (W) p

N · P max

(21)

where P max is the single-WEC optimal power absorption, obtained similarly by Eq. (20). Calculating the q-factor requires the hydrodynamic interaction coefficients (F and B), which can be computationally demanding to obtain with exact solvers. For this reason, the Point Absorber (PA) approximation has been widely used in relevant layout optimization studies [20,71]. PA utilizes a small-body assumption (each body is sufficiently small and widely spaced compared to incident wave lengths) to neglect diffracted waves, resulting in a closed-form q-factor approximation formula [23] that only involves the wave number (k), incident angle (β ) and pairwise distances/angles among all bodies and therefore can be very efficient to evaluate. A comparison between the WEC power absorption approximation accuracy by PA and the one based on Eq. (20) using Hierarchical (i.e., the hierarchical surrogate predicted hydrodynamic coefficients) is first conducted. Inheriting the setup in [20], we consider the wavenumber setting of kr p = 0.4 with the incident wave heading towards the positive x-axis (β = 0). The first body is always fixed at the origin (w1 = [0 0] T ), while the other bodies are restricted in the design domain W = [0 m, 300 m] ⊗ [−300 m, 300 m]. The layout is also required to be symmetric with respect to the x-axis. In Fig. 8–9, we present the contour plots for the q-factor of 2-body and 3-body WEC arrays with the second body being moved in the first quadrant. Note that for the 3-body case, the third body is also moved accordingly to maintain the layout symmetry about the x-axis. In both figures, the surrogate approximated q-factor exhibits a very good correlation with the exact one. PA is also able to capture the rough landscape, but it shows substantial deviations from the exact q-factor even when the second body is far away from the origin (e.g., kl12 > 7.5), despite these correspond to regions for which the small-body assumption of PA is valid and thus PA should facilitate higher accuracy [20]. The numerical setup for the 3-body WEC layout problem is also leveraged to justify our preference on approximating hydrodynamic interaction coefficients rather than the q-factor

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

17

Fig. 9. The q-factor versus the second body’s location for 3-WEC arrays (the third body moves symmetrically with respect to x-axis) estimated by different numerical schemes.

Fig. 10. The damping coefficient b¯ 12 versus the second body’s location for 3-WEC arrays (the third body moves symmetrically).

for the optimization problem, even though the latter quantity is the ultimate QoI. Intuitively, it should be less challenging for surrogate models to achieve an accurate approximation on the former quantity, as the latter involves another layer of nonlinear transformation [essentially Eq. (20)]. To visually demonstrate this, Fig. 10 shows the variation of a particular damping coefficient (b¯ 12 ) with the location of the second body for the 3-body WEC array, similar to the plot for the q-factor in Fig. 9. Comparing the variation trends between Fig. 9 (q-factor) and Fig. 10 (hydrodynamic characteristics), it is evident that the damping coefficient variation follows a simpler pattern, which should be easier to approximate by a surrogate model. This benefit of reduced approximation difficulty level does not come free: as the q-factor is a scalar quantity, working directly with the q-factor will simplify and speed up the surrogate model development and prediction process. An interesting question would be how to make the choice between approximating the hydrodynamic coefficients or the q-factor. Our recommendation is to approximate the hydrodynamic interaction coefficients when using expensive numerical solvers, as the training data efficiency will be a priority and the reduced approximation difficulty is definitely needed; for cases where the solver is relatively cheaper to evaluate, it might be a good practice to start with approximating the q-factor, and observe if the surrogate model can reach the desired accuracy level. If not, one can move to approximating hydrodynamic interaction coefficients. Note that another drawback for working with the q-factor is that MBE will lose its

18

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

Table 8 WEC layout performance for 5-body array. Optimum identified by

PA Hierarchical Exact

Corresponding q-factor estimated by PA

Hierarchical

Exact

2.777 2.754 2.689

2.702 2.729 2.717

2.662 2.732 2.753

Fig. 11. Optimal layouts under the three different q-factor estimation schemes.

physics-based interpretability since the controller is dependent on the entire array, and thus the MBE convergence rate is also expected to be affected (i.e., approximating higher-order clusters might be needed to achieve the same accuracy level). This also provides some additional incentives to work with the hydrodynamic interaction coefficients, as advocated in this work. Lastly, the classical problem to find the optimal layout maximizing the q-factor for a 5-body WEC array in [20] is also considered. WECs are symmetrically placed with respect to the x-axis but beyond that no other constraints are considered. The optimization is conducted by a pattern search algorithm [72], starting from the global optimum for the PA q-factor, which should be a well-educated initial guess for such highly multi-modal problem. Optimization algorithm is set to terminate after 10000 function evaluations. The same procedure is repeated to find the best q-factor using either the hierarchical surrogate-approximated or the exact hydrodynamic coefficients. The results and optimal layouts found by different approaches are presented in Table 8 and Fig. 11. Fig. 11 shows the optimal locations of the WECs identified by the three different numerical schemes, PA, Hierarchical and Exact, whereas Table 8 presents the q-factor for the corresponding designs, estimated by all three numerical techniques. The ground truth performances should always be referred to the ones given by the exact solver (last column of Table 8). The proposed hierarchical surrogate modeling approach leads to better solutions than point absorber approximation: it identifies designs with better performance as evident in Table 8 comparison and reaches closer to the actual optimum in the design space as evident in Fig. 11. Furthermore, for all listed layouts in Table 8, the surrogate approximated q-factors are always more accurate than the PA ones (closer to the exact values contained in the last column at Table 8). This discussion demonstrates the superiority of the proposed surrogate model approximation over PA for the WEC layout design problem. Interestingly, although q-factors can noticeably differ from each other in Table 8, all optimal layouts look very similar in Fig. 11. This is expected for the considered problem, which has been confirmed to have multiple local optima associated with sharp peaks, and so large sensitivity can exist around these optima [20,71]. Although this section is dedicated to showcasing the accuracy and solution quality improvement of employing hierarchical surrogate models over PA, we also want to emphasize another significant, and perhaps more important benefit, that is the expanded application scope. This includes, but is not limited, to working with closely spaced WEC arrays [71], or considering more practical control strategies in power absorption [18], as no assumption constraint exist within the proposed approach for evaluating the optimal power (hydrodynamic characteristics are directly utilized). Surrogate model approximations can enable more thorough design explorations in such scenarios, which might have been limited by computational resources in the past due to the high cost of exact numerical solvers and the inapplicability of the PA approximation. 6. Conclusions In this paper, a novel surrogate-model based efficient solution methodology was established for computing hydrodynamic characteristics of waves interacting with an array of multiple identical bodies (wave diffraction and wave radiation).

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

19

The foundation of the approach is an innovative application of the many-body expansion (MBE) principle to overcome the curse of dimensionality through a hierarchical decomposition. Instead of a single surrogate model to predict the interested characteristics involving the whole multi-body configuration (and therefore a high dimensional input and output), multiple surrogate models for clusters with fewer bodies are established and combined to form the predictions for the whole system. The lower-order surrogate models, involving fewer number of bodies, can be developed at a substantial smaller computational cost due to their lower dimensionality, especially for the first components of the many-body expansion, which are most influential to the total hydrodynamic characteristics and the overall prediction quality. Additional enhancements for the surrogate modeling implementation pertain to the characterization of the surrogate model input and output, exploiting symmetry and invariance principles for the hydrodynamic problem. A data-driven approach, leveraging cross-validation accuracy information, was also discussed for guiding this selection, and it was demonstrated to align well with the physics of the problem. Theoretical comparisons were also offered between the adopted MBE scheme and the multi-scattering (MS) approach, which is a popular numerical solution for estimation of hydrodynamic characteristics. It was argued that convergence of the MBE scheme (with respect to the maximum cluster size considered) is faster than the MS approach (with respect to the scattering order), indicating that the MBE-based surrogate model implementation can be truncated after few terms, based on the well-established convergence properties of MS in the literature for the application considered. The hierarchical interaction decomposition surrogate model was formulated within a general setting, while specific details were provided for applications that consider MBE for up to three bodies. Although all illustrations are focused on the use of Kriging as surrogate model for predicting the hydrodynamic characteristics for vertical axisymmetric bodies, the overall philosophy should be extendable to a wide class of multi-body interaction problems or the use of other surrogate models. Several numerical studies were presented to demonstrate the computation efficiency, prediction accuracy as well as the scalability (to arrays sizes) of the proposed method. For the lower-order body contributions, it was shown that excellent surrogate model prediction accuracy can be established across various wave conditions. The wave diffraction problem was demonstrated to be more challenging than the wave radiation problem, indicating that greater care should be placed at the design of experiments (informing the surrogate model development) for that problem. Selection of the surrogate model input and output is also important and can offer substantial improvement in prediction quality. For multi-body systems, the proposed scheme was shown to yield high quality predictions for all examined hydrodynamic characteristics. In particular, the MBE decomposition was demonstrated to be absolutely indispensable, as the direct surrogate model implementation (without MBE) provides poor predictions even for clusters with small number of bodies. Comparing to MS solvers of lower numerical fidelity, the proposed approach can achieve better or comparable accuracy with orders-of-magnitude saving on the computational effort. The application of the hierarchical interaction decomposition surrogate modeling scheme within the context of layout optimization wave energy converters (WECs) was also examined. Results showed that the proposed method produces nearly identical results to the exact solver, and noticeably outperforms the widely used point absorber approximation. This promising result further suggests that the established surrogate model approximation can enable more thorough design explorations for WEC layouts, which might have been limited by computational resources in the past due to the high cost of exact numerical solvers. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This research effort was supported by the National Science Foundation (NSF) under Grants No. CBET-1235768 and CBET1235732. This support is gratefully acknowledged. The statements in this paper reflect the views of the authors and do not necessarily reflect those of the sponsors. The first author gratefully acknowledges the financial support from the Patrick and Jana Eilers Graduate Student Fellowship for Energy Related Research at the University of Notre Dame. Appendix A. Multi scattering solution MS is an iterative, multi-level scheme to calculate the total potential field, first proposed in [10] for acoustic wave and later extended to water wave problems in [13]. Its main idea is to decompose the total interaction to a series of consecutive scattering events, and then reach the exact solution by iteratively solving this sequence. At each iteration, the diffraction and radiation problems are solved for each isolated body considering added potential in the previous iterations from all remaining bodies. For the diffraction problem, the solution procedure starts by considering each isolated body within the array under the incident wave. In response to the incident wave potential φ I , the pth body emits a wave field 1 p φ S as the first-order scattering potential, obtained by solving the Laplace’s equation for the first order diffraction potential p 1 p φ D = φ I + 1 φ S considering its prescribed boundary conditions. The summation of all first order scattering waves from

N 1 q q = p φ S corresponds to the next-level ‘incident wave’ to the pth body. Solving the diffraction problem  p p N 1 q 2 p for φ D = ( q = p φ S ) + 2 φ S yields the second order scattering wave 2 φ S . Such procedure moves similarly forward to the remaining bodies

20

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

mth order scattering, while the potential is anticipated to decrease with increasing m due to fluid-body interactions, leading to convergence to the exact potential solution around the pth body, ultimately given by p

φD =

∞ 

∞ 

p

(m φ I ) +

m =1

p

(m φ S )

(A.1)

m =1 p

Knowing this potential, forces on the pth body in the ith direction F i can be estimated. Radiation is similarly treated, with the primary difference being how the first-order scattering solution is obtained. In this case, the pth body is initially considered to undergo a forced oscillation along the jth mode of motion. This creates pp the zero-order radiation 0 φ B j , estimated by solving Laplace’s equation with appropriate boundary conditions to account for pq

its motion. The potential radiated on the pth body by the qth body, 1 φ I j , represents a first order excitation, for which the 1

pth body radiates a first order scattering

pq φB j

in response. Laplace’s equation and respective boundary conditions need to pq

pq

pq

be satisfied for the corresponding total first order potential 1 φ R j = 1 φ I j + 1 φ B j around the pth body due to the motion of pp pp the qth (q = p) body. Note that 1 φ I j = 1 φ B j = 0. N pq lq second-order excitation 2 φlj = l = p 2 φ B j for the

Subsequently, the first-order scattering from remaining bodies acts as a pq

pth body, which in response radiates the second-order scattering 2 φ B j . The summation of all order of scattering up to infinity finally yields: pq

pp

φ R j = δ pq · 0 φ B j +

∞ 

m

pq

pq

φI j + m φI j

(A.2)

m =1 pq

pq

where δ pq stands for the Kronecker delta function. This potential determines the added mass ai j and added damping b i j coefficients for the pth body in the ith direction due to the interaction with motion of the qth body in the jth direction, as discussed in the previous section. This MS approach represents an efficient iterative methodology for estimating array hydrodynamic characteristics. The key computational component is the solution of the diffraction and radiation problem for each isolated body at each individual scattering order m. For vertical axisymmetric bodies as examined in the illustrative application, this can be performed using the matched axisymmetric eigenfunction expansion approach as proposed first in [6] for the diffraction problem and later in [47] for the radiation problem. This approach expresses all potentials as a product of scattering coefficients and partial waves involving Hankel and Bessel’s functions of first and second kind, and then uses the boundary conditions to determine unknown scattering coefficients by solving a linear system of equations. The matrix inversion involved in this solution represents ultimately the computational costly part of MS, whose complexity increases with the truncation order (i.e., the number of retained terms) for the eigenfunction expansion. In practical applications, the scattering order for MS is truncated to some finite value (denoted by N m ), which introduces numerical error. The computational cost increases linearly with the max scattering order (N m ) and quadratically with the array size (N). Therefore, albeit more efficient than FEM or BEM based solvers [21], MS still becomes expensive to evaluate when facing a large array or slow convergence over the scattering order. This further hinders its application for analyses that require multiple evaluations of the hydrodynamic characteristics, like the WEC layout optimization problem [26]. Appendix B. Details for surrogate models for vertical axisymmetric bodies considering two and three body clusters This appendix discusses surrogate model implementation details for arrays of vertical axisymmetric bodies with MBE expansion up to the third order. Without loss of generality, we assume that the wave propagates along the positive x-axis. Other configurations can be mapped by rotation. The diffraction potential field is invariant under y-axis translation (i.e., response unchanged if adding a constant to the y coordinate of all bodies), symmetric with respect to x-axis (i.e., response unchanged if reversing the y coordinate sign of all bodies), whereas translation along the x-axis by D x of the array simply creates a phase shift of exp(i · k · D x ) on all diffraction forces. The radiation potential field is invariant under either rotation or translation of the coordinate system. For further convenience, we denote the pairwise distance and angle between the pth body and qth body by

l pq =



x2pq + y 2pq and θ pq = atan2

y pq x pq

(B.1)

where atan2 is the four quadrant inverse tangent function, with xqp = xq − x p and y qp = y q − y p . To avoid overlapping of bodies only configurations with l pq ≥ lmin are considered, with l min representing the body width (equal to the diameter for cylindrical bodies). Fig. 12 reviews the geometric characteristics that will be considered as surrogate model input for the two and three body clusters. For an isolated body, w = [x y ] T , the radiation output is independent of its coordinates:

a(w) = a([0 0] T ), b(w) = b([0 0] T )

(B.2)

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

21

Fig. 12. Geometric characteristics used as surrogate model input for the two and three body clusters.

whereas the diffraction output obeys the following relationship:

F (w) = F ([0 0] T ) exp(i · k · x)

(B.3)

Therefore, the single-body hydrodynamic problem needs to be numerically solved only once for the body at the origin using any desired numerical approach, for example via the matched axisymmetric eigenfunction expansion discussed in Appendix A. In the following surrogate model development for two-body systems (denoted by [w1 w2 ]), for the sake of brevity, we shall frequently drop the dependency of hydrodynamic characteristics (i.e., F, a and b) and their additive effects on such two bodies when there is no ambiguity, such that F  F([w1 w2 ]) and similarly a  a([w1 w2 ]). For the diffraction output di f v2 , we choose to describe the additive effect on the first body’s wave excitation force F 1 , written as:

F 1 = F 1 − F (w1 )

(B.4) 2

This characterization is sufficient because the other required force-related element F can be predicted by flipping the order of bodies. A further phase-shift by exp (−i · k · x1 ) is placed on this diffraction additive effect, in order to make output independent of the absolute position of the two bodies and thus dependent only on the relative position with respect to the first body:

F¯ 1 = F 1 · exp (−i · k · x1 )

(B.5)

This ultimately allows us to reduce the input dimensionality for the surrogate model. Since F¯ 1 is a complex quantity, we need to include its real Re and imaginary Im part as two surrogate output components:



di f

v2 =

Re( F¯ 1 ) Im( F¯ 1 )

 (B.6)

Instead of separately approximating such two components, one can also treat them as a single vectorized output component and formulate a single vectorial Kriging model. This strategy has been shown advantageous in a similar context (interatomic force field prediction) [73], and we plan to investigate its use in the future. As for radiation, the following output definition is exhaustive since a11 = a22 and b11 = b22 :

⎤ ⎡ 11 ⎤ a11 a − a(w1 ) ⎥ ⎢ a12 ⎥ ⎢ a12 ⎥ ⎢ ⎥ ⎢ vrad 2 = ⎣ b 11 ⎦ = ⎣ b 11 − b (w ) ⎦ 1 b12 b12 ⎡

(B.7)

Exploiting existing symmetries and invariances, the two-body system radiation solution depends strictly on the two-body pairwise distance, whereas the diffraction solution depends additionally on the pairwise angle. Based on the symmetry with respect to the x-axis, this angle only needs to be examined in the range of [0, π ) for the surrogate model development since any angle in [−π , π ) can be mapped to that solution. This ultimately means that for the two-body system, one can express the surrogate inputs with respect to the relative placement of the second body (in other words using the local coordinate system around the first body), and that for radiation problem only one half of the domain needs to be considered. The diffraction surrogate model input along with its corresponding transformation is defined as: di f u2

 =

l12

θ12



 ;

di f g2 ([w1

w2 ]) =

 x221 + y 221 atan2

( y 21 )·sgn( y 21 ) x21

 (B.8)

22

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

where sgn(.) is the sign function, introduced to constrain the range of the surrogate model input to [0, π ). The radiation input selection further simplifies to:

u rad 2 = [l12 ] ;

g 2rad ([w1 w2 ]) =

 x221 + y 221

(B.9)

For defining the input domain U 2c , the range of the distance component l12 is lower-bounded by lmin to avoid overlapping bodies, and upper-bounded by the max distances between bodies within W (denoted as lmax ) to avoid extrapolation, whereas the range for θ12 will be [0, π ), as has been discussed above. di f Finally, given the surrogate output v2 , its transformation to the first-body additive force effect for an arbitrary two-body 1 clusters F is described as: di f



T2



di f

di f





di f

di f





di f

v2 , [w1 w2 ]  v 2,1 g2 ([w1 w2 ]) + i · v 2,2 g2 ([w1 w2 ]) di f

· exp (i · k · x1 )

(B.10)

di f

where the notation v 2,r is used to describe the rth component of vector v2 . For radiation coefficients, the transformation is simply identity:









rad Trad vrad grad 2 2 , [w1 w2 ]  v2 2 ([w1 w2 ])

(B.11)

The training output estimation (Step D.3 in Section 3.4.2) for two-body clusters does not entail any output from lower-order clusters as shown in Eq. (B.6) and (B.7). Ultimately, the two-body surrogate models approximate the following input/output mappings:



di f

u2 =



l12

θ12

∈ R+ ⊗ [0, π ) → v2 = di f



Re( F¯ 1 ) Im( F¯ 1 )

⎤ a11 ⎢ a12 ⎥ + rad 4 ⎥ ⎢ u rad 2 = [l12 ] ∈ R → v2 = ⎣ b 11 ⎦ ∈ R b12 ⎡

 ∈ R2

(B.12)

For the three-body surrogate model involving ([w1 w2 w3 ]), we again simplify the dependency of hydrodynamic characteristics on the three bodies, similarly to the treatment in the two-body case. The surrogate outputs need to include one force element for diffraction, plus one diagonal and off-diagonal elements in the added mass and damping matrices for radiation. The former is chosen similarly as the two-body case, i.e., the real/imaginary part of the phase-shifted additive effect on the first body’s wave excitation force:



di f v3

=









Re F 1 · exp (−i · k · x1 )

Im F 1 · exp (−i · k · x1 )

1

1

1

(B.13)

1

F = F − F ([w1 w2 ]) − F ([w1 w3 ]) − F (w1 ) For radiation coefficients, the choice is not unique. The off-diagonal element [.]23 is used as it is most compatible with our selected input based on parametric studies detailed in Section 4. This leads to the following radiation output definition:

⎤ ⎤ ⎡ 11 a11 a − a11 ([w1 w2 ]) − a11 ([w1 w3 ]) − a(w1 ) ⎥ ⎢ a23 ⎥ ⎢ a23 − a12 ([w2 w3 ]) ⎥ ⎥ ⎢ ⎢ vrad 3 = ⎣ b 11 ⎦ = ⎣ b 11 − b 11 ([w w ]) − b 11 ([w w ]) − b (w ) ⎦ 1 2 1 3 1 b23 b23 − b12 ([w2 w3 ]) ⎡

(B.14)

Remaining hydrodynamic characteristics for F, a and b of the three-body cluster are obtained by successively permuting the order of the second, or third body to become the first one. The three-body surrogate input vector is expanded from the two-body input in Eq. (B.8) and Eq. (B.9), with additional information about the third body. For diffraction, the third body’s distance and angle with regard to the first body are added, extending the principles for the two-body cluster. No restriction is adopted for the angle of the third body (all symmetries/invariance have been exploited for the second body), leading to the following ranges for added inputs: l13 ∈ [lmin , lmax ] and θ13 ∈ [−π , π ). The surrogate model input along with its transformation is therefore:

⎡ di f u3

l12

⎡

⎤

⎢ θ12 ⎥ ⎥ =⎢ ⎣ l13 ⎦ ; θ13

di f g3 ([w1



x221 + y 221

⎢ ⎢ atan2 ( y 21 )·sgn( y 21 ) x21 ⎢ w2 w3 ]) = ⎢  ⎢ x231 + y 231 ⎣ atan2

y 31 ·sgn( y 21 ) x31

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

where sgn( y 21 ) is again used to accommodate the restriction for the placement of the second body.

(B.15)

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

23

For radiation, the third body’s distance from the first body, as well as the angle between vectors connecting the first body to the other two bodies (θ213 ) are chosen to expand the previous input urad 2 . Alternative input choices are compared in the parametric study in Section 4. Considering the rotational invariance, only the range of [0, π ) needs to be considered for θ213 in the surrogate model development. Without loss of generality, we can always take the second body as the nearest neighbor to the first when obtaining the input. The surrogate model input along with its transformation is therefore:

⎡

⎡ urad 3

⎤

l12

= ⎣ l13 ⎦ ; θ213

grad 3 ([w1

⎢ ⎢ ⎢ w2 w3 ]) = ⎢ ⎢ ⎣

 arccos

⎤

 x2 + y 221  21 x231 + y 231

y 21 y 31  x21 x31 + x221 + y 221 x231 + y 231

⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(B.16)

with the additional constraint l12 ≤ l13 . Evidently according to Eq. (B.16), this input is invariant with respect to the ordering of the second and third body, which allows the introduction of the additional constraint for l12 ≤ l13 . The output transformations are given as follows, analogously to the two-body ones:



di f

T3



Trad 3



di f

di f





di f

di f





di f

v 3 , [w1 w2 w3 ] = v 3,1 g3 ([w1 w2 w3 ]) + i · v 3,2 g3 ([w1 w2 w3 ])

   rad vrad grad 3 , [w1 w2 w3 ] = v3 3 ([w1 w2 w3 ])

· exp (i · k · x1 )

(B.17) (B.18)

The training output estimation (Step D.3 in Section 3.4.2) now involves output from their two-body subsets as shown in Eq. (B.13) and Eq. (B.14). Providing approximations to such quantities using the developed two-body surrogate models is straightforward (discussed in Step D.3 in Section 3.4.2), since the two-body input vector (i.e., uc2 ) has always been included in uc3 from the advocated expansion philosophy at Step D.1 in Section 3.4.2. Ultimately, the three-body surrogate models approximate the following input/output mappings

⎡

l12

⎤

  ⎢ θ12 ⎥ Re( F¯ 1 ) di f ⎥ ⎢ ∈ R2 =⎣ ∈ [lmin , lmax ] ⊗ [0, π ) ⊗ [lmin , lmax ] ⊗ [−π , π ) → v3 = l13 ⎦ Im( F¯ 1 ) θ13 ⎤ ⎡ ⎤ ⎡ a11 l12 23 ⎢ a ⎥ 4 ⎥ ⎢ ⎣ l13 ⎦ ∈ [lmin , lmax ] ⊗ [lmin , lmax ] ⊗ [0, π ) → vrad urad 3 = 3 = ⎣ b 11 ⎦ ∈ R θ213 b23 di f u3

(B.19)

It should be stressed that none of the above surrogate input or output choices are unique, and that any input modification would have also altered the necessary surrogate input domain to avoid extrapolations when predicting hydrodynamic characteristics for arbitrary multi-body configurations in W . Using MBE expansion up to N b = 3 and the discussed above surrogate models, the prediction of Eq. (11) for the radiation force on body p and the added damping coefficients involving this body (similar equations apply for added mass coefficients) are simplified to:

 F p (W) = F ([0 0] T ) +

N



di f





di f

di f



r =1,r = p N −1 

N −r 

di f





di f

v 2,1 g2 ([w p wr ]) + i · v 2,2 g2 ([w p wr ])



di f

di f



 

di f

v 3,1 g3 ([w p wr ws ]) + i · v 3,2 g3 ([w p wr ws ])

+

  · exp i · k · x p

r =1,r = p s>r

b

pp

T

(W) = b([0 0] ) +

N 

v rad 2,3



grad 2 ([w p

N −r N −1      di f wr ]) + v rad 3,3 g3 ([w p wr ws ])

r =1,r = p

b pq (W) = v rad 2,4



 grad 2 ([w p wq ]) +

(B.20)

r =1,r = p s>r N 



di f



v rad 3,4 g3 ([wr w p wq ])

r =1,r = p ,q

References [1] [2] [3] [4]

J. Dubinski, A parallel tree code, New Astron. 2 (1996) 133–147. A. Leonard, Vortex methods for flow simulation, J. Comput. Phys. 37 (1980) 289–335. T. Lee, C. Yang, Many-body problem in quantum mechanics and quantum statistical mechanics, Phys. Rev. 105 (1957) 1119. K.M. Watson, Multiple scattering and the many-body problem—applications to photomeson production in complex nuclei, Phys. Rev. 89 (1953) 575.

24

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

[5] M. Gordon, D. Fedorov, S. Pruitt, L. Slipchenko, Fragmentation methods: a route to accurate calculations on large systems, Chem. Rev. 112 (2012) 632–672. [6] S. Mavrakos, P. Koumoutsakos, Hydrodynamic interaction among vertical axisymmetric bodies restrained in waves, Appl. Ocean Res. 9 (1987) 128–140. [7] Y.M. Marzouk, A.F. Ghoniem, K-means clustering for optimal partitioning and dynamic load balancing of parallel hierarchical N-body simulations, J. Comput. Phys. 207 (2005) 493–528. [8] L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (1987) 325–348. [9] P.A. Martin, Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles, Cambridge University Press, 2006. [10] V. Twersky, Multiple scattering of radiation by an arbitrary configuration of parallel cylinders, J. Acoust. Soc. Am. 24 (1952) 42–46. [11] Y.-l. Xu, Electromagnetic scattering by an aggregate of spheres, Appl. Opt. 34 (1995) 4573–4588. [12] R.-S. Wu, Multiple scattering and energy transfer of seismic waves—separation of scattering effect from intrinsic attenuation—I. Theoretical modelling, Geophys. J. Int. 82 (1985) 57–80. [13] M. Ohkusu, Wave action on groups of vertical circular cylinders, J. Soc. Nav. Archit. Jpn. 1972 (1972) 53–64. [14] H. Kagemoto, D.K. Yue, Interactions among multiple three-dimensional bodies in water waves: an exact algebraic method, J. Fluid Mech. 166 (1986) 189–209. [15] M.A. Peter, M.H. Meylan, C. Linton, Water-wave scattering by a periodic array of arbitrary bodies, J. Fluid Mech. 548 (2006) 237–256. [16] A. Babarit, J. Hals, M. Muliawan, A. Kurniawan, T. Moan, J. Krokstad, Numerical benchmarking study of a selection of wave energy converters, Renew. Energy 41 (2012) 44–63. [17] B. Czech, P. Bauer, Wave energy converter concepts: design challenges and classification, IEEE Ind. Electron. Mag. 6 (2012) 4–16. [18] J. Scruggs, P. Jacob, Harvesting ocean wave energy, Science 323 (2009) 1176–1178. [19] J. Falnes, Ocean Waves and Oscillating Systems: Linear Interactions Including Wave-Energy Extraction, Cambridge University Press, 2002. [20] C. Fitzgerald, G. Thomas, A preliminary study on the optimal formation of an array of wave power devices, in: The 7th European Wave and Tidal Energy Conference, Porto, Portugal, 2007, pp. 11–14. [21] M. Folley, Numerical Modelling of Wave Energy Converters: State-of-the-Art Techniques for Single Devices and Arrays, Academic Press, 2016. [22] K. Budal, Theory for absorption of wave power by a system of interacting bodies, J. Ship Res. 21 (1977). [23] J. Falnes, Radiation impedance matrix and optimum power absorption for interacting oscillators in surface waves, Appl. Ocean Res. 2 (1980) 75–80. [24] P. McIver, D. Evans, Approximation of wave forces on cylinder arrays, Appl. Ocean Res. 6 (1984) 101–107. [25] M. Simon, Multiple scattering in arrays of axisymmetric wave-energy devices. Part 1. A matrix method using a plane-wave approximation, J. Fluid Mech. 120 (1982) 1–25. [26] B. Child, V. Venugopal, Optimal configurations of wave energy device arrays, Ocean Eng. 37 (2010) 1402–1417. [27] J.-R. Nader, S.-P. Zhu, P. Cooper, B. Stappenbelt, A finite-element study of the efficiency of arrays of oscillating water column wave energy converters, Ocean Eng. 43 (2012) 72–81. [28] A. Babarit, G. Delhommeau, Theoretical and numerical aspects of the open source BEM solver NEMOH, in: 11th European Wave and Tidal Energy Conference (EWTEC2015), 2015. [29] D. Sarkar, E. Contal, N. Vayatis, F. Dias, Prediction and optimization of wave energy converter arrays using a machine learning approach, Renew. Energy 97 (2016) 504–517. [30] E. Darve, The fast multipole method: numerical implementation, J. Comput. Phys. 160 (2000) 195–240. [31] F.A. Amirkulova, A.N. Norris, Acoustic multiple scattering using recursive algorithms, J. Comput. Phys. 299 (2015) 787–803. [32] M. Ganesh, J. Hesthaven, B. Stamm, A reduced basis method for electromagnetic scattering by multiple particles in three dimensions, J. Comput. Phys. 23 (2012) 7756–7779. [33] A. Vion, R.V. Sabariego, C. Geuzaine, A model reduction algorithm for solving multiple scattering problems using iterative methods, IEEE Trans. Magn. 47 (2011) 1470–1473. [34] M. Murai, H. Kagemoto, M. Fujino, On the hydroelastic responses of a very large floating structure in waves, J. Mar. Sci. Technol. 3 (1999) 123–153. [35] M. Kashiwagi, Hydrodynamic interactions among a great number of columns supporting a very large flexible structure, J. Fluids Struct. 14 (2000) 1013–1034. [36] Q. Zhong, R. Yeung, Wave-body interactions among energy absorbers in a wave farm, Appl. Energy 233 (2019) 1051–1064. [37] M. Göteman, J. Engström, M. Eriksson, J. Isberg, Fast modeling of large wave energy farms using interaction distance cut-off, Energies 8 (2015) 13741–13757. [38] M. Eldred, A. Giunta, S. Collis, Second-order corrections for surrogate-based optimization with model hierarchies, in: 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, 2004, p. 4457. [39] G. Jia, A. Taflanidis, J. Scruggs, Layout optimization of wave energy converters in a random sea, in: The Twenty-fifth International Ocean and Polar Engineering Conference, International Society of Offshore and Polar Engineers, 2015. [40] R. Jin, W. Chen, T.W. Simpson, Comparative studies of metamodelling techniques under multiple modelling criteria, Struct. Multidiscip. Optim. 23 (2001) 1–13. [41] K. Yao, J.E. Herr, J. Parkhill, The many-body expansion combined with neural networks, J. Chem. Phys. 146 (2017) 014106. [42] A.P. Bartók, M.C. Payne, R. Kondor, G. Csányi, Gaussian approximation potentials: the accuracy of quantum mechanics, without the electrons, Phys. Rev. Lett. 104 (2010) 136403. [43] A. Glielmo, C. Zeni, A. De Vita, Efficient nonparametric n-body force fields from machine learning, Phys. Rev. B 97 (2018) 184307. [44] J. Ling, R. Jones, J. Templeton, Machine learning strategies for systems with invariance properties, J. Comput. Phys. 318 (2016) 22–35. [45] A.P. Bartók, R. Kondor, G. Csányi, On representing chemical environments, Phys. Rev. B 87 (2013) 184115. [46] J. Sacks, W.J. Welch, T.J. Mitchell, H.P. Wynn, Design and analysis of computer experiments, Stat. Sci. (1989) 409–423. [47] S. Mavrakos, Hydrodynamic coefficients for groups of interacting vertical axisymmetric bodies, Ocean Eng. 18 (1991) 485–515. [48] J. Newman, Wave effects on multiple bodies, in: Hydrodynamics in Ship and Ocean Engineering, RIAM, Kyushu University, Japan, 2001. [49] J.P.C. Kleijnen, Design and Analysis of Simulation Experiments, Springer, 2008. [50] M.C. Kennedy, A. O’Hagan, Bayesian calibration of computer models, J. R. Stat. Soc., Ser. B, Stat. Methodol. 63 (2001) 425–464. [51] S.N. Lophaven, H.B. Nielsen, J. Sondergaard, Aspects of the MATLAB toolbox DACE, in: Informatics and Mathematical Modelling Report IMM-REP2002-13, Technical University of Denmark, 2002. [52] F. Bachoc, Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model misspecification, Comput. Stat. Data Anal. 66 (2013) 55–69. [53] G. Jia, A.A. Taflanidis, Kriging metamodeling for approximation of high-dimensional wave and surge responses in real-time storm/hurricane risk assessment, Comput. Methods Appl. Mech. Eng. 261–262 (2013) 24–38. [54] F.A. Viana, C. Gogu, R.T. Haftka, Making the most out of surrogate models: tricks of the trade, in: ASME 2010 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, 2010, pp. 587–598. [55] X. Ma, N. Zabaras, An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations, J. Comput. Phys. 229 (2010) 3884–3915.

J. Zhang et al. / Journal of Computational Physics 408 (2020) 109298

25

[56] H. Rabitz, Ö.F. Alis, ¸ General foundations of high-dimensional model representations, J. Math. Chem. 25 (1999) 197–233. [57] I.M. Sobol, Sensitivity estimates for nonlinear mathematical models, in: Mathematical Modelling and Computational Experiments, vol. 1, 1993, pp. 407–414. [58] G. Blatman, B. Sudret, An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis, Probab. Eng. Mech. 25 (2010) 183–197. [59] K. Tang, J.M. Wang, J.B. Freund, Adaptive sparse polynomial dimensional decomposition for derivative-based sensitivity, J. Comput. Phys. 391 (2019) 303–321. [60] A.F. Cortesi, G. Jannoun, P.M. Congedo, Kriging-sparse Polynomial Dimensional Decomposition surrogate model with adaptive refinement, J. Comput. Phys. 380 (2019) 212–242. [61] B. Paulus, K. Rosciszewski, N. Gaston, P. Schwerdtfeger, H. Stoll, Convergence of the ab initio many-body expansion for the cohesive energy of solid mercury, Phys. Rev. B 70 (2004) 165106. [62] A. Marrel, B. Iooss, F. Van Dorpe, E. Volkova, An efficient methodology for modeling complex computer codes with Gaussian processes, Comput. Stat. Data Anal. 52 (2008) 4731–4744. [63] P. Angelikopoulos, C. Papadimitriou, P. Koumoutsakos, X-TMCMC: adaptive kriging for Bayesian inverse modeling, Comput. Methods Appl. Mech. Eng. 289 (2015) 409–428. [64] C. Rasmussen, Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian integrals, in: Bayesian Statistics 7: Proceedings of the 7th Valencia International Meeting, Oxford University Press, 2003, pp. 651–659. [65] Y. Bengio, O. Delalleau, N.L. Roux, The curse of highly variable functions for local kernel machines, in: Advances in Neural Information Processing Systems, 2006, pp. 107–114. [66] M. Kadiyala, J. Jones, R. Mylavarapu, Y. Li, M. Reddy, Identifying irrigation and nitrogen best management practices for aerobic rice–maize cropping system for semi-arid tropics using CERES-rice and maize models, Agric. Water Manag. 149 (2015) 23–32. [67] S. Mavrakos, A. Kalofonos, Power absorption by arrays of interacting vertical axisymmetric wave-energy devices, J. Offshore Mech. Arct. Eng. 119 (1997) 244–251. [68] C. Fitzgerald, G. Thomas, A preliminary study on the optimal formation of an array of wave power devices, in: The 7th European Wave and Tidal Energy Conference, Porto, Portugal, 2007. [69] A. Georgiadis, Gain, phase imbalance, and phase noise effects on error vector magnitude, IEEE Trans. Veh. Technol. 53 (2004) 443–449. [70] J.T. Scruggs, S.M. Lattanzio, A.A. Taflanidis, I.L. Cassidy, Optimal causal control of a wave energy converter in a random sea, Appl. Ocean Res. 42 (2013) 1–15. [71] J.P. McGuinness, G. Thomas, Hydrodynamic optimisation of small arrays of heaving point absorbers, J. Ocean Eng. Mar. Energy 2 (2016) 439–457. [72] T.G. Kolda, R.M. Lewis, V. Torczon, Optimization by direct search: new perspectives on some classical and modern methods, SIAM Rev. 45 (2003) 385–482. [73] A. Glielmo, P. Sollich, A. De Vita, Accurate interatomic force fields via machine learning with covariant kernels, Phys. Rev. B 95 (2017) 214302.