HAMS: A Frequency-Domain Preprocessor for Wave-Structure Interactions—Theory, Development, and Application

: This paper presents the theoretical background, the numerical implementation, and the applications of a new software that has been developed in recent years for the analysis of wave-structure interactions. The software is developed in the frequency domain, as a preprocessor of computing the wave excitation force, the added mass, and the wave radiation damping, for the input to a time-domain solver via the Fourier cosine and sine transforms. In addition, it can also predict the motion responses of a marine structure with sufﬁcient accuracy, with or without the presence of a mooring system. Unlike other frequency-domain software, such as WAMIT ® and Hydrostar ® , the present software currently employs the least squares method in association with a partially extended boundary integral equation method to remove the so-called “irregular frequencies”. Calculation of the free-surface Green’s function employs a combination of fast-convergent series expansions in different parametric sub-regions. The solution of the resultant linear algebraic system employs the lower-upper (LU) decomposition method. Symmetry properties can be exploited, and the open multi-processing (OpenMP) parallelization technique can be applied to reduce the computation burden. The accuracy and the efﬁciency of the developed software are ﬁnally conﬁrmed by numerical validations on three benchmark cases of a ﬂoating ellipsoid, a truncated circular cylinder and the OC4 DeepCwind semisubmersible ﬂoating wind turbine. A free executable version of the software is available to the research communities with a hope of facilitating the advancements in the researches that are relevant to ocean engineering and marine renewable energies.


Introduction
Safety and performance are always considered as the first priorities for a marine structure when it is in operation.The sea circumstances can be varied depending on the site location, weather condition, and many other factors.It is preferable for assessments with sufficient reliability to be done before the operation or even before the construction of the marine structure is taken into action.To this end, laboratory experiments are usually performed in advance for model calibrations, wave/load measurements, and system evaluations, etc.However, due to the expensive cost, instead of carrying out the model test for all sea conditions, a reliable numerical tool to assist the assessment in association with a model test in some critical sea conditions are normally employed, which can therefore reduce the overall cost.
So far, to evaluate the overall performance of large marine structures, the boundary integral equation method (BIEM) based on the potential flow theory may still be the first choice over other computational fluid dynamics (CFD) methods based on the N-S equations [1].This is due to the fact Figure 3. Illustration of the local coordinate system of each panel, in triangular shape shape, respectively.The origin o normally locates at the centroid of each panel, and th points to the normal direction of the panel.The sequence of vertices is arranged in a the normal direction that points into the inside of the body.

Solution of the Linear Algebraic System
In order to solve the linear algebraic system Equation ( 8), a direct solv elimination will commonly be adopted.These direct solvers are generally robust computations (N denotes matrix size).For a large-scale computation of complex offshore structures, a direct inversion or inefficient iteration of such a large, den equations with   unknowns for a set of wave frequencies seems prohibitive even with modern computers.On the other hand, although some iterative me GMRES method [11], can reduce the effort to   operations; however, thes methods can easily have a problem of ill preconditions.In addition, considering th the radiation-diffraction problem needs to be solved with multiple wave headin wave frequency, the computation effort of using either a direct method or an itera not acceptable because it needs to successively solve the linear equations for every this reason, the LU decomposition is adopted in the solver, since it only needs to d (N 3 ) operations where N is the dimension of the linear systems.The iterative methods, such as the generalized minimum residual (GMRES) method [11], can reduce the effort to By taking advantage of the above formulation, the singularities encountered in the boundary inte equations can be evaluated successfully without difficulties.

Solution of the Linear Algebraic System
In order to solve the linear algebraic system Equation ( 8), a direct solver such as Ga elimination will commonly be adopted.These direct solvers are generally robust but require  computations (N denotes matrix size).For a large-scale computation of complex three-dimensio offshore structures, a direct inversion or inefficient iteration of such a large, dense system of lin equations with   unknowns for a set of wave frequencies seems prohibitively time consum even with modern computers.On the other hand, although some iterative methods, such as GMRES method [11], can reduce the effort to   operations; however, these kinds of itera methods can easily have a problem of ill preconditions.In addition, considering the usual cases w the radiation-diffraction problem needs to be solved with multiple wave headings for each sin wave frequency, the computation effort of using either a direct method or an iterative method is not acceptable because it needs to successively solve the linear equations for every wave heading.this reason, the LU decomposition is adopted in the solver, since it only needs to decompose the (N 2 ) operations but may have a problem of ill preconditions.Meanwhile, taking into consideration the fact that, normally, the radiation-diffraction problem needs to be solved with multiple wave headings for each single wave frequency, the computation effort of using such an iterative method is still high because it needs to successively solve the linear equations for every wave heading.
In addition, another computational issue that should be given attention is the suppression of the so-called "irregular frequencies".The discrete boundary integral equation is ill-conditioned and not uniquely solvable at these frequencies.The occurrence of such phenomena was firstly discovered by Ref. [12] and later studied by several scholars.There are basically two branches of numerical methods to remove these irregular frequencies: modification of the integral operator or modification of the domain of the integral operator [13,14], and extension of the boundary conditions in Dirichlet or Neumann type [15][16][17].Although the first branch of methods can theoretically find the unique solutions at all frequencies, they have a critical numerical drawback in calculating the double derivatives of free-surface Green's function.The second branch of methods requires the solution of a set of completely extended boundary integral equations for which an additional computational effort of evaluating the logarithmic singularity on the free surface is necessary.
This paper presents a new software that solves the above challenging issues.It has been developed for years and tested repeatedly via a series of research and industrial projects.The calculation of free-surface Green's function employs a combination of series expansions in different parametric subregions with the aid of the epsilon acceleration algorithm, while compromising the accuracy and efficiency of calculations.The solution of the resultant linear algebraic system adopts the lower-upper (LU) decomposition method, which needs to be evaluated only once for the left-hand side matrix over a distribution of wave headings.Removal of the irregular frequencies uses a partially extended boundary integral equation method in association with the least squares method, which reduces the computational effort and avoids evaluation of the logarithmic singularity on the free surface.In addition, symmetry properties are exploited wherever possible for bodies with 1~2 symmetry planes, and the open multi-processing (OpenMP) parallelization technique is applied to the codes in order to reduce the computation time on multi-core machines.The software is therefore named HAMS (Hydrodynamic Analysis of Marine Structures).In the subsequent computations, the computer codes are compiled by Intel ® Fortran Compiler 18.0 (Intel ® Parallel Studio XE 2018) and Microsoft ® Visual Studio 2015.
The remaining part is organized into the following sections: The background theory and the mathematical algorithms of wave-structure interactions are introduced in Section 2. The special techniques of numerical implementation are introduced as separate topics in Section 3. Validation of the software for a simple analytical geometry is given in Section 4.1.Application of the software to marine hydrodynamics of a complex floating structure is carried out in Section 4.2.The conclusions are given in Section 5.

