1The 36th International Workshop on Water Waves and Floating Bodies, Seoul, Korea, 25-28 April 2021

Computational accuracy and efficiency for Diffraction Transfer
Matrix using hybrid source-dipole formulations
Yingyi Liua †, Hui Liangb , Massashi Kashiwagic and Peiwen Congda

Research Institute for Applied Mechanics, Kyushu University, Fukuoka 8168580, Japan
Technology Centre for Offshore and Marine, Singapore (TCOMS ), 118411, SingaporecDepartment of Naval Architecture & Ocean Engineering, Osaka University, Osaka 5650871, Japand
State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, Chinab

1. INTRODUCTION
The interaction theory presented by Kagemoto & Yue (1986) significantly reduces the computational burden in
the wave interaction problem of multiple surface-piercing bodies, particularly arrays of wave energy converters inrecent years. One of the essential operators is the so-called Diffraction Transfer Matrix (DTM). Many subsequent
researchers (Goo & Yoshida 1990; Flavià et al. 2018) have implemented the theory using the source distribution
method in evaluating the two linear operators of a single unique geometry. However, nowadays, a great majorityof boundary element method codes have been written by virtue of the hybrid source-dipole distribution method(i.e., potential formulation) taking account of its high accuracy (Dai & Duan 2008).
There are basically two ways to solve the scattering potential in the potential formulation: (1) using theconventional method and then obtain the diffraction potential, or (2) solve directly the diffraction potential byapplying incident-wave potential as the forcing term on the right-hand side of the boundary integral equation.
The two ways are similar in wave analysis with a single body (Lee & Newman 2005) in accuracy, but the latter
one runs a bit faster than the first because there is no need to evaluate the normal derivative and its integration
over the wetted surface. However, when being applied to a multi-body problem, since the right-hand side of
the boundary integral equation needs to be evaluated for hundreds and thousands of times, the difference ofcomputation time can be increased substantially. In this regard, it is necessary to check how this difference
could be, such that the equivalence of the two alternative approaches may need to be reappraised.

