Multibody operations are routinely performed in offshore activities, for example, the floating liquefied natural gas (FLNG) and liquefied natural gas carrier (LNGC) side-by-side offloading case. To understand the phenomenon occurring inside the gap is of growing interest to the offshore industry. One important issue is the existence of the irregular frequency effect. The effect can be confused with the physical resonance. Thus, it needs to be removed. An extensive survey of the previous approaches to the irregular frequency problem has been undertaken. The matrix formulated in the boundary integral equations will become nearly singular for some frequencies. The existence of numerical round-off errors will make the matrix still solvable by a direct solver, however, it will result in unreasonably large values in some aspects of the solution, namely, the irregular frequency effect. The removal of the irregular effect is important especially for multibody hydrodynamic analysis in identifying the physical resonances caused by the configuration of floaters. This paper will mainly discuss the lid method on the internal free surface. To reach a higher accuracy, the singularity resulting from the Green function needs special care. Each term in the wave Green function will be evaluated using the corresponding analysis methods. Specifically, an analytical integral method is proposed to treat the log singularity. Finally, results with and without irregular frequency removal will be shown to demonstrate the effectiveness of our proposed method.

## Introduction

Multibody operations are routinely performed in offshore activities. This scenario is economically efficient especially in remote offshore locations. One classical example of such kind is the side-by-side case of floating liquefied natural gas (FLNG) and liquefied natural gas carrier (LNGC). The integrated concept has moved the land-based factory offshore. FLNG will produce and process the natural gas and off-load it to the close-by LNGC. This scenario will reduce the cost and the negative environmental impact to the local area. It has great potential in future offshore gas explorations.

However, when two huge floaters are in close proximity, the hydrodynamic interactions will become complex and large. It is necessary to understand the complex physical phenomenon happening between the two vessels, make accurate predictions, and take proper measurement when operating in practice. Therefore, the multibody hydrodynamics analysis will be essential in the design phase.

In such problems, the boundary element method will result in “resonance” when there is a gap if adopting the Green function approach. This “resonance” was over predicted in numerical simulations and thus considered as unphysical. Besides, the resonance resulted from the boundary element method can be divided into two categories: one is the resonance from the irregular frequencies and the other perhaps from the limitation of the method. Sometimes, the resonance caused by the irregular frequencies will confuse the true resonance. Thus, it needs to be removed.

Irregular frequencies result from the ill condition of the linear system in boundary integral problems. In other words, the matrix to be solved is almost singular at some frequencies. The solver can still provide some results because of limited accuracy of numerical techniques and computer round-off errors. The irregular frequency effect in wave–body interaction was first found by John [1]. The harmful effects in hydrodynamic analysis were identified by Frank [2]. The characteristic “sharp spikes” resulted from the irregular frequency effect were observed in added mass and damping, which affected the accuracy of the final results. The undesirable “spikes” will be confused with the physical resonance especially in multibody interaction problems. Therefore, this effect must be removed in order to expand the applications of boundary integral methods in hydrodynamic analysis.

Researchers first sought for the feasible approaches to remove this effect. Afterward, several applicable methods for more frequencies and more general shapes were proposed. The modified-integral method and the extended-boundary-condition one are the two approaches to resolve this problem.

In the study of the modified-integral method, Ursell [3] investigated this problem in an analytical way. He put wave source at the center of the circle and did not observe irregular frequencies for shorter wavelengths. This analytical study indicates that such a technique can be adopted in numerical evaluation of such problems. Schenck [4] applied the combined integral equations at selected internal points, leading to an overdetermined system. Burton and Miller [5] also adopted a modified Green function in acoustic wave scattering problem. Jones [6] added a source at the origin to remove interior eigenmode effects in acoustics. Inspired by Ursell [3], Ogilvie and Shin [7] modified the Green function integral by adding a source or dipole at the center of the internal free surface for the wave–body interaction problem. Sayer [8] also examined the suggestions from Ursell [3] on a symmetric body in finite depth. Later on, Ursell [9] demonstrated that a sequence of singularities can remove all the irregular frequencies, resulting in a method applicable for the general case. Wu and Price [10] extended this method to a twin-hull problem. Lau and Hearn [11] adopted the combined integral equation to study this problem. Lee and Sclavounos [12] further developed the modified integral method. Lee [13] pointed out that the final effect was determined by the choices of linear combination coefficients. He converted this problem into an optimization one and used the condition number at the first irregular frequency as the objective function. Besides, Kress [14] did the similar problem conversion in acoustic and electromagnetic scattering. Zhu [15] also discussed this method and proved its effectiveness.

On the other hand, Martin [16] applied null equations in setting up the equation system. Liapis [17] combined null equations to the original equation set, making it valid for all frequencies.

The extended-boundary-condition method was proposed by Paulling [18]. They enforced a fixed lid condition on the internal free surface. Ohmatsu [19] validated the method in the two-dimensional case. Kleinman [20] revisited this method by strict mathematical derivation. He proved the uniqueness in the potential formulation. Rezayat et al. [21] improved the method from Schenck [4] and applied the “lid” method in elastodynamics. Zhu [15] followed Kleinman [20], validating the effectiveness of his method. Lee et al. [22] further discussed this approach in a more general way, including the second-order effect.

