Higher order panel methods for potential flows with linear of non-linear free surface boundary conditions

TECHNISCHE UNIVERSITEIT

Laboratorium voor Scheepshydromechanica Meketweg 2 - 2628 CD DELFT"

CHALMERS

UNIVERSITY

OF TECHNOLOGY

GOTEBORG

SWEDEN

HIGHER ORDER PANEL METHODS FOR POTENTIAL

FLOWS WITH LINEAR OR NON-LINEAR FREE

SURFACE BOUNDARY CONDITIONS

Shao-Yu Ni

ARCHEF

4

Division of Marine Hydrodynamics

Goteborg, Sweden

1987

el1

CHALMERS TEKNISISA HOGSKOLA

fri

HIGHER ORDER PANEL METHODS FOR POTENTIAL FLOVE WITH LINEAR OR NON-LINEAR FREE

SURFACE BOUNDARY CONDITIONS

Shao -Yu Ni

...

Division of Marine Hydrodynamics

Goteborg, Sweden. 1987

Chalmers University of Technology

Goteborg, Sweden, 1987 ISBN 91-7032-323-2

CONTENTS

Abstract

Keywords

Introduction

The attitude and righting moment of heeled ship in a stationary wave

The higher order panel method 4

The potential flow problem with the free surface 7

Solution of the linear free-surface problem 8

Solution of the nonlinear free-surface problem 10

Reference Acknowledgements Appended Papers Paper A Paper B Paper C Page 1 2 1,

ABSTRACT

In this work potential flows around ships with linear and nonli-near free surface conditions are investigated. The thesis con-sists of three parts as follows:

A computer program package has been developed, which can

gene-rate panels automatically for the calculation of the stability of a ship in a stationary wave and for the hull definition in

the nonlinear free-surface potential flow problem. The ship attitude and the righting moment are determined by the

integra-tion of the hydrostatic pressure including the Smith effect. Results are given for two cases, and in one of them comparisons are made with measurements.

A higher order panel method is investigated here to represent

the boundary surface covered by sources in accurate way. In this method the panel is supposed to be a parabolic quadrilate-ral with a linearly varying source density. The higher order panel method is applied to some cases which the flow is unboun-ded and to the linear free-surface problem. The free surface boundary conditions can be linearized along streamlines or arbitrary smooth body-fitted curves. Improved results as pared with the first order methods are shown in the test com-putations.

In the nonlinear free-surface method developed, the free surfa-ce boundary conditions are exactly satisfied on the wavy sur-face instead of on the undisturbed free sursur-face as in the linear free surface problem. An iterative procedure is applied to solve the source distribution and wave elevation

modifica-tion. To obtain convergence a higher order global algorithm

with a relaxation factor is essential. A single model is propo-sed, where the hull panels are wavy to fit the free surface. They can also change in each iteration. The model has been

extended to the nonlinear free-surface problem in shallow

water also. Convergent results have been obtained for all the test computations.

This thesis is based on the work contained in the following papers:

Ni, S Y: "A method for Calculating the Attitude and Righting Moment of a Heeled Ship in a Stationary Wave"

SSPA Report 2912-4, SSPA Maritime Research and Consulting, Goteborg, Sweden, 1986

Ni, S Y: Higher Order Panel Method for Double Model linea-rized Free Surface Potential Flows"

SSPA Report 2912-5, SSPA Maritime Research and Consulting, Goteborg, Sweden, 1987

Ni, S Y: "A method for Calculating Non-Linear Free Surface Potential Flows Using Higher order Panels"

SSPA Report 2912-6, SSPA Maritime Research and Consulting, Gbiteborg, Sweden, 1987

Keywords:

righting moment, generate panels, potential flow, higher order panel method, Rankine source, free surface, boundary condition,

1. Introduction

1

Numerical hydrodynamics has become a neccessary tool to guide model experiments and ship design in the field of naval architec-ture. The calculation of the potential flow is a very important

part of numerical hydrodynamics.

To perform numerical calculations of the non-lifting potential

flow about arbitrary three-dimensional bodies, the panel method was developed by Hess & Smith in 1962, [1]. This method turned

out to be very successful for flows without a free surface. In 1976 Gadd first introduced simple Rankine Sources on panels of

the ship hull and a part of the undisturbed free sface for determining a flow field, where the body and free surface

bounda-ry conditions are satisfied, [2]. Soon thereafter in 1977, Daw-son presented the free surface boundary condition linearized

along the streamlines of the double-model solution, [3]. In

recent years, Dawson's approach has become very popular.

In Gadd's and Dawson's methods the free surface boundary

condi-tion is not really satisfied on the wavy surface, but later many

investigations have been carried out to solve the nonlinear free-surface problem. Cgiwara took the nonlinear effect into account iteratively by using relaxation factors in 1985, [14]

and Xia proposed a iterative procedure in 1986, [5]. In Xia's method the free surface boundary condition is linearized about the initial wavy surface, using the small pertubation principle and new wave elevations and source distributions are solved for in the next iteration. It was found that the convergence problem

for the iterative procedure is quite severe.

Conceivably the assumption that the panel is flat and covered by

constant sources is not satisfactory for the panel method itself. Hess presented a higher order panel method for two dimensional flows in 1973, [6]. In this method the two-dimensional panel is taken as a parabolic element. with a linear variation of the

sour-ce density. Then, Hess in 1977 [7] investigated the nonlinear two-dimensional free surface problem and concluded that a higher

order global algorithm with an initial flat free surface is essential for the iterative convergence. In 1979 Hess extended

the higher order panel method to the three-dimensional case,

[8].

In the present study a computer program has been developed first.

for the higher order panel method in three-dimensional flow. Then the higher order panel method is applied to the linear

free-surface problem. The free-free-surface condition can be linearized

along streamlines or along arbitrary smooth body-fitted curves.

The momentum approach to calculate wave resitance is investigated also. This work is reported in part B. In part C the higher order global algorithm is applied to the nonlinear three-dimiensional

free-surface problem. The results of the test computations are

always convergent if a relaxation factor is used for the wave elevation modification, a single model with a new paneling to fit the wavy free-surface is used and the vertical derivatives of the

induced velocity are kept in the free-surface boundary condition.

A program package to generate panels on the hull surface automa-tically has also been developed not only for the nonlinear free-surface problem, but also for the calculation of the attitude and righting moment of a ship in a stationary regular wave. This work

is reported in part A.

Future investigations may consider the following:

to search an accurate method to calculate the wave resistance to extend the free surface method to the calculation of the

righting moment.

2. The attitude and righting moment. of heeled ship in a stationary wave (Part A)

2

The discretization of boundary surfaces as panels is frequently

needed for numerical calculations in the field of hydromechanics.

A computer program package has been developed, which can genene-rate panels automatically for the calculation of the stability of

a ship in a stationary wave and for the hull definition in the

nonlinear free-surface potential flow problem.

For a certain wavy free surface the submergence of the stations of the ship hull can be obtained numerically. Then the corner points of the panels are obtained by dividing each station in equal arc lengths from the keel up to the intersection with the

free-surface. Of course, some new added stations and corner

points on the stations are allowed for at some locations of the hull surface

SB where the curvature is large. All. the geometric

quantities of theNB panels can be computed from the corner

points.

For a certain sinkage zo , trim angle e and heel angel (1)

the displacement. is given by

D(zo, 0, (P) =

no

pj ASi

(2-1)

The trim moment and the righting moment are obtained as follows

NB

TE(z,0'60) =

-I n

p. AS. E

+ DEG j 3 3 _rilDi PM = NB nrj pj ASj nnpi - DnG

j

Pj = Pg Cnpj (2-2) (2-3) 3

where the pane] area A Si, the

null

point

npjnpjnpj

)

and the normal vector F (ncj, nnj, nci ) for the j-the panel are

known. The gravity center of the ship is ocated at _{(C G,} CG ) If the free surface is calm the pressure is simply

and if the free surface, is a stationary regular wave with wave

length A , the pressure becomes

'Crijp pg [c. t,

