The development of numerical methods for the compution of unsteady propeller forces

3 SEP. 19!I

ARCHtE

PAPER 23/8 - SESSION 2

(I)

Lab. v. Scheepsbouwkund

Technische Hoyeschool

Deift

SPQNSOR: DET NORSKE .VERITAS SYMPOSIUM ON

"HYDRODYNAMICS OF SHIP AND OFFSHORE PROPULSION SYSTEMS"

HOVIK OUTSIDE OSLO, MARCH 20.-25., 1977

"THE DE'1ELOPMENT OF NUMERICAL METHODS FOR THE COMPUTION OF UNSTEADY PROPELLER FORCES".

o

By

0. Frydenlund, Det norske Veritas, Oslo J.E. Kerwin Dept. of Ocean Engineering,

By

Oddvar Frydénlund1 Justin E. Kerwin2

Prepared for the Symposium on

HYDRODYNAMICS OF SHIP AND OFFSHORE PROPULSION SYSTEMS Oslo, Norway

March 20-25, 1977

Sponsored by: DET NORSKE VERITAS

1Det Norske Veritas, Visiting Engineer, Department of Ocean Engineering, M I T , 1975-1976

2Professor of Naval Architecture, Department o Ocean Engineering, M.I.T.

ABSTRACT

This paper is concerned with the development of a discrete vortex method suitable for application to the problem ofnumerical computation of unsteady propeller forces. The method is developed in detail for two-dimensional flows,andsystematic results are given to demonstrate the convergence of the method and to check agreement with known

analytical results. Extension of the method to three-dimensional flows is briefly discussed, and preliminary results are given.

The authors wish to acknowledge the contribution of Mr. Scott S. Tsao, Research Assistant to section IX and X of this paper. Mr. Tsao had

a major role in the Writing of the two computer programs which are described in these sections.

This research was initiated at M.I.T. under the sponsorship of Det.rrske Veritas and continued under the sponsorship of the U. S. Naval

Sea Systems Command General Hydromechanics Research Program, admin.isterd by. the David Taylor Naval Ship Research and Development Center contract

N0014-76-C-0357,

and the United States Maritime Administration, contract

6-38068.

TABLE OF CONTENTS

Title Page

Abstact

2

Acknowledgements

3

Table of Contents

4

List of Figures

6

Nomenclature

9

INTRODUCTION

. . . 12

II

FORMULATION OF TWO-DIMENSIONAL PROBLEM

15

III

DISCRETE VORTEX ARRANGEMENTS

17

IV

NUMERICAL KUTTA CONDITION

2.3

DISCRETE TIME STEP SOLUTION OF THE

BOUNDARY VALUE PROBLEM

29

VI

FORCE. CALCULATION 31

VII

TWO-DIMENSIONAL RESULTS

. 33

TABLEOF CONTENTS

IX

STEADY-STATE PROPELLER PERFORMANCE

PREDECTION

61.

UNSTEADY PROPELLER BLDE LOADING

63

FIGURE CAPTIONS

Figure 1

Distributions of vortex and controlpoints along chord.

a) Cosine spacing

b) Constant spacing.

Figure 2

Vortex distribution as a function of chord.. Vortidity

con-tinoüs at trailing edge.

Figure 3-6.

Vortex distribution on the foil and in the wake at the

time of maximum oscillatqry lift.

Heaving motion reduced frequencies k

u /2. Number of

ele-ments along chord 10, 20, 30,. 40 respectively.'

Figure 7-10

SinusodaI Gust, reduced frequency k = u/2. Vortex

distri-butiön on the foil and in the wake at the time of

maxi-mum oscillatory lift.

Nu.irther of elements along chord 10,

20,

30, 40 respectively.

Figure 11

Oscillatory lift coefficients Versus reduced frequence.

6

Phase angle (in degrees) as a function of reduced

frequen-c:ies.

Nu.thber Of elements along chord: N = 20.

Figure 13-18

Vorticity distribution along chord and in the wake.

Figure 13-18 shows the effect of reduced frequencies on

vorticiiy.distr. for heaving motionk

u, ir/2, n/3,

iT/4, u/5, it/6.

Figure 19-24

Vorticity distribution along öhord and in the wake.

Figure 19-24 shows theeffect of reduced frequencies on

vorticity distribution for sinusodial gust.

k = ii, Tr/21 u/3, ir/4, ii/5, ir/6.

Figure 25

Vorticity distribution along chord and in the wake in the

case of step change in angle. of attack.

A wake-length corresponding to six chords were used in

the numerical calculatiOns.

8

-Figure 26

Lift coefficients 'versus reduced frequency. for different

äspect-ràtioS.

Figure 27

Computer plot of vortex and controlpoint arrangement on the

camber sUrface of a propeller.

Inthis particular case cosine spacing have been used.

Figure 28à-h

Shedding of radial vortex element in unstady motion.

The motion is stared frOm rest with all blade at mean

loading.

Blade 1 is taken to be 'the reference blade on which the

boundary conditiOn is solved.

Figure 29

Thrust coefficient asa functIon of 'angular position.

For the first 3 revolutions a 30 degree time-step has been

applied..

NOMENCLATURE

A

Projected area of the foil, (fig. 26).

Constants defined by equation 14 and 15.

b1

Weight factcrs defined by equatiOn 17 and 18.

Chordlength of foil.

c1...c15

Algebraic coitthinations of the weighting faOtors b1.

Total lift coef:ficient.

F(K)

_{General function of space and time.}

G(x)

Non-dimensional vortex sheet strength.

Position of lifting surface.

Time step i.ndex

eferi

to present time slep.

Reduced frequency

c/2U.

NOMENCLATURE

Number of elements along chord.

Pressure at infinity (equation 24)

'L

Pressure on upper and lower surface of the Vortex

sheet respectively.

Free stream velocity.

'

L

Pèrtur]ation velocity on upper/lower surface of

vortex sheet.