In numerical evaluation, special care is needed for the integral of the Green function across free surface panels. It is because there will exist a log singularity if the panel is located on z = 0 [23,24]. Newman and Sclavounos [25] proposed one method to evaluate the log singularity, however, the important information is missing about the final expression and assumptions. Based on the idea of converting the integral across the pyramid bottom surface to that on the surrounding four surfaces, we have developed our own method for evaluation, ending up with different expressions but with accuracy up to 10−7 when comparing against Newman's results in the free surface panels only.

From the timeline of the development, the extended-boundary-condition method gradually showed its advantages. It is convenient to use, especially for users without abundant experience with such problems. In this paper, we will briefly review the modified-integral method and extended-boundary-condition method. The incomplete parts in the previous literature are clarified. In the implementation section, we will adopt the latter one. To achieve better accuracy, the integral of the wave Green function at the internal free surface will be discussed. Finally, the irregular frequency removal effect will be evaluated for the single-body and two-body cases.

## Methods

In this section, we will discuss the methodology of irregular frequency removal and the numerical evaluation of the wave Green function terms.

### Formulation of Boundary Value Problem.

The Cartesian coordinate setting is indicated in Fig. 1: V stands for the internal volume of the floater, bounded by the boundary surface Sb and Si; Sf stands for the outside free surface; V is the fluid domain, bounded by free surface Sf, body surface Sb, bottom SB, and infinite control surface Sc. n is the unit normal vector pointing outside fluid domain V, inside the floater body.

Fig. 1
Fig. 1
Herein, we will implicitly consider V and V to be the definition domain. If eliminating the evaluation on $ϕ−$, the remaining part will naturally corresponds to the discussion of only V to be the definition domain. Besides, we adopt the deep water assumption. The Green function based on Ref. [26] is defined as
$G=1r+1r′+2f*[R0(h,v)−iπJ0(h)exp(v)]$
(1)

where, $f*=ω2L/g, ρ=[(x−ξ)2+(y−η)2](1/2)$, $r=[(ρ)2+(z−ζ)2](1/2), r′=[(ρ)2+(z+ζ)2](1/2), h=f*ρ, v=f*(z+ζ)$, J0 is the Bessel function of the first kind of order 0, and R0 is the function defined in Telste and Noblesse [24].

The following case is the basic formulations for boundary integral problems when the field point is located on the body boundary. If the point of interest $x$ approaches the body boundary from V and xSb, for the external potential $ϕ$, we have
$2πϕ(x)=∬Sb∂ϕ(ξ)∂nξG(x;ξ)dSξ−∬Sbϕ(ξ)∂G(x;ξ)∂nξdSξ, x∈Sb$
(2)
For internal potential $ϕi$
$−2πϕi(x)=∬Sb∂ϕi(ξ)∂nξG(x;ξ)dSξ−∬Sbϕi(ξ)∂G(x;ξ)∂nξdSξ, x∈Sb$
(3)
The −2π is a result of consistent normal direction with that defined by the external domain. Taking the subtraction, then
$2π[ϕ(x)+ϕi(x)]=∬Sb[∂ϕ(ξ)∂nξ−∂ϕi(ξ)∂nξ]G(x;ξ)dSξ−∬Sb[ϕ(ξ)−ϕi(ξ)]∂G(x;ξ)∂nξdSξ, x∈Sb$
(4)

In this case of enforcing $ϕ(ξ)=ϕi(ξ)$ on Sb, we will get the source distribution equation; if $∂ϕ(ξ)/∂nξ=∂ϕi(ξ)/∂nξ$, then, we will get the doublet distribution equation.

### Mathematical Background.

The final output is based on the solution of the potential or source strength. To get the potential value from body boundary conditions, we need to solve the following equation:
$2πϕ(x)+∬Sbϕ(ξ)∂G(x,ξ)∂nξdSξ=∬Sb∂ϕ(ξ)∂nξG(x;ξ)dSξ, x∈Sb$
If converting it into a matrix form, we get
$[A][ϕ]=[b]$

From the theory of ordinary differential equations, we assume that the matrix equation has its homogeneous form. Thus, there will be a homogeneous solution and a particular one: $ϕ=ϕh+ϕp$. Based on the property of matrices, the following conclusion can be drawn:

• (a)

If $det[A(ω)]=0, ϕh=0, ϕp$ has only one solution.

• (b)

If $det[A(ω)]=0, ϕh$ has infinite solutions, $ϕp$ has infinite solutions or none, depending on the value of [b].

Therefore, the second case needs to be avoided, and it is necessary to modify the structure of the matrix [A] to become [A*]. The objective is to make the $rank(A*)=rank(ϕ)=rank(A*,b)=n$ to ensure the uniqueness of the solution. In other words, the homogeneous solution must be zero and only zero.

There are two general approaches to achieve this: One is to construct [A*](m × n), where m > n, rank(A*) = n, n is the number of unknowns. This method will result in an overdetermined matrix. In some cases, [A*] has the same rank with [A], but a different structure.

