An Extremely Efficient Boundary Element Method for Wave Interaction with Long Cylindrical Structures Based on Free-Surface Green ’ s Function

Yingyi Liu 1,*, Ying Gou 2, Bin Teng 2 and Shigeo Yoshida 1 1 Renewable Energy Center, Research Institute for Applied Mechanics, Kyushu University, Kasuga 8168580, Japan; yoshidas@riam.kyushu-u.ac.jp 2 State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China; gouying@dlut.edu.cn (Y.G.); bteng@dlut.edu.cn (B.T.) * Correspondence: liuyingyi@riam.kyushu-u.ac.jp; Tel.: +81-80-8565-7934


Introduction
Cylindrical structures have been widely used in the rapidly-developing coastal and offshore engineering industries in recent decades, in the form of such as floating breakwaters, oscillating water columns (OWC) for power generation, etc.These devices are used to either passively avoid the large wave kinematic energy from attacking harbors or actively convert the wave energy into other kinds of energies.Cylindrical structures are important in the industries probably due to their simplicity in geometry and the relatively lower fluid forces they may experience.Extensive efforts have been made on investigation of such kind of important structures, theoretically [1][2][3], numerically [4][5][6][7], and experimentally [8,9].In the numerical approaches, boundary element method should be one of the most popular tools for analysis [5][6][7].However, the present work is very different from the classical research in the following aspects: (1) since the governing equation used herein is the Laplace equation instead of the Helmholtz equation, Wehausen's free-surface Green function [10], constituted by several simple arithmetic functions, can be employed instead of Haskind's Green function [5,11] which may require many evaluations of the modified Bessel functions; (2) since the free-surface Green function is used as the kernel instead of the Rankine Green function, meshing of the geometry could be restricted to only the body surface, such that there is no need to deal with the open boundaries, as has been accomplished by Zheng et al. [7]; (3) application of a three-node higher-order element rather than a traditional constant or linear element guarantees the accuracy of geometrical/physical discretization; and (4) thanks to the exponential integral functions, the free-surface Green function could be written in a simpler form and then evaluated in a faster speed with a precise result.All of the above advantages facilitate numerical investigations for hydrodynamic performances of such cylindrical structures.
These horizontal cylindrical structures are so long in its axis direction that the problem to be solved could be considered in two-dimensional.In comparison to the three-dimensional wave-structure interaction model within the framework of linear potential flow theory, the present simplified two-dimensional model has many fewer unknowns on the body surface, since the surface integrations have been substituted by line integrations.In our HOBEM (higher-order boundary element method) model, for a typical frequency domain problem, about 10~50 elements (or, in other words, less than 100 nodes) are sufficient to represent an arbitrary cross-sectional shape.Therefore, generally about 10 2 ~10 4 evaluations of the Green function are needed for each incident wave period, compared to those ) application of a three-node higher-order element rather ment guarantees the accuracy of geometrical/physical ential integral functions, the free-surface Green function n evaluated in a faster speed with a precise result.All of investigations for hydrodynamic performances of such s are so long in its axis direction that the problem to be mensional.In comparison to the three-dimensional e framework of linear potential flow theory, the present y fewer unknowns on the body surface, since the surface e integrations.In our HOBEM (higher-order boundary ency domain problem, about 10~50 elements (or, in other o represent an arbitrary cross-sectional shape.Therefore, reen function are needed for each incident wave period, he three-dimensional cases (see [12]).Furthermore, by putational technologies, some special technique may be ti-processor machines.
, efficient evaluation of the free-surface Green function is rous studies have been performed in the field since 1980s.ost important contributions for this issue [12][13][14][15][16][17].They eparating the local component from the far-field one and m, or making the singular functions slow-varying by approximate the resulting functions by Chebyshev simplified in the present model since the problem to be finite water depth, in which an extremely convenient Green's kernel in the boundary integral equation.The dity and efficiency of the present method.incompressible, and the motion is assumed irrotational.onic in time with an angular frequency ω, the velocity (10 6 ) evaluations in the three-dimensional cases (see [12]).Furthermore, by taking advantage of the contemporary computational technologies, some special technique may be applied to parallelize the algorithm on multi-processor machines.
Apart from the HOBEM discretization, efficient evaluation of the free-surface Green function is another important issue in this work.Numerous studies have been performed in the field since 1980s.Noblesse and Newman have made the most important contributions for this issue [12][13][14][15][16][17].They developed several popular methods, e.g., separating the local component from the far-field one and then calculate them by tabulation algorithm, or making the singular functions slow-varying by subtracting some component and then approximate the resulting functions by Chebyshev approximation.These methods have been simplified in the present model since the problem to be considered is two-dimensional and in infinite water depth, in which an extremely convenient analytical solution can be found for the Green's kernel in the boundary integral equation.The numerical results in Section 3 show the validity and efficiency of the present method.

Governing Equation and Boundary Conditions
The problem is to consider interactions between linear water waves and a long prismatic rigid structure in arbitrary cross-sectional shape, either floating or submerged in water of infinite depth, as shown in Figure 1.The right-handed Cartesian coordinate system (x, z) is defined, with its origin located at the undisturbed free surface level and the z-axis taken vertically upward.has been accomplished by Zheng et al. [7]; (3) application of a three-node higher-order element rather than a traditional constant or linear element guarantees the accuracy of geometrical/physical discretization; and (4) thanks to the exponential integral functions, the free-surface Green function could be written in a simpler form and then evaluated in a faster speed with a precise result.All of the above advantages facilitate numerical investigations for hydrodynamic performances of such cylindrical structures.
These horizontal cylindrical structures are so long in its axis direction that the problem to be solved could be considered in two-dimensional.In comparison to the three-dimensional wave-structure interaction model within the framework of linear potential flow theory, the present simplified two-dimensional model has many fewer unknowns on the body surface, since the surface integrations have been substituted by line integrations.In our HOBEM (higher-order boundary element method) model, for a typical frequency domain problem, about 10~50 elements (or, in other words, less than 100 nodes) are sufficient to represent an arbitrary cross-sectional shape.Therefore, generally about 10 2 ~10 4 evaluations of the Green function are needed for each incident wave period, compared to those O(10 6 ) evaluations in the three-dimensional cases (see [12]).Furthermore, by taking advantage of the contemporary computational technologies, some special technique may be applied to parallelize the algorithm on multi-processor machines.
Apart from the HOBEM discretization, efficient evaluation of the free-surface Green function is another important issue in this work.Numerous studies have been performed in the field since 1980s.Noblesse and Newman have made the most important contributions for this issue [12][13][14][15][16][17].They developed several popular methods, e.g., separating the local component from the far-field one and then calculate them by tabulation algorithm, or making the singular functions slow-varying by subtracting some component and then approximate the resulting functions by Chebyshev approximation.These methods have been simplified in the present model since the problem to be considered is two-dimensional and in infinite water depth, in which an extremely convenient analytical solution can be found for the Green's kernel in the boundary integral equation.The numerical results in Section 3 show the validity and efficiency of the present method.

Governing Equation and Boundary Conditions
The problem is to consider interactions between linear water waves and a long prismatic rigid structure in arbitrary cross-sectional shape, either floating or submerged in water of infinite depth, as shown in Figure 1.The right-handed Cartesian coordinate system (x, z) is defined, with its origin located at the undisturbed free surface level and the z-axis taken vertically upward.The fluid domain is denoted by Ω, whose boundaries S consists of a free surface boundary SF, an up-side open boundary SU, a lee-side open boundary SL, and a wetted surface boundary SB on the structure.The fluid is assumed to be inviscid and incompressible, and the motion is assumed irrotational.For the linear small amplitude wave harmonic in time with an angular frequency ω, the velocity potential can be expressed by: The fluid is assumed to be inviscid and incompressible, and the motion is assumed irrotational.For the linear small amplitude wave harmonic in time with an angular frequency ω, the velocity potential can be expressed by: where φ is a time-independent complex velocity potential, which can be further decomposed into: where the three components denote incident potential, diffraction potential, and radiation potential, respectively; χ j represents displacement of the body motion in each mode (sway, heave, or roll), and φ j stands for the corresponding radiation potential to each motion mode.
The incoming wave of amplitude A and frequency ω, propagating in the positive x direction in the water of infinite water depth, can be described by the following incident velocity potential: where K is the infinite depth wave number defined by K = ω 2 /g.The four induced wave potentials φ j (j = 1~4) must satisfy the Laplace equation: and be subjected to various boundary conditions in the fluid domain, including the free surface condition on S F : the bottom boundary condition as z→∞: the boundary condition on the surface S B of the structure: and the radiation condition in the far field boundaries S U and S L : where n is the normal direction of the body geometry, with its three components where n x and n z are the x and z components of the unit inward normal, respectively, and (x c , z c ) is the rotation center.The subscripts j = 1, 2, 3 denote the direction of sway, heave, and roll for radiation, respectively, and j = 4 stands for the diffraction.

Numerical Techniques
According to Green's second theorem, by employing Wehausen's free-surface Green function as the kernel, a boundary integral equation can be obtained as: where α is the solid angle.The free-surface Green function is defined to be: where the path of the contour integral passes below the poles at µ = K; coordinate of the source point is x 0 = (ξ, ζ); r is the distance between field point and source point, and r 1 is the distance between field point and the image of source point with respect to the free surface.
On the other aspect, if Rankine Green function were employed as the kernel, the boundary integral equation would be: where the kernel would be simply expressed by: In this paper, we denote the method based on Equations ( 9) and (10) as FSG_BEM, and the method based on Equations ( 11) and ( 12) as RKG_BEM.The former is applied as the present numerical method, while the latter is used as a comparison for computational efficiency.
The three-node isoparametric element is selected to discretize both the geometry of body surface and the physical variables, the shape functions of which being expressed by: where η is the local coordinate (−1 ≤ η ≤ 1).Therefore, the velocity potential and its normal derivative on the boundary surface can be expressed straightforwardly as: φ, ∂φ ∂n Applying the above discretization and the body surface condition Equations ( 5)-( 8) leads to the following discrete form of the boundary integral equation of Equation ( 9): where N B represents the number of total elements along the body surface, and J(η) the Jacobi matrix for local-global coordinate transformation, the determinant value of which is calculated by: By employing a collocation process for Equation ( 18) that the source point is arranged to be put on each grid node on the immersed body surface mesh, a linear algebraic system could be obtained in closed form: where N is the total number of nodes.Solution of the above linear system is sensitive to the diagonal terms of the left-hand side influence matrix [A] which, therefore, needs to be evaluated precisely with caution.However, direct calculation of these diagonal terms is usually inaccurate and troublesome, due to the high singularity of the Green function in the case when the field point and the source point coincides with each other.Fortunately, this weakness can be avoided by considering a constant flux across the fluid (φ = 1); thereafter we obtain: In calculation of each influence coefficient A ij (j = i), OpenMP parallelization technique is employed to distribute the computation burden on multiple processors of a single computer.The parallelization works well since calculation of the influence coefficient on one element is independent from that on another element.After that, the Gauss elimination algorithm is used to solve the linear system, which is extremely robust regardless of arbitrary shape of the structure.
Given solution for the linear system, we can get the wave exciting force, added mass and added damping by directly integrating the corresponding hydrodynamic pressure over the immersed body surface, respectively, i.e.: and:

Direct Calculation of Free-Surface Green's Function
As pointed out in the introduction, accurate calculation of the free-surface Green function is of great importance to the final solution of the problem.At a preliminary step, we may apply the function decomposition method which was proposed by Newman [17].Equation ( 10) can be decomposed into a simplified form: where i denotes the imaginary unit.The singular part of Equation ( 23) is: where PV denotes the Cauchy principle value of the integral.The principle task here is to evaluate the real function F 1 (x, z) for all relevant values of input parameters (x, z) of possible physical interest [17].
Using the identity: Equation ( 24) can be written [18] in a more convenient form from the view point of numerical evaluation: where: In the neighborhood of µ = K, linear approximation may be applied such that: Hence, Equation ( 25) can be evaluated accurately by an adaptive Gauss-Kronrod-type quadrature algorithm (see Appendix A).
Following a similar procedure, the derivatives of the Green function with respect to x and z can be normalized as: where their singular parts are: Equations ( 30) and (31) can be formulated as the same form as Equation (25), where: Computation 2016, 4, 36 7 of 20

Fast Evaluation by the Analytical Method
Although calculation of the free-surface Green function becomes applicable following the method described in Section 2.3, a large amount of computation time would be consumed due to the direct integration by the meticulous adaptive numerical quadrature method.The reason for that is the effort of dealing with singularity in the denominator, as well as the oscillating inherence of the integrand.A possible way for its improvement is to derive an alternative analytical expression which can automatically remove the troublesome singularity, as described below.
Based on McIver [19], we can obtain the following representation for the principle value of the singular integral without too much of difficulty: where the exponential integrals are defined as: Let Z = K (z + ζ) + iKX, the real part of Equation ( 10) can then be written as: while the imaginary part is obtained by applying the residue theorem, after which we find: Based on the following identities of the exponential integral, i.e.,: It is possible to calculate derivatives of the Green function with respect to x and z, in which their real parts are corresponding to: respectively.Their imaginary parts are easily obtained by applying the residue theorem as: Through this method, we are able to calculate the free-surface Green function in a fast manner, since Equations (39), (40), and (43)-( 46) are all in analytical form, which just simply consists of the exponential functions and the trigonometric functions.
In terms of the following two coordinates: The three singular functions in Equations ( 24), ( 30) and ( 31) can be expressed as F 1 (X, Y), F 2 (X, Y), and F 3 (X, Y).Figures 2-4 show a comparison between plots of the three singular functions calculated by the direct integration method and the analytical solution method.In general, the two methods get almost same results which are hard to be distinguished from each other.It is obviously to see that F 1 (X, Y) and F 3 (X, Y) are even functions in symmetric with respect to the Y axis, while F 2 (X, Y) is an odd function which is anti-symmetric about the Y axis.Remarkable variations with a period of π in parallel to the X axis can be observed in all the plots for the region of Y ∈ [0, 3].It is also important to see that the variation becomes slow-varying with the increase of Y in the region of Y ∈ [3, 9+].

Numerical Results and Discussion
Based on the direct calculation method and the analytical solution method described above, values of the free-surface Green function and its derivatives are compared, as shown in Figures 5 and 6, against variation of the physical horizontal distance |x − ξ| between the source and the field points.Both the real part and the imaginary part are compared, showing that they coincide fairly well with

Numerical Results and Discussion
Based on the direct calculation method and the analytical solution method described above, values of the free-surface Green function and its derivatives are compared, as shown in Figures 5 and 6, against variation of the physical horizontal distance |X − ξ| between the source and the field points.Both the real part and the imaginary part are compared, showing that they coincide fairly well with each other.It should be noted that calculation of the imaginary parts is relatively straightforward since the real parts contain troublesome principal values of the singular integrals.Due to different nature of the free-surface Green function when the source point locates at the free surface or not, we need to verify the results for both floating bodies and submerged bodies.In order to compare with some existing analytical results, we select horizontal floating/submerged circular cylinders for the benchmark examples.When the cylinder is submerged, its wet body surface should be considered as a completely immersed circle; when it is floating with its centroid located on exactly the mean water surface, whereas its wet body surface should be treated as a half circle.Their analytical solutions are all obtained based on the so-called multipole expansion method, which were published in [2,20] for a submerged circle in water of infinite depth, and in [1,21] for a half circle in water of infinite depth.The RKG_BEM is also implemented for comparison, in which all of the Due to different nature of the free-surface Green function when the source point locates at the free surface or not, we need to verify the results for both floating bodies and submerged bodies.In order to compare with some existing analytical results, we select horizontal floating/submerged circular cylinders for the benchmark examples.When the cylinder is submerged, its wet body surface should be considered as a completely immersed circle; when it is floating with its centroid located on exactly the mean water surface, whereas its wet body surface should be treated as a half circle.Their analytical solutions are all obtained based on the so-called multipole expansion method, which were published in [2,20] for a submerged circle in water of infinite depth, and in [1,21] for a half circle in water of infinite depth.The RKG_BEM is also implemented for comparison, in which all of the boundaries should be taken into consideration; whereas for the FSG_BEM, only the body surface needs to be meshed.The computations are implemented on a SONY laptop (Sony Corporation, Tokyo, Japan), with an Intel(R) (Intel, Inc., Santa Clara, CA, USA) Core(TM) i7-2670QM CPU of 2.2 GHz, on 64-bit Windows operating system.The OpenMP parallelization technique has been applied in the parallel mode.
Figure 7 shows modulus of complex exciting force, added mass and added damping of a semi-immersed cylinder of radius a in comparison with those computed by the RKG_BEM and the analytical multipole expansion method [20].The corresponding meshes used by the two boundary element methods are specified in Table 1.The computation takes 39.08 s for the RKG_BEM and 2.67 s for the FSG_BEM, both in parallel mode.In Figure 7, the present method based on the analytically evaluated free-surface Green function achieves good agreement with both of the other two methods.Noted that, there are some odd points on the curve calculated by the FSG_BEM.This phenomenon should be attributed to the so-called "irregular frequencies" [10,22].The irregular frequencies are given by the eigenvalues associated with the interior Dirichlet problem, within the fluid domain on which the integral equation is applied.The discrete boundary integral equation is ill-conditioned and not uniquely solvable at these frequencies.
boundaries should be taken into consideration; whereas for the FSG_BEM, only the body surface needs to be meshed.The computations are implemented on a SONY laptop (Sony Corporation, Tokyo, Japan), with an Intel(R) (Intel, Inc., Santa Clara, CA, USA) Core(TM) i7-2670QM CPU of 2.2 GHz, on 64-bit Windows operating system.The OpenMP parallelization technique has been applied in the parallel mode.
Figure 7 shows modulus of complex exciting force, added mass and added damping of a semi-immersed cylinder of radius a in comparison with those computed by the RKG_BEM and the analytical multipole expansion method [20].The corresponding meshes used by the two boundary element methods are specified in Table 1.The computation takes 39.08 s for the RKG_BEM and 2.67 s for the FSG_BEM, both in parallel mode.In Figure 7, the present method based on the analytically evaluated free-surface Green function achieves good agreement with both of the other two methods.Noted that, there are some odd points on the curve calculated by the FSG_BEM.This phenomenon should be attributed to the so-called "irregular frequencies" [10,22].The irregular frequencies are given by the eigenvalues associated with the interior Dirichlet problem, within the fluid domain on which the integral equation is applied.The discrete boundary integral equation is ill-conditioned and not uniquely solvable at these frequencies.
Table 1.Mesh specifications for the case shown in Figure 4. LF, LU, LL and LB denote the length of the boundaries as shown in Figure 1, respectively, and NF, NU, NL, and NB denote the number of elements meshed on the boundaries, respectively.(e) (f) Figure 8 shows modulus of complex exciting forces, added mass and added damping of a submerged cylinder in comparison with those computed by the RKG_BEM and the analytical multipole expansion method [21].The computation takes 49.21 s for the RKG_BEM and 2.65 s for the FSG_BEM, both in parallel mode.In the solution of multipole expansion method, we derive a new exact formulation (see Appendix B) for the multipole expansion coefficient Amn which is troublesome for calculation due to its high singularity in the integrand.In this case, the radius of the cylinder is a, and the submergence (vertical distance from its centroid to the mean free surface) is f/a = 1.5.The meshes used by the two boundary element methods are specified in Table 2.In Figure 8, similar to the case of semi-circle, the present method based on the analytically evaluated free-surface Green function highly agrees with both of the other two methods.In addition, there is no 'irregular frequencies' phenomenon, which proves the knowledge that for submerged bodies, the solution is always unique. Figure 9 shows a comparison of computation time (unit: s) for 60 incident wave periods between the BEMs (boundary element methods) based on the direct integration and the analytical solution of the free-surface Green function, in either sequential mode or parallel mode.In Figure 9, a remarkable trend of reduction in CPU time is shown by using the parallel mode, which tends to be more apparent with increasing number of the total elements, in both Figure 9a,b.On the other hand, the analyticalbased Green function has saved a significant amount of computation time for the BEM analysis.Roughly speaking, it has improved the computation speed for around 27~36 times in the sequential mode, and 12~60 times in the parallel mode, in comparison to the direct-integration Green function method, depending on the number of total elements in the input mesh of geometry.The computation time is further compared in Figure 10, for each calling of the subroutine, involving calculation for both the value of Green's function and its derivatives.It is evidently shown that the direct integration method consumes much more time than the analytical solution method.In addition, the computation time tends to increase proportionally as the increase of the point distance |x − ξ|, which is different from that of the series expansion method as reported in the three-dimensional problem in [17,[23][24][25].Table 1.Mesh specifications for the case shown in Figure 4. L F , L U , L L and L B denote the length of the boundaries as shown in Figure 1, respectively, and N F , N U , N L , and N B denote the number of elements meshed on the boundaries, respectively.Figure 8 shows modulus of complex exciting forces, added mass and added damping of a submerged cylinder in comparison with those computed by the RKG_BEM and the analytical multipole expansion method [21].The computation takes 49.21 s for the RKG_BEM and 2.65 s for the FSG_BEM, both in parallel mode.In the solution of multipole expansion method, we derive a new exact formulation (see Appendix B) for the multipole expansion coefficient A mn which is troublesome for calculation due to its high singularity in the integrand.In this case, the radius of the cylinder is a, and the submergence (vertical distance from its centroid to the mean free surface) is f/a = 1.5.The meshes used by the two boundary element methods are specified in Table 2.In Figure 8, similar to the case of semi-circle, the present method based on the analytically evaluated free-surface Green function highly agrees with both of the other two methods.In addition, there is no "irregular frequencies" phenomenon, which proves the knowledge that for submerged bodies, the solution is always unique. Figure 9 shows a comparison of computation time (unit: s) for 60 incident wave periods between the BEMs (boundary element methods) based on the direct integration and the analytical solution of the free-surface Green function, in either sequential mode or parallel mode.In Figure 9, a remarkable trend of reduction in CPU time is shown by using the parallel mode, which tends to be more apparent with increasing number of the total elements, in both Figure 9a,b.On the other hand, the analytical-based Green function has saved a significant amount of computation time for the BEM analysis.Roughly speaking, it has improved the computation speed for around 27~36 times in the sequential mode, and 12~60 times in the parallel mode, in comparison to the direct-integration Green function method, depending on the number of total elements in the input mesh of geometry.The computation time is further compared in Figure 10, for each calling of the subroutine, involving calculation for both the value of Green's function and its derivatives.It is evidently shown that the direct integration method consumes much more time than the analytical solution method.In addition, the computation time tends to increase proportionally as the increase of the point distance |x − ξ|, which is different from that of the series expansion method as reported in the three-dimensional problem in [17,[23][24][25].Figure 11 demonstrates the convergence rate of the present FSG_BEM using analytical solution Green function.It is shown that with the increase of number of elements, the numerical solution gets close quickly to the exact solution [21], especially from 30 to 60 elements in the implemented case the difference is almost negligible.This may suggest that, for such cross-sections in simple geometries, only a few elements are sufficient to obtain a high accuracy.To further demonstrate the proposed FSG_BEM, wave interactions with a horizontal rectangular cylinder have been examined, for either floating or submerged cases.For simplicity, both of width and length of the rectangular cross-section shape are taken as 2a, respectively.f /a is the submergence of its centroid, as defined in the circular cylinder case.We examine four different submergences (particularly, f /a = 0.0 stands for the floating cylinder), and the results are shown in Figure 12.From the results, it is apparent that the floating cylinder is very different from those submerged in the water, especially in the low frequency region.In general, the peak value of the forces and its corresponding peak frequency tend to decrease as the submergence increases.S represents the cross-section area (note that S of a floating rectangular cylinder is half of a submerged one).For captions of the subplots please refer to Figure 7.

Method
Figure 11 demonstrates the convergence rate of the present FSG_BEM using analytical solution Green function.It is shown that with the increase of number of elements, the numerical solution gets close quickly to the exact solution [21], especially from 30 to 60 elements in the implemented case the difference is almost negligible.This may suggest that, for such cross-sections in simple geometries, only a few elements are sufficient to obtain a high accuracy.
To further demonstrate the proposed FSG_BEM, wave interactions with a horizontal rectangular S represents the cross-section area (note that S of a floating rectangular cylinder is half of a submerged one).For captions of the subplots please refer to Figure 7.
ions s between linear water waves and a long prismatic rigid , either floating or submerged in water of infinite depth, rtesian coordinate system (x, z) is defined, with its origin and the z-axis taken vertically upward.The fluid domain nsists of a free surface boundary SF, an up-side open and a wetted surface boundary SB on the structure.utation domain of the problem.
The fluid domain is denoted by Ω, whose boundaries S consists of a free surface boundary S F , an up-side open boundary S U , a lee-side open boundary S L , and a wetted surface boundary S B on the structure.

Figure 1 .
Figure 1.Computation domain of the problem.

Figure 1 .
Figure 1.Computation domain of the problem.

Figure 2 .Figure 2 .
Figure 2. Comparison of the singular function F1(X, Y) calculated by the two methods: (a) contour plot by direct integration; (b) oblique view by direct integration; (c) contour plot by analytical solution; and (d) oblique view by analytical solution.

Figure 2 .Figure 3 .Figure 3 .
Figure 2. Comparison of the singular function F1(X, Y) calculated by the two methods: (a) contour plot by direct integration; (b) oblique view by direct integration; (c) contour plot by analytical solution; and (d) oblique view by analytical solution.

Figure 4 .
Figure 4. Comparison of the singular function F3(X, Y) calculated by the two methods.For captions of the subplots please refer to Figure 2.

Figure 4 .
Figure 4. Comparison of the singular function F 3 (X, Y) calculated by the two methods.For captions of the subplots please refer to Figure 2.

Figure 5 .
Figure 5.Comparison of two methods for calculating Wehausen's Green function and its derivatives (K = 1.2 m −1 , ζ = −1.0m, z = −1.0m): (a) real part value of G; (b) imaginary part value of G; (c) real part value of Gx; (d) imaginary part value of Gx; (e) real part value of Gz; and (f) imaginary part value of Gz.

Figure 5 .Figure 6 .
Figure 5.Comparison of two methods for calculating Wehausen's Green function and its derivatives (K = 1.2 m −1 , ζ = −1.0m, z = −1.0m): (a) real part value of G; (b) imaginary part value of G; (c) real part value of G x ; (d) imaginary part value of G x ; (e) real part value of G z ; and (f) imaginary part value of G z .

Figure 7 .
Figure 7.Comparison of hydrodynamic characteristics of a floating cylinder, in semi-immersed circle of radius a: (a) sway exciting force; (b) heave exciting force; (c) sway added mass; (d) heave added mass; (e) sway added damping; and (f) heave added damping.

Figure 7 .
Figure 7.Comparison of hydrodynamic characteristics of a floating cylinder, in semi-immersed circle of radius a: (a) sway exciting force; (b) heave exciting force; (c) sway added mass; (d) heave added mass; (e) sway added damping; and (f) heave added damping.

Figure 8 .
Figure 8.Comparison of hydrodynamic characteristics of a submerged cylinder.For captions of the subplots please refer to Figure 7.

Figure 8 .
Figure 8.Comparison of hydrodynamic characteristics of a submerged cylinder.For captions of the subplots please refer to Figure 7.

Figure 9 .
Figure 9. CPU time of the two boundary element methods based on different calculation schemes of the free-surface Green function in sequential or parallel mode: (a) comparison for the direct integration method; and (b) comparison for the analytical solution method.

Figure 10 .Figure 11 .
Figure 10.Computation time (μs) for each implementation of the Green's kernel, by the two methods, respectively, as a function of point distance |x − ξ|.The figure is obtained based on averaged CPU time of 1 million evaluations of the two codes, respectively, for every input of the wave period ω.

Figure 9 .Figure 9 .
Figure 9. CPU time of the two boundary element methods based on different calculation schemes of the free-surface Green function in sequential or parallel mode: (a) comparison for the direct integration method; and (b) comparison for the analytical solution method.

Figure 10 .Figure 11 .Figure 10 .
Figure 10.Computation time (μs) for each implementation of the Green's kernel, by the two methods, respectively, as a function of point distance |x − ξ|.The figure is obtained based on averaged CPU time of 1 million evaluations of the two codes, respectively, for every input of the wave period ω.

Figure 9 .
Figure 9. CPU time of the two boundary element methods based on different calculation schemes of the free-surface Green function in sequential or parallel mode: (a) comparison for the direct integration method; and (b) comparison for the analytical solution method.

Figure 10 .Figure 11 .
Figure 10.Computation time (μs) for each implementation of the Green's kernel, by the two methods, respectively, as a function of point distance |x − ξ|.The figure is obtained based on averaged CPU time of 1 million evaluations of the two codes, respectively, for every input of the wave period ω.

Figure 11 .
Figure 11.Convergence test of the present method to the exact solution[21] with respect to the number of elements: (a) sway exciting force; (b) heave exciting force.

Figure 12 .
Figure 12.Variation of hydrodynamic characteristics of a floating rectangular cylinder (f/a = 0.0) and submerged rectangular cylinder for different submergences (f/a = 2.0, f/a = 4.0 and f/a = 8.0).S represents the cross-section area (note that S of a floating rectangular cylinder is half of a submerged one).For captions of the subplots please refer to Figure7.

Figure 12 .
Figure 12.Variation of hydrodynamic characteristics of a floating rectangular cylinder (f /a = 0.0) and submerged rectangular cylinder for different submergences (f /a = 2.0, f /a = 4.0 and f /a = 8.0).S represents the cross-section area (note that S of a floating rectangular cylinder is half of a submerged one).For captions of the subplots please refer to Figure7.

Table 2 .
Mesh specifications for the case shown in Figure5.Refer to Table1for the notes of the symbols.

Table 2 .
Mesh specifications for the case shown in Figure5.Refer to Table1for the notes of the symbols.