Governing Equation and Boundary Conditions
The flow is assumed to be inviscid, free of separation or lifting effects, irrorational, and incompressible.In this potential flow framework, the flow can be described by a velocity potential φ(x, y, z) which can be further decomposed into three parts: the incident wave potential φ i (x, y, z), the diffracted wave potential φ d (x, y, z), and the radiated wave potential φ r (x, y, z): Mathematically, the above velocity potentials satisfy the Laplace equation in the entire fluid domain subjecting to boundary conditions respectively at the free surface, on the body surface, at the sea bottom, and in the far field, as shown in Figure 1, which can be expressed as where v = ω 2 /g is the wave number in deep water, g is the acceleration of gravity, V n denotes the normal velocity at a point on the immersed body boundary S B , h is the water depth in case of finite depth water, and R denotes the horizontal distance from the body.The last boundary condition in Equation ( 3), also referred to as the Sommerfeld radiation condition, shows that the velocity potential gradually decays with the horizontal distance and eventually vanishes in the far field.It should be noted that the Sommerfeld radiation condition applies only to the scattered potential (including diffracted and radiated wave potentials).
If further using subscripts from 0 to 7 to denote the incident wave potential, the six components of the radiated wave potential, and the diffracted wave potential, then the total velocity potential is expressed as where the constant ξ k (k = 1,2, . . .,6) denotes the complex amplitudes of the body oscillation motion in its six rigid-body degrees of freedom, and φ k (k = 1, 2, . . ., 6) denotes the corresponding unit-amplitude radiation potentials.In the framework of the Airy wave theory, the incident wave potential φ 0 is defined by for infinite depth water and for finite depth water, where β is the angle between the direction of propagation of the incident wave and the positive x-axis, and k is the real positive root of the water wave dispersion equation.
for finite depth water, where  is the angle between the direction of propagation of the incident wave and the positive x-axis, and k is the real positive root of the water wave dispersion equation.

Mixed Source/Dipole Formulation and Discretization of the Integral Equations
By applying Green's theorem, the radiated and diffracted wave velocity potentials on the immersed body surface SB are solved by the mixed source/dipole boundary integral equations.The integral equation satisfied by the radiation velocity potential on the body boundary takes the form of 2 () +  () (; )  d =  , ()(; )d , ( = 1,2, ⋯ ,7), where  denotes the source point,  denotes the field point, and  , denotes the kth component of the body surface boundary condition.The boundary surfaces are discretized into a set of quadrilateral or triangular plane panels to approximate the exact geometry.The vertices of each panel are numbered in the counter-clockwise direction when the panel is viewed from the fluid domain.The radiation and diffraction velocity potentials are also represented by piecewise constant functions over each panel.By applying such a "collection method", as shown in Figure 2, the boundary integral equations (Equation ( 7)) are discretized into where  is the number of panels.The integrations of sources and dipoles over each panel are represented by where  , denotes the jth panel surface.

Mixed Source/Dipole Formulation and Discretization of the Integral Equations
By applying Green's theorem, the radiated and diffracted wave velocity potentials on the immersed body surface S B are solved by the mixed source/dipole boundary integral equations.The integral equation satisfied by the radiation velocity potential on the body boundary takes the form of where ξ denotes the source point, x denotes the field point, and V n,k denotes the kth component of the body surface boundary condition.The boundary surfaces are discretized into a set of quadrilateral or triangular plane panels to approximate the exact geometry.The vertices of each panel are numbered in the counter-clockwise direction when the panel is viewed from the fluid domain.The radiation and diffraction velocity potentials are also represented by piecewise constant functions over each panel.By applying such a "collection method", as shown in Figure 2, the boundary integral equations (Equation ( 7)) are discretized into where N p is the number of panels.The integrations of sources and dipoles over each panel are represented by where S B,j denotes the jth panel surface.

Evaluation of Green's Function and Self-Influences
In the boundary integral equations (Equation ( 7)), numerical evaluation of free-surface Green's function G(x; ξ) becomes the most essential task due to its highly complex nature.It is not trivial, because of the three sticking points that need to be overcome: the singularity in the denominator of the integrand, the oscillation nature of the Bessel function, and the integration range from 0 to infinity.The Green function G(x; ξ), or source potential, is usually defined as the velocity potential at the field point (x, y, z) due to a point source of strength (−4π) located at the source point (ξ, η, ζ).The Green function is defined by for infinite water depth and for finite water depth, where the path of the contour integral passes below the pole at  =  in Equation (11) and at  =  in Equation (12).J0(x) is the Bessel function of zero order.r is the distance between the source point and the field point,  denotes the distance between the field point and the image of the source point with respect to the mean free-surface, and  denotes the distance between the field point and the image of source point with respect to the sea bottom: Currently, the algorithm in Ref. [18] has been applied to evaluate Equation (11), and a newly developed special algorithm [10] has been applied to evaluate Equation (12).Due to the strong singularity of the Rankine source r −1 , the integration of it over each panel is evaluated separately using an analytical algorithm [19].To this end, as shown in Figure 3, a local coordinate system, , needs to be established for each panel by taking the  plane at the panel surface and the  axis along its normal direction.Assuming the field point in the new local system to be (, , ), and the coordinates of vertices of each panel in a counter-clockwise direction to be ( ,  ,  ), where  = 1,2, ⋯ ,  ,  is the number of vertices on each panel, the formulation for calculating the Rankine part integration can be derived as follows:

