A semi-analytical computation of the Kelvin kernel for potential flows with a free surface

A semi-analytical computation of the Kelvin kernel for potential flows with a free surface

Jorge D'elía^*; Laura Battaglia^*; Mario Storti^*

Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC) Instituto de Desarrollo Tecnológico para la Industria Química (INTEC) Universidad Nacional del Litoral - CONICET Güemes 3450, 3000-Santa Fe, Argentina E-mails: [e-mail] / [e-mail] / [e-mail] / web page: http://www.cimec.org.ar

ABSTRACT

A semi-analytical computation of the three dimensional Green function for seakeeping flow problems is proposed. A potential flow model is assumed with an harmonic dependence on time and a linearized free surface boundary condition. The multiplicative Green function is expressed as the product of a time part and a spatial one. The spatial part is known as the Kelvin kernel, which is the sum of two Rankine sources and a wave-like kernel, being the last one written using the Haskind-Havelock representation. Numerical efficiency is improved by an analytical integration of the two Rankine kernels and the use of a singularity subtractive technique for the Haskind-Havelock integral, where a globally adaptive quadrature is performed for the regular part and an analytic integration is used for the singular one. The proposed computation is employed in a low order panel method with flat triangular elements. As a numerical example, an oscillating floating unit hemisphere in heave and surge modes is considered, where analytical and semi-analytical solutions are taken as a reference.

Mathematical subject classification: Primary: 33F05; Secondary: 65N38.

Key words: green function, boundary integral equation, three dimensional potential flow, free surface, computational techniques.

1 Introduction

In seakeeping flow problems for ship hydrodynamics, a rigid body placed on the free surface of an incompressible inviscid fluid can oscillate in any of the six degrees of freedom around its mean position due to a passing front wave [1]. The standard potential flow theory assumes that the motion is relatively small and harmonic in time [2, 3].

The classical analysis with a linearized free surface boundary condition splits the problem into seven parts. First, six radiative modal potentials Φ_k(x, t) have to be determined, for k = 1, 2, ... 6, where the rigid body performs imposed small harmonic oscillations in each degree of freedom, where x is the position vector and t is the time. Next, a diffraction potential Φ₇(x, t), due to a passing harmonic monochromatic wave of small amplitude, has to be found. These modal velocity potentials Φ^k(x, y, z, t), for k = 1, 2, ..., 7, are found by solving seven boundary integral equations, where the left hand sides have the same integral operator and only the independent terms are specific for each mode, e.g. see [4].

As it is well known, boundary element methods, or panel methods [5], are a natural choice for obtaining numerical solutions of boundary integral equations [6] through collocation or Galerkin techniques [7], as well as they are closely related to the Green function theory [8].

The Green function (x, t) for seakeeping is expressed as the product of a time factor T (t) and a spatial G (x) one. Since the incident front wave is assumed to be monochromatic in time, with absolute circular frequency ω, then, the time factor takes the simple form T (t) = e^iωt, and all computations can be performed in the frequency domain. The spatial part of the Green function G(x) is also known as the Kelvin kernel which, in turn, is decomposed into the sum of two Rankine kernels and a wave kernel. Both Kelvin and Rankine kernels are widely used in numerical ship hydrodynamics, although neither of them satisfy the slip boundary condition over the wetted hull surface and, consequently, such condition must be enforced for a numerical computation.

On the one hand, the Rankine kernel has rather simple mathematical properties; however, it does not satisfy either the outgoing radiation or the free surface boundary conditions. Thus, a finite portion of the free surface must be also discretized in order to impose these missing conditions. Aside from these drawbacks, the Rankine kernel gives a great advantage in unsteady potential flow problems with non-linear boundary conditions [9, 10].

On the other hand, the use of the Kelvin kernel avoids the discretization of the free surface, and the outgoing radiation boundary condition is automatically satisfied. However, it involves several rather elaborated mathematical expressions and tends to be ill-conditioned for field points nearby the axisymmetric axis of the local cylindrical frame at each panel, which is a serious numerical drawback, particularly in hull meshes with a relatively high number of panels.