27- e27

)11 3 nP3 where C- is attitude, i. ment D and D (z D =,13 o ' TM (zoi0f(1), . 0

and then. the 'righting moment is determined from Eq i(2L3).

However, the heeled ship

in

the calm water or in the stationary wave will radiate waves. As the first approximation, the Dawson's

linear free-surface could be taken into account. as a dynamic

cor-rection to the righting moment. Furthermore, the nonlinear free-surface solution could be introduced on the stationary wave In

futu-re work.

The paper A gives some computational results for the equations,

mentioned above.

3. The higher order panel method (Part B)

In the first order panel method each element is assumed to be flat and covered by a constant source density. The constant

sour-ce C. on the j-th flat panel A. will indusour-ce the potential ]

ct, at the i-th contrOl point P(x,y,z,1 in the panel

codr-dinate system 1-1c, ,

the wave elevation above the j,-th panel. The ship

e.

its sinkage zo and trim e to a certain displace-heel angle cP1 , is computed from the equations

4

= (e

(1) =

A.2rf

ds = c. °)

ij 3

(3-1)

where the unit potential is defined as

(0

I 1 ds

=

Aj r

( 3-2 )

and the distance from the control point P to any point in the j-th panel is

rf = ibt-E,")

2

+(Y-n)2+z2

However, the panel is supposedly parabolic in the higher order panel method

= zo + AZ + Bn +

K2

+ 2(2n +

Rn2

= D (cfrl) E(c,n) P + F(,n) R

where D,E and F are determined by the four input points of the panel. The arbitrary values of P and R are adjusted so as to make the distances from the parabolic surface to the neighbouring input points of the adjacent panels as smallas possible in the least squares sence.

The source distribution G(

_,n)

_{is considered to be linear on}

the j-th panel

GV,,n) = G. +GC+GT1

_x

5

and is fitted in the least. squares sense to the values of the source densities at the control points of the adjacent panels. The induced potential may then be written, omitting higher order terms,of and n = ds A. A (10) (Q) (El

=o.(1)+0.(1p,f' -

+2W 144, )+0x{(1x) (1y) , 2 , 2 2 where

r =

J

kx-t=1 ty-1-1) (z-C)

,(o)

The potential Y which refers to the first order panel, is defined in Eq (3-2).

The potentials, caused by the curved shape of the panel, are obtained from 7,2 cb(p) =

f===-

dCdn 3 A. rf 3 (1)(Q)

f

3 clZdh A. rf (10 r zn2 =

3ddq

A. rf 6 3

The potentials, caused by the linear variation of the source

density, are given by

(lx) ci) IA A. r f (1y) n A. r f

The induced velocities are gradients of the potentials above. The

solutions of the flow field are more accurate by the higher order

panel method than the first order one, as shown in Paper B.

4. The potential flow problem with the free surface (Part B continued)

The reference coordinate system with axes of x and y in the undisturbed free surface and z axis in the vertically upward direction is used on a ship. It is assumed that the ship is fixed in a uniform onset flow of velocity U. in the direction of the x

axis. The flow field maybedefined by a velocity potential since

the flow is supposed incompressible, inviscid and irrotational. The potential is govered by the Laplace equation.

(4-1)

At infinity the flow approaches the freestream -V =Uoo ddri dc1.71 (4-2) 7

The boundary condition on the wetted hull surface SB is

(4-3)

where the subscript n means a partial derivative along the normal

direction to that surface.

The boundary condition on the free surface SF consists of two

parts:

the pressure on SF is constant

1

C+ (42)X+42)1711)2

2

Z o

and the normal velocity to SF is zero

(4-4)

The free surface boundary conditions (4-4) and (4-5) are

non-linear not only because the equations are nonlinear, but also

because the location ç of the free surface is not known. In the 8

ci),Jx-qyCy-(Pz = 0

(4-5)

where is the wave elevation of the surface SF.

To obtain the numerical solutions of Eqs (4-3), (4-4) and (4-5),

the surface singularity method is employed because the solutions

will satisfy Eqs (4-1) and (4-2) automatically.

5. Solution of the linear free-surface problem (Part B continued)

linear free surface problem the boundary conditions are only app-lied on the plane z=0 and the linearization is made on this pla-ne in terms of the double model solutions (1).

If the linearization is carried out along a streamline of the double model flow, the boundary condition is simplified as

29,(1))

(1) = 2402 (1) Z

(5-1)

where denotes the direction along the streamline.

If the linearization is performed along an arbitrary body-fitted cure the boundary condition can be written as

(1) Co i-(1) c° -(1) - {(I) (2 4)2 2 [(I) +(/) (1) +4)) 2 y 2 (1) [4) cl) +4) cp -(4) +cD

)1

.1=0

xxyy

x

yjy

(5-2) o .

where C is the elevation of Bernoullis wave.

To solve the problem, the hull surface SE and a part of the plane z=0 are divided into NE and NF source panels respecti-vely. The source density can be determined by solving the linear equations (4-3) and (5-1) or (5-2). Then the flow field is

obtai-ned from the known source distribution. The wave resistance is calculated by integrating the pressure on SB as usual, or by employing a momentum approach

NF 47 I U .a.AS. . B3 Cw- .7 NB U- E AS. (5-3)

where AS, is the j-th panel area and is the velocity in

x -direction at the j-th control point on the plane z=0 induced by all the panels on the hull surface SE.

The test computations prove that some improvement can be obtained by using the higher order panel method and the momentum approach as indicated in paper B.

6. Solution of the nonlinear free surface problem (Part C)

In the nonlinear free-surface problem the free surface boundary conditions should be exactly satisfied on the wave surface SF although the elevation C of SF is unknown a priori. An itera-tive procedure is applied here because the surface singularity method requires that the locations of all the boundaries be

known. Starting with an initial estimate of the wave elevation C°

and the velocity V4) on c° the boundary conditions (4-4) and

(4-5) on the unknoc,m, surface c=c°+fc can be expanded as linear

equations with unknown wave elevation modification fc and unknown velocity VO on

c°.

10 c +0 c0 -0 +0 (% +0 dc +(0 c0 +0 c -0 )6c=0

xxyyzxxyy xzxyzy zz

1 1 fc= f [0 0 +0 0 +0 0 _

g xxyy zz

1 1+-(4) 10 +4) 4) +4 4) ) g x xz y yz z zz 2 2 2 2 o, 1 (4) +4) +4) +U ) 4-C x y z (,) (6-1) (6-2) 3 w

and the boundary condition on the hull surface SB is given by

Eq (4-3).

If the solutions of the wave elevation ej+co6c and the

velocity Vq5+07(1)zcSc do not satisfy Eqs (4-4) and (4-5), the solutions may be used as the new initial estimates for next ite-ration. The under relaxation factor w is set simply to 0.5 to

ensure that the pertubat ion is small. Eqs (6-1) and (6-2) are based on the small pertubation assumption.

co and q) --(1) _{on convergent iterations, so finally the}

linearized equations (6-1) and (6-2) becomes the exact free

sur-face boundary conditions (4-4) and (4-5).

The higher order panel method is essential to improve the itera-tive convergence because the higher order panels (covered by a

varying source density) represent the wavy free surface and the

hull surface better than the first order panels.

The hull surface SB can be the surface of a double model such as in the linear wave problem, but it is more reasonable to _{use a}

singe] model, where the hull panels _{are wavy to fit the free}

sur-face, of course, they then have to be changed in each _iteration.

The automatic procedure for generating panels _{on the hull is} nee-ded for this kind of _{model. The singel model can also be}

em-ployed to investigate the nonlinear free-surface _{problem in a} restricted depth of water when the water bottom is assumed to be

a symmetry plane.

Convergent results are obtained for all the test computations

presented in paper C. Better methods to calculate wave resistance should be investigated in further work, since the integration of

pressure is not always satisfactory.

Reference

Hess, J L & Smith, A MO: "Calculation of Non-lifting

Poten-tial flow about Arbitrary Three-Dimensional Bodies"

Douglas Report NO E S 40622, 1962

Gadd, G E: "A MEthod of Computing the flow and Surface Wave

Pattern around Full Dorms"

Transactions of the Royal Intstitution of Naval Architecture,

1961

Dawson, C W: "A practical Camputor method for Solving

Ship-Wave problems

Proceedings of the 2nd International Conference on Numerical Ship Hydrodynamics, 1977

Ogiwara, S & Maruo, H: "A Numerical Method of Non-linear

Solution for steady Waves around Ships"

Journal of the Society of Naval Architects of Japan, Vol 157,

1985

Xia, F: "Numerical Calculation of Ship Flows, with Special

Emphasis on the Free Surface Potential Flow"