V(x,T)

General expressiofi for velocity.

Ve(K)

External vertical velocity at a controlpoint

at the K!th time step.

e

Location of controlpoints along foil. (eq. 7).

Spacing between vortex elements on foil (eq. 8)

6W

Spacing between vortex elements in wake (eq. 9).

Location of discrete vortex elements on foil.

10

-y(x)

Vortex sheet strength (eq. 14 and 15).

TN-i

Refers to vortex strength at pOsition N-i, i = 0

> last element on foil.

TN+i

Vortex elements in wake i

1,

Total circulation around foil at the K"th time

step.

Strength, of consentrted vortex shed into the

wake at

'th time step.

12

-I. 'IN:TRODUCTION

During the past several years considerable attention has been given to the problem of ship vibration caused by transient propeller cavitation. Both model and fUll-scale hull pressure measurements have shown conclusively that this can, in many cases, be the primary source of hull excitation. Theoretical prediction of this phenomenon.reqiiires a knowledge, in particular, of the time history of cavity volume. This, in turn, requires a knowledge of the unsteady pressure distribution on the blades as they traverse the non-uniform ship wake field. This problem is far from trivial even for non-cavitating flows, and becomes increasingly complicated as the interaction between the cavities and the flow field is brought into consideration.

A satisfactory solution to this problem clearly requires an elaborate numerical procedure in which the shape and extent of the cavities are established by an iterative procedure beginning with the

unsteady non-cavitating pressure distribution. The procedure employed to compute blade pressure distributions must be capable of producing accurate results for both mean and fluctuating cOmponents of blade loading, since both directly influence the shape of the cavity.

The development of such a program has been underway at MIT since January, 1975. The first phae of this research has been, the develop-ment of a numerical discrete-vortexmethod for the computation of

This method is an extension o,f the discrete vortex/source _procedures developed for propeller design /1/ and for the analysis o'f pro.el1er performance in steady flow /2/.

This discrete element approach appears to offer _{number of} advantages as a starting point for the cornutation of unsteady, partially cavitating flows:

It is capable of yielding accurate predictions of mean loading both at design and off-design conditions.

Being a numerical procedure, blade geometry _{can be} incorporated exactly so that propellers with large skew, rake and varying pitch distribution can be accommodated. _{This is considered essential, since} it is. through the variation of these parameters _th.t optimum propeller designs can be evolved.

(C) Since the procedure includes all three _components of induced velocity, there is no particular problem in including tangential and radial wake field _components.

(d) Since no loading mode functions are employed, the

modifications ultimately required to include the cavities' would appear to be feasible. _{Source elements presently} included to represent blade thickness can assume the further role, of representing cavity volume.

(e) A discrete element method lends itself naturally _{to a} step-by-step time domain solution, which is also essentia1 for the subsequent inclusion of unsteady

r

14

-Two fundaEnental additions are required in conve.rting from steady to.Unstèady flow. _{Radial vorticity must be shed into the wake at the} trailing edge and. convected downstream, and _{the.contributjon of the} unsteady te.rm.in Bernoulli's equation _{must be included in determining} blade pressure loading.

In order to insure that these Unsteady effects could be properly incorporated in the existing discrete _{vortex scheme, the method was} extensively tested for simple two-dimensional _{flows consisting of a} flat plate undergoing vertical oscillation, _{traversing a sinusoidal} gust and starting suddenly from rest within _{a finite angle.of attack.} This permitted direct comparison with known _{exact solutions unhampered} by the complex geometry present in _{a realistic propeller flow. Due} to the low cost of such computations, calculations _{could be made to} explore convergence of the method with large _{variations in the number} of discrete vortex elements.

As the simplest possible tests of three-dimensional unsteady motions, the procedure was then applied tO the case of rectangular flat plate undergoing simple harmonic oscillation.in _{heave. The method} was then tested for a range of reduced frequendes and _{aspect ratios} comparable to those of a marine propeller.

Finally, the method was applied to a case of a marine propeller of arbitrary form. _{The computer program has been completed and the} essential process of testing and refinement is und.erway. .

In the present paper we will be concerned principally with the development and testing of the numerical _{procedure for'two-dimensional}

unsteady flow. An outline of the general procedure is given for the propeller program, together with a preliminary reult.

II. FORMULATION OF THE TWO-DIMENSIONAL PROBLEM

We will first consider the problem of a two-dimensional thin hydrofoil advancing with constant velocity U , which may either be passing through a. spatially varying velocity field, orundergoing

some form of rigid-body oscillation. The linearized problem will

be solved, to enable comparion of the numerical results wIth known exact solutions;

A. fla.t plate, with chord length, c , is situated on, or close to the interval (.o,c) of the x axis in a flow with speed U Oriented in the pOsitive x direction. The various motions being considered may all be specified in terms of a vertical velocity v(x,t)

prescribed along this interval of the x axis. For simple harmonic heaving motion

v(x,t) _Re

{voet

}

while for 'the sinusoidal gust problem

v(x,t)= Re (2)

Finally., for the Wagner problem,

v(x,t)= 0 for t < 0 (3)

16

-The integral equation for the distribution of. vorticity -y:(t) over the foil may then be, written in the following form

C

v(x,t) +

I

c

y(t.)

d +

I

2rr 27r

0 C

where and r are dummy coordinates on the foil and wake, respectively, and

y1

is the strength of the shed vorticity in the wake. The induced velocity at.a point x on the foil is therefore the sum of the

contributions from the vortiity on the foil and vorticity in the

wake. The upper limit of integration in the wake integral will increase continuously with time as the motion progresses.

According to Kelvin's theoren,, the total circulation of the system must remain zero,

C . Ut j ,t)d +

J ;(t)

d 0 (5) 0 C

so that any change in circulation on the foil must be followed by. an equal change.of circulation in the wake.

In particular the circulation around a differential element of wake d , at position ri behind the trailing edge is the difference

