A boundary integral equation method for free-surface flow problems

MECHANICS RESEARCH COMMUNICATIONS Vol. 20(3), 227-235, 1993. 0093-6413/93 $6.00 + .00

Printed in the U.S.A. Pergamon Press Ltd.

A BOUNDARY INTEGRAL EQUATION METHOD FOR FREE-SURFACE FLOW PROBLEMS

S. Liapis Aerospace and Ocean Engineering, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA (Received 24 August 1992; acceptedfor print 9 December 1992)

INTRODUCTION This study uses a boundary integral equation method to compute water-wave loads on floating structures. A source distribution is used to represent the flow around the body. The source distribution is submerged inside the body, unlike the traditional approach where the source distribution lies on the body surface. The idea of submerging the source distribution is old and dates back to the classical work of Von Karman [1] where an axial source distribution was used to evaluate the flow about axisymmetric bodies. Webster [2] studied infinite flow problems using a piecewise linear source distribution which was submerged below the body surface. In more recent work Jensen [3] has solved nonlinear ship wave resistance problems using Rankine sources located above the free surface. Submerging the source distribution results in general to a smoother and more accurate solution. It is also possible to use simple sources and sinks located inside the body instead of a continuous source distribution therefore eliminating the intermediate integration step [4],[5],[6]. However, it is unclear whether this procedure can produce accurate results for bodies of arbitrary shape especially in areas of high curvature. As is well known for free surface problems, the integral equation breaks down at certain critical frequencies also known as "irregular" frequencies. These frequencies correspond to eigenvalues where the interior Dirichlet problem has a non-trivial solution. This phenomenon was first pointed out by John [7] and verified numerically by Frank [8] for oscillating cylinders in two dimensions. It must be emphasized that the irregular frequencies are not caused by any physical phenomenon but are related to the solution procedure. Large numerical errors occur in a frequency bandwidth close to each irregular frequency and the fact that these frequencies are not a priori known for a given body adds to the severity of the problem. Several methods for removing the irregular frequencies have been proposed some of them motivated from acoustics where a similar non-uniqueness problem occurs. A survey of different methods used may be found in Lau and Hearn[9].

227

228

S. LIAPIS

In the present paper the advantages of using a submerged source distribution for freesurface problems are examined. Results using this method were obtained for a floating hemisphere oscillating in heave and surge and a Wigley hull ship oscillating in heave and compared with exact results and results obtained by the conventional method. As the results suggest, submerging the source distribution results in a smoother solution which is remarkably more accurate at no extra computational cost. Furthermore, submerging the source distribution significantly reduces the errors close to the "irregular" frequencies where the numerical solution breaks down. This remarkable improvement in the solution accuracy is mostly attributed to the fact that employing a piecewise constant source distribution results in a solution which is logarithmically singular at the panel edges. Submerging the source distribution removes the singularities away from the computational domain.

MATHEMATICAL FORMULATION

Consider a rigid body floating on the free surface of deep water. Plane progressive waves interact with the body which undergoes small sinusoidal oscillations about a mean position. An xyz coordinate system is defined with z positive upwards and the xy plane coincident with the calm free surface. Under the usual assumptions of an inviscid, irrotational flow the problem may be defined by a velocity potential ¢(x,y,z,t). It is convenient to decompose into the sum of the incident potential,the diffraction potential and the six radiation potentials as : 6