Similar approaches have also been considered by Telste-Noblesse [11] and by Ponisy et al. [12]. For instance, in [11] eight expressions were proposed in complementary regions of the spatial coordinate d, between the field point and the mirror image of the source point in the mean sea plane, given by: two asymptotic expansions for large values of the distance d, two ascending series for small values of d, two Taylor series around the vertical axis, and two expressions for intermediate values of d.

In this work, a computation of the Kelvin kernel is proposed through a singularity substraction technique, where the boundary integral is split into the sum of a regular term and a singular one. For the regular term, a globally adaptive numerical quadrature is employed, while for the singular one an analytic integration is performed. The proposed computation is performed with a low order panel method where only the wetted surface of the body in hydrostatic state is discretized with flat triangles. As a numerical example, the oscillating floating hemisphere of unit radius in heave and surge modes is considered, for which there are analytical and semi-analytical solutions.

2 Seakeeping flow problem

2.1 Differential formulation

A Cartesian (x, y, z) coordinate system is chosen, where the z = 0 plane matches the still water plane and the z-axis is positive upwards. The complex e^iωt dependency of the time t is implicitly assumed, where ω is the circular frequency of the periodic motion.

An infinitesimal rigid body oscillating in the k mode and placed under the free surface of a fluid without a uniform mean current, is described with the linearized governing equation [13]

where Δ = ∂_xx + ∂_yy + ∂_zz is the three-dimensional Laplacian operator, Φ_k is the k-modal radiation potential, and K = ω²/g is the wave-number for gravity waves in deep water.

2.2 Boundary integral equation

A boundary integral equation for solving Eq. (1) is given by [2]

for x ∈ S, where x and ξ are the field and source points, respectively, and S is the boundary of the flow domain ω. The independent term is

while Φ_k(x), for k = 1, ..., 6, is the k-radiation velocity potential, and σ_k are known fluxes. A standard panel method imposes the integral boundary equation (2) by means of a collocation technique at the panel centroids, obtaining a complex valued linear system AΦ_k= Cσ_k= b_k, where Φ_k is the k-velocity potential vector, and σ_k is the k-flux vector corresponding to the k-mode. The dipolar matrix, which is non-symmetric and regular, is given by

and the monopolar one, symmetric, is

Both the monopolar C and the dipolar A matrices are square and full populated. They include the spatial Green function G (x, ξ) and the normal derivative G _n(x, ξ), respectively.

2.3 Kelvin andRankine kernels in the spatial Green function

The spatial Green function G(x, ξ) in Eq. (5) that satisfies Eq. (1), is known as the Kelvin kernel, which gives the interaction between the field point x = (x, y, z) and the source point ξ = (ξ, η, ζ) [14]. The physical meaning of the Green function is given by the real part Re {Ge^iωt}, which is the disturbed velocity potential measured at the field point x, at time t, caused by a pulsating source ξ of circular frequency ω and unit intensity [1]. It should be noted that the outgoing radiation and free surface boundary conditions are automatically satisfied by the Kelvin kernel.

Due to the local axisymmetry around the source point ξ, it is convenient to introduce the non-dimensional cylindrical coordinates

where X is the radial coordinate and Y the vertical one. Then, the Kelvin kernel for seakeeping is written as

where

are the Euclidean distances between the field point x and the source point ξ, and between the field point x and the image point ξ' = (ξ, η, -ζ), respectively. In Eq. (7) the first two terms, r^-1 and s^-1, are the Rankine kernels, while the term inherits the spatial wave properties of the Kelvin kernel and, then, it is termed the "wave-kernel".

2.4 Haskind-Havelock representation of the wave-kernel

The wave kernel involves several transcendental functions and, consequently, the computational cost can be rather expensive. Moreover, the Haskind-Havelock representation tends to be ill-conditioned for field points located near the axisymmetric axis, as well as for points in far away regions. Thus, a semi-analytical integration strategy is proposed as a compromise solution between numerical cost and complicated mathematical expressions, especially in numerical simulations with non-linear boundary conditions, e.g. when there is a mesh motion and the Jacobian of the system matrix is required.

The Haskind-Havelock representation for the wave part of the Kelvin kernel is written as [14]

where H₀( X) is the Struve function of zero order, J₀( X) and Y₀ (X) are the zero order Bessel functions of first and second kind, respectively [15], and P₀( X, Y) is the Haskind-Havelock integral [14]

