SEA-KEEPING COMPUTATIONS USING THE SHIP MOTION
GREEN'S FUNCTION
BOIN BA M.** et GUILBAUD M.*
Laboratoire d'Etudes Aérodynamiques (UMR CNRS 6609)
*CEAT - Université de Poitiers, 43, rue de l'Aérodrome - 86036 Poitiers Cedex France **ENSp, I rue Clément ADER, BP 40109, 86960 Futuroscope Chasseneuil Cedex France
Key words: Naval Hydrodynamics, Sea-keeping, Diffraction-radiation, Green function,
Frequency domain
Abstract
We present here a seakeeping boundary element method of computations using the diffraction-radiation with forward speed Green function in the frequency domain with an analysis of the
accuracy and computational time needed for the boundary integration. The results of an
analytical surface integration of the Green function and of its first and second derivatives have been compared with a numerical procedure based on a Gauss method and further refinements dependant of the relative location of the source and field points and to their distances to the
free-surface. The results are compared for a single panel and a field point in the whole
domain, which is divided in various zones where different techniques of integration are used.
These surface integration techniques are introduced in a panel method for seakeeping
computations adapted to the case of surface-piercing bodies or ships in forced oscillations.The results are also compared with results available from other calculation methods or from test measurements for the DNV barge or a simple surface-piercing flat plate.
Deift University of
Technology
Ship HydromeChafliCS
laboratory
Library
Mekelweg 2 26282 CD Deift
phone: +31 (0)15 2786873
Introduction
The study of the seakeeping for a ship is a very complex problem and no entirely satisfactory
solution has been yet obtained. Most of the studies are achieved using boundary element method (B.E.M.) assuming inviscid fluid, the main physical phenomena being the wave
production and the RANSE methods requiring very high CPU time for unbounded fluid. To
solve the Laplace equation, the advantage of using the diffraction-radiation with forward speed Green function with respect to the corresponding function in an unbounded domain (Rankine source), is the automatic satisfactory of the free surface boundary and radiation
conditions, leading to lower order system of equations to solve and avoiding any problem with gravity waves reflection on the free surface grid boundaries as in the Rankine methods (for a fine review of Rankine methods, see Sclavounos, 1996). Furthermore, these methods are able
to take into account very small wavelengths (no filtering due to the free surface grid). The
main difficulties are to compute the Green function and its surface integration, but the
progress of the computers and of the computation algorithms, see for example steepest descent method (Ishawita and Ohkusu, 1989 or Brument and Delhommeau, 1997, Brument, 1998), the Super Green functions (Chen and Noblesse, 1998) or an adaptative Simpson method (Ba and
Guilbaud, 1995) have overcame a part of the various difficulties, so it is today possible to
develop such a seakeaping computational code in the frequency domain with a good accuracy and relatively moderate computational time.
We present here the development of such a code. lt is necessary not only to control the accuracy of the computations of the Green function and its first derivatives but also to check the accuracy of the boundary integrals. The results of an analytical surface integration of the function and of its first and second derivatives on planar panels have been compared with a
numerical procedure based on a Gauss method (with different number of Gauss points and
The comparison of the results, based on both the accuracy and the computational time, has led
to the division of the computational domain in various zones where different integration
methods are used. Seakeeping computations using this improved method of boundary inte-gration have been performed and results concerning added mass and damping coefficients and pressure distribution on various models (flat plate, DNV barge) are presented and compared with the results obtained by other methods and test measurements available in the literature.
Problem to solve
A right-handed co-ordinate system translating in the x direction with the steady forward speed
U of the body is chosen, so the undisturbed upstream velocity is U; the xOy plane is the
undisturbed free surface and the z axis is vertically upwards. The water depth is infinite. The
irrotational flow of an incompressible and inviscid fluid bounded only by the free surface around a surface piercing body moving at a constant forward speed is considered, so the
potential flow theory can be used. With convenient assumptions, steady and unsteady
components of the velocity potential can be uncoupled and only the last one will be treated here in the frequency domain, (p: spatial component of the harmonic velocity potential defined by ç1(M, t) = ç(M)e''), w being the circular frequency of encounter. The Froudenumber based on the body length L is F = U /
jjJ
and t!ZJ= w,.jLJj is the dimensionlessfrequency. The potential p satisfies the Laplace equation, the Kelvin linearized free surface condition and a radiation condition at upstream infinity. Finally, the body condition on the body Sfi has to be satisfied:
+ f
V4rcosû-1 27zK(l -i)(eK&+ e)
vrsinûc ff/2 K3[g3(K3) + g3(K3)]- K3[g1(K4)+ gl(K4)]8+ $
l-4rcos9
(2b)where and Ö are the body translation and rotation velocities. By applying the third Green's formula to a domain bounded by the closed surface S=SB C-SL L (SL free surfaces, S control surface surrounding the body, located away from it and below z0), the potential ip at the field
point M(x,y,z) is obtained as a function of the potential and of its derivatives at the source point M'(x',y',z'). The Green function of the problem G, is defmed by Ba and Guilbaud,
1995, using the formulation of Guével and Bougis, 1982:
G(M,M',t) = 9{[G0(M,M')+ G1(M,M')+ Conjg(G2(M,M'))]e}
G(M,M')= K1[g1(K1)+g1(K1)]- K2[g(K2)+ gI(K2')]
dû (2a)
V1+4rcosû
Z-[g,(Z3)+ g,(Z3)]- Z4[g,(Z4)+ 2(Z4)]9
J4rcos8-I
R and are respectively the real and imaginary components. G0 is the Rankine term very
often studied, so no treated here and:
= z+
z'+i[(x
- x ') cosû±(y- y') sinû].
The limits
of integration
for G2 are
if t<O.25,6 =
= 0;if
O.25<t<O.5,=cos1
, 8 =0, et =e and if
0.5<t,=cos11
,0 =cos'--, a
4,
4r1
2T,We have chosen c=1OE6. These formulas involved the modified complex integral functions, gi,
g and g3, defined by:
g1() = ee1() if O<arg()<2ic; g2() = eE1() if -ir<arg()<t;
g3()
= e{e1(4) +2i7r] if t<arg()<2t, with
()=E1() if
()
O,()=E1 ()-2iir if ()
<O. The complex integral function is: E1()=
$_dt if -ic<arg()<7c and 9)> O. The poles
K, and Z of eq.(2b) are given by:
l+2rcos8+(_1)'jl+L4rcos9
,i=l to 4;L1 =L, =1;L3 =L4 =-1 K, = 2F2 cos2 O1 2rcosO +(lyJ4rcosO1
=34.
42
Z-
2F2 cos2 OSimilar equations are obtained for the first and second derivatives (Ba and Guilbaud, 1995). It is then possible to compute either q) from the Green's formula:
(M)=
(M ')G(M, M') (M')
aG(M,M')
4r aflM
or, by applying the grad operator to equation (3a), the following equation is obtained:
¿iço(M) i
/q(M')G(M,M')
o(MI)G(MM)dS
M (3b).aflM 4 aflM
In a first step, the waterline integrals and the surface ones on semi-infinite surface
extending from the trailing edge to model the lifting effects are not studied. By discretizing the body in quadrilateral panels with area S, equations (3a and b) lead to the following integral, two of them having to be computed:
a2G J1 = JJGds; J2 = S'JflM ds; J3 =
if
ds (4). S1 SIThe method developed being a first order one, the potential and its derivatives have been
assumed to be constant on each elementary surface. The study of the accuracy will be
developed for the three integrals of equation (4), so either eq. (3a) or eq. (3b) can be used in a panel method.
Methods of integration of the Green function ant its derivatives
Numerical integration for surface integral
A Gauss method of integration is used following a procedure already used (Ponizy et al., 1998) for the wave resistance Green function. The number of Gauss points chosen is a
function of the ratio of the distance between the field point M and the centre of gravity of the
panel M
to the largest diagonal (called diaL with L=1,2)of
the panel:PU = MM1 /max(dial1,
dia2).
For lo, a unique Gauss point is used; For 5<ø<lO, four Gauss points are used;
For l<p<5, the panel is divided into several triangles in function of the location of M1; on each triangle, the integration is performed by a seven points method, the total number of Gauss points ranging from 14 to 224.
Analytical integration for surface integral
For a polygonal panel with m nodes, the surface integral can be transformed into a contour one by use of the Stokes theorem, see for example Bougis, 1981:
1= Í$-j-fds=C,
f(k+l)-f(Ck)
(5).d k=1
-s-,
To compute the coefficients Ck the following properties for the modified complex exponential
integral have been used, calling () and () the first and second derivatives of the functions
with respect to the parameter :
fg1 ()d = g. ()+ in
with O<arg()<27c, forj=1,3; 5g2()d= g2()+ in2 with-it<arg()<it; $d Jg () d
= g () + (
+ i) in-
with O<arg()<27c, ifj=1 or 3 andf dJ g2()d = g,() + (+ 1) in2
-
with -ir<arg()<ir, in
= in + i arg(e) with O<arg()<2ir,in2 = 1n + iarg() with -ic<arg()<it.
The term G0, often studied in aerodynamic or Rankine hydrodynamic studies, will not be studied here. For G1 and G2, the following expressions are obtained:
1,
{Bn}S
j4rcos8i
1
J2=[Bj J
o i
for n=O,1,2 where °Q/x°=Q, L0 =L2 =1 andL3
=1, Q=Q,Q=Q,Q=Q.
fJ{ ¡x} G,ds =--{(1IL0y}{1 '2 +13+14]
i
Ck [Z Fikz; F2k +C'kr [Ln(Z kZ4 F'2k)1 4rcosO-1 k=i 13 =[(BnK)2it[C
+[L]C'k±ik
,r/2 14 LB J i[; HIkK
T2o1 + C'0 [Ln(K;' K2
1k 4h1'2k)]}d9 .Jl-4rcosOwhere:
= (qirsin0)(xk+l _xk)_(p_irsinO)(y
_yk) Mk(xk,yk,zk)
is the nodek of a panel with IvP'=M'; the outer unit normal to the panel is flext = (p,q,r). So:
[Jj g1
(d11
ix)[ff g2 ()d]
[jf g1()dY
DJk k+l kJ1,2 EJk
j3,4,
kx
f1 f
f+l f
d9 r/2 d8 Jf{"iaX"}G1dS=I{(B/L0)"} J 1, ,= Ck r -i I K, DlkK" D2k L J I+C'k I f Ln(Kin D'kK" L s, 1i + 4rcos 9F2k
[$$g3j]X+l
[JJg3()d]1'
[JJgi()]
;H1 , H2 =-
(+1Z6
-Quantities ( )' are deduced from the () ones by replacing X by .
=i!_[z+zk
i{(x_xk)cosû±(y_yk)sifle]],jl
to 6; TjKj ifjl,2 TjZj ifj3,4
1
} =
fZ+ Zk +i[(x _xk)cosO ±(y- yk)sjne]]
and TJ=K2 ifj=5,6;LoL
Concerning the first derivatives, coefficients (..) are deduced from the previous ones by
replacing the two integrations by a unique one. For the second derivatives, the coefficients (..) involved only the functions gj. Furthermore, we have the following relations: x1=x, B=icosO,
B=icos8 if i=l;x=y, B=isinB, B=iin8 if i=2 and X3=Z, B=B1 if i3.
Results of the integrations on a single panel
b Function )G1 /)x
i
e Function 2G1 !x2
Figure 1 Absolute errors and computational time differences between analytical and numerical integration methods for function G1 on a horizontal plane z=-.Ol (F=O.2; (13=1.4; 'r=O.28)
Figures 1 present the results of the difference between analytical and numerical
integrations of G1 (graph a), G1/x (graph b), and 2G1/ax2 (graph e), on a square panel with unit side parallel to axis x, submerged (upper part) of z=-O.O1, at Froude number F=O.2, non dimensional frequency (13=1.4 (reference length is L= 1), and Brard parameter t=O.28. Other derivatives give similar results and lead to the same conclusions. The left parts of the graphs
show the contours of the absolute values of the difference between the two methods of
integration and the right ones of the CPU time (negative values correspond to faster analytical integration). As the various functions are even with y, only results for yO are presented. The
graphs show a V-shaped zone where it is quite difficult to obtain a good accuracy (error is
about 0.5%) with the numerical integration; the shape of this zone is dependent of the flow
parameters as it will be shown later. Outside of this zone, very good agreement between
numerical and analytical integrations can be observed. The mean computational time of the analytical integration is about identical to the time required for a numerical integration with 4 Gauss points. For field points close to the panel, the agreement between both methods is quite good because, for the numerical integration, the number of the Gauss points increases till244
leading to a CPU time of about O.16s instead of 4.ts (time is for a DEC Alpha 255-233
workstation) for the analytical one. The time contours are almost circular close to the origin;
the O value contour has a radius of about r=lOmm and the numerical integration becomes
more and more efficient as r increases.
a Function G2
b Function G2/ax
cFunction 2G,/*x2
Figure 2 Absolute errors and computational time differences between analytical and numerical integration methods for function G2 on a horizontal plane z=-.Ol (F=O.2; 03=1.4; t=O.28)
Figures 2 show similar results for G2 and its derivatives. These computations are more difficult to perform than the previous ones, with difficulties located nearly in the same zone that for G1. The error, quite higher, is about 10% for G2 but may grow till 150% for the second derivatives, for X=Y- 10mm. Here, only a numerical integration with at least 14 Gauss points can give correct results but it requires higher CPU time than an analytical one. It is worthwhile
to notice here that computations involving G2 require less computational time than the
corresponding ones for G1 whatever is the integration method, numerical or analytical.
Figure 3 Absolute errors between analytical and numerical integration methods on a vertical plane x=y, function G1
a Function G1 b Function
G, lx
-W 43 - -10 0 lo X N c Function 2G, lax2 Doe ODO Dol 000 Doe 02) Doe Dois I oiFigure 4 Absolute errors between analytical and numerical integrations on a vertical
-Doe Doe 52) Dolo Doe Gol 02) .3 Doe -2-Doe 02) 52) Doe N DolDols N 000 N .3 -3 .4 .4
.5 tllllJl
oe - .oe - -io o i .2) .05 .2) - -2) - O 10 X Xplane x=y, function G2
Figures 3 and 4 show the results of the integration on a vertical plane oriented with 450
with respect to the x-axis (where the main difficulties have been observed previously). The
error is quite weak except close to the free surface (z and z' are here very close to O and the error may become quite high) for G1 and G2. The numerical integration is more efficient than the analytical one, the Gauss points being never too close of the free surface, contrary to the analytical integration where a part of the contour is exactly on the free surface.
Results presented on figures 1 and 2 show important numerical errors in a dihedral zone limited by angles c and ci with respect to x axis and which are functions of the
Brard parameter 'r. In this zone, it is not possible to accurately compute the various integrals involved, particularly for G2, in these computations except for high depths of immersion.
Deg
tau
Figure 5 Characteristic angles a and C(.)( versus r.
We have tried to determine these limiting angles and their dependence on the
pa-rameter
t
for values ranging from 0.27 to 1.5 by various regression, but we succeed only fora, c0 being also function of the Froude number. For a constant value of 'r, c,( increases
with the forward speed, but always
a53
1600. Figure 5 plots a and c versus 'r; symbols170 160 =1= 4
___.
4 153 :4= =1= + - -140 + 120 120A'
110 153 - L, AI p6anifl c13(1atiCI1aIx,pta6
60-(
- - - aIphninf(tu) -I 80 0.5 1.5correspond to numerical tests and curves to the limits obtained by regression. For t<0.27, the
whole domain has to be taken into account, giving a=O° and a=l8O°. On the
free-surface, z=O, Gauss integration gives better accuracy than analytical one because the Gauss points are always under the free surface but a branch of the panel contour is always on the free surface (z0), as already explained previously.We want to control the accuracy of the various integrations of the Green function and its derivatives on a single panel, but it is not necessary to have a better accuracy than the one of the function (1 Os). The error of the analytical boundary integration is also checked by the Simpson Adaptative method for the Fourier integration and are similar to the one of the Green function. From numerical tests, we have determined zones where the analytical integration is
more efficient than the numerical one (better accuracy andlor weaker computational time). Figure 6 is a schematic representation showing the various zones using the parameter R0 (equal to 10 times the source panel diagonal), the characteristic angles being given by the
following expressions:
(V-o.52)
if'r< 0.27, ct= 0°, Um= 1800;
if t>0.27, amin =160-72e
0.2,
a,= 160°.
(2) (3)
So the integration methods chosen for the integration of the terms and G2 of the Green function and of their first or second derivatives on the source panel are summarized in table 1:
CMII
CAl I
Table I Integration techniques in the various domain zones
Numerical results of seakeeping computations
Comparison of with simple Gauss method of integration
Tm
Figure 7 Added-mass and damping coefficients for the DNV barge
Zone G1 and derivatives - G2 and derivatives
O i Gauss point i Gauss point
Z=O numerical Z<O analytical
ZO numerical Z<O analytical
2 numerical Z=O numerical
Z<O analytical
The integration method defined in table 1 has been used in a sea keeping code using only a source distribution in equation (2b) without waterline integral. Figure 7 plots
added-mass (CM1 i and CM33) and damping (CAL L and CA33) coefficients for the DNV barge at F=O.34. The results obtained show only few differences between a fully numerical method or the present method (table 1), probably due to the improvement of the calculations of surface integrals. The main differences can be seen for heave for period T between 10 and 20s, but the agreement is excellent. The results are in good agreement with those obtained by Brument and Delhommeau, 1997, taking into account that our method neglects the waterline integral which must have a certain importance for a body as this barge.
14.9k CPU
7.7 h CPU pI'9''&n1 onct/iod
jul/v numerical method
1.51642 1.34693 1.17744 100795 0.838461 o 86997 o 496479 0,329988 o 160498 .000599313 .0 179404 .0347975 .0517486 .0886956 .0956447 -1 02594 119543
Figure 8 Wave amplitude around a flat plate in sway motion (top: fully numerical method, bottom present method) (AR=0.5; F0.32; =2; t=0.64)
Figure 8 shows the wave amplitudes around an aspect ratio AR=0.5 surface piercing
flat plate at F0.32 and 3=2 in forced sway motion with 128 panels. The graph upper part
corresponds to calculation with numerical Gauss method of integration and the lower one to the present technique to compute the integrals (table 1). Differences may be observed only in
surface). Finally figure 9 presents the relative error on respectively the wave pattern (fig.9a) and the pressure (fig.9b) amplitudes around the same flat plate. The errors are located in the near field, close to the free surface and to the trailing edge of the body.
Wa
L
.._
_.I
I..
_.I
I..__ii
____-fi
____
Figure 9 Influence of the integration method on the accuracy of the computations (same conditions as figure 8)
Conclusion
We have presented here a study of the accuracy of the boundary integrals dealing with the
diffraction-radiation with forward speed Green function and its first and second
derivatives. This work is a first step to develop a sea keeping method of computation using a panel method with control of the accuracy of the surface integration. An integration technique mixing both analytical and numerical Gauss methods (with variable Gauss point numbers) has been defined. The analytical method uses a Stokes theorem. The whole computational domain has been divided into four zones where the method of integration has been chosen taken into account both the accuracy and the computational time. The first zone is upstream of the body and limited downstream by the angle cu1 with x-axis out of a vertical cylinder zone which isthe second one (with radius 10 times the panel diagonal and is of infinite extent in the z
direction) around the centre of gravity of the source panel; the third one is downstream, in a a) wave pattern b) pressure distribution
dihedral zone limited by angles a, and OEmax and the last one downstream of this last angles. The values of these angles have been defined as a function of the Brard parameter t. Finally, close to the free-surface, a numerical integration method has been chosen because no Green function computations are required on the free surface, contrary to the analytical method for
which a branch of the contour is on this surface. This integration technique has been
introduced into a velocity based panel method using only source distribution (no lifting effects) and neglecting the waterline integral. Results for DNV barge in surge and pitchmotions or a surface piercing flat plate in sway motion have shown good agreement with other
results showing very few discrepancies with respect to a fully numerical Gauss integration
method for global coefficients, but a better representation of the near wave pattern and of the pressure distribution close to the free-surface and to the body trailing edge. Furthermore, the
computational time to calculate the wave pattern has been divided by a factor of about 2.
Work is on progress to compute flows around other ship models as Wigley or Series-60 hulls, to calculate accurately the waterline integral and to introduce lifting effects in the code.
References
Ba, M, and Guilbaud, M (1995). "A fast method of evaluation for the translating and pulsating Green's function', Ship Technology Research, Vol. 42, pp. 68-80.
Bougis, J (1981). "Etude de la diffraction-radiation dans le cas d'un flotteur indéformable
animé d'une vitesse moyenne constante et sollicité par une houle sinusoïdale de faible
amplitude", Thèse de doctorat, Université de Nantes.
Brument, A, and Dethommeau, G (1997). "Evaluation numérique de la fonction de Green de la tenue à la mer avec vitesse d'avance", Proc. of the6èm Jour, de l'Hydrodynamique, Nantes (France), pp. 147-160.
Brument, A (1998). "Evaluation numérique de la fonction de Green de la tenue à la mer",
Thèse de doctorat, Université de Nantes.
Chen, XB, and Noblesse, F (1998). "Super Green functions", 22nd Symposium on Naval Hydrodynamics, Washington (USA), Thursday/Friday session, pp. 139-153.
Guével, P, and Bougis, J (1982). "Ship motions with forward speed infinite depth", mt. Ship. Progress, Vol. 29, pp. 103-117.
Iwashita, H, and Ohkusu, M (1989). "Hydrodynamic Forces on a Ship Moving at Forward
Speed in Waves",J.S.N.A. Japan, Vol. 166, pp. 87-109.
Ponizy, B, Guilbaud, M and Ba, M (1998). "Numerical computations and integrations of the wave resistance Green's function", Theoretical and Computational Fluid Dynamics, Vol.12, N°3, pp. 179-194.
Sciavounos, PD (1996). "Computations of wave ship interaction", Advances in Marine