in circulation on the fOil at time t + and t - /U , t

U

being measured from the start of the motion.

We complete the formulation of the problem with a statement of the. Kutta condition,which requires that y (x,t) be bounded at the trailing edge.

Ut

(4) d = 0

Exact analytical solutiOns to the linearized problems thus formulated can be obtained by several classical tehniques and are

therefore available as indisputable mèasurésof accuracy of the, numerical procedure which we will next develop.

III. DISCRETE VORTEX ARRANGEMENTS

A discrete vortex representation of a two dimensional hydrofoil in unsteady flow can be characterized by the followingfour _elements.:

The choice and nuniber of circulation node functions. The number and distribution of the discrete bound vortices along the chord.

(c). _{The. number and distribution of the discrete shed vorticies} in the wake.

.

(d) _{The number and location of the control points When the}

boundary condition is to be satisfied.

While these elements.are to some extent interrelated,, many possible ornbinations exist., and in very fw cases is there any analytical

justification behind a particular choice. The fact that the accuracy of the solution of a numerical lifting-surface problem is _critically dependent on. this choice is, well known to most workers in this field.

One of the few, analytical justifications in existence has been qiven by R. James /3/, who showed that the discrete vortex solution

converages in the limit of infinite number of vorticies to the analytical solution for both a flat plate and a parabolic camber line for one

particular arrangement of vortices and control points. _This arrange-ment, which is illustrated in figure i consists of dividing the chord

a)

CONSTANT. SPACING

1. E

L. .E CONTROL POINTS CONTROL POINTS F VORTE X ELEME N IS

/VORIEX.

ELEMENTS

T.E

1. E

FIG. 1

b)

COSINE SPACING

into N _{equal intervals, and piacinu the vortices, at the. quarter chord}

and the control points at the three-quarter chordof _{each interval. The} solution. for the N _{discrete vortex strengths is otained directly by}

collocation at the N _{control points, without the use of chordwise} mode functions. _{This arrangement has: the additional advantage of} yielding the exact value of lift in these two cases for all values àf

N . The location of control points mjdway between adjacent

vortex elements is essential for.convegence to the Cauchy _{principal va.lue'of}

the underlying integral, which is a _{necessary condition for all} discrete vortex arrangements.

Mother possible arrangement is obtained, by. subdividing _the chord into equal intervals of the transformed Glauer't _coordinate

0 = cos1 (1

_{- ?}

(6)

which results in a concentration of elements near the leading and

trail'inq edges.

Disretevortice

are ther'i located .at the mid-point of each sub-interval. _{Control points must again be located midway between} adjacent vortex points, which .res.ults in their relative _{position in} each sub-interval not being a constant. _{This arrangëmênt is illustrated}

'in figure lb.

The latter spacing is a natural choice if the Glauert chordwise mode function as employed, and it can be demonstrated that the flat-plate term the Glauert series becomes linearly independent of _the remaining terms of the series if this spacing is used in _{a collocation} procedure. . However, we have not been .successful in developing a stable

20

-the addition of a variety of trial circulation mode _{functions '.iith} finite values at the trailing edge. _{This led us to the decision to} try a discrete vortex approach in which eac.h vortex element was an unknown to be determined. _{The, latter approach ias found to be} inherently stable, and was considered to be _{clearly superior when.} planning for future extention of the method to partially _cavitating

flows. The principal reason fOr avoiding _{a hmodelessu approach in} the past was the fear of the excessive _{computing cost associated with} the inversion of large matrices. _{However, computer economics has} changed significantly in recent years with the result that this is no lonoer a valid objection, even in the case of three-dimensional

propeller lifting-surface theory.

It was still not obvious at the outset whether a constant spacing or.a constant Glauert angle spacino would work best for' unsteady motion. In this case it is necessary to consider the _{influence of the discrete} approximation of the shed vorticity in the wake. _{For a discrete time'} solution, vortices are shed at uniform time, _{and 'hence, uniform}

spatial intervals. _{An important parameter to be optimized is the} position of the discrete shed vortex _{within each advance interval.}

U6t.

This proved to be a more sensitive parameter when _Glauert

spacing was used, due to the abrupt change in _{spacing when passing into} the wake frorn"th trailing edge. _{est results were obtained with}

constant spacing on the foil, with the _{same spacing maintained in the}

wake. _{However! it was found that even a doubling of the wake spacing}

The _{remaining consideration in choosing between these two} spacing options is the question of accuracy in the prediction _{of moan} loading due to camber. While this is riot of immediate concern, it is a matter which must. be ultimately faced in dealina with. _the propeller problem.

In particular, James /3/ showed that convergence was very slow using constant spacing to approximate an MACA a = 1.0 constant load camber line. We have verified this conclusion and have found that Glauert spacing is more accurate in this particular case. This result is not surprising in view of the concentration of vortices _{near the} ends in the latter case.

However, we have made calculations for the _{NACA a= .8 mean line} using both spacing options, and find them to be of comparable

accuracy, with errors substantially less than for the NACA a = 1.0 mean line. _{Evidently the logarithmic singulaity at the trailing edge} causes the most trouble, and this is fortunaly absent in the

.NACA a= .8 camber line.

We have, therefore, concluded that constant chordwise spacing is optimum, being "exactt1 for steady flat plate and parabolic camber lines, better for unsteady loading, and of the same _{accuracy as Glauert} spacing for .NACA a= .8 camber lines. The inaccuracy of constant spacing in the case of a = 1.0 camber linesis not considered of

-practical importance since real fluid effects, for very similar reasons, substantially degrade their performance in actual use.

22

-We can now establish a nornenclatui' for the various quantities needed to construct a discrete vortex computational procedure. We will non-dimensionalize all lengths with respect to the chord length, c

so that the locations of the discrete vortices, , and the control

points, x0 , can be written

=l,2...N.

The position of the first vortex shed fn the wake is x

while all successive vortices tn the wake. are uniformly spaced with a

ëommon interval

6 UcSt.

w C

The strength of each concentrated vortex at a particular time step _tK

= (K-

1) ót is

fl= 1,2, N (10)

while the total circulation arOund the hydrofoil t the K'th time step

is N

r(K)

(8) (9) (7) (, -0.25) 2. = 1,2 ... N

We will adopt the iotation that k _{denotes the present time step,}

whfle. K S: a dummy index used to sum oyer the past history of the motion. The time index is omitted ii (10) with the understanding that the, distribution of circulation over the chord pertains to the present time step. . .

The strength of the concentrated vortex shed into the wake at

the K 'th trne step is

r(K)

₍₁₂₎

whose position, relative to a coordinate system moving with the foil

is

x+ (k--),

(13)

Finally, the averae vortex sheet strength in each wake interval can be obttined by dividing (12) by -,

IV. . NUMERICAL KUTTA CONDITION

It is necessary tO impose 'a suitable form of the Kutta condition at the trailing edge to. insure that the numerical solution will converge to the desired analytical result'. The proper Kutta condition for the anal ytical solution is that the vorticity be bounded at the, trailing edge. This is not astrong enough condition to impose. in anumerical scheme since the discrete vortex strength must obviously be bounded,

even at, the leadi.ig. edewhëré.t expect a square-root singularity. Yet

we must be careful not to impose a cOndition which i.sunnecessarily

24

-There are many possible approaches which could be adopted. We

have tried several, and find thatoractically anything will work in the sense 'that the results converqe to the right answer given a

sufficiently larqe number of elerents. Obviously, the best numerical Kutta condition is the one which converges to within the desired 'accuracy with the minimu,m number of elements.

A scheme which works well in this respect is illustrated in figure (2). Near the trailing edge, the vortex sheet strength on the

foil is assumed to be a parabola of the form y(x) = a0 + a1 Cx - c) -i-' a2(x - c)2 x < c

while in, the Wake immediately behind the trailing edge, it is a linear function

y(x) = a0 + a(x - c) x >c . (15)

Thus the vorticity is assumed to be continuous at the trailing edge with a value of a0 , but need not have a continuous slope. The four

unknown coefficients 'a0 - a3 may be found by requiring that the approximating parabola match the discrete vortex sheet strengths

1N-3

and

1N-2

on the foil and and y'

N+2

in the wake.

'We then have two dependent vortices

- and N which, are required to lieU on the parabola. We, therefo-e, need to express

and

1N' as a linearcoribination of

N3

''N-2'

N+l

and

N+2

For the present development we have used the notation

KUTTA CONDITION

to approxirate the vortex sheet strengths _{in the first two wake} intrva1s behind the traiiing edge.

To do this, we first construct the _{coefficient matrix of the} set of sinultaneous equations '.ihich require that the approximating curves pass throuoh the base points locatedat

N13

N- 2 4 A a ₌ , for i 1, 4 j=1 - C. - C 2.

2-c

- c)2 0 ₍₁₆₎ 0 0

x-c

0 0

x+_c

where the index i _{refers to the four base points.}

We can next solve

for.

N - 1 in terms of the base points by inverting (16) and substitultinq

the results in (14). . . .

26

-1 _k

- c)