PHD thesis, Chalmers University of Technology, 1986

Hess, J L: "Higher order numerical Solution of the Integral Equation for the Two-Dimensional Neumann Problem".

Computer Methods in Applied Mechanics and Engineering 2, 1973

Hess, J L: "Proyless in the Calculation of Nonlinear Free Sur-face problems by SurSur-face Singularity Techniques"

Proceedings of the 2nd International Conference on Numerical Ship Hydrodynamics, 1977

Hess, J L: "A Higher Order Panel Method for Three-Dimensional

Potential Flow" Douglas Report. N 62269-77-C-0437, 1979 12 3., 4. 8.

Acknowledgements

I wish to express my genuine respect and sincere thanks to

pro-fessor Tars Larsson for his erudite knowledge and effective super-vision.

Grateful thanks are also given to Dr C Kallstrom, Dr N Norrbin (SSPA), professor C Falkemo, professor G Dyne, professor J Lunde (Marine Hydrodynamics), professor B Qvarnstrom (Control Enginee-ring), Dr S Holm (Mathematics) and Dr N Fjdllbrant (Main Library)

for their excellent courses and helpful discussions.

I would like to thank Miss E Samuelsson, Mr G Forssberg, Mr L Broberg, Mrs B Karsberg (SSPA), Mrs S Bernander, Mir C-0 Tarsson

(Marine Hydrodynamics) as well as other colleagues and friends at SSPA and CTH for their valued help in my studies.

The financial support from the SwedishBoard for Technical Development and the Defence Material Administration of Sweden is gratefully acknowledged.

A

SSPA Report 2912-4

A METHOD FOR CALCULATING THE ATTITUDE AND RIGHTING MOMENT OF A HEELED SHIP IN A STATIONARY WAVE

by

Shao-Yu Ni

1. Introduction

At SSPA and Chalmers University of Technology work _has been under way for some years to develop numerical methods

for potential flow calculations. _{Recent results are}

pre-sented in [1] and [2]. The methods developed _{are all of} the panel type, i e the hull and part of the free water

surface are approximated by a set of quadrilateral panels.

The generation of these panels takes _{considerable effort,}

and, particularly in _{cases when the hull is free to sink}

and trim during the calculations, _{it would be very useful}

to have an automatic procedure for generating _{the panels} on the hull. Also, such a procedure might be _{required in}

the further development of Xia's _{non-linear method, where}

the free surface panels are adjusted to the wavy surface.

To fit this surface the hull panels _{also have to be wavy,}

and, furthermore, they have to _{change in each iteration.}

The development of an automatic _{panelization method was}

the primary objective when starting this project. However,

as pointed out by Soderberg [3], _{once this method has been}

developed it could be used for _{an entirely different} pur-pose, namely the calculation of the stability _{of a ship in} a following wave, stationary with respect _{to the ship. The} basic idea is just to add the _{contributions to the}

right-ing moment by the pressure forces _{on all panels. In calm}

water and zero speed the calculation _{of this (hydrostatic)}

force is trivial, while in general the problem is quite

complex. As a first approximation, _{however, the}

hydro-static pressure might be corrected by the so called Smith

effect, taking into account the dynamics _{of the}

undis-turbed wave. A method based _{on this idea is presented in}

this report together with a description of the

paneliza-tion procedure.

It should be pointed out, that the following wave case

is very important to the ship broaching problem.

2. Mathematical Statement of the Problem

2.1 Coordinate systems

The 0-xyz system fixed on the hull and the

z0-/-1

system fixed on the still water surface are shown in Fig 1.

Fig 1

The hull surface is defined in the 0-xyz system as

f.(y,z) = 0

for each station x., defined by discrete measured points.

The still water surface is given in the z0-7-1C system as

= 0

and the wave surface can be described in z as

0 X = CA sin (-27 + 2

where (i is the wave initial phase referred to the origin 3

The attitude of the ship in water is specified by _three parameters: position of origin

zo,

pitch angle e and heel

angle (4).

2.2 Submergence of each _station

The expression of the hull surface fi(y,z) = 0 for each

station x. in the

Z0 system is

g. (E,1-1,) = 0

and the transformation between the two systems is

= A z-Z 0

[case

sinesimp

sinecost.2

w

where A = 0 co -sing)

-sine cosesimp cosecosq)

The intersection point _{of the i-th station}

1 1 1 1

with a still water surface is obtained by solving the

following equations

IC = 0

Igi(

n,c)

= 0

(1)

Ixi

-1

y. =A

ni

T.-z -Ci

1 0

To obtain the intersection point of each station with the wave surface a further transformation should be performed.

= B

rn]

and the expression for the hull surface at station i may be written

h.(C',1-11,c) = 0

1

The intersection point

Pi(q,ni,Ci)

of the i-th station

in the

z0 system is the solution of the equations

X

lc

= CA sin (1 +13)

h.(',T1',c) = 0 1 Fig 3

where A-1 =AT =

where B = [ cosy sinyi

cose 0 -sine

sinesimp COSW cosesinw

sinecoscp -sinw cosecosw

-siny cosy

i-th station

4

This intersection point in the original system 0-xyz, P.(x.,y.,T.), can be obtained by inverse transformation

(3)

and

The girthline between the keel and the intersection _point P. may now be divided in some regular way to generate the₁

quadrilateral panels, see Fig 4.

Fig 4

2.3 Displacement, trim moment and righting _{moment in} still water

To calculate the displacement, trim moment and righting moment it is necessary to transform the _{centre of gravity}

G and all the input points into the _{zo-CnE, system again}

by using equation (1). = _{A-1 11.} 1 1 keel 5

The solution may be transformed to the 0-xyz system as

follows

cosy -siny

= B-1

rl

where B-1 = BT = I

i

siny cosy

For the j-th panel the area S., the null point

npj

C .) and the normal vector

npj npj

= n .T + n + n

I n3 CJ

can be obtained. (See Hess & Smith, [4].)

The displacement for the present attitude of the ship

zo,

0 and 4) in still water is: (Summations in the

fol-lowing are taken over all panels. Pressures and moments

are divided by Pg.)

V(z0,0,(4)) =

The trim moment about the y-axis and righting moment

about the x-axis are as follows:

TM(z0,6,4)) = -E n .c .S. . +

VG

cj npj j npj