Evaluation of Green's Function and Self-Influences
In the boundary integral equations (Equation ( 7)), numerical evaluation of free-surface Green's function G(x; ξ) becomes the most essential task due to its highly complex nature.It is not trivial, because of the three sticking points that need to be overcome: the singularity in the denominator of the integrand, the oscillation nature of the Bessel function, and the integration range from 0 to infinity.The Green function G(x; ξ), or source potential, is usually defined as the velocity potential at the field point (x, y, z) due to a point source of strength (−4π) located at the source point (ξ, η, ζ).The Green function is defined by

Evaluation of Green's Function and Self-Influences
In the boundary integral equations (Equation ( 7)), numerical evaluation of free-surface Green's function G(x; ξ) becomes the most essential task due to its highly complex nature.It is not trivial, because of the three sticking points that need to be overcome: the singularity in the denominator of the integrand, the oscillation nature of the Bessel function, and the integration range from 0 to infinity.The Green function G(x; ξ), or source potential, is usually defined as the velocity potential at the field point (x, y, z) due to a point source of strength (−4π) located at the source point (ξ, η, ζ).The Green function is defined by for infinite water depth and for finite water depth, where the path of the contour integral passes below the pole at  =  in Equation (11) and at  =  in Equation (12).J0(x) is the Bessel function of zero order.r is the distance between the source point and the field point,  denotes the distance between the field point and the image of the source point with respect to the mean free-surface, and  denotes the distance between the field point and the image of source point with respect to the sea bottom: Currently, the algorithm in Ref. [18] has been applied to evaluate Equation ( 11), and a newly developed special algorithm [10] has been applied to evaluate Equation (12).Due to the strong singularity of the Rankine source r −1 , the integration of it over each panel is evaluated separately using an analytical algorithm [19].To this end, as shown in Figure 3, a local coordinate system, , needs to be established for each panel by taking the  plane at the panel surface and the  axis along its normal direction.Assuming the field point in the new local system to be (, , ), and the coordinates of vertices of each panel in a counter-clockwise direction to be ( ,  ,  ), where  = 1,2, ⋯ ,  ,  is the number of vertices on each panel, the formulation for calculating the Rankine part integration can be derived as follows: for infinite water depth and Illustration of the collection method: Collection points are located at the centroid of each panel, and the influence coefficients Sij and Dij between each two of these points are evaluated subsequently.