O(x,y,z,O = Re [(~o + qb7 + ~.,dPt(x,Y,Z) ~t )e-i''q

(1)

j-I

where ~i is the complex amplitude of the body motion in each of the six degrees of freedord of the body and ~o "

i g A e VZ_ivxcosl~_ivysinf~ f.o

is the potential of the incident waves of amplitude A, frequency co, wavenumber v and direction of propagation ~. In the fluid domain each of the potentials must satisfy the Laplace equation:

v2¢j - o

(2)

229

BOUNDARY INTEGRALS WITH FREE-SURFACE FLOW = 4

The boundary condition on the free surface is in its linear form: -v~j +

&

- 0

on

z - 0

(3)

On the body surface OD each of the potentials must satisfy:

~7

~o

On

an

~J

vg~ (4)

io)+ = v p 9

On

j = 1,2....6

where n i are the components of the generalized normal directed out of the fluid domain "defined by the relations : (nl,n2,n3) - g, with

(n4,ns,n6) - ?' x g, r - (x,y,z)

in addition to the above, appropriate radiation conditions at infinity are required to make the solution unique. The Green function for this problem which is the potential of a submerged pulsating source is given by Wehausen and Laitone [10,13.17] as Re[G(P,Q)e -~*t] with G(P,Q) ffi 1 where

÷

+ 2v /0

1

P - (x,y,z)

ek(Z+OJo(k~dk

is the field point

r2 _ (x_~)2 + (y_n)2 + ( z _ 0 2 , R2

+ 2~ive v¢z+OJo(vR)

Q " (~Jl,0

is the source poim

r~ - (x-~) 2 + (y-n)~ + ( z + O ~

(x_~)2 + (y_~)2,

v - - - is the wavenumber g Jo is the Bessel function of the first kind denotes principal value integral

(5)

230

S. LIAPIS

The potential may be represented in terms of a source distribution. We write:

%(I")- fsf

o(Q)G(P,Q)dSQ

(6)

where the sources are distributed on a closed surface S which, unlike the traditional approach, is submerged below the actual body surface OD. Following Cao et al [4] the normal distance L d between the two surfaces is given by the formula : -

.

(7)

where D m is a measure of the panel size (the panel length in two dimensions and the square root of the area in three dimensions). The optimum value of the constant depends on the problem but in general a is less than 1. Special care must be taken so that the submergence does not exceed the local radius of curvature due to the well known analytical result that a singularity sheet distribution submerged beneath the centers of curvature cannot represent the flow past the body. An integral equation for the source strength may be obtained by differentiating equation (6) in the normal direction and using the body boundary condition to yield:

fj

a(Q) OG(P,Q)d$0 . Vj(P)

Onp

(s)

From the values of the potential on the body surface, the added mass and damping coefficients may be obtained. Following the standard definition given in Newman [11] the added mass and damping coefficients are:

Aij - Re

[( - p

f a J ~J(Q)ni(Q)dSco/i(o]

Bij - Re tp o f-of

(9) i, j = 1,2....6

The integral equation (8) is solved numerically using a panel method. The surface S is approximated by an ensemble of flat quadrilateral panels of constant source strength. This reduces the problem of finding a continuous source distribution to determining a finite number of unknown source strengths. These source strengths are determined by collocation where the integral equation (8) is satisfied at one point for each panel.

BOUNDARY INTEGRALS WITH FREE-SURFACE FLOW

NUMERICAL RESULTS To illustrate the performance of the numerical method, the added mass and damping coefficients for a floating hemisphere oscillating in heave and surge and a Wigiey hull ship oscillating in heave have been computed. Since the sphere has two planes of symmetry, only one-quarter of the underwater surface needs to be discretized. The special case of an infinite fluid flow past a sphere was considered as the first test problem. This problem corresponds to the limiting case of heave motion for v--~,. The numerical results obtained are compared to the exact solution which, as is well known, is the flow generated by a dipole. Figure 1 presents the RMS error in the computed values of the potential as a function of the total number of panels (panels used *8). The source panels are submerged into the body their submergence distance given by equation (7). Results were obtained for several different submergence distances corresponding to different values of the constant a. As shown in figure 1, submerging the panels leads to a substantial improvement of the accuracy and convergence of the source method. The accuracy of the numerical solution is greatly improved even for small values of the submergence distance. Turning now to the free surface flow problem case, figures 2 and 3 show the nondimensional added mass and damping coefficients for a sphere in heave and surge versus the nondimensional parameter vR where R is the radius of the sphere. 65 plane quadrilateral elements were used to discretize one quarter of the underwater surface. Because of the simple geometry of the sphere there exist published results using simplified approaches [12],[13] which can be used as a benchmark for checking the solution accuracy. The asterisks in figures 2 and 3 are in the case of heave, results obtained by Baracat [12] using the method of multipoles and in the case of surge, results obtained by Hulme [13] using ring sources and a one-dimensional integral equation. The dashed line curves correspond to results obtained by using the traditional approach where the source panels are on the body surface. In the vicinity of the irregular frequencies substantial numerical errors occur. The solid line curves are results obtained by using source panels submerged into the body their normal distance from the body surface given by equation (7) with a--0.4. Figures 2 and 3 verify that submerging the source distribution produces a solution which is significantly more accurate. This is particularly true in the high frequency range of the added mass results. It was found that when using the traditional method the number of panels must be increased to at least 250 to achieve an accuracy comparable to that of the submerged source approach. The influence of the irregular frequencies is present in all cases. However, as figures 2 and 3 indicate, submerging the source panels reduces considerably the discontinuity close to the irregular frequencies both in magnitude and frequency bandwidth. Given that irregular frequencies correspond to the eigenfrequencies of the Dirichlet problem for the domain interior to the source surface S the location of the irregular frequencies will change with the submergence distance. For the case of a sphere the irregular wavenumbers are given by the relation vR S = const where R S is the radius of the surface S. Therefore by decreasing the radius R S of the source surface S the irregular

231

232

s. LIAPIS

frequencies will increase. It must be noted that the dependence of the critical wavenumber location entirely on the dimensions of the surface S may lead to a practical way of circumventing the problem for the case of a body of arbitrary geometry where the location of the critical frequencies is not known. One can always distribute sources on an interior surface S of simple geometrical shape where the critical wavenumbers are known a priori and avoid these wavenumbers in the numerical computations. The last example considered was a Wigley hull ship. The Wigley hull ship has parabolic sections and waterlines, a length to beam ratio of 10 and a beam to draft ratio of 1.6. Figure 4 presents the nondimensional added mass and damping coefficients for the ship oscillating in heave. The asterisks are values found by extrapolating results obtained by the conventional method with 320 and 500 panels assuming an error inversely proportional to the number of panels and will be used as a benchmark for numerical accuracy. The dashed line curves are results using the traditional approach with 80 panels on one-quarter of the underwater surface while the solid line results correspond to submerging the source panels according to equation (7) with ~, = 0.4. As figure 4 indicates, similar to the sphere case, submerging the source panels significantly improves the solution accuracy at all frequencies as well as reducing the solution misbehaviour close to the irregular frequencies.

CONCLUSIONS An investigation of the merits of using a submerged source distribution to evaluate water-wave structure interactions was undertaken. It was found that submerging the source distribution results in a much more accurate solution at no extra computational cost. The proposed modification can be easily incorporated in existing boundary integral codes resulting in much more accurate numerical predictions.

REFERENCES

[1] T. Von Karman, NACA Technical Memorandum No. 574, (1930). [2] W.C. Webster, J.Ship Res.,19, 206 (1975). [3] J. Jensen, Proc. 3d Int. Workshop on Water Waves and Floating Bodies, 83 (1988). [4] Y.Cao, W.W.Schultz and R.F.Beck,Intl.J.Num.Meth.Fluids,12, 785 (1991). [5] P.S.Han and M.D.Olson, Int.J.Num.Meth.Engrg.,24, 1187 (1987). [6] R.L. Jonhnston and G.Fairweather,Appl.Math.Model.,8,265 (1984). [7] F.John,Comm. Pure Appl. Math.,2,13 (1950). [8] W.Frank,Rep. 2375, David Taylor Research Center (1967). [9] S.M.Lau and G.Hearn, Intl.J.Num.Meth.Fluids,9,763 (1989).

BOUNDARY INTEGRALS WITH FREE-SURFACE FLOW

233

[10] J.V.Wehausen and E.V.Laitone,Handbuch der Physic,9,446 (1960). [111 J.N.Newman, Marine Hydrodynamics,3d Ed.,MIT Press, Cambridge Mass.(1980). [12] R. Baraeat,J. Fluid Mech.,13, 540 (1962). Also corrections in an unpublished report entitled : Forced periodic heaving of a semi-immersed sphere. [13] A. Huime, J. Fluid Mech. ,121, 443 (1982). ~'~ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 .I

i:

~ ...........

~-. . . . . . . . . . .

t I I

I I I

a----O

"

~

.....

¢~ = 0.1 a=0.2 ...............

°:o.4 i1

......

,

. . . . . . . . . . .

°.t

.~ . . . . . . . . . . . .

.

.

.

.

.

,

1

I

z" . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

-~

i

I- - ~

,~ ............

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

l I

I i,

0 ............

ZOO

I

I

I

I

I

I

"~----------'~300

- - - -.-.- ; - - - , - r

4.00

- - -,- - - ' ~ 7 -.-,° - ~ - - . - -

500

600

N u m b e r o f panels

Figure

I

Efl'eet o f the source distribution submergence for a sphere moving in an infinite fluid

~"": -.-.- ~",

700

234

S. LIAPIS

0401

0 9¸ Baracat [12] * * * * * Present calculations a=0 ......

0 8~ A33

p*V

o

055 ]

B33

p*V*~

0 ~0

a =0.4

7-

025 0 20

06

0 15 0 5-

0 I0 t 0 05

04

LI

"-.

i

0 O0 03

-o 05 -o~o

02

I.

0

2.

3

Figure 2

0.7[ ,~ 0 6 ],~ AIi

0 5 I~

4

5 vR

6

2.

3

4

5

G

7

Heave added mass and damping coefficients for a sphere

o51

P~s~"nOC~'lculati°ns

~

I

vR

Hulme [13] . . . . .

l

,

0

7

0.4 BIi

~ = 0 . 4 ~

p*V*~ 04"

03

0 :

(3.2

0 1

'd 0. ~

"

0

0.

0.

2.

4

Figure 3

6.

uR

8'.

10.

/ 2

4

6

~:R

Surge added mass and damping coefficients for a sphere

8

10

B O U N D A R Y INTEGRALS W I T H FREE-SURFACE F L O W

235

0s

Li

A33

Exrapolated r e s u l t s Present calculations a=

L,

i.|

* * * * * B33

,

p*V*~.,

0.4

is

p *V ,~ iz i.i i,i

,i

os

oo

Qa

a~

I.s e,4

Io

Figure 4

],

)o

vL

vL

Heave added mass and damping coefficients for a Wigley hull ship