The other is to construct a square matrix [A**](m × m), where rank(A**) = m, the number of unknowns equal to m. This method needs more unknowns and needs to extend the boundary conditions, ensuring it is still a square matrix.

### Methods to Remove Irregular Frequency

#### Method I: Modified Green Function Method.

The essence of this method is to find more equations which the unknowns could satisfy, leading to an overdetermined system but guaranteeing a unique solution. More exactly, it is to choose some points inside the domain V that the unknowns could still satisfy some conditions. For the external potential $ϕ$, we have
$2πϕ(x)=∬Sb∂ϕ(ξ)∂nξG(x;ξ)dSξ−∬Sbϕ(ξ)∂G(x;ξ)∂nξdSξ, x∈Sb$
(5)
$0=∬Sb∂ϕ(ξ)∂nξG(x;ξ)dSξ−∬Sbϕ(ξ)∂G(x;ξ)∂nξdSξ, x∈V−$
(6)
If the uniqueness of the solution can be proven, then, the irregular frequencies effects can be removed. The proof is similar to the procedure in Schenck [4]. Under the assumption that $ϕx$ satisfies Eqs. (5) and (6) and $ϕy=ϕx+Aϕz$. The key equation in the proof is
$A∬Sbϕz∂G∂nξ=0 x∈V−$
(7)

The conclusion is that $∬Sbϕz(∂G/∂nξ)$ is not necessarily zero for every $x∈V−$. Then A = 0, ensuring that the solution is unique. So if the chosen points are not the node points of the homogeneous solution for the Dirichlet internal potential problem, the solution will be unique. To ensure it, a sufficiently large number of interior points might be needed.

This method will result in an overdetermined problem, which can be solved using a least square approach. However, special treatment will be needed in selecting the interior points and setting up parameters to remove the irregular frequency effect for an arbitrary shape. To make it convenient to users, the fixed lid method is discussed by Kleinman [20], Zhu [15], and Lee et al. [22].

#### Method II: Extended Boundary Condition Method.

The motivation of this method is to convert the overdetermined linear system into a square matrix. Then, a direct matrix solver can be utilized. Meanwhile, the irregular frequency corresponds to the sloshing mode of interior space. Therefore, it would be natural to place a lid on the interior free surface to suppress it.

The final equations to be solved are listed here. For a detailed derivation, please refer to Liu and Falzarano [27]
$4π∂ϕ(x)∂nx=2πσ(x)+PV∬Sb σ(ξ)∂G(x;ξ)∂nxdSξ+∬Si σ′(ξ)∂G(x;ξ)∂nxdSξ, x∈Sb$
(8)
$4π∂ϕ−(x)∂nx=4πσ′(x)+∬Sb σ(ξ)∂G(x;ξ)∂nxdSξ+PV∬Si σ′(ξ)∂G(x;ξ)∂nxdSξ, x∈Si$
(9)

In both the potential formula and the source formula, the derivative of the Green function needs to be evaluated at the free surface. Based on Newman [23] and Telste and Noblesse [24], the expressions for the Green function both contain a log term, which will result in a singular value when the field point is infinitesimally close to the source point and both are near the free surface.

Newman:
$G(x;ξ)=1r+1r′−2kek(z+ζ)[log(r′+|z+ζ|)+(γ−log 2)+r′]+O(r′2 log r′)$
(10)

where $x$ is the position of field point, and $ξ$ is the source point. $r2=(x−ξ)2+(y−η)2+(z−ζ)2, r′2=(x−ξ)2+(y−η)2+(z+ζ)2, γ=0.577…$ is the Euler constant.

Noblesse:
$4πG=−(1r+1r′)−2f*[R0(h,v)+iπJ0(h)exp(v)]$
(11)

where $f*=ω2L/g, ρ=[(x−ξ)2+(y−η)2]1/2, r=[(ρ)2+(z−ζ)2]1/2, r′=[(ρ)2+(z+ζ)2]1/2, h=f*ρ, v=f*(z+ζ)$, J0 is Bessel function of the first kind at order 0, and R0 is the function defined in Telste and Noblesse [24]. In the limiting case when $d→0, R0=−ln(d−v)+ln(2)−γ$ with error $O(d ln(d))$.

When the field point and the image source point become infinitesimally close, the log term inside both Green functions will become singular. However, the integral value of the log term across the panel is still finite, depending on the panel size. Based on the numerical evaluation in Sec. 4, it is necessary to consider the log singular effect.

The Green function will have different behaviors when the field point is infinitesimally close to the source point and both are located at internal free surface. The relations are as below
$∂G∂z=kG, when z=0, ζ=0, x→ξ or x→ξ$
(12)
$∂G∂z=kG, when z→0−, ζ=0, x→ξ$
(13)
$∂G∂z=−2zr3+kG, when z→0−, ζ=0, x→ξ$
(14)

In the constant panel method, each panel has a uniform source distribution. All the internal lid panels are exactly on the z = 0 surface, while the field points are inside the floater body. Thus, Eqs. (13) and (14) are better descriptions of the practical model.

In calculating the panel effect on itself, the integral of –2z/r3 will result in 4π. Also please note that the above three equations are derived based on Noblesse's Green function. Noblesse uses a source point as the singular part, while Newman uses a sink point.