The asymptotic behavior of the Kelvin kernel given by Eq. (7), at very low and very high frequencies, will be dealt with in Secs. 4.2 and 4.3.

3 Evaluation of the Kelvin kernel

3.1 Rankine kernels

The Rankine kernels r ^-1 and s ^-1 can be evaluated in several ways. Onepossibil-ity is a numerical integration, which has the advantage that high order distributions can be considered without further complications. However, the numerical integration is rather sensitive to the mesh quality and, moreover, the diagonal terms would deserve a special treatment. Another alternative is an analytic integration, where the surface integral over each panel is replaced by its closed contour integration and a side local reference frame is used for each side contri-bution[16, 17, 18].

3.2 Normal derivative of the Haskind-Havelock kernel

The normal derivative of the Haskind-Havelock kernel is found from ,_n = (∇_ξ,n_ξ), where n = (n_ξ,n_η,n_ζ) is the unit normal of dS_ξ and is the gradient of G, both evaluated on the source point ξ = (ξ,η,ζ). By the chain rule in Eqs. (6) and (9)

where

Note that the gradients of the wave-kernel of the Green function, evaluated on the field point x = (x, y, z) and the source point ξ = (ξ, η, ζ) are linked as

The complex kernel is

where the real part, Re {..} = (..)', and the imaginary one, Im {..} = (..)", are given by

The partial derivatives of are

and the corresponding ones of , with J₀ = d J₀(X)/dX,

3.3 Ill-conditioning of the Haskind-Havelock kernel

The Haskind-Havelock finite integral is given by

and its partial derivatives are

The Haskind-Havelock finite integral P₀ evaluated at t = 0 tends to be ill-conditioned when X << 1, that is, for field points near the axisymmetric axis. This is a serious numerical drawback, in particular in hull meshes with a high number of panels. For overcoming this disadvantage, a singularity subtraction technique is proposed, where the integral is split into the sum of a regular term and a singular one. For the regular term, a globally adaptive numerical quadrature is employed, while for the singular one an analytic integration is performed. On the other hand, a direct computation of the Struve functions H₀ and J₀ can be performed through their definitions and asymptotic expansions.

4 Semi-analytical computation of the Haskind-Havelock kernel

4.1 Singularity subtraction technique

The Haskind-Havelock integral given by Eq. (18) is split into the sum of a regular term and a singular one. For the regular term, a globally adaptive numerical integration can be used, while for the singular one an analytic integration is performed. Thus, Eq. (18) is rewritten as

where