RM(z0,(3,(4)) = n .S.

11

-

V

G

cj npj j npj

For a fixed value of displacement Vo only one set of

z''

60, w should be chosen to satisfy

0 V(z1 ' ,61 (p)

_{- Vo}

= 0 1 0 1TM (q,01,(.0 = 0 2.4 Wave effects

The dynamic effect on the pressure of the orbital motions in the waves may be taken into account by the so called Smith effect, see e g [5]. Assume that the static pressure

is

P. =ç.-;

. 17 1 nPj 6 (4) j

where C, is the wave height .above the j-th panel and the

total pressure is

P P,. + AP.

where AP.

is

the correction of the Smith effect, which.

can be computed from (et the Appendix)

A X _

C -7- C

AP. =27

nPj

- e2.7

Accordingly, for the wave case the displacement, trim moment and righting moment become

1V(zo,,e,(0 'I= E nciPljSj 'TM(zotel,,(p) =-:En ..P 1 C3 13

3 3

1-, \/'G R1vI(Z0Je,W) = E nciyliSjni Vn-G

if the Smith effect as omitted and

1 ,V4z0,0,(p) = ,E ntiP2i,Sili .1,.TM(z0"; 0 4) =-_-.E n .P23.S-E. 3 3 3 1

.RM(zo,,81,0

,=-. E'' n AD .S.T1, C3 23 "3 J

+ VG

G

if the Smith effect is taken _{into account.} 2.5 Analysis of the mathematical statement

The total velocity potential around a ship in, Waves is

usually resolved into three parts: incident potential, scattering potential and radiant potential.

gbr

7

=

-3.1 Intersection points of the water surface and some

curves

The curves can be body plan, profile and corner line, which are known as discrete points.

a) Body plan

The i-th station in the z0 or the z0 system is

a space curve. For the still water case this curve is in equation (2) regarded as

f (C )

(5)

Ti = f2 (C )

and the solution of equation (2) is obtained by

inter-polating equations (5) at = 0, using a cubic spline

subroutine, INT. For the wave case this space curve is

in equation (3)

8

For convenience only following sea is considered. It is assumed for the mathematical statement that the ship speed is equal to the wave speed, so that the wave profile is

fixed to the ship hull, but relative motion between the

ship and the water still exists. Therefore cl)s = 0 and

cp. has to be considered as a Smith effect correction. Some

error may be caused if cihr is neglected. As the first

approximation (Pr may be the solution of the Dawson method

[6] for an asymmetric case. This is an interesting exten-sion of the present work, where only the effect of ¢j is

considered.

= _(C)

= f2 (C)

In the coordinate plane zo--E,'C the intersection _{point (see} Fig 3) _{is obtained by calling subroutine POINT, which is}

designed to resolve the intersection of two plane curves

= fl

cAsin (-27T E

Then n' _{= f2(C) is interpolated.}

b) Corner line

Some high speed vessels have a corner line on the hull

surface with or without step, which should _{be taken as}

a boundary of sections.

Fig 5

The intersection x3 can be obtained in a similar way as

the treatment of the body plan. _{The maximum number of}

intersection points is three in the _{program, so the wave}

length cannot be too short.

9

(6)

c) Profile

It is also easy to obtain the draughts at the bow post and the stern post. In accordance with xl and x2 the length of the submerged hull and the numbers of stations and sections

under water are determined.

Fig 6

The stations in front of

x1 and after x2 are discarded,

but should be included in some special cases, as Fig 7,

when the ship is heeling heavily.

Fig 7

d) Deck

When the water surface is higher than the deck, the added wet surface on the deck will be considered automatically.

x2

The intersection of the straight line AB and the water surface is obtained by the same approach as mentioned

above,

Fig 8

i-th

station,

3.2 Panel formation

Each station is divided

in

equal arc lengthS from the keel up to the intersection with the surface. Since the Subroutine INT

can

only be applied to single valued

functions the following different shapes of stations

should be treated individually.

E

sine wave

11

y f (z) I,: one panel _I: z = cp(y)

v(y)

II: y = f(z) II: y = f(z) (a) _{(b) (c)} (d) Fig 9 = _z =

Fig 10

In each section new stations and new rows of panels can

be added for the location where the surface curvature is large.

Added new stations

5 1,2, 3,4 2 4 3 12

For case (c) the subroutine INT is called twice, and the number of measuring points should be more than two

for each piece. For case (b) the flat keel must be a single panel without calling subroutine INT, so two measuring points are enough for part I.

In each section the number of panels in one column is constant and to avoid too thin panels in some cases the special features of Fig 10 are used. Triangular elements can thus be generated and the width of a panel may be

reduced to zero.

Section I Section II

Fig 11

Added new

13

Finally it should be mentioned that the keel of a sailing

boat and other appendages are taken _{as a single section.}

3.3 Iteration to obtain the ship attitude

Equation (4) _{is non-linear and implicit. To solve it the}

subroutine FM (zo,

e,

(I), V, TM, RN, IWV), which generates

panels and calculates the displacement V, trim moment TM and righting moment RN according to ship attitude zo, e

and (4) _{are called according to the following procedure. The}

parameter IWV = 0 means the still water case, otherwise

is known only Call

.FM(Z.,0,0,V.,TM.,RM.,0)

1 1 1 1

(i> 3)

Interpolating Z0lV=V0 from Zi=f(Vi)

Call FM(Z ,0,0 V ,TM

RM 0) 0

'

1

l'

1,

= 0 and TM1=RM1=0? f

(

V3=V0 and TM3=0?

(V4=V0

and TM4=0? i=1 Yes No is known onl Yes Yes_ Stop Yes RM3 as well as Z0 and 6

= 0

RM4 as well as Z4 and 8 = 0

RM. as

well as Z., 1

and Z are known

Call

FM(Z ,0,0,V

,

Call

FM(Z.,0,4),V.,TM.,RM.,1)

where

Z.(i<3)

around

Zo

1 1

Interpolating Z IV=V from

Z.=f(Vi)

Call

FM(Z.,6..,(4),V..,TM..,RM..,1)

13 13 13 13

where

8..=f(TMij) is

_Ij

obtained for each Z. > 3)

Interpolating

eilTmi=o

from

eii=f(Tmii)

Call

FM(Z.,0.,(1),Vi,TMi,RMi,1)

'TM2=RM2=0? No Stop --Stop (OK) --Stop (OK) --Stop (OK) 14 !Call

FM(Z ,0,(4),V3,TM

)

Ki>3?

No i=i+1 Yes Call FM(Z 10,41)1V f , 1 1. 2.

,RM ,1)

\ I. V ,TM ,0)1

For sequence

ze.,v.(Tm.=o),

i>3

a. 1 1 1

interpolating e*Iv=v0 from Oi=f(Vi) and interpolating Z*!V=V0 from Zi=f(V.)

Cali FMAZ*, * (;),V*,TM*,RM*,1)

(V*=Vo and No Stop (Failure)

TM*=0? Yes

RM* as well as Z* and e*

Stop (OK)

16

4. Calculated Results

The method has been applied to two ships a high-speed

vessel and.

a

_{Ro-Ro ship. Unfortunately measured data}

are very rare for the wave case and the only ones avail= able to the author are the GI curves for the high speed

ship.

4.1 The. high speed ship

For proprietary reasons the particulars of this ship Cannot be released. It is known that the metacenter is

1.30 M.

In Fig 12 results of measurements and calculations are shown for the following sea case- Measurements _were obtained at two speeds by keeping the model restrained

It roll, yaw, surge and, sway. It is. seen that. there is

a considerable difference between the measured hogging and sagging cases, but

in

the calculations the differ

ence is smaller. Considering the variations with speed

for the, measured sagging Case and extrapolating to zero

(as in the calculations) the numerical prediction seems

good. The correspondence between calculations and

measurements does, howevef, teem considerably worse for the hogging case. As a matter of fact the measured

righting moment for this case seems surprisingly low, Considering the fact that the actual GM yields

_a

still,

water curve close to, the measured sagging case. This

curve is quite well predicted by the program and IS about midway between the two wave cases, as expected.

4.2 The Ro-Ro ship

Table 1. Particulars of the Ro-Ro ship

There are no measurements available for this ship; but calculated results are given for two waves; at two

directions in Figs 13-20. The difference between. the

sagging and hogging cases is, substantial in some instances. 17 Length Breadth Draft Displacement Metacenter Block coefficient Slenderness coefficient V GM CB, L/V13 iM) (m) .613 (rn) 29 180.0 27. 9.102 190 0.81 0.6525 5.992 L B T (m)

5. References

18

Xia, F: "Numerical Calculations of Ship Flows, with Special Emphasis on the Free Surface Potential Flow". PhD thesis of Chalmers University of Technology, 1986

Xia, F & Larsson, L: "A Calculation Method for the Lifting Potential Flow around Yawed Surface-Piercing

3-D Bodies". 16th Symposium on Naval Hydrodynamics,

1986

Soderberg, PI Private communication

Hess, J & Smith, A M 0: "Calculation of Non-Lifting Potential Flow about Arbitrary Three Dimensional Bodies". Douglas Report No ES 40622, 1962

Bhattacharyya, R: "Dynamics of Marine Vehicles". A Wiley-Interscience publication, USA, 1978

Dawson, C W: "A Practical Computer Method for Solving Ship-Wave Problems". Proceedings of 2nd International Conference on Numerical Ship Hydrodynamics, 1977

[2]

Appendix

Derivation of the Smith correction.

Bernoulli's equation for unsteady flow may be written as

P. + + zpU1 2 ₌

-pn-

60

Assume small waves,i.e small orbital velocity U2=-'0, and

wcA Kz

(1) = K e cos(Kx = wt)

A

where K = and w2 = Kg,

Therefore ,we have

P -pgz +

pgA

el<z sin(Kx - wt)

Kz

= -2,g2

4 pg.

e

Especially,the pressure On the free, surface is

Kc.

P0= pgC_. e- 3

Eq( 1 ) minus ( 2 ) yields

P - P pg(cj - pgc.(eKz - e

K.

31)

19

( 2 )

( 3 )

When the field point is located at _{the i.-th null point,} the Smith correction is obtained by _{inserting z=}

into Eq( 3 ). pgz = ( 1 - - z) + -npj

SSPA

Righting Moment Arm

OZ [ml

10

0.8 0.6 0.4

02

0

GZ Curve for a High Speed Ship

Regular Wave. Heading 0 Deg. X/L=1.25 X/H=22 FIG. 12 GM=1.30 M Calculation Calm Water Sagging Hogging Measurement 15 Knots 25 Knots Heel Angle 1 A II lb L II Sagging a a

A/

illHogging I

W

rit

irlig

A

A

/01

AllAl

A

/

,i,

,

/

IA Lfl SO Lk.) 0 - -I

-SSPA

GZ (METER)2 0

1.5

1.00

0.50

GZ CURVE

SSPA MODEL 2062-A REGULAR WAVE HEADING 0 DEG

H = 8 METER T = 10 SEC HEEL ANGLE

mm

FIG* 13 0

0

0 CALM WATER A-

-

-Li SAGGING

-I- - -+ HOGGING

A

.

.

.

.

.

-.

.

.

.

A

.

.

/

/

-3(

/

/

/'I

.4

/

.

A/ ..

/

/

...---.-+

.

-V

V

'0

10 00

20 00

30.00

40.1 -A

fz,

2

0.0 10.0 SINKAGE

00

O.

1.

P 0.0 TRIM ANGLE

mm

Trim and Sinkage Curves

SSPA MODEL 2062-A REGULAR WAVE HEADING 0 DEG

H = 8 METER T = 10 SEC'

10.0

20.0

14

20.0

HEEL ANGLE (DEG)

CALM WATER fr-- A SAGGING + HOGGING 30.0 F 6: 40.0 --m .m-4.... ... --.. . -,.-.-.. -._

,_

..._e_____. III I: 1 , 1 i[ Ai-

- -

-

-.. I. ... .... , 1. -_

SSPA

0

- - - A

30.0

40.0

2

.0

GZ (METER))

1 5

1.0

0.50

00

0.

61 CURVE

SSPA MODEL 2062-A

REGULAR NAVE HEADING 0 DEG

H=5 METER T 8 SEC

HEEL ANGLE (DEG)

0 , 1 , 1 C) CALM WATER SAGGING HOGGING,

ta-

C)--- A

+

1--

-, , 1 . . . 7 I rr , 1 1 I 1 11 1 1 , 1

.

..

/

A ./ t.'

I

/

. ( 1 i 1

,

i

/

/

*

.

1 , , A' , , I , , ,

10.00

20.00

30.00

40 _I

SSPA

= 15 A ,,

SSPA

Trim and Sinkage Curves

SGPA MODEL 2082-A

REGULAR WAVE HEADING 0 DEG

H 5 METER T 8 SEC FIG: 16

0.

§ a

0

0 2 8 PQ

00

10.0

0

TRIM ANGLE (DES)

-SINKAGE 00

10.0

HEEL ANGLE (DEG)

CALM WATER

tr

-

- A SAGGING

-F- -

HOGGING

---

A. -0 ---

_,

-1. , -k i _ i i i 114.

20.0

30.0

40.0

+

SSPA

1.5

1.00

0.5

2 0

. (METER)GZ.ETER) GZ CURVE SSPA MODEL 2062-A

REGULAR WAVE HEADING 30 DEG, H = 8 WIER = 10 HEEL ANGLE

mm

0 II

/

A

/

/

CALM WATER SAGGING HOGGING I

A-

A--- A

11

-1 I

/

i

/

/

/

1

/

/

/

A

/

/

i 1

/

/

/I-I

/

A/

/

1

/

/

I I 1 , H

.V

'le

--sO 10.00

20.00

30.00

40 1 GZ

0.

FIG: 17

0

-SSPA

Trim and Sinkage Curves

SSPA MODEL 2062-A

REGULAR WAVE HEADING 30 DEG

H = 8 METER T = 10 SEC

FIG: 18

4r1

₀

cn

TRIM ANGLE (DEG)

40.0

o

o

C ALIN WATEFI At- A SAGGING + HOGGING §.-

i .0

0 Q 0 m 0 2.0 Ti $.4 ro 04

1.0

0. (12, A ik._:=1._=...1r , , 0 10.0 SINKAGE 00 20.0 30.0 _ ,

---415---4---,0 ---

--- +

A --- + .0 rcl --4 cc 3 3 0

--.-.-e

`-0

0

10.0

20.0

30.0

HEEL ANGLE (DEG)

40.0

-SSPA

GZ (METER)

2.

1.

1.0

05

0 0

0 GZ CURVE

SSPA MODEL 2062-A

REGULAR WAVE HEADING 30 DEG

R = 5 METER ,T = 8 SEC

HEEL ANGLE (DES)

19 0, - 0 CALM WATER SAGGING HOGGING

fr

-

-

ti

-

---i-1 1 A

.

ri 1

i

I ,

.

1 .#

A

/

il !

.

' to

/

/

i-/

Z-i

/

A

/

.

/

,

/

../°. I

./

! ...-t,

.--....--....-+ , 1 .6

10.00

20_00

30.00

40 1 FIG:

-i-,

I I

/

I I

//

I

SSPA

§

o

0 0 2.0 0.0 10.0 SINKAGE (M) 1.0 fj 0 2 -2.

00

TRIM ANGLE

mm

10.0

Trim and Sinkage Curves

SSPA NOBEL 2062-A

REGULAR WAVE HEADING 30 DEG H - 5 METER T = 6 SEC 20.0 30.0 HEEL ANGLE

mm

FIG: 2 -0 CALM WATER fr - - A SAGGING + HOGGING 40.0 .0

__ ..

_.

"I.'''. ... _{..., ,...} ...1. ... .0. 4.

SSPA Report 2912-5

A HIGHER ORDER PANEL METHOD FOR DOUBLE MODEL LINEARIZED

FREE SURFACE POTENTIAL FLOWS

by

Shao-Yu Ni

On (2-2)

where n denotes' the outward normal to the hull _{surface- when the}

INTRODUCTION

In the present report a method for calculating the potential

flow

about ships, piercing the free water surface is described. The special feature of the method is that the problem is discretized more accurately than in other methods of the same type.

A brief mathematical formulation of the theory behind methods of the panel type is given in Chapter 2 and in Chapter 3 the first order discretization is explained. Better accuracy can be

obtained with a higher order discretization, as described in Chapter 4, and in Chapter 5 the application to ships, based on the so called double model approximation is reviewed.

The free surface boundary condition has to be imposed when the wave pattern and the wave resistance are to be determined. The boundary condition can be linearized along streamlines as is described in chapter 6 or along arbitrary smooth curves as in chapter 7. Once the velocity distribution on the hull has been computed, the wave resistance can be obtained by using the usual pressure integration or the momentum approach presented in chapter 8. Finally, some calculated results and discussions are shown to prove the improvement obtained by the higher order panel

method in chapter 9.

THE FUNDAMENTAL PROBLEM

A steady inviscid flow past a body can be described by a velocity potential 0 which is generated by a certain distribution of

singularities on a control surface S and by the uniform onset flow in the x-direction (see Fig 1). In this report only sources are employed. The potential 0 is given by

0(x,y,z), = fj(a.(0)/(r(p,q))dS

Ux

(2-1)

where Um is the speed of the onset flow and r is the distance from the integration point q(E,T1,) on S to the field point p(x,y,z) where the potential is being evaluated.

The potential 0. in Eq (2-1) is governed by Laplace's equation

v20 =

and satisfies the regularity condition at infinity

= rim

The source density a should he selected so that certain boundary conditions may be met. On the wetted hull surface the solid

boundary condition is

+

S

wave elevation and wave resistance Cw are to be predicted the free surface conditions should be imposed, i e the pressure

should be constant

1/2(01,2 .4. 0172

oz2 _ u.2)

= 0 _(2-3)

and the normal velocity on the free surface

c/Sxwx 9517`-'wy Oz = (2-4)

The numerical solution begins with the discretization of the

integration domain S, which is expressed by quadrilateral panels. Thus the integral in Eq (2-1) is replaced by summation.

Of

=ij +

where

Oi

is the total potential at the i-th field point and Oii

is the perturbation potential at the i-th field point induced by the j-th panel Aj on which source density a is distributed.

O.. = f(a/r)dS (2-5)

Aj

3. THE FIRST ORDER PANEL METHOD

In this method each panel is assumed to be flat and quadrilateral

with constant source density [1]. It is characterized by the null

point, transformation matrix T and various geometrical

properties. The projections of the orthogonal vectors 1, j and E.

of the panel coordinate system T-1c on the reference coordinate

system xyz are

1 = (ix, iy,

iz)

7

= (ix, iy,

iz)

= (nx, ny, nz)

The transformation matrix T is given by

T=

Oij

=

paj/rf)dS

aj.0(°)

A

(3-1)

The control point at which the boundary condition is to be applied is taken as null point, i e the point where the panel

itself induces no velocity in its own plane. Therefore the

constant source density

aj

on the j-th flat panel Aj which will induce the potential

Oij

at the i-th null point p(x,y,z) in the panel system is 2 . . 1 (3-2) =

where

0(0) . _{f(i/rf) d01-1}

and

I(

X -

) 2

-

)1 2 z2

Accordingly, the velocity induced by a unit source density on the j-th panel in the panel coordinate system is the gradient of the unit potential (PO)

V(0) , v0(0) or (°) (V3x.) f (1/rf)iddri A. (0), .

(a/ay)

_f(l/rf)dEdn

(-5)

Vc (a4

z)

. (1/rOddri

Vc(0), Vn(°) ana3Vc0) can be calculated numerically for each panel. The component Vc(0)

0)

- f (Urf3)dCdn

Ai

requires special handling when z -4 0. Vc('°) 0 as z 0 if p is

approaching a point in the plane outside the boundaries of the panel. If p approaches a point within the panel Vc(°) 2n(sign

z) as z -4 _{O. Here sign z = + 1 because the field point p always}

approaches the hull surface from the exterior flow.

The components Of V(0) in the reference coordinate system _are

called influence coefficients are. computed as follows

Xiji

17V .'( °-)7

1

;

Y.. i_ij = Tim l Vil(0, _(3-6)

, . 1 V 0)', Zi_; ,, L c _J j

4. THE HIGHER ORDER. PANEL, METHOD

The approximation that the real hull surface is represented by flat panels with constant source density unavoidably _causes

errors. If more panels are designed in order to reduce the _errors it would take much more computer time so that the cost might not be acceptable. Sometimes the calculation results are not

acceptable even though a large number of panels is used. The interior flow is an example. Therefore a higher order panel

(3-3) (3-4) 3 rf = + (0) = A. =

M1

= E D E(CkTik)P F(Ek,iik)Rl}2

k=1

Once P and R are known, all other parameters of Eq (4-1) can be

obtained from Eq (4-2).

In this method the control point is the nearest point on the curved surface to the origin of the first order panel coordinate system. Then the panel coordinate plane is moved to a new

projected flat plane which is tangential to the curved surface, passing through the control point as the origin. In the higher order panel coordinate system (Fig 2a) the curved panel of Eq

(4-1) becomes

C = PC2 + 2Qcn + RT12 (4-4)

The outward normal unit vector n,n) is

n(1-1) =

= (nr?, nn, nc) = (-(C/F),-(cn/F), (l/F))

= min

4

method has been developed. This method, which essentially follows a theory presented by Hess [2], is more accurate than the first order ones since the panels are curved and have a varying source density. Of course, the calculation is also more complicated. The curved panel is supposedly parabolic

C = Zo + AC + Bin + K2 + 2QCT-1 + RT-12

(4-1)

in the first order panel coordinate system defined as inchapter 3.Let

the curved surface pass through the four input points of this panel exactly. The coefficients Zo, A, B and Q can be solved as linear combinations of unknown P and R

A = Ao + ApP + ARR

B = Bo + BpP + BRR

(4-2)

O = O0 QPP ORR

Zo = ZpP + ZQQ + ZRR

where the coefficients Ao, Bo, Q0, Ap, Bp, Qp, Zp, AR, BR,

QR, ZR

and ZQ are determined by the coordinates of the four points.

The quantities A, B, Q and Zo are replaced in Eq (4-1) by their representations in Eq (4-2), the result being

= D(1-1) + E(,n)P + F(,T)R (4-3)

The arbitrary values of P and R are adjusted so as to keep the distances from the surface denoted by Eq (4-3) to M1 neighbouring input points of adjacent panels as small as possible in the least squares sence

-[ _(Ek,11k) +

where

= 2(Pr, + Qn)

= 2(Qc + Rn)

and

F = /1

+-c 2

i:-C 2

Some simpler expressions are obtained as follows by dropping

higher order terms of _{,which is small everywhere on the panel in} question. The distance r from a field point p(x,y,z) to any point

(TI,c)

on the panel in question

r2 =

(x

-E)2

(y - n)2 +

(z -

L12 = rf2

- 2zr,

hr = (1/rf)(1 + (zc/rf2)) + . (4-5)

where

rf

is defined as Eq (3-4). The area element dS on the

curved surface is

dS = (1/n

_{)(1Ccin = 11 +}

c 2 4. c

2 don =

_{(1 + 2t1)2)cidn (4-6)}

where

= (p2

Q2)2

2(pQ

QR)Efl 4. (Q2 4. R2),12

On the other hand, the source distribution

a(E,n)

is considered

to be linear on the j-th panel in question

cr(,n) = aj

+

ax

+ ayn

(4-7)

This is fitted in the least squares sense to the values of the

source densities at the control points of the M2 adjacent panels (Fig 2.b) where

2M2 = M1

By minimizing the function J2

M2

J2 =

E [ak

-(aj

axic

aynk)]2=

min

k=0

the derivatives of the source density

ax

and a are solved.

accordingly M2

ax = E

CkkxJak

k=0

(4-8)

M2 f 0

= E

CickY)ok k=0 5 + + +

IL

where ao =

ai

and the derivative coefficients Ck(x) and Ck(17) are

purely geometrical. Finally, the moments of the area of the panel are required

Inm = If

cnrimdCdn ft = 0, 1,...4, rn = 0, 1, (4-9)

Aj

Specially

100

is the projected area of the panel.

Now a more complicated formula can be expanded from Eq (2-5) rather than Eq (3-2), omitting higher order terms of and n

1

Oij = f(a/r)ds

-A.

= A(I (1/tfrMaj-t-axi-ayri)1[1+(z(p2+2Q,En +12T12)/rf2)](1-1-211)2)dEdn

=, a]0(°)+ai(P0(13)+2095M+RO(R))+ax0(1x)+ay00Y) (4-10)