## Evaluation of Log Singularity

When the source point is close to field point, we will have a natural log singularity in the Green function. This happens when studying the internal free surface panels. The log term can be written as
$f=ln(d−v)$
(15)
where $d=(x−ξ)2+(y−η)2+(z+ζ)2, v=z+ζ$. To solve this log singular integral, the approach of Newman and Sclavounos [25] is adopted, which is to form a pyramid first and convert the bottom surface integral into that along the four triangular facets. The derivation herein is based on the panel coordinates, and the final results are different from Newman's. The function f is harmonic in the pyramid region. Green's second identity can be used, and then, we can get
$∫S(ζ∇f−f∇ζ)·n dS+∫∑i=14Ti(ζ∇f−f∇ζ)·n dS=0L+∑i=14Li=0$
On this bottom panel, ζ = 0, $(∂/∂n)=(∂/∂z)$, then
$L=−∫Sf dξdηLi=∫Ti(ζ∇f−f∇ζ)·n dS$
If the field point is at a finite distance from the panel, ζ in Li is not zero. However, the log term will become a singularity only if the field point is infinitesimally close to the panel, which makes it necessary to evaluate the log integral in an alternative way. Thus, in this case, we need to assume ζ → 0 even on the triangular facets. Thus, we will have
$Li=∫Tif dS$
(16)

Please note that on the triangular surface, the approximate normal vector n is pointing along negative z-axis, i.e., $(∂/∂n)=−(∂/∂z)$. From the above, we need to integrate f along the triangular surface. It will be convenient to complete it in the polar coordinate on the triangular surface. Thus, the coordinate is described in Fig. 2.

Fig. 2
Fig. 2

We are more interested in the internal free surface panels. To have a better view of these angles, the panel is flipped over. That is the reason why z-axis is pointing downward.

Assume O-uvw is the Cartesian coordinate on the triangular surface. The w-axis is perpendicular to the surface BCE in Fig. 2. We determine u such that u will point to the edge $BC¯$, and the surface u–O–w is perpendicular to the edge $BC¯$. This has ensured that the u-axis is perpendicular to edge $BC¯$, and v-axis is, therefore, determined. The angle between the w-axis and negative z-axis is defined as $φ$. For an arbitrary vector a on the surface BCE, α is the angle between the vector a, and the u-axis, and θ is the angle between the vector a and the positive z-axis. Naturally, the u-axis becomes the reference axis when converting to the polar coordinate on the surface BCE. In the O-uvw coordinate system, we define B (u1, v1), C(u2, v2). The line equation can be A0u + B0v = 1. If in polar coordinate and $u=ρ cos α, v=ρ sin α$, it is easy to get
$ρ(α)=1A0 cos α+B0 sin α=1C0 cos β$
where $C0=A02+B02, β=α−δ, δ=arctan(B0/A0)$. After solving for the following two equations:
$A0u1+B0v1=1A0u2+B0v2=1$
Solve for A0, B0 as:
$A0=v2−v1u1v2−u2v1B0=u1−u2u1v2−u2v1$
(17)

Based on the geometric relation, we have $cos θ=cos((π/2)−φ)cos α=sin φ cos α$.

Therefore, on the triangular surface, this equation can be transformed to
$f=ln[d(1−vd)]=ln[d(1+cos α sin φ)]$
(18)
In the discussion later, we distinguish the source point on the quadrilateral panel and those on the triangular facet. We will use r instead of d. Then, the integral equation will become
$∫α1α2dα∫0r(α)ln[r(1+cos α sin φ)]$
(19)
If enforcing the assumption that v → 0, then $φ→0,sin φ≈φ, cos φ≈1$. The final result will be
$∫Tif dS=∫α1α2dα∫0r(α)ln[r(1+cos α sin φ)]rdr={− tan β2C2[lnC+ln(cos β)+32]−φ2C2[cos δln[tan(π4+β2)]− sin δ cos β]+β2C2}|β1β2$
Finally, we can get the integral of log singularity on the panel by this equation
$L=∫Sf dξdη=∑i=14∫Tif dS$

## Results and Discussion

In this section, we will first verify the accuracy of our method to evaluate the log singularity, then, justify the method of subtracting the log singularity from the wavy Green function, and finally, we evaluate the irregular removal effects.

### Log Singularity.

As mentioned in Sec. 1, Newman and Sclavounos [25] did not provide complete details about the assumptions, and some other information for the final expression was also not included. To be conservative, we choose some special cases in which Newman's final expression might be valid and applied a systematic trial-and-error approach to tune the parameters until we get similar results with Maple. Afterward, based on the idea, we have developed our own method for evaluation, ending up with somewhat different expressions. When comparing against Newman's results, not surprisingly, they are very close. The comparison validates our assumptions about Newman's expression and also proves our approach is accurate.

Table 1 contains the comparison results for the method in this paper and modified Newman's. As can be seen from the table, the difference is 10−7. We can draw the conclusion that our method will be accurate enough in evaluating the log singularity. Please note that Table 2 contains the node position vector in panel coordinate for different cases. The log analytical evaluation should be used for panels on and near the free surface. Based on our numerical testing, when $|v|<0.05$, separate evaluation of the log singularity will be needed. The alternative approaches are discussed in Liu and Falzarano [28].