A2'

+

c)2A37

j!l {

A1 ± N-1 - c)A2 'N -1

All

4 b1 'y ₍₁₇₎ j=l Similarly,

1N can be expressed as a linear combination of. the base point circulations:

=

h2

y (18)

The weightingfactors, b , depend only on the ratio of vo'tex spacingon

t.he foil and in the wake, and and on the position of the first shed vortex, and may therefore be computed in advance. _{In particular, the} results

for

= and .x, = 1 + .25w are =

046667N3 + 1.27273 ''N-2 +

0.24.242

-

0.04848

1N+2

(19) 1N -0.40000 + O.8l8l8y 2 + 0.72727

_1N+l

014545N+2

The strengths of the two wake vortices appearing, in (19) are, of course, not independent, but related to the.past _{and present} circulation as indicated in (12).. _{Since the present circulation is..}

unknovn, these shed vortices

_{are runctions of the N -2 unknown}

_vortex strengths, _y , and on the known total circulation around _{the foil} durina thepreceding two time steps. _{This can he expressed as follows}

U -4 1 = c1 _{+ C2YN} 3 + C3 n=l - 28 N - 2

+ F+

C5 F2)

5w

c

N-4 1N =

'fl + C'y

3 + C8Y -2 + - 1) + C1 2) (20) -n=l N- 4 + 1 ll n C12YN ₃ + C1SYN

-+ Cl4 r(k)

₊ C15 (k-2)

when the coefficients _{C1, c2, ... C15 are simple alrebraic combinations} of the weighting coefficients b1 For a cOtflPon uniform spacing _on the, foil and in the wake, these coefficients _are

C1, ... C5 - 0.12308 - 0.48308 0.89232 0:12308 - 0.02461

c2, ... c10 - 0.36923- 0.44923-0.32308 _{0.36923 - 0.07385 (21)}

C

...

_{c15. -0.S0769 - 0.06769-1.56924 0.50769 0.09846}

The linear dependence of

N - 1 N + 1 on the strengths of

the remaining vortices On the foil,

N

N-2

and the total

circula-tion around the foil during the preceding two time steps expressed in (?O) insures that both the Kutta _{cOndition and Kelvin's circulation}

theorem are satisfied at each time step.

In executing a particular time step, we have N - 2 unknown circulation strength

-' 2 - 2 to be determined by satisfying

strenqths of the last two vortices on the foil, as well as the vortex shed into the wake are unknown, but in acco-dáne

with (20),

_{these are} dependent on he 'N - 2 unknown votex strenths and on the known

total circulations obtained from the last two time steps.

.V. DISCRETE TIME STEP SOLUTION OF THE BOUNDARY VALUE PROBLEM For, the discrete vortex model eipioyed here, the governing integral equation (4) is reduced to a system of linear algebraic equations: N L _x -n=l 2 n k K 1 F(K_1 ) -x 2. w (k - K)O (22) where (k)

is the external vertical velocity at the cont'ol point located at x _{at the k'th time step. We can now rearrange (22) in}