wher-e O(,0) is defined as Eq (3-3)

0(P) f (zC2/rf3)dCd1-1 Aj =

_{5(zWrf3)d1Wni}

A-3 0(R) =

f

(zn2/rd)dain Aj 0(1x) =

f

(/rf)dal-rl Ai O(1Y) =

f

(n/rf)dEdil Aj

Then the induced velocities in the higher order panel coordinat

system are gradients, of the potentials above

i/7(0) , v0(0)

V(P)

Vo(P) V(Q) -

vo(o)

(4-11)

V(R)

vo(R)

i7(1x)

vo(lx)

7i7(1y)

vo(ly)

All the induced velocities above in Eq (4-11) can be calculated numerically and they are explained as follows: V(0) corresponds to a flat panel with a constant source density as Eq (3-5), V(P)

V(Q) and V(R) are caused by a parabolic panel shape, y(lx) and

V(117) come from a linear variation of source density.

11

6

...4

=

7 Especially in Eq (4-11)

V() = (/az)

f

(Urf)dOn

C =

-

f

(Ez/rf3)dEdn = - i [(zx/rf3) + (z(-x)/rf3)]dCdn = x V(()) - tV(°) C C

In the same way we have V (1Y) = vv(0) _ zv(0)

Y

n

In chapter 3 it is shown that as z -0- 0