Table 1

Comparison of log singularity in different methods

CaseLiu et al.NewmanAbs error
1−1.472113−1.4721133.48 × 10−7
290.62547590.625475−4.12 × 10−7
3−1.268191−1.268191−2.48 × 10−7
480.64736880.647368−5.39 × 10−7
5−1.188854−1.188854−2.33 × 10−7
CaseLiu et al.NewmanAbs error
1−1.472113−1.4721133.48 × 10−7
290.62547590.625475−4.12 × 10−7
3−1.268191−1.268191−2.48 × 10−7
480.64736880.647368−5.39 × 10−7
5−1.188854−1.188854−2.33 × 10−7
Table 2

Explanation on cases setting

Casexy
2 × 2 square[−1, 1, 1, −1][−1, −1, 1, 1]
Arbitrary quadrilateral[−1, 10, 1, −1][−1, −1, 10, 1]
Trapezoid[−1, 1, 0.1, −0.1][−1, −1, 1, 1]
Long rectangular[−1, 1, 1, −1][−1, −1, 20, 20]
Triangle[−1, 1, 0, 0][−1, −1, 1, 1]
Casexy
2 × 2 square[−1, 1, 1, −1][−1, −1, 1, 1]
Arbitrary quadrilateral[−1, 10, 1, −1][−1, −1, 10, 1]
Trapezoid[−1, 1, 0.1, −0.1][−1, −1, 1, 1]
Long rectangular[−1, 1, 1, −1][−1, −1, 20, 20]
Triangle[−1, 1, 0, 0][−1, −1, 1, 1]

### Integral of R0.

Since the log singularity was evaluated analytically, one may be curious about the function shape of R0 after subtracting the log integral from it. This section will illustrate the shape and justify a proper numerical method to evaluate R0. The numerical evaluation method for R0 is from Telste and Noblesse [24].

Figures 36 illustrate the function value of $R0−ln(d−v)$ across a unit square panel.

Fig. 3
Fig. 3
Fig. 4
Fig. 4
Fig. 5
Fig. 5
Fig. 6
Fig. 6

As shown in Figs. 36, if the panel size is small, a four-node Gaussian quadrature method can still give accurate results. However, if the panel is a larger size, the wavy behavior will nullify the four-node Gaussian quadrature. To balance accuracy and efficiency, it is preferable to construct smaller panels for the internal lid surface. Moreover, please note that in implementing the Gaussian quadrature method, the input is dimensionless. The variable related to length is multiplied by f = ω2L/g, where ω is the wave frequency, L is the wavelength, and g is the gravity constant. For a floater with a large L, when wave frequency ω is relatively higher, the node position could be amplified. Therefore, the four-node quadrature method may not produce accurate enough results. Nevertheless, the shorter wave may not lie in the range of interest for such floaters in sea-keeping analysis.

Fig. 7
Fig. 7
Fig. 8
Fig. 8
Fig. 9
Fig. 9

### Irregular Frequency Removal Effect.

The accurate numerical evaluation of the log singularity is important to eliminate the irregular frequency effect. This section will demonstrate the irregular frequency removal effects for single-body and multibody case. The results are generated by MDL-MultiDYN, an in-house program developed by Marine Dynamic Laboratory, Texas A&M University. This program is a redesigned program based on MDL HydroD by Guha [29]. It is able to conduct multibody analysis with improved log singularity evaluation and an irregular frequency removal module [30]. It also has the capability to analyze the forward speed effect of multiple floaters discussed in Liu and Falzarano [31,32] and the drift forces or added resistances addressed in Liu and Falzarano [33].

From John [1], the irregular frequency effect is more likely to happen in shorter waves. Thus, we adopted a miniboxbarge as the test case. The miniboxbarge has 500 panels. If we add a lid at the internal free surface, the panel number will be 700. All the cases are evaluated under head sea condition. For two-body cases, the separation distance is 10 m. Besides the miniboxbarge with the most significant irregular frequency removal effects, we have also validated our approach using boxbarges side-by-side (10 m apart), a cylindrical dock, the U.S. navy ships BOBO and Bob Hope. The particulars are listed in Table 3 (unit: meter). Please note that all the cases are conducted in the head sea condition. We use 1 6 to denote the six degrees-of-freedom (DOFs: surge, sway, heave, roll, pitch, and yaw) for the body 1, 7 12 for the body 2 in the multiple-body cases.

Table 3

Floater particulars

LBTkxxkyykzz
Miniboxbarge20105111
Boxbarge80201020520
Cylinder Dock404020111
BOBO187.32532.1566.68512.8624146.831346.8313
Bob Hope269.4532.2588.79511.4461066.5838066.58380
LBTkxxkyykzz
Miniboxbarge20105111
Boxbarge80201020520
Cylinder Dock404020111
BOBO187.32532.1566.68512.8624146.831346.8313
Bob Hope269.4532.2588.79511.4461066.5838066.58380