is a regular integral which can be evaluated accurately by a globally adaptive integration, for example, the qag routines of the Netlib Repository (http://www.netlib.org). The remaining integral

contains a logarithmic singularity when X = 0, and it is ill-conditioned when X → 0. Then, it is evaluated in a closed form by performing the following variable changes

for which

then

The partial X-derivative of Eq. (21) is similarly decoupled as

where

is a regular integral, whereas

is the integral that contains the singularity and it is computed in closed form. The variable change α = X sinh θ is introduced again and

As cosh²θ - sinh²θ = 1, then

The A term is given by

with the variable change

Replacing

and then

Next, the B term is given by

introducing the variable changes

for which

it results

Finally, the trivial C term is

However, when the field point x = (x, y, z) is on the axisymmetric axis of the source point ξ = (ξ, η, ζ), then X = 0 and these expressions are not applicable. In such case, the asymptotic representation [13]

can be used when X << 1, where

with

where Ei (Y) is the exponential integral [15]. Then

that is

which is valid for X << 1. As W_k (X, Y) and its derivatives tend uniformly to zero in a small neighborhood of X = 0, then, Eq. (45) can be written at X = 0 as

In summary, seven expressions are employed for the complementary regions of the coordinates X and Y given by: (i) two regular integrals given by Eqs. (22) and (28), (ii) two analytic integrals given by Eqs. (26) and (31), for X >> 1, with the constants A, B and C given by Eqs. (35), (39) and (40), and (iii) three approximate expressions given by Eq. (46) for X << 1.

4.2 Kelvin kernel at very low frequencies or near the vertical axis

For very low frequencies K << 1, or in the neighborhood of X = 0, it can be shown that the terms H₀, Y₀ and P₀ in Eq. (9) tend to cancel each other out and, consequently, the square bracket results bounded [...] < D, with D as a constant independent of K. Then, when the field point x = (x, y, z) does not match the source one ξ = (ξ, η, ζ), the Kelvin kernel has the asymptotic form G →r^-1 + s^-1 for K →0 (low frequencies) or X →0 (near the vertical axis).

4.3 Kelvin kernel at very high frequencies or points far away

For very high frequencies K >> 1 or points far away from the origin, it is verified that X >> 1 and, then, the following expansion can be used [1]

which is valid for X >> 1. From this,

Moreover (see Abramowitz-Stegun [15])

therefore J_n, Y_n, H₀ ~ O (X^-1/2) and

Thus, for very high frequencies K >> 1 or points far away from the origin, it is verified that, for X >> 1, the wave kernel has the asymptotic form . Then, the Kelvin one has the asymptotic form G →r^-1 - s^-1 for K >> 1 (high frequencies) or X >> 1 (far away).

4.4 Direct computation of the special functions

The Struve differential equation is (see chap. 9, Abramowitz-Stegun [15])

where z = x + iy is the complex variable and γ is the factorial function. Its general solution is

where J_v and Y_v are the Bessel functions, of first and second kind, respectively, H_v is the Struve one, all these of integer order Ν, and a, b are constants. A direct computation involves ascending and descending series. The ascending series for the Bessel function of first class J_n (X) and order n = 0, 1 are

The corresponding ones for the Bessel function of second kind Y_n (X) are

and

where

with ψ_k+3/2 = ψ_k+1 + ψ_k+2 and

where γ = 0.5772 ... is the Euler constant and ψ₁ = -γ. On the other hand, for the Struve functions H_n (X) of order n = 0, 1 (see chap. 9, Abramowitz-Stegun [15])

and

When the abscissa X is far from the origin, these series show slow rate of convergence and numerical instability. Therefore, they are replaced by the asymptotic expansions

for X > 25, while for the Struve ones the expressions adopted are

where

for X > 30. The derivatives of Eqs. (55) and (65) with respect to X are, respectively,

Finally, an asymptotic expansion for the exponential integral Ei( Y) (see chap. 5 [15]) is

Plots in Fig. 2 show: (i) the Bessel functions of the first kind J₀ and J₁ (top-left), (ii) the Bessel functions of the second kind Y0 and Y1 (bottom-left), (iii) the Struve functions H0 and H1 (top-right), and (iv) the exponential integral Ei ( x) (bottom-right).

5 Numerical examples

5.1 Free surface test

The Kelvin kernel computation is validated through a numerical test, where the free surface boundary condition

is explicitly computed on a grid over the plane z = 0, inside a finite region |x, y|/L < 200. The source is a square panel of side length L = 0.1, submerged at depth H and harmonically pulsating at frequency Ω. Thus, Eq. (71) is verified at machine precision. In Fig. 3 (top), a three-dimensional view of the wave pattern on the plane z = 0 produced by the submerged source panel is shown.

5.2 Body coordinates and motions

In seakeeping, the body coordinate system (X, Y, Z) is fixed to the hull, i.e. it moves together with it. The Z-axis is upwards, the X-axis is to bow and, when there is not motion, the plane Z = 0 matches the still water plane, as represented in Fig. 4 (left). The harmonic body motion is given by the instantaneous position of the body coordinate system ( X, Y, Z) with respect to the moving-frame ( x, y, z) and it is decomposed into surge, sway and heave oscillating translations along the body-axes, and into roll, pitch and yaw oscillating rotations around the body-axes, see Fig. 4 (left).

5.3 Oscillating floating unit hemisphere

An oscillating hemisphere of radius R = 1 in surge (k = 1) and heave (k = 3) modes is considered. Due to the symmetry between the surge and sway (k = 2) modes, it is not necessary to determine the sway mode, nevertheless, for a code validation, it was also verified at machine precision, in a perfectly symmetric mesh. The added mass coefficients are computed as A'_kk = A_kk/(ρ V), and the damping ones as D'_kk = D_kk/(ρ Vω), where V = (2/3)πR³ is the hemisphere volume, ρ is fluid density and ω is the imposed circular frequency. Figure 3 (bottom) shows a convergence plot for the surge added mass at very low frequencies (ω → 0), obtained with a lineal regression analysis (exact value A'_n = 1/2). The mesh is shown in Figure 4 (right) and has 3 000 Kelvin panels over the wetted body surface. Plots of the added mass A'_kk and damping coefficients D'_kk, as a function of the wave number coefficient KR, are shown for the surge and sway modes in Figure 5 (left), and for the heave mode in Figure 5 (right). All of them are in good agreement with the literature results [19]. The asymptotic values of these coefficients, for very low and very high frequencies, can be obtained analytically, e.g. by variable separation or image methods. For the surge/sway mode at very low frequency, the boundary condition φ_,z = 0 is equivalent to a symmetry operation with respect to the plane z = 0 and, then, corresponds to the solution of a sphere oscillating in an infinite medium. The added mass for the last case is half of the displaced volume, then, the surge/sway added mass coefficient is A'₁₁ = 1/2 with respect to the true displaced mass (2/3)πR³ρ, where the half factor is due to the analytic prolongation. On the other hand, the asymptotic values of the added mass in heave mode are nΠot easy to obtain and could be computed with spherical harmonics (e.g. see Storti-D'Elia [20]). Bounds for the surge A'_n and heave A3₃ added mass coefficients of the oscillating unit hemisphere at very low and very high frequencies are summarized in Table 1. The first column corresponds to those found in Storti-D'Elia [20]. The values for the surge/sway mode in the second column correspond to those found in Sierevogel [21] and Prins [22], while the corresponding ones to the heave mode are taken from Korsmeyer [23] and Liapis [24]. It should be noted that only the intervals [0.25, 1.50] and [0.6, 1.5] were considered in Prins [22], respectively, and, then, the extrapolations are rather doubtful. The third column corresponds to the results found in Hulme [25].

Korsmeyer [23] used a panel method with Fourier transform and complex impedance extended to very low frequencies, while Hulme [25] used spherical harmonics. The Sierevogel [21], Prins [22] and Liapis [24] results were obtained with other panel methods with Kelvin kernels. In general, the concordance among the present results and the literature is good. Another validation test may include a comparison with related approaches found in the literature of ship engineering, such as those based on Chebyshev polynomials, e.g. the WAMIT software (http://www.wamit.com).

6 Conclusions

A semi-numerical scheme for computing the Kelvin kernel for seakeeping flow problems has been proposed. The Kelvin kernel is decomposed as the sum of two Rankine sources and a wave one. The Rankine sources are the standard Green functions for the Laplacian equation, one due to the generic panel on the body surface, placed below the plane z = 0, and the other one due to the mirror image with respect to the same plane. The wave kernel (i) tends to be ill-conditioned for field points near or over the local axisymmetric axis; and (ii) involves a rather heavy computation, due to the Haskind-Havelock integral which, in turn, involves the computation of Bessel and Strouve functions. The Haskind-Havelock integral was accurately computed with a singularity subtraction technique that involves a regular closed term and a numerical adaptive quadrature, while the Bessel and Strouve functions were calculated with asymptotic expansions. The proposed semi-numerical scheme was validated with analytical and semi-analytical solutions for the unit hemisphere in surge and heave motions, without showing numerical instabilities nor precision loss.

Acknowledgements. This work has received financial support from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina, grant PIP 5271-05), Universidad Nacional del Litoral (UNL, Argentina, grant CAI+D 2009-III-4-2), Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, Argentina, grants PICT 1506-06, PICT 1141-07 and PAE 22592-04 nodo 22961) and was performed with the Free Software Foundation/GNU-Project resources such as GNU-Linux OS, GNU-Gfortran and GNU-Octave, as well as other Open Source resources as Scilab, TGif, Xfig and ETf-X. The authors thank the referees for their constructive suggestions and careful reading.

Received: 25/XI/09.

Accepted: 08/II/10.

#CAM-155/09.

• [I] J.N. Newman, The theory of ship motions. Advances in Applied Mechanics, 18 (1978), 221-285.
• [2] J. Nossen, J. Grue and E. Pato, Wave forces on three-dimensional floating bodies with small forward speed. J. of Fluid Mechanics, 227 (1991), 135-160.
• [3] S. Finne and J. Grue, On the complete radiation-difraction problem and wavedrift damping ofmarine bodies in the yaw mode of motion. J. of Fluid Mechanics, 357 (1998), 289-320.
• [4] R.F. Beck and W.C. Webster, Seakeeping and controllability. In: E.V. Lewis, editor, Principles of Naval Architecture, volume III. SNAME (1989).
• [5] F. París and J. Canas, Boundary Element Method. Fundamentals and applications. Oxford Press (1997).
• [6] W. Hackbusch, Integral equations. Birkhäuser (1995).
• [7] J. D'Elía, L. Battaglia, A. Cardona and M. Storti, Full numerical quadrature ofweakly singular double surface integrals in Galerkin boundary element methods. Int. J. for Num. Meth. in Biomedical Engng., 27(2) (2011), 314-334. doi:10.1002/cnm.1309.
• [8] F. Hartmann, Introduction to Boundary Elements. Springer-Verlag (1989).
• [9] Y. Huang and P.D. Sclavounos, Nonlinear ship motions. J. of Ship Research, 42(2) (1998), 120-130.
• [10] D.C. Kring, Ship seakeeping through the Τ = 1/4 critical frequency. J. of Ship Research, 42(2) (1998), 113-119.
• [II] J.G. Telste and F. Noblesse, Numerical evaluation of the Green function of water-wave radiation and difraction. J. of Ship Research, 30(2): 69-84, June (1986).
• [12] B. Ponizy, F. Noblesse, M. Ba and M. Guilbaud, Numerical evaluation ofthe free-surface Green functions. J. of Ship Research, 38 (1994), 193-202.
• [13] F. Noblesse, The Green function in the theory of radiation and difraction of regular water waves by a body. J. of Engineering Mathematics, 16 (1982), 137169.
• [14] J.N.Newman, Algorithms for the free-surface Green function. J. of Engineering Mathematics, 19 (1985), 57-67.
• [15] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions. Dover (1972).
• [16] D.E. Medina and J.A. Liggett, Three-dimensional boundaryelement computation of potential flow in fractured rock. Int. J. for Num. Meth. in Eng., 26 (1988), 2319-2330.
• [17] J. D'Elía, MA. Storti and S.R. Idelsohn, A closed form for low order panel methods. Advances in Engineering Software, 31(5) (2000), 335-341.
• [18] J. D'Elía, M.A. Storti and S.R. Idelsohn, Iterative solution of panel discretizations for potential flows. The modal/multipolar preconditioning. Int. J. Num. Meth. Fluids, 32(1) (2000), 1-22.
• [19] A. Papanikolau, On the integral-equation-methods for the evaluation ofmotions and loads of arbitrary bodies in waves. Ingenieur-Archiv, 55 (1985), 17-29.
• [20] M.A. Storti and J. D'Elía, Added mass ofan oscillatinghemisphere at very-low and very-high frequencies. ASME, Journal of Fluids Engineering, 126 (6), 10481053, November (2004).
• [21] L. Sierevogel, Time-domain Calculations of Ship Motions. PhD thesis, Technische Universiteit Delft (1998).
• [22] .J. Prins, Time-domain Calculations of Drift Forces and Moments. PhD thesis, Technische Universiteit Delft (1995).
• [23] F.T. Korsmeyer and P.D. Sclavounos, The large-time asymptotic expansion of the impulse response function for a floating body. Applied Ocean Research, 11(2) (1989), 75-88.
• [24] S.J. Liapis, Time-domain Analysis of Ship Motions. PhD thesis, University of Michigan (1986).
• [25] A. Hulme, The wave forces acting on a floating hemisphere undergoing forced periodic oscillations. J. Fluid Mechanics, 121 (1982), 443-463.

*

Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET): PIP 5271-05. Universidad Nacional del Litoral (UNL): CAI+D 2009-III-4-2. Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT): PICT 1506-06, PICT 1141-07 and PAE 22592-04 nodo 22961.

Publication Dates

• Publication in this collection
  27 July 2011
• Date of issue
  2011

History

• Accepted
  08 Feb 2010
• Received
  25 Nov 2009