Evaluation of Green's Function and Self-Influences
In the boundary integral equations (Equation ( 7)), numerical evaluation of free-surface Green's function G(x; ξ) becomes the most essential task due to its highly complex nature.It is not trivial, because of the three sticking points that need to be overcome: the singularity in the denominator of the integrand, the oscillation nature of the Bessel function, and the integration range from 0 to infinity.The Green function G(x; ξ), or source potential, is usually defined as the velocity potential at the field point (x, y, z) due to a point source of strength (−4π) located at the source point (ξ, η, ζ).The Green function is defined by for infinite water depth and for finite water depth, where the path of the contour integral passes below the pole at  =  in Equation (11) and at  =  in Equation (12).J0(x) is the Bessel function of zero order.r is the distance between the source point and the field point,  denotes the distance between the field point and the image of the source point with respect to the mean free-surface, and  denotes the distance between the field point and the image of source point with respect to the sea bottom: Currently, the algorithm in Ref. [18] has been applied to evaluate Equation (11), and a newly developed special algorithm [10] has been applied to evaluate Equation (12).
Due to the strong singularity of the Rankine source r −1 , the integration of it over each panel is evaluated separately using an analytical algorithm [19].To this end, as shown in Figure 3, a local coordinate system, , needs to be established for each panel by taking the  plane at the panel surface and the  axis along its normal direction.Assuming the field point in the new local system to be (, , ), and the coordinates of vertices of each panel in a counter-clockwise direction to be ( ,  ,  ), ,  is the number of vertices on each panel, the formulation for calculating the Rankine part integration can be derived as follows: for finite water depth, where the path of the contour integral passes below the pole at µ = v in Equation (11) and at µ = k in Equation (12).J 0 (x) is the Bessel function of zero order.r is the distance between the source point and the field point, r 1 denotes the distance between the field point and the image of the source point with respect to the mean free-surface, and r 2 denotes the distance between the field point and the image of source point with respect to the sea bottom: Currently, the algorithm in Ref. [18] has been applied to evaluate Equation (11), and a newly developed special algorithm [10] has been applied to evaluate Equation (12).Due to the strong singularity of the Rankine source r −1 , the integration of it over each panel is evaluated separately using an analytical algorithm [19].To this end, as shown in Figure 3, a local coordinate system, oξηζ, needs to be established for each panel by taking the ξoη plane at the panel surface and the ζ axis along its normal direction.Assuming the field point in the new local system to be (x, y, z), and the coordinates of vertices of each panel in a counter-clockwise direction to be (ξ is the number of vertices on each panel, the formulation for calculating the Rankine part integration can be derived as follows: where By taking advantage of the above formulation, the singularities encountered in the boundary integral equations can be evaluated successfully without difficulties. where By taking advantage of the above formulation, the singularities encountered in the boundary integral equations can be evaluated successfully without difficulties.

Solution of the Linear Algebraic System
In order to solve the linear algebraic system Equation ( 8), a direct solver such as Gauss elimination will commonly be adopted.These direct solvers are generally robust but require ( ) computations (N denotes matrix size).For a large-scale computation of complex three-dimensional offshore structures, a direct inversion or inefficient iteration of such a large, dense system of linear equations with ( ) unknowns for a set of wave frequencies seems prohibitively time consuming even with modern computers.On the other hand, although some iterative methods, such as the GMRES method [11], can reduce the effort to ( ) operations; however, these kinds of iterative methods can easily have a problem of ill preconditions.In addition, considering the usual cases when the radiation-diffraction problem needs to be solved with multiple wave headings for each single wave frequency, the computation effort of using either a direct method or an iterative method is still not acceptable because it needs to successively solve the linear equations for every wave heading.For this reason, the LU decomposition is adopted in the solver, since it only needs to decompose the left-

Solution of the Linear Algebraic System
In order to solve the linear algebraic system Equation ( 8), a direct solver such as Gauss elimination will commonly be adopted.These direct solvers are generally robust but require O N 3 computations (N denotes matrix size).For a large-scale computation of complex three-dimensional offshore structures, a direct inversion or inefficient iteration of such a large, dense system of linear equations with O N 4 unknowns for a set of wave frequencies seems prohibitively time consuming even with modern computers.On the other hand, although some iterative methods, such as the GMRES method [11], can reduce the effort to O N 2 operations; however, these kinds of iterative methods can easily have a problem of ill preconditions.In addition, considering the usual cases when the radiation-diffraction problem needs to be solved with multiple wave headings for each single wave frequency, the computation effort of using either a direct method or an iterative method is still not acceptable because it needs to successively solve the linear equations for every wave heading.For this reason, the LU decomposition is adopted in the solver, since it only needs to decompose the left-hand side matrix once for each wave frequency.This decomposition can be applied to calculate the wave forces for a distribution of wave headings via a forward and backward substitution for triangular matrices (L and U), which can be solved directly without using the Gaussian elimination process.

Removal of Irregular Frequencies
Numerical solutions of Equation ( 7) possess substantial errors in the neighborhood of the so-called "irregular frequencies".This phenomenon is caused by the water-plane section of the members of floating bodies that intersects the free water surface.The irregular frequencies actually coincide with the eigenfrequencies of the corresponding sloshing modes of the interior tank (assuming flow filling inside the tank).
In order to suppress these frequencies, a kind of partially extended boundary integral equation can be developed, which assumes that the potentials on the interior water plane are zero.This method originates from Ref. [20].Therefore, by applying Green's theorem in the interior domain of the floating body, an additional boundary integral equation is introduced in a combined application with Equation ( 7): where S WP denotes the interior water-plane area.By discretizing S WP into M elements, a set of over-determined linear algebraic equations can be obtained as follows: Equation ( 24) can be solved by the least squares method.Defining the square error function as and optimizing Equation ( 25) to be of minimum square error a new set of equations which are free of irregular frequencies can be developed instead of Equation ( 8): which is not over-determined and thus can be solved directly.The present method for removing irregular frequencies avoids evaluation of the logarithmic singularity [21] of free-surface Green's function occurring in the limiting case when that panel is on the free surface.This can be viewed as an advantage which makes the programming simpler.Besides, since the size of the left-hand side matrix in Equation ( 27) is N × N, which is less than that of (M + N) × (M + N) in Lee's method [21], the additional numerical work on resolving a larger linear algebraic system can be avoided.

Exploitation of Symmetrical Properties
The efficiency of the solution method for the boundary integral equations may dramatically increase for bodies having one or two planes of symmetry, because the computation time and storage are saved.In the analysis below, symmetry is therefore exploited wherever possible [22,23].
The discretized form of the boundary integral equation, i.e., Equation (8) or Equation (27), can be expressed in the following matrix form: It is possible to partition the matrix A, the vector φ, and B according to the number of symmetry planes.
For the body having two planes of symmetry, the following partitions can be found where A 11 is the submatrix of coefficients associated with the discretization of boundary region 1, while A 12 , A 13 , and A 14 are corresponding submatrices arising from the coupling between the boundaries 1 and 2, 3, and 4, respectively.φ 1 , φ 2 , φ 3 , and φ 4 are the vectors of components of the velocity potential corresponding to elements in boundaries 1, 2, 3 and 4 respectively.Similarly, B 1 , B 2 , B 3 , and B 4 are the right-hand side vectors associated with boundaries 1, 2, 3, and 4 respectively.Additionally, the following relations are necessary to be employed in the calculation: An orthogonal transformation matrix [R] can be introduced to reduce the size of the linear system such that where N s denotes the number of symmetry planes.This transformation enables the linear system resulting into a simplified diagonal form: Combining Equation (37) with Equations ( 33)-( 35), the solution of {φ} can be finally obtained with sufficient accuracy.The partition process is also illustrated in Figure 4.  3. The legend "Kudo" in Figure 8 has been substituted by the consistent form "Kudou (1978)" with Ref. [25] in the following new Figure 8. Correspondingly, in the text paragraph before Figure 8, "Kudo (1977)" should better be changed with "Kudou (1978)", as the follows.
A comparison of the hydrodynamic coefficients against the normalized length  is given in Figure 8.The rotation center of the torque is defined at the origin.The present numerical results show good coincidence with the analytical results obtained by Kudoh (19771978) [25],

OpenMP Parallelization on Multi-core Machines
The problems to be solved are often of a very large size such that resolving the resultant linear systems requires huge computational resources.Nowadays, with the facility of a fast multi-core computer, it is natural to maximize the advantages of the current hardware technology in our computations.A large portion of programs that people write and run daily are serial programs, especially in the marine hydrodynamic field as far as the author knows.The serial program runs on a single computer, typically on a single processor, which might be considered as a waste of computational resources.For this reason, it is better to code the program in parallel mode, which will enable the program to run simultaneously on multiple processors.Taking into consideration that in a typical case, the number of panels involved in the hydrodynamic computation is usually below 10~30 thousand, and that off-the-shelf computation machines contain multiple processors, the OpenMP technique is evidently the most suitable for the work in the present study.
The OpenMP standard for Fortran was first released in 1997.It is a standard application programming interface (API) for writing shared memory parallel applications in C, C++, and Fortran.OpenMP has the advantage of being very easy to implement on currently existing serial codes and allowing incremental parallelization [24].It also has the advantage of being widely used, highly portable, and ideally suited to multi-core architectures, which are becoming increasingly popular in up-to-date desktop computers.
As shown in Figure 5, an OpenMP program usually begins with a single process called the master thread.The master thread executes sequentially until a parallel region is encountered.At this point, the master thread "forks" into a number of parallel worker threads, each of which shares a portion of the computation task.The instructions in the parallel region are then executed by the team of worker threads.At the end of the parallel region, the threads synchronize and join to become the single master thread again [24]. of worker threads.At the end of the parallel region, the threads synchronize and join to become the single master thread again [24].Based on the Amdahl's Law, if we suppose p (0 < p < 1) part of the program can be run in parallel, on an Nt-thread computational platform, the speed-up ratio   can be expressed as This means that the maximum speed-up ratio will be Nt by letting p approach unity.However, since the sequential part of work in typical cases is not always zero, this maximum speed-up ratio will generally not be easily reached.Nevertheless, enhancing the serial code to minimize the sequential work is quite necessary to improve the computational speed.

Computation of an Analytical Geometry for Verification
A numerical benchmark case is performed on a floating ellipsoid (since the geometry of which can be precisely expressed) to demonstrate the validity of the developed solver.The geometry is firstly meshed by 20 × 20 panels (20 in longitude and 20 in latitude directions) on the immersed body surface, and additional 10 × 20 panels (10 in longitude and 20 in latitude directions) at the interior waterplane to suppress the so-called "irregular frequencies".To check its symmetrical property, the other two meshes that have one symmetry plane and two symmetry planes are also employed for comparison in altogether 300 panels and 150 panels, respectively.All the meshes are shown in Figure 6, where lengths of the major axis and the minor axis of the ellipsoid are L/A = 1.0 and B/A = 0.5, respectively, and normalized by the unit amplitude of a regular incident wave coming towards the positive x-direction.Based on the Amdahl's Law, if we suppose p (0 < p < 1) part of the program can be run in parallel, on an N t -thread computational platform, the speed-up ratio R s can be expressed as This means that the maximum speed-up ratio will be N t by letting p approach unity.However, since the sequential part of work in typical cases is not always zero, this maximum speed-up ratio will generally not be easily reached.Nevertheless, enhancing the serial code to minimize the sequential work is quite necessary to improve the computational speed.

Computation of an Analytical Geometry for Verification
A numerical benchmark case is performed on a floating ellipsoid (since the geometry of which can be precisely expressed) to demonstrate the validity of the developed solver.The geometry is firstly meshed by 20 × 20 panels (20 in longitude and 20 in latitude directions) on the immersed body surface, and additional 10 × 20 panels (10 in longitude and 20 in latitude directions) at the interior waterplane to suppress the so-called "irregular frequencies".To check its symmetrical property, the other two meshes that have one symmetry plane and two symmetry planes are also employed for comparison in altogether 300 panels and 150 panels, respectively.All the meshes are shown in Figure 6, where lengths of the major axis and the minor axis of the ellipsoid are L/A = 1.0 and B/A = 0.5, respectively, and normalized by the unit amplitude of a regular incident wave coming towards the positive x-direction.
Figure 7 shows the effect of the irregular frequencies.In the vicinity of the irregular frequencies, un-removal of them results in sharp leaps and declines within small intervals, leading to substantial numerical errors.In contrast, the "Removed" method based on the modified boundary integral equation obtains smooth results in the entire neighborhood.It is interesting to see that the irregular frequencies are different in the horizontal and the vertical exciting forces since they correspond to solutions of the interior sloshing problem in different modes.
A comparison of the hydrodynamic coefficients against the normalized length vL is given in Figure 8.The rotation center of the torque is defined at the origin.The present numerical results show good coincidence with the analytical results obtained by Kudoh (1978) [25], in both the translational mode and rotational mode.Meanwhile, the comparison also shows that although utilization of the symmetrical properties decreases the computation burden, it does not affect the accuracy of results.Figure 7 shows the effect of the irregular frequencies.In the vicinity of the irregular frequencies, un-removal of them results in sharp leaps and declines within small intervals, leading to substantial numerical errors.In contrast, the "Removed" method based on the modified boundary integral equation obtains smooth results in the entire neighborhood.It is interesting to see that the irregular frequencies are different in the horizontal and the vertical exciting forces since they correspond to solutions of the interior sloshing problem in different modes.
A comparison of the hydrodynamic coefficients against the normalized length  is given in Figure 8.The rotation center of the torque is defined at the origin.The present numerical results show good coincidence with the analytical results obtained by Kudo (1977) [25], in both the translational mode and rotational mode.Meanwhile, the comparison also shows that although utilization of the symmetrical properties decreases the computation burden, it does not affect the accuracy of results.    in both the translational mode and rotational mode.Meanwhile, the comparison also shows that although utilization of the symmetrical properties decreases the computation burden, it does not affect the accuracy of results.4. In the Concluding Remarks, a case statement is missed.It has been corrected as the follows using the Word "Track Changes" functionality:

Computation of a Truncated Circular Cylinder
Circular cylinder is another common element in the offshore structures.In this test, a truncated circular cylinder with a radius of 1.0 m and a draft of 0.5 m is chosen as a basic example, as displayed in Figure 9. Thanks to its two symmetry planes, a quarter of the cylinder is discretized into 16 × 8 × 8 constant panels (16 in circumferential, 8 in radial, and 8 in vertical directions) of a quadrilateral or triangular element shape.To remove the irregular frequencies, the water-plane area is discretized into 16 × 8 constant panels (16 in circumferential, and 8 in radial directions).Computations are performed using the HAMS solver as well as the commercial software WAMIT ® , and Hydrostar ® , for comparison.
Hydrodynamic quantities involving added mass and wave damping are compared between the computation results derived from the three software.Dimensional results are shown against the wave angular frequency ω in Tables 1 and 2. In Table 1, the computations are carried out in the deep water, while in Table 2, the water is shallow with a depth of h = 1.0 m.The tables illustrate that no matter if the water depth is infinite or finite, the results from HAMS and WAMIT ® are much closer between each other against the wave frequency, while those from Hydrostar ® are a bit shifted away.It may possibly be due to Hydrostar ® being based on the source distribution method while the other two are based on the mixed source/dipole method in calculating the wave potentials.
to surge motion and (b) cross terms of the surge added mass and radiation damping due to pitch motion.

Computation of a Truncated Circular Cylinder
Circular cylinder is another common element in the offshore structures.In this test, a truncated circular cylinder with a radius of 1.0 m and a draft of 0.5 m is chosen as a basic example, as displayed in Figure 9. Thanks to its two symmetry planes, a quarter of the cylinder is discretized into 16 × 8 × 8 constant panels (16 in circumferential, 8 in radial, and 8 in vertical directions) of a quadrilateral or triangular element shape.To remove the irregular frequencies, the water-plane area is discretized into 16 × 8 constant panels (16 in circumferential, and 8 in radial directions).Computations are performed using the HAMS solver as well as the commercial software WAMIT ® , and Hydrostar ® , for comparison.Hydrodynamic quantities involving added mass and wave damping are compared between the computation results derived from the three software.Dimensional results are shown against the wave angular frequency ω in Tables 1 and 2. In Table 1, the computations are carried out in the deep water, while in Table 2, the water is shallow with a depth of h = 1.0 m.The tables illustrate that no

Computation of a Complex Marine Structure
To illustrate the capabilities of the HAMS preprocessor in the wave-structure interaction analysis, the floating foundation of the complex OC4 DeepCwind semisubmersible floating wind turbine is presented herein as a benchmark test in practice.The platform consists of a central column, three outer offset columns, and a host of slender bracings to connect between the columns and make the floating structure sufficiently stiff.Geometrical specifications of the DeepCwind semi-submersible are found in Ref. [26].
A hydrodynamic mesh is only generated for the wetted part of the floater, and the waterplane at the cross-sections of the columns intersecting the mean sea surface.Around 3000 panels are chosen as an appropriate number of panels for the following computation, as displayed in Figure 10.Each of the three footings is discretized into 38 × 8 × 4 constant panels (38 in circumferential, 8 in radial, and 4 in vertical directions) of a quadrilateral or triangular element shape.Each of the three outer offset columns and the central column are discretized into 20 × 4 × 8 panels and 12 × 2 × 12 panels, respectively.The bracings are also meshed by sufficiently dense panels.A symmetry plane x-z is applied to remarkably reduce the computation time.Computation results using the same hydrodynamic mesh from the commercial software Hydrostar ® are also given to validate the present numerical results.
of the three footings is discretized into 38 × 8 × 4 constant panels (38 in circumferential, 8 in radial, and 4 in vertical directions) of a quadrilateral or triangular element shape.Each of the three outer offset columns and the central column are discretized into 20 × 4 × 8 panels and 12 × 2 × 12 panels, respectively.The bracings are also meshed by sufficiently dense panels.A symmetry plane x-z is applied to remarkably reduce the computation time.Computation results using the same hydrodynamic mesh from the commercial software Hydrostar ® are also given to validate the present numerical results.Hydrostatic properties including displacement, center of buoyancy, water plane, restoring, etc., of the OC4 DeepCwind floater have been calculated by HAMS and Hydrostar ® , as displayed in Table 3.These properties are evaluated based on the input hydrodynamic mesh of the immersed part of the floater.The comparison shows that the results generated by the two software are very close to each other, with quite small relative errors.A sound accuracy of the hydrostatic properties ensures a good prediction of the motion responses in subsequent calculations.Figure 11 shows distributions of the wave excitation force in the x-direction on the OC4 DeepCwind floater with respect to the wave angular frequency and wave headings.Firstly, it is seen that the distributions computed by respectively HAMS and Hydrostar ® are in good agreement with Hydrostatic properties including displacement, center of buoyancy, water plane, restoring, etc., of the OC4 DeepCwind floater have been calculated by HAMS and Hydrostar ® , as displayed in Table 3.These properties are evaluated based on the input hydrodynamic mesh of the immersed part of the floater.The comparison shows that the results generated by the two software are very close to each other, with quite small relative errors.A sound accuracy of the hydrostatic properties ensures a good prediction of the motion responses in subsequent calculations.Figure 12 shows the analysis of response amplitude operators (RAOs) on the motions of the OC4 DeepCwind floating wind turbine in no-wind and parking condition, and a comparison between the present results and those computed by the commercial software Hydrostar ® .The floating wind turbine is moored with three catenary lines at the near bottom of the outer columns spread symmetrically about the platform z-axis [26].The mooring system is considered in the calculation, and the mooring stiffness matrix is obtained via linearization of the mooring system.The motion RAOs are measured by the motion response amplitude over the incident wave amplitude.The RAOs are plotted in the logarithmic coordinate to clearly show the comparison.It is seen that the two results basically coincide well with each other, except for some tiny differences below the scale of the 10 −4 degree.The peaks at the resonance region (0.0 < ω < 0.5) predicted by the two software packages are Figure 12 shows the analysis of response amplitude operators (RAOs) on the motions of the OC4 DeepCwind floating wind turbine in no-wind and parking condition, and a comparison between the present results and those computed by the commercial software Hydrostar ® .The floating wind turbine is moored with three catenary lines at the near bottom of the outer columns spread symmetrically about the platform z-axis [26].The mooring system is considered in the calculation, and the mooring stiffness matrix is obtained via linearization of the mooring system.The motion RAOs are measured by the motion response amplitude over the incident wave amplitude.The RAOs are plotted in the logarithmic coordinate to clearly show the comparison.It is seen that the two results basically coincide well with each other, except for some tiny differences below the scale of the 10 −4 degree.The peaks at the resonance region (0.0 < ω < 0.5) predicted by the two software packages are also in quite good agreement.The motion RAOs generally decreases with the increase of the wave angular frequency to a negligible level when the wave angular frequency exceeds 2.0.The computation time and the speed-up ratio of the present case against the number of OpenMP threads for each wave frequency are shown in Figure 13.The elapsed time is an average of the computation time for 60 wave frequencies.The speedup ratio is calculated by the ratio of the computation time of a single thread over that of N t OpenMP threads.It is found that the computation time in either finite depth or infinite depth decreases with the increase of the number of threads.The computation time in finite depth is a bit less than two times of that in infinite depth.For instance, when eight threads are used, the computation time for every frequency is 9.10 s for the finite depth and 4.86 s for the infinite depth due to the fact that the calculation of the finite depth Green's function is much more time-consuming than that of infinite depth.The speedup ratio generally increases proportionally with the number of threads but seems to drop a bit when the number of threads exceeds six.Moreover, it is worth noting that the speed-up ratio is less than the number of threads because not all of the codes in the solver can be run in parallel, as indicated in Equation (38).

Concluding Remarks
A new software is presented herein for analysis of wave-structure interactions in the frequency domain.To give a clear map of how the software has been developed, the background theory as well as its numerical methodologies and techniques have been introduced in detail.The irregular frequencies have been removed based on the partially extended boundary integral equation method, and symmetry properties can be exploited to reduce the computation burden significantly.The freesurface Green's function is calculated by a combination of series expansions in both infinite depth and finite depth.OpenMP parallelization is employed to speed up computations on multi-core machines.The accuracy and the efficiency of the developed software were confirmed by numerical validations on a floating ellipsoid benchmark case and the OC4 DeepCwind semisubmersible floating wind turbine in practice.
Author Contributions: Y.L. designed the code structure, developed the software, performed the computation, and drafted and polished the paper.Acknowledgments: The author appreciates many helpful discussions with Bin Teng, Hidetsugu Iwashita, and Masashi Kashiwagi regarding the issues of calculating free-surface Green's function, removing irregular frequencies, exploiting symmetrical properties, and integrating Rankine terms, etc.The author appreciates the support from Changhong Hu and Makoto Sueyoshi on many other issues when he was undertaking his Ph.D.

Concluding Remarks
A new software is presented herein for analysis of wave-structure interactions in the frequency domain.To give a clear map of how the software has been developed, the background theory as well as its numerical methodologies and techniques have been introduced in detail.The irregular frequencies have been removed based on the partially extended boundary integral equation method, and symmetry properties can be exploited to reduce the computation burden significantly.The free-surface Green's function is calculated by a combination of series expansions in both infinite depth and finite depth.OpenMP parallelization is employed to speed up computations on multi-core machines.The accuracy and the efficiency of the developed software were confirmed by numerical validations on three benchmark cases of a floating ellipsoid, a truncated circular cylinder and the OC4 DeepCwind semisubmersible floating wind turbine.
Author Contributions: Y.L. designed the code structure, developed the software, performed the computation, and drafted and polished the paper.

Figure 3 .
Figure 3. Illustration of the local coordinate system of each panel, in triangular shape or quadrilateral shape, respectively.The origin o normally locates at the centroid of each panel, and the positive ζ-axis points to the normal direction of the panel.The sequence of vertices is arranged in accordance with the normal direction that points into the inside of the body.

Figure 1 .
Figure 1.Definition of the coordinate system in the three-dimensional space.Q denotes the source point on the immersed body boundary SB, and P denotes the field point anywhere in the fluid domain.

Figure 1 .
Figure 1.Definition of the coordinate system in the three-dimensional space.Q denotes the source point on the immersed body boundary S B , and P denotes the field point anywhere in the fluid domain.

Figure 2 .
Figure 2. Illustration of the collection method: Collection points are located at the centroid of each panel, and the influence coefficients Sij and Dij between each two of these points are evaluated subsequently.

Figure 2 .
Figure 2. Illustration of the collection method: Collection points are located at the centroid of each panel, and the influence coefficients S ij and D ij between each two of these points are evaluated subsequently.

Figure 2 .
Figure 2. Illustration of the collection method: Collection points are located at the centroid of each panel, and the influence coefficients Sij and Dij between each two of these points are evaluated subsequently.

Figure 3 .
Figure 3. Illustration of the local coordinate system of each panel, in triangular shape or quadrilateral shape, respectively.The origin o normally locates at the centroid of each panel, and the positive ζ-axis points to the normal direction of the panel.The sequence of vertices is arranged in accordance with the normal direction that points into the inside of the body.

Figure 3 .
Figure 3. Illustration of the local coordinate system of each panel, in triangular shape (a) or quadrilateral shape (b).The origin o normally locates at the centroid of each panel, and the positive ζ-axis points to the normal direction of the panel.The sequence of vertices is arranged in accordance with the normal direction that points into the inside of the body.

2 .
The space between the subfigure (a) and (b) of Figure4is too large.And the space between the subfigure (b) and the Figure4title is too small.I have adjusted the spaces in Figure4.Please substitute with the below new Figure4if possible.

Figure 4 .
Figure 4. Partitions of the computation domain (a) when the body has unique symmetry and (b) when the body has two symmetries.

Figure 4 .
Figure 4. Partitions of the computation domain (a) when the body has unique symmetry and (b) when the body has two symmetries.

Figure 5 .
Figure 5. Schematic of the two main OpenMP parallel processes in the HAMS (Hydrodynamic Analysis of Marine Structures) program.

2 :Figure 5 .
Figure 5. Schematic of the two main OpenMP parallel processes in the HAMS (Hydrodynamic Analysis of Marine Structures) program.

Figure 6 .
Figure 6.Mesh of the floating ellipsoid (a) without symmetry in 600 panels; (b) with one symmetry in 300 panels; and (c) with two symmetries in 150 panels.

Figure 6 .
Figure 6.Mesh of the floating ellipsoid (a) without symmetry in 600 panels; (b) with one symmetry in 300 panels; and (c) with two symmetries in 150 panels.

Figure 7 .
Figure 7. Modulus of the complex exciting wave forces of a floating ellipsoid, as a function of the normalized length .(a) Horizontal exciting force.(b) Vertical exciting force.

Figure 7 .
Figure 7. Modulus of the complex exciting wave forces of a floating ellipsoid, as a function of the angular wave frequency ω.(a) Horizontal exciting force.(b) Vertical exciting force.

Figure 8 .
Figure 8. Added mass and radiation damping coefficients of a floating ellipsoid as a function of the angular wave frequency : (a) diagonal terms of the surge added mass and radiation damping due to surge motion and (b) cross terms of the surge added mass and radiation damping due to pitch motion.

Figure 8 .
Figure 8. Added mass and radiation damping coefficients of a floating ellipsoid as a function of the angular wave frequency ω: (a) diagonal terms of the surge added mass and radiation damping due to surge motion and (b) cross terms of the surge added mass and radiation damping due to pitch motion.

Figure 9 .
Figure 9. Hydrodynamic mesh of the immersed part of a truncated circular cylinder for use in HAMS, WAMIT ® , and Hydrostar ®.

Figure 9 .
Figure 9. Hydrodynamic mesh of the immersed part of a truncated circular cylinder for use in HAMS, WAMIT ® , and Hydrostar ®.

Figure 10 .
Figure 10.Hydrodynamic mesh of the immersed part of the OC4 DeepCwind floater, for use in HAMS and Hydrostar ®.

Figure 10 .
Figure 10.Hydrodynamic mesh of the immersed part of the OC4 DeepCwind floater, for use in HAMS and Hydrostar ®.

Figure 11 21 Figure 11 .
Figure11shows distributions of the wave excitation force in the x-direction on the OC4 DeepCwind floater with respect to the wave angular frequency and wave headings.Firstly, it is seen that the distributions computed by respectively HAMS and Hydrostar ® are in good agreement with each other.Secondly, due to the symmetry of the OC4 DeepCwind floater, the distributions are symmetric with respect to the line of the wave heading β = 180• .It should be noted that there are several major regions where the floater is acted on significantly by the wave force.The maximum value occurs at the region of 170• < β < 190• and 1.0 < ω < 1.2.Outputting such distributions of wave excitation forces may be of great importance to the design of such floaters in practice[27].J. Mar.Sci.Eng.2019, 7, x FOR PEER REVIEW 16 of 21

Figure 11 .
Figure 11.Modulus of the wave excitation force in x-direction acting on the OC4 DeepCwind floater as a function of the wave angular frequency ω and wave headings β, computed by (a) HAMS and (b) Hydrostar ® , respectively.

Figure 12 .
Figure 12.Motion response of the OC4 DeepCwind floating wind turbine in no-wind and parking conditions as a function of the wave angular frequency ω: (a) surge response; (b) heave response; and (c) pitch response.The results denoted by symbols are computed by the commercial software Hydrostar ® .

J 21 Figure 13 .
Figure 13.Computation efficiency of the HAMS solver for each wave frequency against the number of threads in finite and infinite depth, respectively.

Funding:
Development of the software was based on the author's personal interest, but received no specific external funding.A part of the APC (395 CHF) for publishing this paper is funded by the Q-PIT (Kyushu University Platform of Inter/Transdisciplinary Energy Research) Support Program for Young Researchers (Grant Number 18112).

Figure 13 .
Figure 13.Computation efficiency of the HAMS solver for each wave frequency against the number of threads in finite and infinite depth, respectively.

Funding:
Development of the software was based on the author's personal interest, but received no specific external funding.A part of the APC (395 CHF) for publishing this paper is funded by the Q-PIT (Kyushu University Platform of Inter/Transdisciplinary Energy Research) Support Program for Young Researchers (Grant Number 18112).

Table 1 .
Hydrostatic quantities of the truncated circular cylinder calculated by HAMS, WAMIT ® , and Hydrostar ® , under the infinite water depth. A

Table 2 .
Hydrostatic quantities of the truncated circular cylinder calculated by HAMS, WAMIT ® , and Hydrostar ® , under a finite water depth h = 1.0 m.

Table 3 .
Hydrostatic properties of the OC4 DeepCwind floater calculated by HAMS and Hydrostar ®.

Table 3 .
Hydrostatic properties of the OC4 DeepCwind floater calculated by HAMS and Hydrostar ®.