The irregular frequency effect is more significant in added mass and damping, phases of forces (Froude-Krylov and Diffraction (FKD)) and motion, less apparent in amplitude of forces and motion response amplitude operator. Herein, we benchmark our results against WAMIT version 6. To follow the notation in WAMIT, we use “IRR” to denote whether the irregular frequency module is activated or not. “IRR1” means that it is in effect, while “IRR0” means it is not.

Figures 714 are the results for the single miniboxbarge (Fig. 15). We may find the significant irregular frequency effect that appear in the higher frequency range. After implementing the lid method, the irregular frequency effects are removed. However, we may observe a tiny discrepancy appearing in the results for added mass A15, B11, B55 near the irregular frequency. The relative error is very small. The discrepancy may be caused by the numerical error in the discretization of the analytical equation. For example, it may be affected by the method we define the collocation point at which the potential is calculated or the way we calculate the part of wavy Green function. The methods of our in-house program were discussed by Guha [29]. In the calculations of drift forces, we observe an interesting phenomenon. In the drift force along the surge direction, we have obtained totally different results compared against WAMIT. The results from MDL multi-DYN are in dashed line, while those from WAMIT are in dotted line. We find the dotted line has obvious discrepancy against the results when the irregular frequency module is not invoked. The dashed line is closer to the results when not removing the irregular frequencies. All of the spikes almost disappear on the dashed line. We believe that our results are more reasonable. However, it is still not clear about the causes of the discrepancy.

Fig. 10
Fig. 10
Fig. 11
Fig. 11
Fig. 12
Fig. 12
Fig. 13
Fig. 13
Fig. 14
Fig. 14
Fig. 15
Fig. 15
Fig. 16
Fig. 16

Figures 1621 are the results for the two miniboxbarges (Fig. 22), side-by-side. From the results, we may find that the irregular frequencies are more likely to happen in the higher frequency range and are successfully removed. From the results for the damping term B57, we may find that one spike is due to the irregular frequency, another spike nearby is due to the physical resonance. This has validated our assumption that the irregular frequency confuses with the true resonance. Thus, in the cases of the box-shape floaters, we must invoke the irregular frequency module.

Fig. 17
Fig. 17
Fig. 18
Fig. 18
Fig. 19
Fig. 19
Fig. 20
Fig. 20
Fig. 21
Fig. 21
Fig. 22
Fig. 22
Fig. 23
Fig. 23

Figures 2326 are the results for case of the larger boxbarge (Fig. 27). For the larger-size boxbarge, the irregular frequency effect happens in the relatively higher frequency range as well. It is less significant compared with cases of the miniboxbarge. By comparison against the miniboxbarge, it also shows that the significance of the irregular frequency effect may be influenced by the sizes of the floater. It is recommended to run the cases without invoking the module and with the module activated.

Fig. 24
Fig. 24
Fig. 25
Fig. 25
Fig. 26
Fig. 26
Fig. 27
Fig. 27
Fig. 28
Fig. 28

Figures 2831 show the irregular frequency effect of the cylindrical dock (Fig. 32). It is also successfully removed. Similar to the large box barge, the irregular frequency effect is not significant. The nondimensional irregular frequency in the figures is around 2. It shows that the irregular frequency is an intrinsic property of the shape of the floater. It should be removed to obtain more accurate results.

Fig. 29
Fig. 29
Fig. 30
Fig. 30
Fig. 31
Fig. 31
Fig. 32
Fig. 32
Fig. 33
Fig. 33

When coming to the verification of the ships, we also output the results from WAMIT (Lee [34]) to demonstrate the effectiveness of our irregular frequency removal module. Figures 3336 are the results for the ship Bob Hope (Fig. 37). The irregular frequency effect is not significant. However, it is still removed. In the results of damping term B15, we find that the discrepancy between the results with irregular frequency inactive and active is more obvious compared against the cases of the miniboxbarge. The reason is not clear yet. We will investigate it more in the future.

Fig. 34
Fig. 34
Fig. 35
Fig. 35
Fig. 36
Fig. 36
Fig. 37
Fig. 37
Fig. 38
Fig. 38

Figures 3843 show the results for ship BOBO (Fig. 44). In the single body case, the irregular frequency effect is very significant in some components in the added mass matrix, for example, A13, A33, etc. For this case, we also find that the irregular frequency effects exist in the drift forces of the floater. The results in the heave drift force and pitch drift force are also presented here to demonstrate the effectiveness of our method.

Fig. 39
Fig. 39
Fig. 40
Fig. 40
Fig. 41
Fig. 41
Fig. 42
Fig. 42
Fig. 43
Fig. 43
Fig. 44
Fig. 44
Fig. 45
Fig. 45

Figures 4550 are the results when ship BOBO is next to the ship Bob Hope (Fig. 51). The separation distance is 3 m. Similar to the miniboxbarge case, the results prove that the irregular frequency effect will confuse the resonance in the relatively higher frequency range. For this case, the irregular frequency effect may exist in a small range of the frequency. For example, in the results of added mass terms A11 and A55, we observe that the discrepancy exists when nondimensional frequency varies from 6 to 7. In the results of the added mass term A33 and the damping term B11, the irregular frequency effect results in a single spike. For the interactions of multiple bodies, it will be necessary to remove the irregular frequency effect.

