Suppression of Irregular Frequency Effect in Hydrodynamic Problems and Free-Surface Singularity Treatment

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⁻⁷ 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.

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

The Cartesian coordinate setting is indicated in Fig. 1: V_– stands for the internal volume of the floater, bounded by the boundary surface S_b and S_i; S_f stands for the outside free surface; V is the fluid domain, bounded by free surface S_f, body surface S_b, bottom S_B, and infinite control surface S_c. n is the unit normal vector pointing outside fluid domain V, inside the floater body.

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+ζ)$⁠, J₀ is the Bessel function of the first kind of order 0, and R₀ is the function defined in Telste and Noblesse [24].

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

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.

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].

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.

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.

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.

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+ζ)$⁠, J₀ is Bessel function of the first kind at order 0, and R₀ 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.

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/r³ 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.

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.

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.

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

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.

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⁻⁷. 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].

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

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

As shown in Figs. 3–6, 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 = ω²L/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.

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.

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 7–14 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.

Figures 16–21 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.

Figures 23–26 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.

Figures 28–31 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.

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 33–36 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.

Figures 38–43 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.

Figures 45–50 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.

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.

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).