A surface panel method for analysis of open and ducted propellers

Jan SZANTYR *

A SURFACE PANEL METHOD FOR ANALYSIS OF OPEN AND DUCTED PROPELLERS

The paper presents a newly developed computer method for analysis of open and ducted propellers. This method is based on surface distribution of dipoles on propeller blades, on the propeller hub and on the duct. The method employs unsteady potential formulation of the boundary condition on the propulsor. For the propulsor of given geometry operating in given three-dimensional non-uniform inflow velocity field the method predicts unsteady pressure distribution and hydrodynamic forces on impeller and duct, describes sheet, bubble and vortex cavitation phenomena. and evaluates pressure pulses generated by the cavitating propulsor in the surrounding space. The resulting computer program may be used on UNIX system workstations. The paper discusses the results of initial experimental verification of the method, which include the cases of an isolated and ducted

propeller.

i. INTRODUCTION

The surface panel method for analysis of unsteady operation of open and ducted

propellers described in this paper has been developed in order to supersede several earlier methods based on lifting surface theory [2,7]. Those methods have been introduced about

15 years ago and they proved to be very efficient and relatively reliable tools for ship designers and hydrodynamicists. However, despite continuous development it became obvious that their potential for improvement is nearly exhausted. Simultaneously, wide

application of more powerftil computers has created the possibility for practical application

of more complicated algorithms and more accurate discrete representation of propeller

geometry and flow.

The main advantage of the surface panel methods over classical lifting surface

approach lies in locating mathematical singularities representing the propeller on the actual

surface of the blades and other elements of the propulsor. This enables more accurate modelling of flow, especially in such crucial regions as leading edges and blade tips. The

importance of this feature for prediction of cavitation phenomena is unquestionable. On the other hand surface panel approach requires twice the number of elements to achieve

the same effective density of elements as the lifting surface. This means quadrupling the

size of matrices and increasing time of their inversion by the factor of 8. One _{must also}

remember that both methods calculate potential flow solutions, which require meaningful

empirical corrections for the effect of viscosity. The magnitude and nature of these

corrections near the leading edge remains uncertain and they may drown the advantage of the surface panel approach even before it comes to the analysis of cavitation.

The subsequent paragraphs of this paper present the theoretical model of the

method, details of the discrete representation of the propulsor, structure of the algorithm

and of the computer program and results of initial analyses of_{test propellers confronted}

with the results of model experiments. The entire research work described in the paper is summarised in conclusions.

TECHNISCHE UN1VERS1TE1T

ScheepahYc1r1leci ca

Archief

Mekelweg 2, 2628 CD

Delf t

Tel: Ol52786873/Fa278l836

Key words: Ship hydromechanics. Propulsors. Unsteady flow.. Cavitation. Numerical analysis.

2. THEORETICAL MODEL

Surface panel methods are in existence for several decades as tools for solving flows around three-dimensional bodies of complicated geometry [1,3]. They may be formulated in several different ways, namely as velocity or potential methods. Velocity methods may be based on distribution of vortices, sources, dipoles or combination of

dipoles and sources. Potential approach may be realised either as perturbation potential or as total potential method. Numerical implementation may in turn be realised as low order or high order method, using selected panel geometry (quadrilateral, triangular etc.). Each

specific formulation has its advantages and disadvantages. The method selected for the

analysis of propulsor flow in this paper may be defined as first order unsteady perturbation potential method employing dipoles located on quadrilateral panels.

The main assumptions of the theoretical model [6] are as follows: -flow around the propulsor is incompressible, inviscid and irrotational;

-flow may be described by the perturbation velocity potential which fulfils the Laplace

equation;

-kinematic boundary condition on the surface of the propulsor S, postulates no flow

through the rigid surface:

n (1)

where: n - unit length normal vector at the propulsor surface U- velocity of the undisturbed flow.

-the wake surface S behind the propulsor is infinitely thin with zero jump of pressure and

normal velocity, but allowing non-zero jump in potential; for lifting unsteady flows this

jump of potential follows the time variation of dipole intensity on the propulsor;

-the Kutta condition, formulated in the most general way, requires the flow velocity to be finite at the trailing edges of the propulsor lifting elements.

Now the flow around the propulsor may be described by the following integral

equation, assuming that potential of the fictitious flow inside the propulsor is equal zero on the propulsor surface:

2ïrq(p)= _j sp i ênR -(2)

where:p - location of the control point q - location of the dipole

R - distance between p and q.

As the values of /â2 are known from the kinematic boundary condition on the propulsor