Fig. 46
Fig. 46
Fig. 47
Fig. 47
Fig. 48
Fig. 48
Fig. 49
Fig. 49
Fig. 50
Fig. 50
Fig. 51
Fig. 51

## Conclusion

This paper reviews the reason of the irregular frequency effects, evaluates the modified integral Green function approach and extended boundary condition approach for removing the irregular frequency effects, thoroughly discusses the analytical method to evaluate the log singularity term inside wavy Green function, and finally, validates the irregular frequency removal effect of in-house program “MDL MultiDYN.”

The comparison of the log singularity ensures the high accuracy of our method. Finally, the irregular frequency removal module is effective based on the figures in the results section. Therefore, the incorporation of this module will enhance the capability of MDL MultiDYN on multibody problems. For the interactions of multiple bodies, the irregular frequency module must be invoked to ensure a higher accuracy in capturing the true resonance. Moreover, the irregular frequency effects in the drift forces still require more investigation in the future.

## Acknowledgment

The authors would like to thank Dr. Paul Hess for supporting this work. The authors also acknowledge the partial funding provided by Society of Naval Architects and Marine Engineers (SNAME) to support this work. Finally, the authors would thank Dr. Francis Noblesse for his suggestions on Green function evaluations.

Funding Data

• Office of Naval Research (N00014-16-1-2281).

## References