2. THEORY AND METHODOLOGY
Wave diffraction of arbitrary-geometric three-dimensional floating structures can be solved by the following
hybrid source-dipole boundary integral equation in frequency domain within the framework of potential flow
theory:ZZZZ∂ϕIj (ξ)∂G(x, ξ)SS
ϕj (ξ)
dS = −G(x, ξ)dS with
(1)
2πϕj (x) +
∂nξ
∂nξ
SB j
SB j
where x = (x, y, z) and ξ = (ξ, η, ζ) represent the field point and the source point, respectively; G(x, ξ)is the free-surface Green function, and can be expressed in polar coordinates in terms of x = (x, y, z) =(r cos θ, r sin θ, z) and ξ = (ξ, η, ζ) = (R cos Θ, R sin Θ, ζ) following Fenton (1978), in the form of
)
(
(1)
∞
X
Jm (kR)Hm (kr)eim(θ−Θ)G =2πiC0 cosh k(z + h) cosh k(ζ + h)
(1)
(kR)J
(kr)
H
m
m
m=−∞
(2)


∞
∞
X
X
Jm (kn R)Km (kn r) im(θ−Θ)
Cn cos kn (z + h) cos kn (ζ + h)e.
+4
Km (kn R)Jm (kn r)m=−∞n=1
The upper terms in the curly braces are used when r > R (the region outside of a circular cylinder thatcircumscribes the body or bodies) and the lower terms when r < R.
2.1. Partial waves in cylindrical harmonics
In a finite-sized array of floating bodies, it is convenient to express the velocity potentials as the scalar productbetween a vector of complex coefficients and a vector of partial cylindrical wave component (Flavià et al. 2018):
φIj = {AIj }T {ϕIj },φSj = {ASj }T {ϕSj },(3)

where the superscript T represents the matrix transpose operator, the curly brace {·} stands for a vector andthe subscript j denotes the jth body. {AIj } and {ASj } are the complex incident and scattered vectors of partial
wave coefficients. Indexes (l, q) are associated with incident waves and (n, m) with outgoing waves. The vectorsof the incident and scattered cylindrical functions are respectively expressed as cosh k(zj + h)
Jq (krj )eiqθj l = 0,
{ϕIj }lq =(4)cosh khcos k (z + h)I (k r )eiqθj l > 1,
l

j

q

l j

† Email address for correspondence: liuyingyi@riam.kyushu-u.ac.jp

2
 cosh k(zj + h) (1)Hm (krj )eimθj{ϕSj }nm =cosh khcos k (z + h)K (k r )eimθj
n j
m n j

n = 0,(5)n > 1.

2.2. Method I
Derivation of the DTM of a specific floating body can be started from considering the wave diffraction by asingle body. The scattering potential of a single floating body in a partial incident wave of mode (l, q) without
the presence of other bodies can be expressed ascosh k(zj + h)[ϕSj (rj , θj , zj )]lq =cosh kh

+

∞
X

∞
X

j,lq (1)Hm (krj )eimθj
D0m

m=−∞
∞
X

cos kn (zj + h)

(6)
j,lq
Dnm
Km (kn rj )eimθj ,

m=−∞

n=1

j,lq
j,lq
and Dnmare scattered complex coefficients. On the other hand, the scattering potential at a field
where D0m

point in the fluid domain (other than the body surface) can be determined by the following equation(Z Z
1∂G(rj , θj , zj ; Rj , Θj , ζj )S
[ϕj (rj , θj , zj )]lq = −[ϕSj (Rj , Θj , ζj )]lq
dS
j
4π
∂nξ
SB
)
ZZ
∂[ϕIj (Rj , Θj , ζj )]lq+dS
G(rj , θj , zj ; Rj , Θj , ζj )j∂nξSBSubstituting Eq. (2) into Eq. (7) and comparing Eq. (6) with Eq. (7) yields
#
ZZ "
I
∂(ϕ
)
i
∂
lq
j
j,lq
(ϕSj )lq+[Jm (kRj ) cosh k(ζj + h)e−imΘj ]dS,
D0m
= − C0 cosh kh
j
2
∂n
∂n
SB
1
j,lq
Dnm
= − Cn

"

ZZ

π

j
SB

(ϕSj )lq#∂(ϕIj )lq∂+
[Im (kn Rj ) cosh kn (ζj + h)e−imΘj ]dS,
∂n
∂n

(7)

(8)

(9)

j,lq
j,lq
where D0m
are exactly the elements of DTM. The only unknown in Eqs. (8) and (9) is (ϕSj )lq , which
and Dnmcan be solved by the following boundary integral equation
ZZ∂G(rj , θj , zj ; Rj , Θj , ζj )[ϕSj (Rj , Θj , ζj )]lq
2π[ϕSj (rj , θj , zj )]lq +
dS
j
∂nξ
SB
(10)
ZZ
∂[ϕIj (Rj , Θj , ζj )]lq=−
G(rj , θj , zj ; Rj , Θj , ζj )dS.
j
∂nξ
SB

2.3. Method II
A second way of evaluating the DTM elements is to consider the following boundary integral equation, which
was first proposed by (Kashiwagi & Kohjo 1995):
ZZ∂G(rj , θj , zj ; Rj , Θj , ζj )D
[ϕDC(rj , θj , zj )[ϕj (rj , θj , zj )]lq +dS = 4π[ϕIj (rj , θj , zj )]lq , (11)j (Rj , Θj , ζj )]lqj∂nξ
SB
where (ϕDj )lq is the total diffraction potential in correspondence to the partial incident wave of mode (l, q).
Following a similar process, the DTM elements can be derived as
ZZ
∂
i
j,lq
[(ϕSj )lq + (ϕIj )lq ] [Jm (kRj ) cosh k(ζj + h)e−imΘj ]dS,
D0m
= − C0 cosh kh
(12)
j
2
∂n
SB
1
j,lq
Dnm
= − Cnπ

ZZj
SB[(ϕSj )lq + (ϕIj )lq ]∂[Im (kn Rj ) cosh kn (ζj + h)e−imΘj ]dS,
∂n

(13)

3. RESULTS AND DISCUSSIONSIn order to verify the present method based on the hybrid source-dipole boundary integral equation (Liu2019; Liang et al. 2020), numerical computations are performed against two benchmark problems given in
Flavià et al. (2018). In Figs. 1, ‘Present’ stands for the numerical results generated using ‘Method I’, as themegascopic difference between the results of ‘Method I’ and ‘Method II’ are indistinguishable.
The real and imaginary parts of the DTM terms of the truncated vertical cylinder are shown in Figs. 1.
In general, a good agreement is found between results computed by the methods of McNatt et al. (2015) and
Flavià et al. (2018), respectively. Exceptions occur with Flavià et al. (2018) at some specific frequencies, e.g.,ka = 2.39, 2.57, 2.75. This should be attributed to the ‘irregular frequencies’ phenomenon, as Flavià et al. (2018)
was using the open-source BEM code Nemoh (Babarit & Delhommeau 2015), which at that time did not have

3

Figure 1: DTM progressive terms for a truncated vertical cylinder of 3 m radius, 6 m draft in a 10 m water
depth: (a) real part; (b) imaginary part.

ka

Re q = 0, m = 0
Method I Method II Error

0.6 −0.049621
1.2 −0.3895
1.8 −0.85964
2.4 −0.94447
3.0 −0.50805−0.04971−0.39203−0.87157−0.96528−0.52609

Im q = 0, m = 0
Method I Method II Error−0.18% −0.21695 −0.21741 −0.21%−0.65% −0.48501 −0.48828 −0.67%−1.39% −0.33147 −0.33618 −1.42%−2.20% 0.17974
0.18373 −2.22%−3.55% 0.4822
0.49947 −3.58%

Table 1: Comparison of the first two diagonal terms of the DTM between Method I and Method II, with anincident wave heading angle 0 degree.

the ‘irregular frequencies removal’ functionality (Penalba et al. 2017). Nevertheless, the agreement between thepresent results and those of McNatt et al. (2015) is pretty good with respect to all the frequencies, as the latter
was using WAMIT with the ‘irregular frequencies’ being removed.
The difference of numerical results between the two methods in computation of the DTM of the truncatedvertical cylinder is assesed in Table 1. Only the first two diagonal terms with q = 0, m = 0 and q = 1, m = 1
00
01(i.e., D00
and D01) are compared since they are the largest terms of the DTM. In general, the absolute value of
Method I’s result is less than that of Method II. The absolute relative error increases against the wave number
ka but is still acceptable around ka = 3.0 (ω = 3.132 rad/s). As it is known that the upper limit (at which thespectrum approaches zero) of a typical wave spectrum can normally be below ω = 2.5 rad/s, an absolute relative
error of the DTM term below 5.0% means that using either of the methods can be fine for the calculation.
Fig. 2 shows the CPU time of per evaluation of the DTM of the floating square box at different stages, on adesktop machine, with the aid of OpenMP parallelism on eight threads. In Fig. 2, Eq. (6) is truncated with amaximum depth mode of NMax = 5 and a maximum angular mode of MMax = 5, indicating that the numberof diffraction problems to be solved at each frequency is (NMax + 1)(2MMax + 1)NWaveHeadings = 2376 and thenumber of the diffraction transfer matrices to be calculated is NWaveHeadings = 36. Assembling the right-handside matrix takes up most of the time (93.90% of the overall computation time) in Method I, while solving
the linear algebraic system costs the most in Method II. In Table 2, the total CPU time and the percentageof ‘assembling the right-hand side matrix’ increase rapidly in Method I with respect to the truncation mode.
While the total CPU time of Method I is twice that of Method II when NMax = 0, it increases significantly to
15.75 times that of Method II when NMax = 10. For Method I, the percentage of ‘assembling the right-handside matrix’ quickly increases to above 90% as the truncation mode increases merely to 3. Method I is much
more computationally expensive due to the more complex right-hand side of the boundary integral equation.

®

4. CONCLUSIONS
Evaluating the DTM using hybrid source-dipole formulations is presented and good agreement is foundbetween the source-distribution method. Two alternative ways of computing the DTM are compared, showingthat while the accuracy does not differ much, the computation cost can be substantially saved by using theapproach proposed in Kashiwagi & Kohjo (1995). Moreover, it is found that when more terms of DTM areneeded, the increase of the truncation mode can lead to a significant increase of the overall CPU time, for whichissue applying the conventional approach may not be acceptable anymore.

4

Figure 2: Percentages of the time cost in the per-frequency diffraction computation using the two methods.
Assembling Assembling
Calculating Calculating
NMax , CPU Time CPU Time right-hand right-hand
DTM
DTM
MMax Method I Method II side matrix side matrix
Method I Method II
Method I Method II
0
1
2
3
4
5
10

23.12 (s)12.31 (s)70.75 (s)15.10 (s)156.48 (s) 20.12 (s)280.32 (s) 27.37 (s)442.26 (s) 36.85 (s)642.31 (s) 48.57 (s)2214.10 (s) 140.60 (s)39.53%
77.50%
87.60%
91.28%
92.99%
93.90%
95.35%

1.39%
6.78%
12.71%
17.44%
20.82%
23.17%
28.01%

15.29%
5.00%
2.26%
1.26%
0.80%
0.55%
0.16%

13.75%
11.21%
8.41%
6.18%
4.59%
3.48%
1.20%

Table 2: Increase of the time cost of per-frequency (e.g., ka = 0.921) DTM computation against the maximum
truncation mode: the 2nd and 3rd columns give the total CPU times using Method I and II; the rest 4th–7th
columns give the percentage of the time cost in the total CPU times at different stages using the two methods.

ACKNOWLEDGEMENTS
This work was supported by Grant-in-Aid for Early-Career Scientists (JSPS Grant Number JP18K13939).

REFERENCES
Babarit, A. & Delhommeau, G. 2015 Theoretical and numerical aspects of the open source BEM solver NEMOH. In
Proceeding of the 11th European Wave and Tidal Energy Conference, Nantes, France.
Dai, Y. S. & Duan, W. Y. 2008 Potential flow theory of ship motions in waves. National Defense Industry Publication.
Fenton, J. D. 1978 Wave forces on vertical bodies of revolution. Journal of Fluid Mechanics 85 (2), 241–255.
Flavià, F. F., McNatt, C., Rongère, F., Babarit, A. & Clément, Al. H. 2018 A numerical tool for the frequency
domain simulation of large arrays of identical floating bodies in waves. Ocean Engineering 148, 299–311.
Goo, J.-S. & Yoshida, K. 1990 A numerical method for huge semisubmersible responses in waves. Transactions-Societyof Naval Architects and Marine Engineers 98, 365–387.
Kagemoto, H. & Yue, D. K. P. 1986 Interactions among multiple three-dimensional bodies in water waves: an exact
algebraic method. Journal of Fluid Mechanics 166, 189–209.
Kashiwagi, M. & Kohjo, T. 1995 A calculation method for hydrodynamic interactions of multiple bodies supportinga huge floating body (in japanese) pp. 247–254.
Lee, C. H. & Newman, J. N. 2005 Computation of wave effects using the panel method. Numerical Models in Fluid
Structure Interaction 42, 211–251.
Liang, H., Ouled Housseine, C., Chen, X. B. & Shao, Y. 2020 Efficient methods free of irregular frequencies inwave and solid/porous structure interactions. Journal of Fluids and Structures 98, 103130.
Liu, Y. 2019 HAMS: a frequency-domain preprocessor for wave-structure interactions – theory, development, andapplication. Journal of Marine Science and Engineering 7 (3), 81.
McNatt, J. C., Venugopal, V. & Forehand, D. 2015 A novel method for deriving the diffraction transfer matrix
and its application to multi-body interactions in water waves. Ocean Engineering 94, 173–185.
Penalba, M., Kelly, T. & Ringwood, J. 2017 Using NEMOH for modelling wave energy converters: A comparative
study with WAMIT. In Proceeding of the 12th European Wave and Tidal Energy Conference (EWTEC2017), Cork,
Ireland .