order to put the unknown

l 2. N 2

on the left hand side, and all remaining known quantities on the right hand side. In doing this, it is necessary to use (20) to eliminate the unknown dependent vortices. The result is the following set of N-2 _simultaneous equations which must be satisfied at each of _{the control points,}

N-4

1'

i

h=1 -+ _.1 N

-2

30

-C1 C C.,1 + +

+.

-

-

-

1

2 -

'N

2 -

w.J 1 CT C12 I

_XNl Xj_NX_XwJ

-

(k) ,-(k-l) -27rv - N -+ C C, C +

J---

+.

U ₊

-X9-1

X-C9 C13 N )

I

c5 C10 C14 - 1 2. - +

-K=l x - -- (k-K) w

k'

-

F 1)

(23)

The computation is performed by starting from rest with the circulation set to zero. Since the simultaneOus equations in (23) must be solved for each time step, it is efficient to invert the coefficient matrix, which depends only on the numberof elements N.

The solution for the N-2 unknown vortex strengths may then be

obtained for each time step by a matrix multiplication with the current value of the right hand side. For the first time tep, only the first term on the right hand side, is present since no wake vortices have yet been shed. The'second term appears in the second time step, and

finally

the

_{remainino tio terms appear starting with the third time} step.

It should be noted that the _{number of eleuents in the last.} summation on the rinht-hand side continuously qrows as the motThn proGresses. However, _{s the shed vortices pass far downstrem}

their influence becomes neglibile, _{so that it is not necessary to}

remember" the ent.ire wake back to

_K=O

. This can be accomplished by

setting the, lower limit.to a suitable _{valu of K > 1}

Finally, it should be noted _{that the particular problem to}

be