IV

0)

--, 0 for points outside the boundaries of the panel

C

V (C))

-4 2n

for points within the boundaries of the panel

C

Therefore V (1x) , 0, V (1Y)c 0 as z-9 0 if the field point P

approaches i point in the r1 plane outside the quadrilateral of

the panel. V() , 0. VE(1Y) -, 0 as z -, 0 even though P is approaching the origin of the panel in question. Because the control point is taken as the origin and boundary conditions are only satisfied at the control point it is concluded that Vc (1x)

0, V(1Y) --. 0 as z -4- 0.

The influence coefficients

Xij, Yij

and Zi mean the velocity components in the reference coordinate system XYZ at the i-th control point induced by unit source density on the control point of the j-th panel where linearly varying sources are in fact

located. The influence coefficients are handled according to the

following procedures.

First the induced velocity is written as

.4_[12 ck(x)iylx) + ck(y)T7(1y)]ok ij k=1 (4-13) where

vE* =

vEo)

v * = v(p)

n Nyr = V(0)

Thena transformation similar to Eq (3-6) is carried out

+

pvV)

+ 204Q) +

Rvr)

4. c0(x),q1x)

+ co(y)x(iy)

+ pv(P) I-)

+ 2w(Q)

fl

+ Rv(R)

fl 4. co(x)v(ix) ri + c (y)v(117) 0 rl + PV(P) + 2OV(Q) + RV(R) + C0(x)V1x) + co(Y)V117) C C =