surface (1), formula (2) is the Fredholm equation of the second kind for the unknown

intensity of dipoles modelling the flow. Through appropriate discretization of the propulsor this equation may be transformed into the system of linear equations with unknown dipole

intensities, which are equal to the values of potential in selected points on the surface of

modelling the propulsor using the Kutta condition. Special care must be taken in numerical evaluation of the integrals in equation (2), in order to ensure acceptable accuracy [5].

After solving the above mentioned system of linear equations the pressure

distribution on the propulsor may be obtained from the formula based on the Bernoulli

equation for unsteady flows:

Cp

_{U 2}

_Uâ

(3)

2

where V is the resultant flow velocity according to the formula:

= U

-

(n. v)n +

-

(.

] +

-

(.

(4)

pxr

where: _p - pressure in the undisturbed flow

- unit length tangential vectors Q - angular velocity of propeller rotation.

Pressure distribution obtained from (3) is the basis for cavitation analysis and

subsequently for calculation of hydrodynamic forces on the propulsor. The method can

detect and analyse three forms of cavitation, namely: -sheet cavitation

-bubble cavitation -tip vortex cavitation.

This analysis is conducted along the lines similar to earlier methods and described in detail in [2,7]. It should only be remembered that sheet cavity results in apparent distortion of the

blade and duct geometry, thus influencing the kinematic boundary condition on the

propulsor. Hydrodynamic forces on the propulsor are calculated through integration of

pressure over propulsor surface. The effect of sheet cavitation on the pressure distribution is taken into account in this calculation.

Apart from that the fluctuating pressure induced by the propulsor in the surrounding space is also calculated, including the effects of:

-variable hydrodynamic loading of the propulsor -finite thickness of rotating elements of the propulsor -steady and unsteady cavitation phenomena.

3. DISCRETIZATION OF THE PROPULSOR

Marine propulsor under consideration may be divided into three basic elements:

propeller blades (usually 3-7), propeller hub and duct. Only the blades are always present

in the analysis. two remaining elements may be regarded as optional. All elements are

represented in numerical calculations by a number of flat quadrilateral panels.

Correct discretization of the propulsor is the result of a compromise between the

Fig. i Discrete model of the propeller blade Fig. 2 Discrete model of the hub

-maximum number of panels results from the size and speed of available computers and

from the required computation time

-appropriate representation of the details of propulsor geometry

-limited allowable difference in size of adjacent panels (usually no more than 25%) in

order to ensure numerical stability of calculations.

-

---GLADE

After considering all requirements, the following density of discretization was suggested as standard: 640 panels on each propeller blade, 320 panels on the section of

propeller hub located between two adjacent blades, 400 panels on the section of the duct

corresponding to one propeller blade. In modelling the blades uniform distribution of

panels both in radial and chordwise directions was assumed (cf. Fig.1). There are 16 rows

along the radius with each row having 40 panels in the chordwise direction (20 on each side of the blade). The section of the hub is modelled by 8 rows in the direction of flow

with 40 panels in each row (cf. Fig.2). The hub may be alternatively modelled as a cylinder with half-ellipsoids on both ends or as a cylinder with ellipsoid on one end and one end left

open. The section of the duct corresponding to one blade is modelled by 10 rows in the

direction of flow with each row having 40 panels (20 on each side of the duct), what

shown schematically in Fig.3. For each panel the following geometrical quantities are defined:

-co-ordinates of the central point, in which the boundary condition is controlled

-components of three unit length orthogonal vectors, one normal and two tangent to the

surface -panel area

Apart from panels modelling the propulsor, another system of panels is employed

for modelling the wake behind the blades and the duct. The wake behind propeller blades

is constructed of regular helical surfaces (cf. Fig.4), divided into variable zone reflecting changes in hydrodynamic loading of the blades during operation in non-uniform inflow velocity field and steady zone defined by the circumferentially average flow condition.

Altogether variable zone contains 160 panels covering the angle of 360/Z [deg], while the

steady zone covers four complete revolutions of the helix with 768 panels. The wake

behind the duct consists of the variable zone only (see Fig. 5), containing 100 panels. The distribution of the dipole intensity in both variable zones is a function of time.

DIREC11OW OF ROTATION

Fig. 5 Free vortex system behind duct sector

For the sake of illustration it should be pointed out that a complete four bladed ducted propeller with hub would be modelled by 5440 panels on the propulsor and 4112

panels in the wake. Such number of panels would be prohibitive from the point of view of computation time. This problem is circumvented by conducting numerical solution for one

blade only, called the key blade. Other blades are included only through their steady.

average flow effect. In this way the maximum size of the system of linear equations

reflecting the kinematic boundary condition on the key blade is 1360, which is acceptable for a medium power workstation..

4. SEQUENCE OF CALCULATIONS

The method analyses one complete revolution of the propeller, during which a

number of specified angular positions of the blade in the nonuniform inflow velocity field

are taken into consideration. The analysis starts with solution of the kinematic boundary

condition on the propulsor for the circumferentially average inflow field. This produces the

distribution of dipole intensity which is applied to all blades except the key blade and to

steady zone of the wake. Initially it also describes the variable zones of the wake.

DUCT SECTOR.

The unsteady analysis starts with so called ,,warming up run", which covers half of

the propeller revolution. The main purpose of this run is to fill the unsteady zones of the

propeller

and duct wakes with

an appropriate distribution

of dipole

intensity. corresponding to the changing hydrodynamic loading of the key blade and duct sector

during operation in the non-uniform inflow field. Initially, both variable zones are

described by the dipole intensity distribution corresponding to the circumferentially

average inflow, but it is gradually shed downstream and superseded with distributions

corresponding to local inflow to the key blade. In the middle of the ,,warming up run" the

cavitation analysis is switched on in order to achieve converged cavitation and wake

picture at the start of the fuliy fledged analysis of one complete revolution of the propeller.

This analysis covers one complete revolution, usually with up to 72 uniformly spaced angular positions of the key blade. Analysis of one angular position includes the

following computations:

-interpolation of the inflow velocity field

-solution of the system of equations reflecting kinematic boundary condition -calculation of the pressure distribution on the propulsor

-updating of the dipole intensity distribution in the variable zone -cavitation inception analysis (sheet, bubble and vortex cavitation)

-solution of the system of equations describing dynamics of the sheet cavity -calculation of the apparent distortion of the blade and duct due to sheet cavity -calculation of the components of hydrodynamic loading on the key blade and duct -calculation of the pressure pulses induced in the surrounding space

When all

prescribed blade positions are analysed, the

final computation is

performed, which includes:

-calculation of the components of hydrodynamic loading on the whole propulsor -calculation of the pressure pulses induced by the whole propulsor

-harmonic analysis of the variable hydrodynamic forces and pressure pulses.

5. STRUCTURE OF THE COMPUTER PROGRAM

The computer program has been developed on the Silicon Graphics Indy R5000

workstation having 128 MB of memory. It is written in simple FORTRAN77 and it

consists of six relatively independent programs which communicate with one another through data files and which are executed under the command of a batch file. The block

diagram of the program is shown in Fig. 6. The six programs perform the following

functions:

PP 1- reads three input data files containing descriptions of the propeller geometry, duct

geometry and non-uniform velocity field; converts the propulsor geometry into a system of discrete panels (files PROP and DUCT); interpolates the inflow velocity field and converts it into Fourier series (file WAKE)

calculates induction factor matrices, which relate intensity of dipoles with potential generated in control points (file MATIND)

executes the solution of propulsor flow for the circumferentially average inflow executes the solution of propulsor flow for a selected key blade position in the

non-uniform inflow velocity field, including calculation of the blade and duct pressure

distribution (file PRESDT) together with pressure pulses generated by blade loading and thickness (file PULSDT)

performs the

analysis of cavitation phenomena (file CAVLIM); calculates

hydrodynamic forces on the key blade and duct (file FORMOM); calculates pressure

pulses induced by the cavitating key blade and duct (added to PULSDT)

calculates hydrodynamic forces on the whole propulsor and pressure pulses induced by the whole propulsor; performs harmonic analysis of these two quantities (file FINRES).

START

PP6

STOP

Fig. 6 Block diagram of the computer program

FORM 0M

FINRES

The typical computation time for a heavily cavitating four-bladed ducted propeller

including hub, analysed in 36 angular positions, is 13735 seconds (3.8 hours) on Iris Indigo

DATA I PROP

---£

DATA 2 PPI DUCT

-DATA 3 WAKE MATIN D PP2 P P3 K= I K=K-4-I PP4 _CAy PRESDT PP5 PULSDT

2 with R8000 processor, which is considered to be on the limit of acceptability from the

point of view of the efficient commercial application of the program.

6. RESULTS OF THE EXPERIMENTAL VERIFICATION

The program was initially tested on two examples of model propellers, for which

the experimental results and results of calculations using other methods were known. One

of them was an isolated propeller, while the other was a ducted system. The main

particulars of the isolated propeller Pl are given in Table 1, while similar information for

the ducted propeller P2 is given in Table 2. The duct was of simplified geometry,

developed in the Institute of Fluid Flow Machinery.

Table 1: Main particulars of the isolated propeller Pl

where: nR - relative radius

PfD - nondimensional pitch c - chord length

s - skewback ordinate

r - rake ordinate

t - maximum blade thickness m - maximum mean line camber

Table 2: Main particulars of the ducted propeller P2

Propeller diameter D [mm] 250.0

Hub diameter hD [mm] 67.5

Number of blades z [-] 4

Expanded blade area ratio AE/AO 0.548

nR P/D c[mm] s[mm] rk [mm] t [miri] m [mm] 0.30 0.7636 45.5 -2.5 0.0 7.30 1.91 0.40 0.8065 66.0 -5.0 0.0 5.75 1.73 0.50 0.8473 81.0 -3.9 0.0 4.55 1.40 0.60 0.8849 89.5 0.0 0.0 3.60 1.21 0.70 0.9000 91.5 6.5 0.0 2.82 1.09 0.80 0.9116 88.0 15.0 0.0 2.15 1.03 0.90 0.9162 74.0 22.0 0.0 1.55 0.99 0.95 0.9155 59.6 29.0 0.0 1.28 0.94 1.00 0.9141 1.0 34.8 0.0 1.00 0.0 Impeller diameter D [mm] 220.4 Hub diameter hD [mm] 53.5 Number of blades z[-] 3

Expanded area ratio AE/Ao 0.85

r/R P/D c[mm] s[mm] r, [mm] _t[mml m[mm] 0.25 0.9186 75.97 2.40 0.0 9.04 2.35 0.30 0.9966 87.08 2.88 0.5 8.20 2.49 0.40 1.0984 107.62 3.84 1.5 6.67 2.49 0.50 1.1463 125.86 4.81 1.9 5.31 2.37 0.60 1.1633 140.44 5.77 2.1 4.13 2.05 0.70 1.1613 150.62 6.73 1.9 3.11 1.54 0.80 1.1542 151.35 7.67 1.5 2.27 1.03 0.90 1.1264 136.77 8.65 0.9 1.60 0.39 0.95 1.1134 109.72 9.13 0.5 1.33 0.15 1.00 1.0984 0.0 9.61 0.0 1.10 0.0

The previously executed experimental program included open water tests in the model basin and tests in non-uniform velocity field in the cavitation tunnel. The latter consisted of cavitation observations and measurements of the pressure pulses induced in

the selected points on the dummy afterbody installed in the tunnel. The calculations

conducted with the new computer program aimed at reproducing selected experimental results. Simultaneously, analogical calculations were performed using earlier computer

programs based on the lifting surface theory, namely [7] in case of the open propeller and

[2] in case of the ducted propeller. The results of comparison of experimental (EXP). calculated by surface panel method (SPM) and by lifting surface (LS) are shown in the

following figures. 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 EXP LS SPM

Fig. 8 Comparison of observed and calculated cavitation on propeller Pl (S0, o1.92, KTO.275)

H11

.LS

V

sP IOK -0 01 02 03 04 05 06 07 0.8 09 11

Fig. 7 Calculated and measured hydrodynamic characteristics of the isolated propeller fl

Figure 7 shows the comparison and measured hydrodynamic characteristics of the

isolated propeller Pl, while in Fig.8 the comparison of calculated and experimentally

observed cavitation pattern for a selected key blade position of this propeller in

non-uniform inflow velocity field is shown. Fig.9 presents the calculated and measured

pressure pulses generated by this propeller in the non-uniform velocity field in the point

located directly above it on the dummy afterbody. In case of the propeller Pl calculations were performed without hub.

0.7 0.6 0.5 0.4 0.3 02 0.1 O -0.1 700 600 500 400

300

-200 100

-

0-hi

I

h2 J D-EXP. -LS -SRM

I

_h3

Fig. 9 Comparison of calculated and measured pressure pulses generated by the propeller Pl

Fig. 10 Calculated and measured hydrodynamic characteristics of the ducted propeller P2

In case of the ducted propeller P2 the calculations using new method include the

hub effect. Figure 10 shows the comparison of computed and measured open water

diaram, while Fig. 11 presents the calculated and experimentally observed picture of

cavitation phenomena on the impeller for one selected angular position of the key blade in the non-uniform inflow velocity field. Finally, Fig.12 compares the unsteady shaft forces on the ducted impeller P2 calculated by the lifting surface and surface panel methods in the same non-uniform field.

----S..

-. 1OKQ XP.

Iíì

_

X 5PM - K1

.uUTD___

J 01 02 03 04 05 06 07 08 09

3 2.5 2 0.5 -o xp t-s

-*

-i' SPM

Fig. 11 Comparison of calculated and observed cavitation extent on the propeller P2 (9=6O, o=l.49, K1-=O.247)

O-LS -SPM

EI EI

Fx Fy Fz Mx My Mz

Fig. 12 Comparison of shaft forces on propeller P2 calculated by the lifting surface and surface panel methods 7. CONCLUSIONS

The following conclusions may be drawn from the analysis of the above presented first

results of experimental verification of the new method:

-in general the accuracy of the new surface panel method seems to be comparable to the earlier lifting surface methods; this may of course improve in the process of program

development and further verification;

-the new method can relatively well reproduce the steady hydrodynamic characteristics of the open and ducted propeller, although the duct thrust is less accurately predicted; visible

differences between calculations and measurements in both lifting surface and surface panel approach may be attributed to the simplifications in the wake model; rather poorer

than expected prediction of the duct thrust by the surface panel method results from

unfavourable interaction between discrete grids modelling duct and impeller;

-surface panel method has demonstrated certain improvement over lifting surface in

prediction of cavitation extent and pressure pulses; this results probably from better

prediction of pressure distribution near the leading edge of propeller blades;

-surface panel method requires about 10 times more CPU time than

a lifting surface procedure for an equivalent calculation; further tests with the new method will show

whether this is a cost-effective increase.

ACKNOWLEDGEMENT

The work described in this paper has been supported financially by the Polish Scientific Research Committee Grant

No. 9 T12C 025 08 executed in the Ship Design and Research Centre (CTO) in Gdansk.

REFERENCES

[11 CAPONNETTO M., BRJZZOLARA S., Theory and Experimental Validation of a Surface Panel Method for the Analysis of Cavitating Propellers in Steady Flow, Proc. PROPCAV'95, Newcastle-upon-Tyne, May 1995,_pp.

239-252

GLOVER E.J., SZANTYR J., The Analysis of Unsteady Propeller Cavitation and f-lull Surface Pressures for Ducred Propellers, Trans. RINA Vol. 132, 1990, pp. 65-78

KERWIN J.E. et al, A Surface Panel Method for the Hydrodynamic Analysis of Ducted Propellers, Trans.

SNAME, Vo!. 95, 1987, pp. 92-122

KOYAMA K.. Comparative Calculations of Propellers by Surface Panel Methods. Papers of the Ship

Research Institute, No.15/93, Tokyo 1993

NE WMAN J .N.. Distribution of Sources and Normal Dipoles over a Quadrilateral Panel, burn, of Engineering

Mathematics, Vo!. 20, 1986. pp. 113-126

SZANTYR J., Podstav.y teoretyc:ne i algorytm programu do wyznac:ania :miennvch sil i n1omentów

lozyskowych, kawitacji i pulsacji cdnieñ dia sruby i zespolu dysza-.ruba. .., (in Polish), Report of the IFFM No.

45 1/96

SZANTYR J., A Method for Analysis of Cavitating Marine Propellers in Non-uniform Flow, Intern. Shipbuilding Progress, Vol. 41, No. 427, 1994, pp. 223-242

METODA PANELOWA DO ANALIZY SRUB ODKRYTYCH I ZESPOLÓW DYSZA-RUBA

Referat przedstawia niedawno opracowan4. metod analizy rub odkrytych i zespoów dysza-ruba. Metoda ta jest oparta na powierzchniowym rozkkadzie dipoli na skrzydlach ruby, najej piacie i na dyszy. Metoda ykorzystuje niestacjonarne sformulowanie warunku brzegowego na pçdniku, oparte na potencjale prçdkoci. DIa pçdnika o znanej

geometrii pracujcego w znanym niejednorodnym trójwymiarowym polu prdkoci metoda dostarcza prognozv

niestacjonarnego rozkladu cinienia oraz sil hydrodynamicznych na rubie i na dyszy, opisuje kawitacj 1aminarna,

pcherzykow i wirow oraz wyzrìacza ¿mienne pole cinienia generowane przez kawitujcy pçdnik w otaczajacej przestrzeni. Oparty na metodzie program komputerowy moze by± uzytkowany na stacjach roboczych pod systemem Unix. Referat prezentuje ponadto wyniki wstçpnej weryfikacji eksperymentalnej metody, które obejmuj przyklad