solveø is defined entirely by _{, which can them be made to} represent any transient or oscillatory _{motion. In the later case,} since the, mot{on must start _{from rest, it is necesa'y to}

carry out the solution for a few cycles _{in order to allow the starting}

transient to die out

VI.. _{FORCE CALCULATION}

For the. unsteady motion of _{an ideal fluid the linearized} Bernoulli equation can be, written _{in the form}

pU2ppuup

.

(24)

where p and U _{are the uniform conditions at infinity,}

u is the perturbation velocity in the _{x direction and is the velocity}

potential.

Applying (24) to one point on the upper surface and _{one point} on the lowersurface of the hydrofoil _{and subtracting gives}

an

32

-:n(x)

- Pu : U (u - UL) +

at - (25)

The jump in horizontal _{velocity is the}

negative of the local _vortex sheet strength, _y(x)

; while the jump in potential is the _inteQral from the leading _{edge to the point}

x _{of the vortex sheet}

strength AP(x)

_-pUy(x)

while the total lift _{is obtained by}

integratingthe _{pressure jump} over the chord

fAP(x)dx

-pur-

Jf

y()ddx

0

₀₀

C _cx

(27)

The first term in (27) _{is the same as the lift due to}

cirulatjon

in steady flow, while

the. second term appears only for unsteady flows. It is unfortunately

impossible to give _{unian'ibigous names to these} two

components, since, in _{particu1a, the separation}

of total lift into a circulatory and added _{nass component is riot always}

the same. _{It is} therefore safer to _{consider only the total}

lift, which, hopefully, must be the same _{regardless of the formulation}

of the theory! The discrete vortex

representation of (27) is _{,e.latively simple.} The first term is already known since

is computed at each time step.. _{The double' integral}

in the second _{term is replaced by} a

double summation of _{discrete Vortex elements}

(28) m=i n=l

first-order differentiation formula

dF F(k) -

F<

(29)

It has been found that (29) works quite well, except that it results in. a phase error of approximately one-half of a time step. A pQssible explanation is that the derivative should be evaluated _at the end of the discrete time interval, rather than at the middle, which is what (29) implies. This suggests that the use of a higher order differentiation formula set to evaluate the derivative at the.end of the interval. The following five point formula has been tried,

dt 12&t

48F'

+ 36

F2,

and has been. found to work very well. Obviously, it is necessary t use (29) fOr the first four time steps, at which point sufficient /

past history is available to permit switchinci to (30).

VII. TWO.0IMENSIONAL. RESULTS

The numerical procedure. was tested for two _{cases whose solutions} are well known; simple harmonic hea,ving motion and the sinusoidal gust. The first is the simplest possible unsteady motion, and was therefore used in. the first trials of the program. The second is

somewhat closer to the problem of a propeller operating in _{an non-uniform} wake and was therefore considered worthwhile to include as a further

-

314

-test of .the program. _{In ordr- to compare the numerical}

resuLt, with the krown analytical solution, _{it is necessary, to obtain}

the solution for a sufficient number of tire _{steps for a steady-state oscillation} to.&e reached. _{It was found that three cycles, of oscillation was} sufficient for this purpose. _{The final amplitude and}

phase for the numerical solution was then found _{from a discrete harmonic}

analysis of the third cycle of the _{motion starting from rest,.}

The exact analytical solution for both problems was programmed to permit accurate calcula-tion of the errors in amplitude _{and phase of the numerical}

solution. Table (1) shows the results of a convergence test with the

numbe, _{of elements over the chord, N}

, varying from 10 to 40, for

a

constant reduced frequency, k , of

rr/2.

The spacing in the wake

was

set equal to the spacing of the _{foil, Sb that the number of}

time steps per cycle, in this _{case, was 2N . The errors' in. amplitude}

ad

phase decrease smoothly with increasing N , and it appears that

sat-isfactory accuracy can be achieved _{with 20 elements. This is} approximately the same as the number of elements needed for _accurate computation of the steady loading due _{to camber for an NACA a= .8} mean line. ,

Figures (3) _{- (10) are 'plots of the distribution of vorticity} over the foil and wake at _{one particular time step of the last} cycle, chosen t,o correspond _{to the time of maximum oscillatory}

'lift coefficien,t. _{It is interesting 'to note tha't the sharp discontinuity.}

in slope of y(x) _{at the tYailing edae for the case Of heaving moti".} and the nearly continuous' _{slope 'in the case of thesinusoided}

-J

NUMBER OF

AMP LI TUDE

PHASE

E R!ROR AMPL IT UDE

AMPLITUDE

ELEMENTS

(ANALYT)

(NUM.)

(ANA LVI.)

(NUM.)

C 0 N V E R&EN CE TE ST

RE CUC ED FRE QU E N CV k

1. 571

10.

20

30

40

5.551

5.266

5.551

5.423

5.551

5 . 468 5 .

551

5 4910

L JO

1 . 5%

1.1%

125.0

126 .0

126 . 0

126.0

127.8

127.0

126 .7

126 .5

.1 8°

1.0°

0.7°

0.5°

1.978

1.978

1.978

1.978

2.074

2.028

2.013

2.005

4

ERROR

A 00 (.) 0 ) 0 1

L.U0

0.0 1 fl9.

.0

In

PHASE

(A NA LVI )

-13.9.4

-139 .4

-139.4

-139 . 4

Zj

PHASE

_(NUM.)

-133.3

-136.9

-137.7

-138. 1

ERROR

5.12

2.50.

1.7°

1 . 3

(N!UM)

A;MPL I TUD:E

ERROR.

PHASE

(A NA LVI.)

VORTICII.TY DISTRIBUTION ALONG CHORD AND

IN WAKE

HEAVING MOTION

0 C

REDUCED FREQUENCY

k Tt12 NUMBER OF ELEMENTS. ALONG CHORD N 10

FIG.

I

REDUCED FREQUENCY

k = it /2

NUMBER OF ELEMENTS N

2,0

ALONG CHORD

-1100

REDUCED FREQUENCY k

itI2

'NUMBER OF ELEMENTS N = 30 ALONG CHORD -1.25 -1.50

HEAVING MOTION

'FIG. 5

-Q 25

-O5O.

REDUCED FREQUENCY

k tI2

NWMBER OF.ELEMENTS _N 40

ALONG CHORD

0.50 0:25

000

0

-0.25

- 0.50 -0.75 -1.00 - 1.25

-1.50

VORTICITY. DISTRIBUTION ALONG CHORD AND IN WAKE

SINUSODIAL GUST

_:

REDUCED, FREQUENCY

k = TtI2

NUMBER OF, ELEMENTS N 10

ALONG CHORD

0

/

0

Q5:O

0.25

o,o o

-075

-1.00.

-1.25

-1.50

REDUCED FREQ(!.LENCY

NUMBER OF ELEMENTS

ALONG CHORD k

N =

Tt t 2

20

0.50

-

-0.75--1.00

-1.25

-1.50

SINUSODIAL GUST

REDUCED FREQUENCY

NUMBER OF ELEMENTS

ALONG CHORD

k :Tt/2

N

30

FIG. 9

0.50

0.25

0..o0

-025

x

-0.5Ô

.

a U. - a a 'a a: a U a U U

.

_. U - a U a - U. U am.

I.

U 'a a. N.

FIG. 10

-0.75-L

REDUCED FREQUENCY

k z

it /2

NUMBER OF ELEMENTS - N

40

-1.25

ALONG CH:ORD-- 15o.

-

114

-result f!orn the identical numerical procedures. _{Convergence of the} shape of y(x) over the chord with increasing N is evident from the plots.

Calculations are also shown for reduced frequencies varying from

to

ir/6

, with a fixed valued of N= 20. The spacing.in the wake and on the foil was again et to be equal, so that as a result, the

number, of time steps per cycle increase as k . A summary of the numerical results in comparison with the analytical solution is given in table (2). _{As one might expect, the accuracy i.s best at the lower} reduced frequency. However, the results are sufficiently accurate over the complete range of principal interest for marine propellers. These. results are plotted in Figure (Il) and (12).

Figures (13)-(24) shows the vorticity distribution 'as a function of reduced frequency. As must be the case, the heave and gust results become the same as the reduced frequency is decreased.

As a final test, the transient problem of a sudden imposition of an angle of attack at t=0 was calculated to enable comparison with the Wagner Function. _{The resulting vorticity distribution i.s} plotted in Figure. (25) after six chord lengths of travel,. _{Here one} can see the distribution of vorticity.over the chord to be almost the steady state "flat plate" distribution with a small remnant of

vorticity in the wake. In this case,.the lift obtained numerically has reached.85% of 2rr , which is within 3% of the Wagner results..

12

11

10

3

NUMBER ELEMENTS. ALONG CHORD.JN 20

ö

HEAVE

NUMERICAL.

GUST 5

SOLUTION

ANALYTIC SOLUTION

H .:

HEAVEA

A

Ar

-U

S.

1.

.2

180 170 160 150 a-120 110 -150 140 -130

- 146

-NUMBER ELEMENTS ALONG CHORDJN 20

o HEAVE

NUMERICAL GUST SOLUTION GU ANALYTIC SOLUTION ST

wO

-J

4

h

AMPLITUDE (.ANALYT.)

AMPLITUDE (NUM..).

.

ERROR

.

PHASE

(A.NALYT.)

PHASE

(NUM.)

ERROR .

AMPLITUDE (ANALYT.)

AMPLITUDE (NUM.)

. E.RRO.R

PHASE

.

(ANALYT.)

PHASE

(NU.M:) ERROR ,

AMPLITUDE

AND.

PHASE. ANGLE

VERSUS REDUCED FREQUENCY

3.14.2

1.571

.1.04.7

.785

.

.628

.524

10.140

:555i

4.308

3.894

3.782

-

3.793

9.751

5.423

4.247

3.85.3

3.744

3.75.3

4.0%

2.3%

1.4%

1.1%

i0%

1.1%

10.8.3

.

126.0

1.41.5

153.7

162.6

169.1

108.9

127.0

142.3

154.1

16:2.8.

169.1

0.6°

1.0°

0.8°

0.4°

0.39

1.410

1.9.78

.

2.396

2.731

3.009

3.247

.

1.346

. ..

2.028

2.476.

2.795

.

3.038

3.254

.4.7%

.

2.5%

.

3.3%

2.0%

1.0%

0.2%

-137.2

-139.4

-141.3

.

-143.1

-144.7

-146.2

-141.8

-136.9

-141.3

-143.9

-145.8

-147.4

4,6°

2.5

O°

0.8°

1.1°

_1.2°

VORTICITY DISTRIBUTION ALONG CHORD AND IN THE. WAKE

HE AV.I NG

MOlt ON

REDUCED FREQUENCY

k.=Tt

NUMBER OF ELEMENTS N = 20

REDUCED FREQUENCY k

it/ 2

NUMBER OF ELEMENTS N 20

ALONG CHORD

HE.AV ING M:OT JO NI

REDUCED FREQUENCY

k =it/3

NUMBER OF ELEMENTS

N

2O

ALONG CHORD

-0.25 I

-

_I

P

0

I

I REDUCED FREQUENCY

k

it 14

NUMBER OF ELEMENTS N = 20 ALONG CHORD

FIG. 16

0.50.-'O.:25. 000.

HEAVING: MOTION

- 0.50.

I

-0.75--1.00- REDUCED. FREQUENCY

k =tt/5

NUMBER OF ELEMENTS N 20 ALONG CHORD

-1.25,

-1.50

_{FIG: 17}

O50 0.25 0.00 - 0.25 -0.50. -0.75 -1.90 -1.25 - 1.5.0 RED.UCEEFREQUENCY NUMBER OF ELEMENTS ALONG CHORD

k =

N = 20

F 1G. 18

0751oo

--t

25--1.50

VORTICITYDISTRIBU.TIONALONG CHORD AND IN THEWAKE

SINUSOOIAL GUST

REDUCED FREQUENCY

k = It

NUMBER OF ELEMENTS N

20

ALONG C;HORD

O50

0.25

-0.50

-0.7.5

-i00

-1.2.5

-1.50

REDU:C.ED FREQUENCY k

It /2

NUMBER OF ELEMENTS N = 20

ALONG CHORD

FIG. .20

SINUSODIAL GUST

REDUCED 'FREQUENCY

k.:

_{It /3}

NUMBER OF ELEMENTS N;

20

ALONG CHORD

0.25

0.50

o.Qo

-a 25

REDUCED FREQUENCY NUMBER OF ELEMENTS ALONG 'CHORD

k = lt/4

N = 20

FIG.. 22

0.50-0.25

0.00

SINUSODIAL GUST

REDUCED FREQUENCY k

it I 5

NUMBER OF ELEMENTS N

20

ALONG CHORD

FIG. 23

0.50

>(

0

:0;25

0.00

-0.25

-0.50

-0.75

-i.00

-1.25

-1.50-REDUCED FREQUENCY

k =

_{It I 6}

NUMBER OF ELEMENTS N

20

ALONG CHORD

FIG. 21.

0.50-'

025

too.

--1.:25

I

INST AlA NE OUS C:HAN:G E

IN

A NGLE OF ATTAC K

NUMBER OF ELEMENTS ALONG CHORD.N

20

VIII. THREE

_DIMSIONAL

_{RECTANGULAR HYDROFQtL} The present method developed

for two-dimensional heaving motion has been applied to _{a rectangular hydrofoil.}

This introduces the three-dimensionality into the prOblem _{but keeps the difficulties}

in

geometry at a reasonable level

The same basic vortex _{arrangement has been applied as in the.} two-dimensional case, with the additIon _{of discrete trail ing}

vortex elements over the.foil and wake.

This work was carried out at M.I.T. as part of a Master of Science thesis by C. S. Lee /4/, _{who made calculations for aspect ratios of} 2,4 and 20.; His results _{are shown in Figure (26), and are in good} areernent with several earlier _{published results.}

IX.

_{STEADY-STATE PROPELLER PERFORMANCE PREDICTION}

A series of studies _{on the numerical modeling of non-linear}

iiftihg surface theory conducated at M.I.T. _{serves as the basis for the}

calculation of unsteady propeller _{characteristics.} The mode collocationmethod _{was fjrs}

introduced for the steady-state calculatIon, as described _{in references /,' and}

/5/. However, there was a need to improve _{the. accuracy of the}

program in order to permit computation of detailed _{pressure distribjtion in steady flow,} as well as to make possible the planned

extension to unsteady flow. Attempts to do this by increasing the number of modes _{were Unsuccessful.}

Encouraged by the stable behavior _{of the "modeless" approach}

in

our study of two-dimensional

15

62

-FIG. 26

4

AR00aTHE0DQRSEN

PRESENT PRESENT

METHOD

METHOD

£

AR :20

- - -

AR:

LAWRENCE GERBER(EXPERIMENTS)

AP: 2

Jo

_{PRESENT METHOD}

LAWRENCEGERBER (EXPERIMENTS)

,,>

.7

/

0.5

tO

2.0

as PI'NV.-4 was developed to calculate the steady-state prooeller perfOrm,anc.e in a more eact manner. The new procedure employs, a discrete vortex/source arrancienent as shown in Figure (27).

The spanwise and chordwise sinciularities can be spaced either with a constant physical spacinq, ora constant Glauert anale. All the.

singularities are situated on the exact blade camber surface, thus accounting for the geometrical influence of pitch distribution, skew and rake, control points are iocated at the geometric center of the vortex panels where the Cauchy-principle value integration is closely approximated. Termination of circulations at. the tip, root, and the trailing edge causes t.he strength of vortices at'. those edges to be functions of their adjacent vortices. The slip stream is modeled by the transition and ultimate wake geometry as presented in. reference

/6/. . .

The vortices are arranged by horseshoes with unknown strength tb bedeterrnined by satisfying the boundary condition at all control

points.. The.boundary condition, requires that theorthogOnality exists between the, local velocity vector and the local surface normal vectOr, thus including the influence of radial velOcities. Forces are

calculated. by suniming.up the forces acting on each discrete vortex Where the Kutta-Joukowsky law is applied.

UNSTEADY PROPELLER BLADE LOADING

The development ofa non-linear unsteady lifting surface program, designated as.PUF-1, isnow underway at MIT in order: to predict the

614

VORTEX ARRANGEMENT

(COS-SPACING)

(COMPUTER PLOT)

i0

1/

\iiuii_

VORTEX SHEDDING AND UPDATING PROCESS

66

-unsteady iQading on propeller blades. _{This program is essentially}

a

combined effort of the PINV-4 numerical _{lifting surface procedure} adapted with non-unifoym i.nf1o: and the unsteady two-dimensional numerical scheme previously described. _{In order to facilitate the} bookkeeping process of the vortices, we link a pair of spanwise vortices with the corresponding pair of chordwise vortices to form a vor.tex quadrilateral which is then equivalent to a three-dimensional dipole of uniform density withi'n the element. Thisarrangement

has the addditional advantage that. the _{Kelvin's.conditión is} auto-matically satisfied by itself.

The, procedure consists of _{a succession of solutions of the} propeller boundary value 'problem with the blade rotation, angle

incremented by a prescribed amount. _{Since the loading on.each blade} is different in the unsteady case, _{some modification in the steady} flow.procedure is required. _{Rather than solving the eqUation.s in}

all

blades simultaneously, we _{use an iterative scheme,, as illustrated in} Figure (28).

Initially all blades are assumed unloaded. _{The boundary value} problem is then solved for the unknown loading _{only on the first blade} without the effect of the other blades, _{as shown in Figure .(28á)}

The propeller is then rotated to the _{next angular interval, resulting in} a newloading on the first blade, _{as shown in FigUre (28h).}

' The vortex

quadrilateral shed immediately behind _{the trailing edge at thi.s time} step has a known strength equal to the blade circulation at the previous time step.... The vortex quadrilaterals _{shed in the wake during earlier} time steps will be convected _{downstream in the transition wake .region,}

as Shown in Figures.

(28c) through _(28e)

, _{and ignored in}

the ultimate

wake regfp. _{As the}

So1utio _{progresse, other blades}

will begin to

occupy angular

positions _{Previously occupied}

by the first _{blade, as} s.hown in _{Figure (28f)}

throuch (28h).

Their loading _{will then be} automatically _{updated, so that}

afer a few

complete revolutions _a Solution will be

obtained with all blades _{having the}

correct loading. Thj _{scheme shows}

rapid

convergence from our test _programs. For a three_blade

propeller, the

convergence occurs _{at the third} revolution, figure 29. _{The overall}

computation can further _be economized by _{selecting a}

course spaced time _{step for the}

first few revolutions until

convergence, and then _{using a}

interpolation technique the loading _{at intermediate}

finer time

steps can be obtained. Thi provides a new starting _{point for}

additional _{revolution with}

finer spacing. _{Forces on each}

vortex are

Computed as 'in the steady _progrm, with the

addition of the _term

proportional to the _{time rate of}

change of the.velocity

potential. Summation over _{the entire}

blade provides all six _{components of,} force onari

individual blade. _{These may then}

be resolved, _{if desired,} by discrete _harmonic

analysis to _{provide mean' and}

CONVERGING PROCESS FOR UNSTEADY CALCULATIONS FOR A THREEBLADED PROPELLER WITH PRE

-DOMINANTLY THIRD HARMONIC WAKE

(NUMBER. OF REVOLUTIONS)

1 2 I

3

I

.09

I'

11

(1

I\

1

08

Ii

It

It

II

I

I'

it

.O7 1

Il

Ii

it

lIt

.1

I

I

1

It

-j

I

I

I

I

/

I

I

I

i

'

IL

w.06

It

I I

I

I

I

I

Il

Lii I

I

I I

t

I

1 I

I

I-i

04

III

III!

till

lit1

03111

I

111111!

Ill

11

I

I'

fJ

IJ

ti

I I . I

.02

I

.1

I I

.01

I I I

0

120

240

0

120

240

0

120

240

0

120

240

360

ANGULAR POSITION

MEAN LOADING

0.55

FIG. 29

P.E FE RE E S

1Kerwin, J. E.,. 'tComputer Techniques for Propeller Blade _{Section Design,"}

International Shipbl4ildin

'grsss, Volume 20, No. 227, July 1973. 6Tsao, S. K., "Documentation of Programs for the Analysis of Performance

and Spindle Torque of Controllable Pitch Propellers, MIT, Department of Ocean Engineering, Report 75-8, _{May 1975.}

3James _{R. M., "On the Remarkable Accuracy of the Vortex Lattice Method,"}

Corrrputer

_{Methods in Applied !'echanics and Engineering 1, 1972,}

pp. 59-79.

4Lee, _{C. S., "A Numerical Method for the Solution of the Unsteady Lifting} Problem of Rectangular and Elliptic Hydrofoils," _{S.M. Thesis, M.I.T.} Department of Ocean Engineering, February _1977.

5Cummings, _{0. E., "Numerical Prediction of Propeller Characteristics,"} Journal of Ship Research, Volume 17, No. 1, March _1973.

6Kerwin, _{J. E., "A Deformed Wake Model for Marine Propellers," M.I.T.} Department of Ocean Engineering, Report _{76-6 , October 1976.}