Consequently the influence coefficients are Xj = \Tx* + r'

12[cit(x)v31(1x)

cil(Y)v21(1Y)] k=1 M2

= v

* E [q(x)v*(1x)

_C(Y)v(-Y)]

k=1 M2 vz* E [q(x)v1(1x) c1t(y)v*(11,)3 k=1 (4-15)

where 4(x), q(Y) and

v*(

lx) v*(1x), vI(lx), v*(1Y), v*(1Y),

VI(1Y) correspond to the source Lrivative coefficients and induced velocity components of the K-th panel in question surrounding the j-th panel. The latter terms in Eq (4-15)

multiplied by a. are also the contribution of the j-th panel to

the total indu d velocity at the i-th control point.

It is obvious that only three integrations in Eq (3-5) are computed for one induced velocity in the first order panel method, while in the higher order panel method further 15 integrations in Eq (4-11) should be calculated.

5. THE DOUBLE MODEL SOLUTION

Only the boundary condition on the wetted hull surface (2-2) is applied when the potential flow without free surface is

investigated. Since the influence coefficients Xij. Yil and Zij are obtained from the first order panel method, Eq (3-4) or from

the higher order panel method, Eq (4-15), the velocity components at the i-th panel induced by all the NE panels on the hull may be written as 8 [ix* VY* V,* -vx(1x1 si (1x)

v

(1x) z (1Y )-V (1Y) V (1Y)_z = TT = TT = TT * Vr* Vi;*

v(lx)

v(lx)

v(lx)

-

-V(1/7) V(1/0 C -(4-14) + z

NB Oxi = E Xijaj + Um NB Oyi = E Y..a. (5-1) NB Ozi = E Zijaj

The normal unit vector ñ1 on the i-th panel is

= (nxi, nyi, nzi)

The boundary condition, Eq (2-2) becomes

On = Oxinxi Oyinyi Ozinzi =

or

NB

E (Xijnxi + Yijnyi + Zijnzi)aj = - nxiUm

i =

1,2 ... NB

(5-2)

The unknowns a. (j = 1, 2 _{... NB) can be obtained by solving the}

J

linear equation (5-2). Iterative procedures may be used because the matrix of the coefficients is diagonally dominant. The

velocity distribution Vi is computed by inserting the ad-values

into Eq (5-1)

2

24-V* = (0_{xi 0a '2ri}

(5-3)

To illustrate the improvements that the higher order panel method has been able to obtain compared to the first order panel method, some cases have been run for which highly accurate solutions are

available. Results are presented in chapter 9.

The hydrodynamic force acting on the hull is usually most

interesting. The pressure coefficient Cp is calculated from the

velocity distribution Vi of Eq (5-3)

Cpi = 1 - (Vi/U.)2 _(5-4)

The dynamic force in the x-direction is thus computed and

norma-lized as

9

NB NB

Cx = - (E Cni nxiASi)/(E AS)

i

(5-5)

NB

cj=

+ E xijaoj

NB

= E Yijaoj

Y1 .

where aoj is the known source density at the j-th panel from the double model solution and the derivatives of 0 are

NB+NF

Oxi = + E Xijaj

NB+NF

Oyi = E Yijaj (6-3)

0,i =-2nai

where a. is to be solved. Because the free surface z = 0 is the

symmetr4 plane the velocity in the z-direction is zero for the double model case. For the first order panel method =

-27o. as is known from chapter 3 and the negative sign denotes the hwnward direction of the normal vector on the free surface panels. For the higher order panel method

ozi = ozi(0) ozi(C) 0i(lx) ozi(117)

where

0,(C)

= IDO,i(P) + 2Q0zi(Q) + Ozi(R)

lz)

= 0 because all the panels

on z = 0 are flat and

o(x)

oi(ly

= 0 as is also mentioned

(6-2)

10

6. THE FREE SURFACE BOUNDARY CONDITION LINEARIZED ALONG STREAMLINES Before the free surface is involved the double model solution

without the surface must be obtained, which means that the double model potential denoted by (I) is known. Now NF panels on the free

surface, z = 0 along the streamlines of (I) are added (Fig 7). Thus

the influence coefficients Xi, Y. and Z. can be calculated for

NB + NF panels by using the E/rst drder pariel method, Eq (3-6) or

tHe higher order panel method, Eq (4-15). Applying the free sur-face conditions (Eq (2-3) and Eq (2-4)) on the symmetry plane z =

0 neglecting higher order terms, the linear free surface condition

is given by Dawson [3] or Xia [4]

(4'12°1)1 gOz = 24'124'11 (6-1)

where the subscript 1 denotes differentiation along a streamline

of the known potential 4, and

0

is an unknown potential. At the

i-th panel on i-the free surface i-the derivatives of 4, are

in chapter 4. Therefore

Ozi = Ozi(°) =

211a-Now let li be the tangent unit vector to the double model streamline at the i-th panel on the free surface

= (lxi, lyi) where

lxi = (4,xi)/(i2 + 4,yi2)12

lyi = (Pyi)/(4,xi2 + Syi2)1/2

The terms in Eq (6-1) are expressed as

= = (4'xi2 4)1Ti

NB+NF

011 = Ox-xi Oylyi = U.lxi E (Xiilxi + Yiilyi)ai

and

NB+NF

(4)129514 =

cD1i2 (xiilxi Yiilyi)aj

To introduce Eq (6-1) into a set of linear equations the coefficients of the four point, upstream, finite difference

operator CAi, CBi, CC i and CDi are employed as follows

11i = CA111 +

+ CCioD1i_2 +

and (1.120.1)11 is obtained in the same way as in Eq (6-5).

For numerical calculation Eq (6-1) becomes

NB+NF

E + Yiilyi) +

j=1

CBi (Pli-12(xi-llxi-1 _Y1-1jiy1-1)

+

+ Yi_3j1yi_3)]ai - 21tga1 =

24.112(CA1t,112 +CBjcIl_l+ +

li xi+ CBi .1i-12U.lxi-1 + +

CD. 2U1

(6-6) 11 (6-4) (6-5) -.xilxi .yilyi

4

CB 1.1i-1 CCiolp + + CD1cP11_3)

-The boundary condition will be satisfied at each panel on the

hull NB+NF

E (Xijnxi +

Yijnyi

+

Zijn,i)aj

= - n

xi

NB linear equations are offered by Eq (6-7) and NF linear

equations by Eq (6-6) to solve unknowns

ai

(j = 1, 2 ... NB +

NF). Gaussian elimination is used since the coefficient matrix of the complete sets of Eqs (6-7) and (6-6) is not diagonally

dominant.

Oncethesourcedensities

aj

(j = 1, 2 ... NB + NF) are

determined, the full flow velocities

Vi

and the pressure

caeffits at the i-th panel on the hull are obtained. The

Cpi.

wave resistance

Is

calculated by summing the X-components of the pressure forces acting on the hull panels. The wave resistance coefficient Cw could be written

NB NB

Cw = -

(E CpinxiASi)/(E ASi)

or Cw is obtained by the momentum approach given in chapter 8. 7. THE FREE SURFACE BOUNDARY CONDITION LINEARIZED ALONG ARBITRARY

SMOOTH CURVES

More panels are expected to be distributed around the bow and the stern, because the panels on the hull and free surface in these regions play an important role for generating waves. However, the shape of some hulls (see Fig 15) at bow or stern is rather blunt and there can only be a few panels on the free surface when the

longitudinal strips of panels are along streamlines. Therefore a set of smooth arbitrary body-fitted curves are generated (Fig 16) and the boundary condition (6-1) is rewritten as follows.

Unknown sources u on the hull and plane z = 0 will induce a

poten-tial 4) and a wave elevation where the subscript w is omitted

for simplification, which satisfy the boundary conditions (2-3)

and (2-4) 12 (6-7) (6-8) 1 D1

4)xCx + (PyCy -

(/)z = 0 2 D2 cpx2 _{+ qDy2} + (Pz + 2gc -

U,2 =

0

But now onlg the double model potential (I) and Bernoulli wave

elevation C are known

o o 1 D10 = cipxCx + c1)17C17 ° D -4) 2 + (1) 2 +

2e

u 2

0 20 x _Y (7-1) (7-2) -10

The increments of D1 and D2 would be found to satisfy

1

DI -# D1 0 + 6D1 ,=, 0 II

r)2 = D20 t 6D2 = 0

where 6D1 and 6D2 are determined by small perturbations (Sq), and (Sc

based on 0 and OD

6D1 and 61312 can be expanded as

SD1 = 6cp

x 'x

o

+ 4 c

o _.,_, 6(P + 0 (5

Y Y x x + _0YcSCY

6D2 =k 2(0 (SO

xx

t 41

yy

6 cp ) + 2g6c

Therefore the boundary conditions are Linearized by inserting

(7-2), (7-4), and (7-5) into (7-3)

Ic° = axict" - ax2c1:

I

Cy,9 -ay/C.0

ay21:

and i1512t '''' axl6Ct

-

aX215CL 45C17 ' ''41716Ct + 41726L where 4171 # 1/(f'1

-

fi

)il

-1-' FY

(7-4), (7-51 OXCX° OZ -F° 4)x08, x (1)178Cy 0 (7-6) 6C = - (1/g)[4,x0x + (1) 0_{y y}

-

_x2 IF 4,/72)1 _(7-7)

Further treatment is given Eqs (7-6) and (7-7) for the

derivatives of and SC. The longitudinal body-fitted curve is y = fL(x) and the transverse grid curve is y = ft(x). In the

reference coordinate system the following relationships can be found

(7-8)

(7-9)

Especially when the transverse grid curves are parallel to the 11-axis these transformation coefficients become

ay2 f't)/1 _f'L2 axl -=" aylf"L ax2 = ay2f't 13 (7-3) = + = + = -= - +

ayl = 1

ay2 0

1 = F'L

Ax2

-117+

F'L2

The derivatives along the L- and t-directions Can, be calculated

using 4-point finite differences A

+ CBLC °i-NL + CCLC ei-2NL + CD1C °1-3NL

1C°ti = GALC°L

4

GBiC°i4.1 + GCiC°i+2 + CDiCei.1.3

1

oCIA = CAOCi + CBioc

I

i-NL 41- CCOCi-2NL + Cpiei-3NL 6cti = GAi6ci + GBo 4.1 + GCOCi+2 + GDOC1+3

where NL is the number of longitudinal Strips on the free surface

and CAL, CB, CC i and CD i are the coefficients of 4-point

backward difference, the same as the coefficients in Eq (6-5):

However, GAL, GBL, GC i and GDL are the coefficients of the

forward difference operator, used in the transverse direction.

Consequently c°x,c°v in Eq (7-8) and 8, 6cy in Eq (7-9) can be

expressed respectively as linear combinations of several wave 1

elevations (7-10) and several changes of wave elevations (7-11)1 The free surface boundary condition can be simplified as a set of

linear equations by inserting Eqs (7-8) and (7-9) into (7-6) and eliminating terms of SC from (7-7)

Y. Inc = = 14 (7-10) (7-11)

NB+NF

E [(ki + k3GAicl'xi + k4CAi4xi)Xii +

(k2 + k3GAicDy1 + k4CArDy1)Yii +

k3GBi( _{(Dyi+1117i+11)}

k3GCi( 4'xi+2Xi+2j (Dyi+2Y1+2j) k3Gpi(4.x1+3X1+3i 4.y.i.+3Yi4.3j)

k4CBi( xi-NLXi-NLj 4)yi-NLYi-NLj)

k4CCi( 49xi-2NLXi-2NLj .171-2NLY1-2NLj)

k4CDi( "Dxi-3NLXi-3NLj cpyi-3NLYi-3NLP]0j 271c71 =

- k3[GA1(V012 -

1.xit40

GB1(V014.12

GCi(VOi+22 Gp1(V01+32 -

4,x1+314)]-1-k4[CAi(Voi2 - (1'xiU.) _{CBi(VOi -NL2}

2 2

CCi(V Oi-2NL Sxi-2NLU.) + CDi(Voi -3NL

where

kl = axl(GAiCi* + GBiCi+lo + + GDiCi+30)

ax2(CAiC1* + -i-NL* + CCir_{-i-2NL CpiCi-3NL)}

k2 =-ayi(GAiCi° + GBiCi+1° + GC1C1+2 + GDiC1+30) +

ay2(CA1C1° + CBj. -i-NL° + CCiE-i-2NL° + CDi Ci-3NL)

k3 = - 1/g(ax14,xi - ay141)

k4 = 1/g(- ax2x1 ay24'yi)

and

xi-3NL14.)] 4)xi+11J.) +

V0i2 = 4,xi2 4,yi2

In Eq (7-12) there are NF linear equations which are similar to Eq (6-6). NB + NF unknown values of ai are solved from Eq (7-12) and (6-7). Eq (7-12) is a little more complicated than (6-7) but the free surface grids are not limited to be streamlines.

8. THE MOMENTUM APPROACH

Besides pressure integration to calculate wave resistance there are some other ways, among which the momentum approach is simple

and effective.

The momentum conservation equation in the x-direction is applied

15 (7-12) + -(1)xi+2140 + -GCiC1+2

to the control volume within the surface S shown in Fig 7%

Note that the interior of the hull is outside the control volume.

if

pcp.T7 FidS ₍₁₃ + pgz)nxdS = 0 (841 S. 5 where NB+NF E X1. + Um

N+N

0 v B 13 3 NB NF w E

Z..a.J

1J and (0x, 9517, Oz)

The mass conservation equation for the control, volume is

If V !lids = 0 (3-2)

Eq (8-1) minus Eq (8-2) multiplied by pU. yields

If

PO

*RdS + ff. (p pg2)nxdS = 0 (8-3)'

S IS

The free surface SF consists of two parts: SFc which refers to the part covered by source panels and SFu which refers to all the other free surface outside SFc. It can be assumed that the front plane, aft plane, both side planes and the lower plane arel

located infinitely far from the hull. The velocity in these planes is equal to the onset flow U. Since

ffPgznxdS - 0, and

f1 pnxdS51

pnxIS =

-Dw S SB, Eq (8-3) Can be simplified. as 'Dw pnxdS PI I = SB SF pff uw dS SF Furthermore, (8-4) 1 6 = S + =

-.rj

u-V*EdS

=

-= + = = =

NB+NF NB NF

u = E= E

X..

a. +

E

X..0.

_J

= + uF j=1 j=1 J j=NB+1 J J NB+NF NF

w=

Z..a.

=B

Z..o. + E 3 Z..u. = wB

+w

F j=1 17 j=1 13 3 . =NB+1 13

where index B refers to the hull sources and their images, F to the free surface sources. Because of the double model

lineariza-tion

wB = 0 on free surface SF. Eq (8-4) now becomes

Dw

=jJ

uBwFdS

-p1

fUW'dS

-

1 uFwFdS (8-5)

SFc SFu SF

If it is imaged that the fictitious flow is generated by the

source distribution only on the free surface SF, the

conservation

of momentum holds also

-P11

uFwFdS -p11 uFwFdS - 0 (8-6)

SF SFi

where SFi is the part of the plane Z _{= 0 inside the hull. Then it}

is deduced that

wF = 0 on SFu and SFi

because no source is distributed on the two parts SFu and SFi. Finally, Eq (8-5) is reduced to the following expression

Dw =

-P[f

uBwFdS (8-7)

SFc

In the practical calculation the wave resistance coefficient can

be expressed as

NF

Cw -

(2 E

_{ugiwFiASi)/(U2EBAS1)}

If it is noted that wFi =-2ncr1 as in chapter 6 Eq (8-8) can be

written NF N Cw = (4n E

uBicriASi)/(U2

EBLS.) _(8-9) co

i

i

i

Comparing with Eq (6-8) NB 2 NB Cw

- (E (1 -

Vi2)nxiASi)/(U. T AS,)

i

(8-8)

17 (8-10)