1.
John
,
F.
,
1950
, “
On the Motion of Floating Bodies II
,”
Commun. Pure Appl. Math.
,
3
(
1
), pp.
45
101
.
2.
Frank
,
W.
,
1967
, “
Oscillation of Cylinders in or Below the Free Surface of Deep Fluids
,” Naval Ship Research and Development Center, Bethesda, MD, Technical Report No.
2375
.http://www.dtic.mil/dtic/tr/fulltext/u2/668338.pdf
3.
Ursell
,
F.
,
1953
, “
Short Surface Waves Due to an Oscillating Immersed Body
,”
Proc. R. Soc. A: Math., Phys. Eng. Sci.
,
220
(
1140
), pp.
90
103
.
4.
Schenck
,
H. A.
,
1968
, “
Improved Integral Formulation for Acoustic Radiation Problems
,”
J. Acoust. Soc. Am.
,
44
(
1
), pp.
41
58
.
5.
Burton
,
A. J.
, and
Miller
,
G. F.
,
1971
, “
The Application of Integral Equation Methods to the Numerical Solution of Some Exterior Boundary-Value Problems
,”
Proc. R. Soc. A: Math., Phys. Eng. Sci.
,
323
(
1553
), pp.
201
210
.
6.
Jones
,
D.
,
1974
, “
Integral Equations for the Exterior Acoustic Problem
,”
Q. J. Mech. Appl. Math.
,
27
(
1
), pp.
129
142
.
7.
Ogilvie
,
T. F.
, and
Shin
,
Y. S.
,
1978
, “
Integral Equation Solutions for Time Dependent Free Surface Problems
,”
J. Soc. Nav. Arch. Jpn.
,
1978
(
143
), pp.
41
51
.
8.
Sayer
,
P.
,
1980
, “
An Integral-Equation Method for Determining the Fluid Motion Due to a Cylinder Heaving on Water of Finite Depth
,”
Proc. R. Soc. London, Ser. A
,
372
(
1748
), pp.
93
110
.
9.
Ursell
,
F.
,
1981
, “
Irregular Frequencies and the Motion of Floating Bodies
,”
J. Fluid Mech.
,
105
, pp.
143
156
.
10.
Wu
,
X.-J.
, and
Price
,
W.
,
1987
, “
A Multiple Green's Function Expression for the Hydrodynamic Analysis of Multi-Hull Structures
,”
Appl. Ocean Res.
,
9
(
2
), pp.
58
66
.
11.
Lau
,
S. M.
, and
Hearn
,
G. E.
,
1989
, “
Suppression of Irregular Frequency Effects in Fluid-Structure Interaction Problems Using a Combined Boundary Integral Equation Method
,”
Int. J. Numer. Methods Fluids
,
9
(
7
), pp.
763
782
.
12.
Lee
,
C.-H.
, and
Sclavounos
,
P. D.
,
1989
, “
Removing the Irregular Frequencies From Integral Equations in Wave-Body Interactions
,”
J. Fluid Mech.
,
207
, pp.
393
418
.
13.
Lee
,
C.-H.
,
1988
, “
Numerical Methods for Boundary Integral Equations in Wave Body Interactions
,” Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA.
14.
Kress
,
R.
,
1985
, “
Minimizing the Condition Number of Boundary Integral Operators in Acoustic and Electromagnetic Scattering
,”
Q. J. Mech. Appl. Math.
,
38
(
2
), pp.
323
341
.
15.
Zhu
,
X.
,
1994
, “
Irregular Frequency Removal From the Boundary Integral Equation for the Wave-Body Problem
,”
Ph.D. thesis
, Massachusetts Institute of Technology, Cambridge, MA.https://dspace.mit.edu/bitstream/handle/1721.1/11691/32279180-MIT.pdf?sequence=2
16.
Martin
,
P.
,
1981
, “
On the Null-Field Equations for Water-Wave Radiation Problems
,”
J. Fluid Mech
,
113
, pp.
315
332
.
17.
Liapis
,
S.
,
1993
, “
A Method for Suppressing the Irregular Frequencies From Integral Equations in Water Wave-Structure Interaction Problems
,”
Comput. Mech.
,
12
(
1
), pp.
59
68
.
18.
Paulling
,
J.
,
1970
, “
Stability and Ship Motion in a Seaway
,” Defense Technical Information Center, Fort Belvoir, VA, Report No.
19.
Ohmatsu
,
S.
,
1975
, “
On the Irregular Frequencies in the Theory of Oscillating Bodies in a Free Surface
,” Ship Research Institute, Tokyo, Japan, Paper No.
48
.http://mararchief.tudelft.nl/catalogue/entries/12180/
20.
Kleinman
,
R. E.
,
1982
, “
On the Mathematical Theory of the Motion of Floating Bodies: An Update
,” David W. Taylor Naval Ship Research and Development Center, Bethesda, MD, Report No.
DTNSRDC-82/074
.http://www.dtic.mil/dtic/tr/fulltext/u2/a121433.pdf
21.
Rezayat
,
M.
,
Shippy
,
D.
, and
Rizzo
,
F.
,
1986
, “
On Time-Harmonic Elastic-Wave Analysis by the Boundary Element Method for Moderate to High Frequencies
,”
Comput. Methods Appl. Mech. Eng.
,
55
(
3
), pp.
349
367
.
22.
Lee
,
C.
,
Newman
,
J. N.
, and
Zhu
,
X.
,
1996
, “
An Extended Boundary Integral Equation Method for the Removal of Irregular Frequency Effects
,”
Int. J. Numer. Methods Fluids
,
23
(
7
), pp.
637
660
.
23.
Newman
,
J. N.
,
1985
, “
Algorithms for the Free-Surface Green Function
,”
J. Eng. Math.
,
19
(
1
), pp.
57
67
.
24.
Telste
,
J.
, and
Noblesse
,
F.
,
1986
, “
Numerical Evaluation of the Green Function of Water-Wave Radiation and Diffraction
,”
J. Ship Res.
,
30
(
2
), pp.
69
84
.
25.
Newman
,
J.
, and
Sclavounos
,
P.
,
1988
, “
The Computation of Wave Loads on Large Offshore Structures
,”
International Conference on Behaviour of Offshore Structures
(
BOSS
),
Trondheim, Norway
,
June
, pp.
1
18
.http://salsahpc.indiana.edu/dlib/articles/00000714/3/
26.
Noblesse
,
F.
,
1982
, “
The Green Function in the Theory of Radiation and Diffraction of Regular Water Waves by a Body
,”
J. Eng. Math.
,
16
(
2
), pp.
137
169
.
27.
Liu
,
Y.
, and
Falzarano
,
J. M.
,
2017
, “
Irregular Frequency Removal Methods: Theory and Applications in Hydrodynamics
,”
J. Mar. Syst. Ocean Technol.
,
12
(
2
), pp.
49
64
.
28.
Liu
,
Y.
, and
Falzarano
,
J. M.
,
2017
, “
A Method to Remove Irregular Frequencies and Log Singularity Evaluation in Wave-Body Interaction Problems
,”
J. Ocean Eng. Mar. Energy
,
3
(
2
), pp.
161
189
.
29.
Guha
,
A.
,
2012
, “
Development of a Computer Program for Three Dimensional Frequency Domain Analysis of Zero Speed First Order Wave Body Interaction
,”
Ph.D. thesis
, Texas A&M University, College Station, TX.https://www.researchgate.net/profile/Amitava_Guha2/publication/272682470_Development_Of_A_Computer_Program_For_Three_Dimensional_Frequency_Domain_Analysis_of_Zero_Speed_First_Order_Wave_Body_Interaction/links/55e4e9bc08aede0b57358484.pdf
30.
Liu
,
Y.
, and
Falzarano
,
J. M.
,
2016
, “
Suppression of Irregular Frequency in Multi-Body Problem and Free-Surface Singularity Treatment
,”
ASME
Paper No. OMAE2016-54957.
31.
Liu
,
Y.
, and
Falzarano
,
J. M.
,
2017
, “
Frequency Domain Analysis of the Interactions Between Multiple Ships With Nonzero Speed in Waves or Current-Wave Interactions
,” ASME Paper No. OMAE2017-62322.
32.
Liu
,
Y.
, and
Falzarano
,
J. M.
,
2017
, “
A Note on the Conclusion Based on the Generalized Stokes Theorem
,”
J. Offshore Eng. Technol.
, in press.
33.
Liu
,
Y.
, and
Falzarano
,
J. M.
,
2017
, “
Improvement on the Accuracy of Mean Drift Force Calculation
,”
ASME
Paper No. OMAE2017-62321.
34.
Lee
,
C.
,
1995
, “
WAMIT Theory Manual
,” Department of Ocean Engineering, Massachusetts Institute of Technology, Cambridge, MA.