Numerical Time Dependent Sheet Cavitation Simulations using a Higher Order Panel Method

Numerical Time Dependent

Sheet. Cavitat ion Simulations

using

a Higher Order Panel Method

Numerical Time Dependent

Sheet Cavitation Simulations

using

a Higher Order Panel Method

Numerical Time Dependent

Sheet Cavitation Simulations

using

a Higher Order Panel Method

P ROE FS C H RIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Deift,

op gezag van de Rector Magnificus Prof.ir. K.F. Wakker, in het openbaar te verdedigen ten overstaan van een commissie

door het College van Dekanen aangewezen,

op dinsdag 15 maart 1994, te 14:00 uur

door

Hendrik Jan de KONING GANS

scheepsbouwkundig ingenieur geboren te Amstelveen

Dit proefschrift is goedgekeurd door de promotor:

Prof. dr. ir. G. Kuiper

Published and distributed by: Deift University Press Stevinweg i 2628 CN Delft The Netherlands Telephone +31 15 783254 Fax +31. 15 781661 ISBN 90-6275-965-3 / CIP

Copyright ©1994 by H.J. de Koning Gans' All right reserved

No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without permission from the publisher: Delft University Press, Stevinweg 1, 2628 CN Deift, the Netherlands.

Acknowledgment.

This research is carried out at the Ship Hydromechanics Laboratory of Delft

University of Technology, Faculty of Mechanical Engineering and Marine

Tech-nology, with support by MARIN (Maritime Research Institute Netherlands). The investigations were supervised by prof. dr. ir. G. Kuiper (MARIN). I would like to thank him, for his support, guidance and encouragement. The

discussions we had, were always interesting, instructive and constructive. He let me free how to carry out this reseach. Furthermore I admire his patience for

correcting the earlier reports and this thesis. I owe him a great deal. The author is also indebted to the following persons:

P.F. van Terwisga. (Student in f990 and now a member of Ship

Hy-dromechanics Laboratory) The discussions concerning cavitation Were

al-ways interesting and instructive. He has also carried out the tests in the cavitation tunnel. From (t)his observations a better understanding of the cavitation mechanism was obtained.

C.P. Poot. (member of Ship ilydromechanics Laboratory), who made the

figures.

The other members of the Ship Hydromechanics Laboratory for their sup-port, during the four years I haved worked on this thesis.

Abstract.

This. thesis déals with sheet cavitatión. The.investigation is aimed at profile design with respect to cavitation control. At present it is possible to predict the. shape of.cavities on an arbitrary two-dimensional profile. in stationary flows. To improve the modeling of cavitation, a new method is developed, which allows three-dimensional flows an d instationary cavitation

To compute the flow around an arbitrary profile, a higher ordèr

three-dimensiönal pane! method program has been developed. The main algorithm used in this program is based' on a special casèofGreen's. theorem, called "de Morino formulation" This computer program (flow program) can calculate the

potential on the body and the velocities at the surface of the body or in the

flow field. The calculations of the influence coefficients has been reconsidered and improved.

A theoretical method is developed for time simulation of unsteady sheet cavitation. The concept of this theory can be summarized as, follows. The

potential is calcúlated at the body. The potential and velocity at the cavity

sheet is calculated with a normal Taylor expansion. This expansion makes use of the divergence and rotation theorems of potential flows. Then the pressure at the sheet can be determined at every time step From this pressure the normal velocity of .the sheet can be determined and by integrating this velocity the new cavity thickness and geometry can be computed. The characteristics of cavity grow.th and stability have been derived theoretically, explaining the behaviour of unsteady cavity.

Numerical simulations of the flow around profiles and of cavitation have been

carried out. The numerical results of the panel methods have been compared

with other calculations of the two-dimensiönal flow around profiles and of

three-dimensional flow around asphere and a wing. Sirnulationsof the growth of sheet cavitation .on a foil have also been carried oùt.

The conclusion is that higher order panel methods are more accurate than the zero order methods. Further refinement of the Kutta condition is required,

however.

The conclusion of thecavitation calculations is that the pressure is calculated very accurately and that the pressure is also correct at the closure of the cavaty

The cavity.growth is unstable, which is in agreement with the physics of a cavity

in potential flow However, the cavity grows predominantly in normal direction

Sarnenvatt ing..

Deze.dissertatie behandelt hetgedrag van viiescavitatie. Bet onderzoek heeft tot doel ontwepmethoden te ontwikkelen orn cavitatié te beheersen. In de huidige stand van de technologie is het mogelijk orn de vorm van viiescavitatie in een tweedimesionale 8troming te bepalen. Verbeteringen in het huidige model zijn mstationaire simulaties in een driedimensionale stroming.

Orn de omstroming van een wilekeurige lichaam te bepalen is een drie di-mensionale panelenprogramrnaontwikkéld. Het programmais gebaseerd op een speciaal geval van de stelling van Green, namelijk "De Morino vergelijkingen". Dit omstrorningsprogramma bepaalt de potentialen en de sneiheden zowel aan het opperviak van het te ornstromen lichaarn als in het veld. In deze dissertatie zijn ook nieuwe theoretische overwegingen en verbeteringen van de berekening

van de invloedscoeflIcienten opgenomen.

Een theoretische method is ontwikkeld voor tijdsafhankelijke simulatie van dynamische viiescavitatie. Het concept van deze methode is als voigt. De po-tentiaal wordt berekend: op het profielopperviak. Vanuit de potentiaai aan het profielóppervlak wordt de potentiaal en de sneiheden aan het vliesoppervlak door een Taylor-reeksontwikkeling loodrecht op het profiel opperviak Hierbij

is gebruik gemaakt van het feit dat potentiaalstromingen zowel divergentié- als rotatievrij zijn. De druk aan bet vliesoppervlak kan worden bepaaid uit de dan bekende grootheden. Bet verschil van de berekende druk met de

damp-druk geeft een maat voor de normaalsnelheden van het vliesoppervlak. Door deze normaalsnellieden te integreren verkrijgt men de nieuwe positie van het vliesoppervlak en de nieuwe dikte van bet vlies.

Ook is er theoretisch onderzoek verricht aan de karakteristieken van de groei van vliescavitatie.

Er zijn berekeningen aan de panelenmethode en cavitatiesimulatie uitge voerd. De numerieke resultaten zijn vergeleken met een driedimensionaie orn-stroming orn een bol en met tweedimensionale omorn-stromingen om een profiel. Tevens zijn er tijdssimuiaties van groeiende vliescaviatie uitgevoerd.

De conclusiè is dat hogere orde paneien methoden nauwkeuriger resultaten geven dan de nulde orde panelen methode. Verdere verfijningen aan de Kutta voorwaarde zijn echter vereist.

De conclusie van de cavitatie berekeningen zijn dat de drukken aan bet

vliesoppervlak uitstekend worden berekend. Tevens wordt de druk ann de

achterkant van de caviteit ook goed berekend De cavitatiegroei is instabiel, hetgeen ook fysisch in een potentiaal strorning optreedt. Echter de caviteit groeit hoofdzakelijk in normaal richting van het profiel, hetgeen fysish minder realistisch is.

Contents

i

Introduction,

_i

1.1 Outline of the Investigation 2

1.1.1 Panel Methods. 2

1.1.2 Calculations of Sheet Cavitation 3

1.2 Outline of the. Thesis. 4

1.3 References. 5

I

Panel methods

7

2 The integral equations.

10

2.1 Introduction. Io

2.2 The Fredhoim equations.

Il

2.2.1 A more general formulation of the integral equations. 12

2.2.2 The principle value integrals and jump relations. 13

2.2.3 The integral equation for a profile in a flow. 14

2.3

The surface at infinity ().

.

. Ï6

2.3.1 _{The fiat surface at infinity (íT) (Trefftz plane). 16}

2.4 The surface at the profile (ap) . 17

2.4.1 The Kutta boundary condition Ï8

2.4.2 Implicit Kutta (boundary) condition. 19

2.4.3 Explicit Kutta (boundary) condition. 20 2.4.4 Velocity on the profile surface. 20

2.5 Conclusion. 21

3

Discretization of the integral equations.

22

3.1 Vectorizing the source and dipole Íorrnulation 24

3.1.1 Implementation Kutta condition 25

4 Higher order panel methods based on calculations in the

refe-rence plane.

27

4.2 The order of the expansion 28

4.2.1 Truncation errors of a power series expansion. 29

4.3 The geometry description ₃₀

4.4 _{From curvilinear coordinates to reference plane coordinates. 30}

4.5 _{The transformation from global space to local space 33}

4.6 The panel surface. 34

4.7 Thè Jacobian 35

4.7.1 The Jacobian for the surface integral. . . ₃₅

4.7.2 The Jacobian for the development of curve in the tangent

plane. 36

4.8 The normal vector expansion ₃₇

4.9 The multipole expansion. ₃₈

4.10 The basic functions of the singularities. 39

4.11 Solution of the integrals 42

4.11.1 Far field integrals. 42

4.11.2 Near field integrals. 43

4.11.3 The velocity contribution of a straight line vortex with a

constant strength. 49

5

Consistence, convergence and stability.

50

5.1 Definitions. ₅₀

5.1.1 Stability. 51

5.1.2 Condition. 52

5.2 Consistence. 53

5.2.1 Some properties of the dipole potential influence matrix. 57

5.2.2 Determination of the maximum eigenvalues . . . 62

5.3 Referezices. ₆₇

6 Curviinearpanel method

70

6.1 Introduction. 70

6.2 The singularities integrals ₇₀

6.2.1 The disadvantages of reference plane method 72

6.2.2 Curvilinear coordinates 73

6.2.3 Choice of the curvilinear coordinates. 73

6.3 Surface elements calculations 76

6.3.1 Singularity distribution 76

6.3.2 Geometry description. 77

6.3.3 Metric tensor. 78

6.3.4 Normal vector. 79

6.3.5 The Jacobian expansion. 79

6.3.6 The radius vector expansion. 79

6.3.7 The multipole expansion. 80

6.3.8 The far and intermediate field integrals. 82

II

Cavitation

99

8

_{Qualitative analysis of steady (sheet) cavitation on a profile.}

103

8.1 Introduction. 103

8.2 2D cavitation. 105

8.3 3D cavitation 105

8.4 Conclusion. 106

9 Time dependent cavitation mechanism.

107

9.1 Assumptions and boundary condition for modeling cavitation. 107

9.1.1 Assumptions 107

9.1.2 Nomenclature. 108

9.1.3 Solution of the flow by use of the potential theory (panel

method). . 109

9.1.4 The Bernoulli equation for a cavitating profile in a

uni-form flow 109

9.1.5

The boundary condition on the

profile and the

sheet. 110

9.2 Introduction of the cavitation model 111

9.2.1 Calculation of the tangent velocity vector at the profile(1).112

9.2.2 Normal expansion into the field. 114 9.2.3 The Fredhoim equation for acceleration potentials 118

9.3 The unsteady Kutta condition. 119

9.3.1 The pressure calculation of a profile in an unsteady

con-dition 121

7

6.4 Vortex calculations.

6.4.1 Choice of the curvilinear coordinate 64.2 Vortex distribution.

6.4.3 Geometry description.

6.4.4 The Jacobian expansion. 6.4.5 The radius vector expansion.

6.4.6 The multipole expansion.

6.4.7 The far and intermediate field line integrals.

64.8 The near field line integrals 6.5 Consistence of curvilinear panels

6.6 Conclùsions.

6.7 References.

Grid generation.

7.1 Grid generation in span direction for simulating two dimensional flow

7.2 Description of the local geometry.

86 93 86 87 88 88 88 89 91 91 92 94 96 98 98

10 A heuristic derivation of the expected solutions of sheet cavity.141

B Equivalence between the velocity induced by a dipole

tion and the velocity induced by a surface vorticity

distribu-tion.

190

B.i References. 192

C A straight line vortex integral with a constant strength.

193

D To' transform a surface integral to a line integral.

196

III

Results

147

11 Results.

149

'11.1 Results of the panel method. 149

11.1.1 Flow around sphere 149

11.1.2 Flow around cylinder. 151

11.1.3 Flow around two dimensional profile. 154

11.1.4 Flow around three dimensional profile 163

11 2 Conditionmg of the mixed Fredholm equation 167

1L3 Simulation of the cavitation. 1.0

1.3.1 Growth of the cavity. 174

11.3.2 Pressure distribution. 174

11.4 Conclusions. 175

12 Recommendations.

176

1:2.1 Recornrnendaions for future research of the panel method. . . 176 12.2 Recommendations for future research on sheet cavitation. . . . 177

1'23 References. 178

13 List of symbols.

180

1:31 Notations 186

A Green's theorem and distribution of singularities.

187 9.4

9.3.2 Acceleration components.

9.3.3 Relation between source strength and normal velocity.

934 The mixed' Fredhoim equation.

9.3.5 Some adaptations of the source influence. . . .

9.3.6' Determination of the cavitating area 9.3.7 Time integration (ODE solvers).

938 An extension of time derivatives

9.3.9 Outline of the proposed method

References. -122 129 123 123 125 130 132 138 140

E Computations of basic integrals.

198

F The basic integrals.

202

G Calculus for curvilinear coordinates.

206

G.1 Curvilinear coordinates. 206

G.2 Base vectors. '206

G.3 Reciprocal space 207

G.4 Fundamental tensor of the first kind, or metric tensor. 207

G.5 Christoffel symbols. 209

G.6 Calculation of the Christoffel symbols 211

G.7 Gradient, divergence, and rotation. 212

G.8 The geometry of a scale and of a surface. 213.

G.8.1 Definition of a scale 213

G.8.2 The metric tensor of the scale and surface. 213

G.8.3 Christoffel symbols of the scale with an index 3. 214

G.9 References. . . 215

H Differential molecules.

216

11.1 Surface differential operator. 216

11.2 Line differential operator. i 217

11.3 Presentation of the molecules . . 218

H.4 Line differential mo1ecules 218

H.5 Surface differential molecules 222

I

Calculation of the normal series expansion coefficients.

230

Chapter 1

Introduction.

The aim of the work undertaken in this thesis is the description of the behaviour of sheet cavitation. Cavitation in generaloccurs when vapor is formed in afluid dúe to local low pressures. These low pressures are caused by high velocities in the fluid in nearly isothermal conditions. As a global pressure limit för the occurrence of cavitation the vapor pressÙre can be used.

Cavitation occurs in many forms. A common form in fluid machinery and on ship propellers is sheet cavitation, where a closed region of vapor is attached to the wall or to the propeller blade. The' volume of a sheet cavitaty cañ change very rapidly due to variations in the external pressure These pressure varia-tions occur when a ship propeller moves through the wake field behind a ship. The dynamic behaviour' of sheet cavitation leads to pressure 'fluctuations in the fluid, causing vibrations in the ship hull. Cavitation is a major cause of

vibra-tions' 'in' ships. The wake distribution of a ship and the requirements of optimum

efficiency makes it impossible to avoid the occurrence of cavitation The present investigation is made in the context of efforts to control the dynanuc behaviour of sheet cavitatión in order to reduce: the cavitation induced pressure fluctua-tions

The 'theoretical description of sheet cavitation uses steady potential theory,

which means that the viscosity of the fluid' is neglected. In most 'regions of

the flow this is 'applicable, but in the steady case it leads to two conflicting boundary conditions at the end of the cavity. The first condition is that the

pressure at the cavity surface isequal to the vapor pressure. The other condition is that the normal velocity at the'end' of the cavity is zero These two conditions can only be satisfied when the end of the cavity has 'the same derivative as the wall contour and this' is contrary to experimental experience. In potential cal-culations of steady sheet cavitation one of the. boundary conditions. is therefore

neglected 'in the region of'the closure'of'the cavity. [15,16,13,14] In he 'unsteady

is to investigate if an unsteady potential flow formulation leads. to .a realistic description óf: the unsteady cavity behaviour.

Physically the conditions at the end of the cavity are much more complex. In the steady case the pressure in the fluid at the cavity closure will rise, causing a

re-entrantjet in the upstream direction into the cavity. This results shedding

of vapor regions and to a periodical behaviour of the cavity. The clOsure of

a steady sheet can also be quite smooth, depending on the three dimensional nature of the cavity closure When the re entrant jet deviates from the flow

direction a strong cross flow at the cavity closure is generated, which makes the closure stable. Three-dirnensional effects are therefore important. The

in-vestigations in this thesis should also form the basis for the investigation of

three-dimensional cavity shapes.

Sheet cavitation also causes entrainment of vapor clouds. These clouds collaps in the fluid or on the surface, causing noiseand erosion. The formation of cloud

cavitation is very complex and has not been dealt with in this thesis. The

phenomena of cloud formation and collective collaps of bubble clouds is being investigated iii Twente University by van Wijngaarden [17,16] The dynamic be havior of the sheet as a single volume is also the basis for the calculation of the cloud formation, so there is a regular cooperation between Twente University and Delft University on these issues, in combination with MARIN.

1.1

Outline of the Investigation.

Iú the mechanism of entrainment of vapor clouds, viscosity may play a rôle. For the description of an unsteady three-dimensional sheet cavity the viscosity

can still be neglected and potential theory can be used. A panel method was

considered the most appropriate way describing the flow

Although the cavity behaviour on a propeller is the ultimate aim, this inves-tigation has been restricted to a three-dimensional profile. This eliminates the problem of the structure of the wake, which has been taken in the plane of the

profile.

1.1.1

Panel Methods.

The first tool which is necessary is the accurate prediction of the pressure dis tribution around non-cavitating foils. Many panel methods are still based on a constant or linear singularity distribution and flat panels,as first developed in 1962 by [1]. In these panel programs lift is generated through a dipole

distri-bution on an auxihary surface in the interior of the body or on the surface of

the body itself [2] . In these programs the von Neumann boundary condition is

used.

has curved panels or flat subpanel& Several second generation panel methods use the Dirichlet boundary condition.

One of the first higher order panel methods is the program 'BAC'(Ï975),

made by [3]. This panel method uses a second order geometry description and

source distribution. The dipole is still a mode-funtion distribution. The first

programs using the Dirichiet boundary condition were introduced in Ï975. Such

programs are e.g.'PAN AIR', 1975, 'MCAERO',11980.

Code's developed in the Netherlands are PDAERO/AEROPAN [4] at the NLR and the program of [10].

Finally panel methods have been developed with an advanced lower order met-hods. These programs have flat panels and a constant source and dipole. dis-tribution. These programs claime that their cakulation are of a higher

or-der accuracy. These programs solve the Dirichlet condition. Examples are

'VSAERO",1981 [5] and 'QUADPAN',1983 [7]

The blade sections of propellers are very thin.. In the region of maximum lOad-ing the thickness to chord ratio of the blade sections is generalily less than five percent. The strong variations in curvature at the leading edge of profiles cause difficulties in the calculation of these pressures, ùnless very small elements are used. Similar problems were expected near the closure of the cavity. Very small panels or a sudden decrease of .the panel size at the leading edge of the profile or at the cavity closure, causes numerical inaccuracies and instabilities. Since a program for a panel method still had to be developed, it was decided to develop a panel method with curved panels and with non-linear distributions of the sin-gularities over the. panels was developed, to act as 'a tool in the investigation of cavitatiOn. The program which would be developed should be similar as the higher order method program BAC, PDAERO/AEROPAN and the program developed by [10] and using 'the. Dirichlet boundary condition as in program such as 'VSAERO" 1981 and QTJADPAN',1983. Now the high order expansion (up to three) and the Dirichiet boundary condition are build in one program.

The order of the method can be chosen (0 to 3) and this makes the program

more flexible and the results can be compared with other programs.

1.1.2

Calculations of Sheet Cavitation

One of the first predictions of the. cavity shape of sheet cavitation was the

method of Geurst [11]. He used conformal mapping in the complex plane. This approach has been implemented for variable camber ratios [12].

Calculations of the steady cavity shape with a panel method have been made by [15]. In this model the panel code of [10] is used. The boundary conditions at the surface ofthe cavity were determined from a linear extrapolation from the profile. To .predict the cavity shape the kinematic and dynamic boundary condition were combined in one equation and linearized At the cavity closure this method exhibits a pressure dip due to the stagnation point.

The first investigationson the unsteady cavity shape based on panel methods

were done in [18] arid' [19]. 'The panel method used is a low order method., In the

cavity model the dynamical boundary condition at the free surface is used tó solve the dipole singularity at the panel fixed on the profile. The cavity shape is determined by integrating the kinematic boundary. This cavitation model albo has a pressure dip at the closure of the cavity.

In this thesis thesingularities describing the cavity arestil positioned on the non-cavitating profile surface, but a second order expansion is used to calculate

the pressures at the cavity surface. This' expansion uses the divergence and rotation theorems of potential flow. Also the cavity is assumed' to be unsteady, so the cavitaty surface is always in motion. At every time step the pressure at the cavity can be determined Because the pressure in the cavity is equal to the vapor pressure the motion of the cavity contour has to be responsible for the difference. From this velocity the cavity thickness after the next time step can be calculated.

1.2

Outline of the Thesis.

The thesis is divided into three parts. The first part descibes the higher order panel method without cavitation. The second part describes the work on sheet cavitation. The third part contains calculation results

In the first part chapters 2 to 6 describe the panel method. At the trailing

edge of the profile the implicit formulation of the Kutta conditión is used. The

calculations of the hydrodynamic influence coefficients is described in chapter 4.

The calculation of these coefficients and the expansions to the proper order re-quires careful bookkeeping. The stability and other characteristics of the panel method are diccussed' in chapter 5

During the development of the panel method it became clear that a curvilinear description could be advantageous. In chapter 6 this method is described. The

advantage is that the càlculations are carried out in the computational space

and therefore no transformations' to the reference plane have to be made. This formulation has not been implemented in the program, however.

The second part consisting of chapter 8 to 10 describes the cavitation mecha-nism First it is shown that two dimensional stationary models are not sufficient to describe cavitation. Therefore an unsteady cavitation model is developed. The singularities describing the sheet cavitation are positioned on the profile surface, but second order Taylor expansions in normal direction are made to calculate the potential, the velocity and the pressure at the surface of the sheet. During the calculations' of the cavity growth it became apparent that the calcu-lations became unstable after some time In chapter 10: an analitycal solution of thefree surface condition is therefore given, which proves that the unsteady model alone is mdeed not sufficient for the description of sheet cavitation in all

circumstances. An extension of the potential description of the cavity is still

necessary.

In the third part the numerical results of the panel methods as well as of the

cavitatiOn simulation are presented and discussed. From these results the re-commendations and conclusions are obtained as given in chapter 12.

1.3

References.

['1] Hess, J.L. & Smith A.M.O.

"Calculation of non-lifting potential flow about arbitrary three-dimensional bodies",

Report No. E.S. 40622, Douglas' Aircraft Co., Long Beach, California, 1962 Woodward, F.A.

"Analysis and design of wing-body combinations at subsonic and supersonic

speeds",

Journal of Aircrafts, vol' 5, No 6, pp.22-44, 1968

Roberts, A., Rundle, K.

"Computation of incompressible flow about bodies and thick wings using the spline-mode system",

BAC(CÄD) Rep. Aero Ma' 19, 1972

Sytsma, H.S., Hewitt, B.L. Rubbert P.E.

"A comparison of panel methods for subsonic flow computation",

AGARDograph :No. 241, Neuilly sur Seine, 1979 Maskew, B.

"Program VSAERO,, Theory document",

NASA CR-4023, 1987

lloeijmakers, H.W.M.

"Panel methods for aerodynamic analysis and design",

AGARD' report 783, Neuilly sur Seine, also N.L.R. report: TP91404U

Yougren, 11.11., Bouchard, E.E., Coopersmith, R.M. Miranda, L.R.

"Comparison of panel method formulations and its influence on the

devel-opment of QUADPAN., an advanced low-order method", AIAA Paper 83-1827, 1983

Katz, J. & Plotkin, A.

"Low-speed aerodynamics from wing theory to panel methods",

Moran, J.

"An introduction to theoretical and computational aerodynamics",

ISBN 0-471-87491-4, John Wiley & Sons, New York, 1984

Romate, J.E.

"The numerical simulatiOn of nonlinear gravity waves in three dimensions using 'a higher ordér panel method",

ISBN, 90-9002667-3,T.0 .Twenthe, Enschede 1989

[11} Geurst, J.A.

"Linearized theory of 'two-dimensional cavity flows", Excelsior, Rijswijk Z.H., 1961

Noordij, L.

Pressure field induced by a cavitating propeller",

International Shipbuilding Progress, Vol 23, No. 260, 1976 also Publication No 512 of N.S.M.B.

Lemonnier, H., Rowe A.

"Another approach in modelling cavitating flows",

Journal of fluid mechanics,Vol. 195,page 557-580, 1988 Pellone, C., Rowe A.

"Effect of separation on partial cavitation",

Transaction of ASME, Vol 110, page 182-189,, june 1988

Buist, L.

"On the origin and acoustical behavior of cloud cavitation",

ISBN 90-9004317-9, Enschede, 1991

Buist, J.,Raven, H.C.,

"A consistence linearized approach for the calculation of partial sheet

cav-itation",

ASM'E 5/88, 1988

Wijngaarden L. van,Buist,

"The emission of sound by statistically homogeneous bubble layer",

Journal of Eúgeneering Mathematics 26: page '195-210, 1992

Kinnas, S.,Fine N.E.,

"A nonlinear boundary element method for analysis of unsteady

pro-peller"sheet cavitation",

Nineteenth Symposium on Naval Hydrodynamics, Seoul Korea, 1992

['19] Fine, N.E.

"Nonlinear analysis of cavitating propellers in non uniform flow",

Report 92-5, Massachusetts. Institute of Technology, Department of Ocean

Part I

Summary: The panel

method.

This part gives a description of the three dimensional panel method. With the

panel method it is possible to compute the velocities at the surface around a

three dimensional profile or body. For these flow calculations it is assumed that the liquid is inviscid and incompressible. So the flow can be calculated by using the potential flow theory.

The surface of the body is divided into panels and is covered with a source and dipole layer. The sourcesand dipoles are distributed on the panel surfaces. The singularity distribution is taken to be of arbitrary high order (up to 3).The panels are curved and also approximated by a higher order approach.

The hydrodynamic influence coefficients are calculated on three ways. When

the collocation points are far away from the panel, the radius is kept constant. When the collocation points which are at the intermediate field from the panel,

the radius can also be expanded with a higher order method. For near fields the influence of the radius is determined exactly for the reference plane and

corrected' with a small curvature expansion. With these descriptions (source-and ipole-strength, radius (source-and curvature expansions) a set of equations can be formulated containing the unknown source- and dipole-strength. By using the implicit Kutta condition and the inner Dirichlet condition the source and dipole strength can be solved. The velocities can be determined by taking the first tangent partial derivative of the dipole strength on the body or profile.

Chapter 2

The integral equations..

2.1

Introduction.

Thepanel method is a numerical method 'to calculate the (potential) flow around a body. This 'body can also be a wing or a foil. With. this method the velocities (in two directions)'on the surface of the body can be calculated.

The' panel 'method is based on the principle of the Green's integral theorem. According to this theorem it is possible to transform a three-dimensional linear homogenic 'differential equation to a two-dimensional integral equation. So, the three-dimensional Laplace (potential) equation can be transformed' to a surface integral equation. This surface integral equation is called Green's identity. The integral equation represents a distribution of sources (or sinks)' and dipoleson

'the surface. In a special case when the source strength is unknown and has to

be solved, the integral equation of Fredholm of the first kind appears When the dipole strength is unknown and has to be solved, the Fredholm equation

of the second kind appears. Both of these equations are also called' the kernel equations.

'For solving the integral equation the body surface will' be divided in several panels. The integral equation can be discretizised for each panel by

substitut-ing an unknown strength of the source and dipole distribution. By ussubstitut-ing the

boundary condition of tangential flow the unknown strength of the sources and dipoles can be solved. If the strength of each dipole or' .source is known the velocities can be determined on the surface.

The advantage of panel methods is that the problem is reduced to a two-dimensional (= surface')' problem instead of a three-two-dimensional (= volume) problém. The grid generation is also reduced to a tw,o dimensional problem.

So the grid has to be generated on the surface of the body only. This is in

contrast to' the three. dimensional 'grid generation, where a 'lot of points have

equations, instead of n3 equations, have to be solved for determining the velocity.

2.2

The Fredhoirn equations.

An incompressible and inviscid flow around a profile or body can be calculated by a potential flow. A potential function (, which describes the potential flow)

has some important properties. One of the properties is that the curl of the

velocity vector is equal to O. _{Further the potential flow is a solution of the} Laplace equation. This is a homogenic and linear differential 'equation andso

the superimposing property can be used. With this property, it is possible to

transform the 'Laplace equation to the Fredholm's equation, with help of the Green's theorem. (See appendix A)

If there is a distribution of sources and dipoles on the surface of the profile or body, the following integral or kernel equation can be described:

where:

the potential at point P on' the surface Q.

p = I'(Q) the strength of the dipole or potential at point Q on

the surface Q.

= On the strength of the source or the normal derivative of the potential at point Q on the surface' Q.

This equation has to be solved. It can be done in two ways. First the

prob-lem is considered that on Q' is known. (Dirichlet boundary coñdition.) An

interpretation of equation (2.1) can be madè as an integral 'equation with an

unknown source strength o (The known part will be placed in the right hand side. and the unknown part will :be placed in the left hand side):

J!_.JdQ=

(P)+-ci

Equation (2.2) is called the Fredhoim equation of the first kind. The unknown source strength can be discretizised and solved numerically. 'But according to

reference [4] the numerical approach of theseequations generally gives problems.

But fortunately equation (2.1) can be handled in another way.

The second way of considering equation (2.1) is described as follows. If on Q is known (von' Neumann boundary condition), anew integEal equation can be formed, with the known paEt at the'right hand sidè. So equation (2.1) yields:

+

_JJ_±dQ

(2.3)

(Q)---(-)dQ

_{ôflQ r}

(2 I)

Equation (2.3) is called the Fredholm equation Of the second kind. In general

this kind of euations can be solved very well by a numerical approach. But

solving this equation has to be done carefully, because the solution of 4) is not unique

A description of the dipole distribution is recognized in the left part of the

equation (23),, with a strength of

= 4). It's well known that a uniform

dipole distribution on a closed surface gives a constant potential (=0) inside

that surface. Thus the (normal )velocities are O at the inner ide of the surface.. An arbitrary constant can be added to the potential, which does not change the right handside of the equation. So the nonuniqueness of the solútion is proved.

The above method is known in the literature as the direct formulation.

2.2.1

_{A more generai formulation of the. integral}

equa-tions.

The method described in appendix A, only gives a description of the potential

of one side óf a volume V. But in the same way it can be dèscribed on the. other side. So the whole space will be described by taking both descriptions

together If a source point inside V is considered, then p is equal to i for the inside description and O for the outside. description (See A1.8a). Superimposing of the both descriptions yields:

4)(P)

_Jf

(±(!

- !-)(4)

-

4)I)P(J))

_(2.4)

A source and a dipole.term is recognized. The.strength of this term are

respec-tively:

.04) 04)'

and On On (2.5)

p = (4) 4)')

In the same way a source point outside V can be considered. In this case the

following equation is' obtained:

4)(P')

ifl ((

_- ) -(4)'

_4))(±))

dQ (2.6)

In this case the strength of the source and dipole are respectively: 84)' 84)

(2.7)

The equations (2.4) and (2.6) describes a moregeneral formulation. Especially for a determined choice of '1'(= O) the direct formulation appears.

The uniqueness of this general integral equation is not always guaranteed,

because it's posalble to substitute.a von Neumann boundary conditión only.

(By example, substituting equation (2.5) ô = - and = O , equation (2.3)

appears.) But by making a good description between the sources and dipoles and a good conditioned boundary condition, the uniqueness of the solution of the integral equation will appear.

2.2.2

The principle value integrals and jump relations.

The principle value appears at the place, where the radius becomes O. The contribution of this principle valúe can be obtained by solving the Cauchy prin-ciple value integral. The prinprin-ciple value integral can be determined by moving

the point P (, which is inside V) to the surface fI. In the case that this point

touches the surface, the hmit of the mtegral has to be determined, which is valid on fi. The result of this limit is given by the following equation (The term of

appears as the result of the potential jump by crossing the dipole layer. The value of this potential jump can be determined with the Cauchy integral theorem. The integral, which remains, is called the principle value integral. See reference [1']):

=

Jf

+

(Q)_(±))

dQ -

(2.8)

PeQ

Taking the normal derivative of (2.8) the following equation appears. (In the same way as the potential jump in the dipole distribution, a jump of o(P) due to source layer appears. See also reference[1]):

=

JJ

(q(Q)-(--)

+

dfi + cr(P)

n

PEfl

(2.9)

Similar the formulation for the potential and normal velocity for a point P', which is approaching the surface at the outside of V (= inside V', see figure

can be given:

cV(P')

=

:-Jj

p(P')

_(2.10)

and

J_J! (Q')__()+

p Q')

_(L!))

dQ

-P'E

(2.11)

These expressions çan be used to solve:, the boundary condition problem for

Figure 2.1: Relation between the profile and the wake surface, including the

inner and outer space. the outer volume of V'.

22.3

The integral equation for a profile, in a flow.

For the 'determinatiön of the flow around a profile asurface has to be generated, which encloses the space of the potential flow. The whole enclosed surface is divided into fourdifferent surfàces (see figure 2.1). The surfaces are the surface of the profile (2.,,), the wake surface

_(2w)

a flat surface at, infinity (') called

Trefftz plane, which is perpendicular to the wake surface or flow directión, and a surface at infinity which closes the space of the potential

00

11

_{4J{i$T j(--}

I ro

i

in which:

t: coordinate in flow direction.,

s: coordinate in' span direction direction.

ST: coordinate of trailing edge in span direction.

) di

(2.15)

¿p(P) ±

_j(T(Q)::;::

_{+ (2.12)}

where Ç = ìp +ìw +

+1oo

Note that the potential has been split into a known part: the incoming

potential

_{o (= U}

) and an

nknown part: the perturbation potential

ça. ( 'o + ça). Because the potential is linear, the incoming potential can

be subtracted. So only the perturbation potential remains. This potential

also satisfies the Laplace equation. To get the exact potential the perturbation potential has to be added to the potential for parallel flow.

For a good description of this integral equation, the integration over each surface shall be discussed. Further it's assumed that ça' O. ((inner) Dirichlet boundary condition.)

The wake surface

The thickness of the wake surface or plane is taken infinitesimally small and the velocity, perpendicular to this surface, wiff therefore be continuous, there is no velocity jump over that surface in that direction So the first boundary conditión requires the following jump relation:

LIôn (2.13)

The potential jump is equal. to the dipolestrength according to (2.8) and (2.10). Also according to this boundary condition the source strength is equal to 'O. So the contribution of the integral of fw yields:

-jj

[ça]-(-.)dQ

(2.14)

û

It's assumed that the dipole strength at the wake surface does not't vary in flow direction and that the dipole strength in the wake varies over the span [tL..tR]

with a strength equal to the dipóle strength at the trailing edge. So the right

2..3

The surface at. infinity (í).

By making

'

O in the inner space (see figure 2.1 and 22), the boundary

condition has to obey to the following equation:

¡='=O

at (2.16)

Apart from the condition that the perturbation potential is equal to O, the'

Figure 2.2: The positions 'of the 'surfaces for solving the flow around a profile, perturbation vefocity is zero too;

at1

(2.17)

Thus the sourceas well as the dipole strength has no contribution for the integral equation and so there is no contribution to the integral from

2.3.1

The flat surface at infinity (fT) (Trefftz plane).

By making çd O in the inner space, the source strength at the flat surface at infinity satisfies is also zero. So:

atQT (2.18) On

The situation is more còmplicated for the dipole strength. The dipole strength at 1' in the neighborhood of the wake will not be 0. There is no contribution to ço'(P E &) from QOO , as shown before and also Ilp gives no contribution,

because of the infinite distance. So the integral equation in this neighborhood

can be described as follows:

ça'(P E T)

=

JJ

+

JJi4..(T)dQ

(219)

fw

Because the contribution of the integral, which is valid on 1.2w and which

influ-ences The potential on 2T, is not 0, the dipole strength of ÇZT is not 0. So the contribution to a point P in V, due to 1T, will be:

ço(P) = (2.20)

However the contribution of (2.20) vanishes, when P is in the neighborhood of the profile because 1T is considered to be at infinity.

2.4

The surface at the profile (ap).

This is the most interesting case, because the velocity distribution on. thissurface

has to be determined. Above it has been shown that there is only a contribution of the dipole of the wake surface. So the integral equation of equation (2.12)

remains:

=

JJ

+

Q()) d

+

II

-PEV

In case that P is. a point of Qp the equation yields:

1111-1

8 1\

+

-(-)) dQ +

11+

(2.21)

The normal derivative of equation (2.22) can be written as:

892 i

f((

ô i

0

0 i\

_(-.-)+iT_-_-(__))dc +

i

i

fr

ô

ô 1

+r+ij [çaI(--)dì

Equation (2.22) represents an incomplete integral equation for the unknown ¡, o and [ç°l. To complete this equation, the Kutta condition will be introduced and

a choice has to made for p and o on Qp. Equation (2.23) has to satisfy to a

Von Neumann boundary condition. But there is an equivalent integral equation for the inner solution with a Dirichlet boundary condition. In that case o can be obtained, by making use of ça' O inside ('p:

(2.24) Substituting this value in the integral equation, which is valid for the inner side

ofp yields:

=

ll

< . >

_±

_JJ/j

(±)

+ (2.25) çw

for P' E p

It can be seen from equation (2.7b), knowing that ça' = O for the inner space,

that a potential jump on lp appears of ça

-

ça' = p. This means that ça at P

on p has the following value:

ça(P E ìp) = p(P)

(2.26)

This relation can be used for the determinatiön of the velocities at the surface of the profile. (See section 2.4.4. The velocities on the profile surface)

So two integral equations (2.23) (with a outer von Neumann condition) and

(225) (with an inner Dirichlet condition) are described, which guarantee a

proper descr.iption of the problem.

2.4.1

The Kutta boundary condition'.

To formulate the Kutta condition in this integrai equations, the potential flow will be studied at the trailing edge of the profile. At the upper side as well as at the lower side the (perturbation) potential will be continuous at the trailing

edge. This means that the potential of a point P1 at the upper side of the profile will approach the same value as the potential of a point P1' at the upper side of the wake surface, when these poiñts are approaching each other. (See figure 2.3) The same situation is valid for the potentials at the lower side ofthe profile

Figure 2.3: The situation near the trailing edge.

(for the points respectively P2 and P ). Further it's known, thät the potential

jump over the wake surface is equal to p (see equation (2.16)) at point P'. So the difference between the potential at the upper side and the lower side of the

profile will be:

W(P2) ça(Pi) p(P2)

(2.27)

Using the jump conditions of o- and p on p (see also equation (2.26)), it is

found that:

,u(P2) - p(P1) - p(P2)

(2.28)

2.4.2

Implicit Kutta (boindary). condition.

In case that the source strength is known and the dipole strength has to be

determined by using the boundary condition on 2p (e.g. the Dirichiet condition

= O) and equation (2.28), there. is enough information to determine the

unknown dipole strength. The above, description is called the implicit Kutta (boundary) condition.

2.4.3

Explicit Kutta (boundary) condition.

If the sòurce strength is not prescribed, but the dipole strength on p is

pre-scribed, equation (2.28') has to be satisfiéd. But this is impossible, without an explicit relation for the dipole strength at the wake in the plane of the trailing

edge This explicit relation is called the explicit Kutta (boundary) condition

For the explicit condition a description can be used, which can be formulated in terms of a direction form or in term of a: pressure form.

Direction form.

The integral equation (2.23) can be taken, for a point P on the wake surface near the trailing edge. This point P approaches the trailing edge This leads

to: - '

- --

11

(_r

____±'\

_{d12 +}

P-.T -

4irjj

r "10fl.P OflQ' r np

1f!

00-1

-. E'1

--(--)d = -

U,.,o

_{np >)P..T}

(2.29)

Where T is the point at the trailing edge

Pressuré form.

There is no pressure jump when the wake surface has been crossed. Acéording to Bernoulli's equation the sum of the squares of the velocity components has

to be equal, when the pressure is constant. So an expression in the pressure

form appears:

(V3(U00 . z + çt)), - (V3(U

. z + ça)), - 0

(2 30)

For P1,P2*T

p(P) can be calculated from equation (2.22) using the Dirichiet boundary

con-dition. But in the case, that the source strength on Ip is chosen in that way, where inside Ip ' = 0, an easier equation can be used (derived from 3.5b)):

y(P) = (P)

2.4.4

Velocity on the profilè surface.

By making use of the potential jump relation of a dipole layer and of the choice

that the inner potential ' = 0, the potential can be determined at the sürfàce according to equation (2.26):

ça(P E

p) = p(P)

According to this equation the potential on Qp is known. The first (partial)

derivative of the potential in any direction delivers the velocity in that direction. So by taking the first partial derivatives of the potential in the tangent plane of the dipole strength the velocities are obtained on the surface of the profile.

2.5

Conclusion..

For a flow calculation around a profile it's possible to make a well conditioned

surface integral equation with a unique solution. Especially for a profile in a parallel flow the integral equation (2.22) (with an inner Dirichlet boundary

condition ça' = 0) supplemented with the implicit Kutta condition (2.28) or

the integral equation (2.23) (with an outer von Neumann boundary condition) supplemented with the explicit Kutta condition (2.26) and (2.29) or (2.30) gives a satisfying solution.

When the solutions of the dipole strength are known, the velocities on the profile surface can be determined by taking the two tangential first derivatives of the dipole strength. (In this condition the inner potential ça' has to be equal to 0.)

Chapter 3

Discretizat:ion of the

integrai

eq.

uations..

The integral equation (2.25) with an inner Dirichiet boundary condition and the implicit Kutta condition is easier to work out than the integral equation

(2.23) with an outer von Neumann boundary condition and an explicit Kutta condition. If there is a stationary flow around a profile, integral equation (2.25) gives a satisfying description of this problem. So only this description is worked

out. The integral equation with an inner Dirichiet condition (' = O with

assumption that the dipole strength at the wake surface is constant in flow

direction) can be written as

+

Jf <Ü flp> r(Q)

+ "P

Ii

Ç(,Q)5 dt(Q) TL Law(1) (3.1)

Note: In this equation the implicit Kutta condition is not used yet, It is

assumed that the dipole strength at the wake surface is known.

In this equation several terms appear, which are dependimt of These

kernel function which have influence on the velocity): and r(P,Q) o

i 1

OflQ "r(P,Q) 8

i 1

r(P,Q) O

8 1 1

Onp Onp'r(P,Q) potential source.

potential dipole or doublet.

normai velocity source.

normal velocity dipole or dOublet.

The influence of these kernel functions, which are distributed on the profile

surface and the wake, on an arbitrary location in the flow, has to be solved.

To solve integral equation, this integral equation has to be discretizised.

The discretization is carried out in the following way The surface of the profile and the wake is divided into a large number of small areas (or elements) called panels A choice has to be made for the collocation points, where the boundary conditionis satisfied. Also the Kutta conditiOn has to be satisfied. A good

choice of the position of the collocation point is the center of every panel. The next step is to approximate the integral for the whole surface by summoning up

the integrals for every panel:

d =

NPXpNPYp

1

(3.2) I 1

JJr(Pjcjricj,Q)

u'1

2P1 _ptipIIipI I

[f

_{__}

1

)d1=

jj

z(Q) ôflq r(Ft ic][jc], Q) lip NPXPNPYP ô

1

(1

-ip=l jp1

-'[ipI ti pj

s!:;; ((L1

r(PÌCjj,

)dQ

-t'pIIpI

J p(t(Q))

_f

₇

Q))ds.(Q)] dt(Q) = Lw()

1

NPXw

¡

I

,j

(r(j]],:Q):

ip=1 IipI (1ÌI]Lf,Q) (3.3) )ds(Q) dt(Qr) (3.4)

In which:

NPXp =. Number of panels in (flow or) x-direction on the profile. NPYp =.Number of panels in (cross flow or) y-direction on the prófile. NPXw = Number of panels in (flòw or ) x-direction on the wake surface. ip = index of panel in x-direction.

jp = index ofpanel in y-direction.

je. = index of cóllocation point in x-direction. je = index of collocation point in y-direction.

Note: If the collocation point lies at the saine panel (Ip = je and jp = je),

the. Cauchy integral has to be determined.

The following stage is to approximate the integral for every panel and every collocation point. The computations of these integrals are based on Taylor series expansions The advantage of this method is that the running time will decrease

in relation with other methods, (like the Gauss or the Simpson integration).

Further an analytical description of the Taylor expansion on the panel surface is known. So the influence of some part of the Taylor expansion (such as the curvature , multipole or singúlarity layer expansion) can be determined.

The functions are.expanded from the panel center point in coordinates of the reference plane or tangent plane of the panel. (For a more detailed description of these procedures see the following chapters) An expansion of the source and dipole strength is written as a function of the (unit) strength in the panel centers of the own and of the neighboring panels. Further an expansion is obtained for the multipole and the curvature. Also an expansion of the Jacobian has to be determined, because the integrating process is worked out 'in the tangent plane, instead of the coordinate system of the profile. The expansions of each part, except the curvature expansion, have to be multiplied with each other and so

a new expansion appears. By integrating this expansion, the contribution of

one panel can be written as function 'of the unknown unit strength of the own

and neighboring panel. The actual source strength can be replaced from the

boundary condition (°ìplFjp]

_{= -}

. Ti >).

3.1

Vectòrizing the, source and dipole

formula-tion.

The double panel and collocation point indices are converted into one single array. A panel or collocation point or can be indicated by a single index:

Ip'

NPY * (ip - 1) + jp

ici

_{NPYp*(ic 1)+jc}

This process will be done for every panel and for every collocation point. So

there is for every strength a relation to a collocation point. (For a higher

order process several source and dipole strengths are used over more panels. The coefficients of one source, (or dipole), which appears of this approach, may summoned together.) So respectively the source and dipole potential 'influence

becomes: and i

1f

np

fJ

r(Pji,

Q))d1

+

tP1)

=

dí

_{= >}

HICSP['][Ii]WcIi

P1',Q) _ip'=l (3.6)

JP(t(Q))n.

_J

ôflQ

r(pi1,

)ds(Q) dt(Q) = TL Lnw(l) NPW

HICDPW['][l') Wp»

ip'1 =i where:

NP

= NPXp * NPYp number of panels at the profile.

NPW

= NPXW * NPYp number of panels at the wake surface.

HICSFI][lJ: hydrodynamic influence component or coefficient of the

source potential of panel ip' on control point ié

HICDP[II][I]: hydrodynamic influence component or coefficient of the dipóle potential of panel ip' on control point ici and

HICDPW[l][l]: hydrodynamic influence component or coefficient of the dipole potential of wake panel ip" on control point ic'.

3.1.1

implementation Kutta condition.

The influence of the wake surface can be added toequation (3.7),, by substituting

the implicit Kutta condition. With this description the dipole strength at the

wake can be determined as the differences of dipole strength at the upper and

(3.7)

lower side at the trailing edge and thus as a function of p' So equations (3.1) remains: NP NP -. >2

- >2

HiCSF1cij1jpij (< >) p =1 $pi=1 (1 < ici < NP) Where:

HICDP($CIJ[IPI] = HICDFÌjcs]Ljpl] - and ip' and ip" are corresponding iñdices with the same dipole strength of the implicit Kutta

condition. PLIp'] E

So there is an unknown vector jii1 which corresponds with the NPD un-known components of the dipole strength. This vector can be solved by the

L.U.decornposition or Conjugated Gradient Square Method (C.G.S.Method). By taking the first tangential derivatives of the dipole strength the velocities components on the wing surface can be determined.

Chapter 4

HIgher

order panel

methods based

on

calculations in. the

reference plane.

4.1

Introduction

Of each panel the influence, of the distribution of the singularity on a collocation

point has to be computed. The influence of the distribution is given by the

following integrals:

=

_JJJdd7)

(4.1)

source potential influence

Qd

=

LI!

4<

>dí

=

L!!

Jded1 (4.2)

dipole potential influence

t7 = L J!

= (4.3)

Vd =

i fHøi

f 7®&

r5 i

ffj'®i'

i

_[

®Jft

-jj

Jddr

4ir

J

IV 8W

with: ji=il®Vp

and ÔQ is the edge of Q. (=curve)

(4.4)

dipole velocity influence

These integrals have to be solved for every panel and collocation point. In this chapter the integrals are solved with the so called "reference plane method" or "physical space method". This method is based on all functions on the panel

being transformed into a tangent plane of the panel This is m contrast with

chapter 6 "Curvilinear panel method", where the solution of the above integrals is based on the curvilinear coordinates. (This is also called the computational space.) in that case the functions on the panels are described by two curvilinear coordinates, Which are laid down at the surface of the panel.

The solution of these integrals can be obtained as follows. First these inte-grals are transformed into a local orthogonal space (e,, ij, ), where the coordi-nates of the panel-center coincide with the natural coordicoordi-nates. In these local coordinates an expansion of several parts of these integrals can be made. Thus

an expansion of the Jacobian (this expansion is due to transformation of the

coordinate system), the radius, the normal vector, the source and dipole distri-bution is made. These expansions have to be multiplied with each other and a new expansion is obtained. Substituting the boundary conditions, determined by the edges of the domain and integrating this function, gives the solution of the integral equation.

4.2

The order of the expansion.

The order of the expansion can be chosen arbitrarily. The advantage is, that the influence of the order of the expansion can be investigated. Further the

algorithm is more general, which decreases the possibility of code errors. Un-fortunately the coefficients cannot be determined directly, which increases the

computing time. Four different kinds of order are distinguished: The order of the source distribution (NS).

The order of Che dipole distribution (ND). The order of the curvature expansion (NC).

4. The order of the (power of the) multipole expansion (NR).

The order of the expansions may be. different. The order of a new compound expansion, which appears after multiplying two expansions, will be equal to the maximum order of the original' expansions

42.l

Tùuncation errors of a power series expansion.

For a single power series 'expansion is obtained:

a = >ajE +O(ëv+i))

(4.5)

in which:

N the order of the expansion. O(), the Landau order symbol

The truncation error will be of the orderof N+ 1, when the highest component of theexpansion is equal to N. The truncation error is ¿N+i. The local truncation error will' be equal to:

o(e(N+l))

f(N+i)(.)

1)

- (N+1)!)

where:

is a value between O and

and f(")() the N + 1 derivative in

The truncation error can be estimated by taking the absolute maximum of f(N+1)() in the interval between O and' ¿. Thus

(N1)!

For an expansion with two variables the power series expansion is:

N N-i

a =

_{a1e'i'n +}

Q(/j(N+1))

1=0 j=0

where:

h: variable with a length approximately equal to ¿ and rj.

(4.6)

(4.7)

(4.8')

The truncatiOn error of two variables can be treated similarly as the truncation error with one variable. 'So the. truncation error will be 'of the order N +1, when the highest component of the expansion is equal' to N.

Finally after a multiplication the components of two series expansions are written as:

Cmn =

a1 bm_i nj

(4.9)

i=O j=O

4o3

The geometry description.

A practical method is to use quadrilateral panel elements. The advantage of

quadrilateral panels is that the index-registering will be easier. For a element five coordinates have to be known., The first coordinate is the panelcenter. The

other four coordinates are the córner points. It is obvious that two corner points of one panel have the same coordinate as the corresponding corner points of the

eighboring panel.

By using quadrilateral elements a wing can divided into a number of rows (NPY = number of panels span direction) and a number of columns (NPX = number of panels in chord direction). An arbitrary panel has the index ip for the chord direction and jp 'for span direction. The geometry of the surface is known at the panel centers. These panel centers are indicated with an index in chord direction ip, an index in span direction jp and an index c, which gives the coordinates x, y, z. (see also figure 4.1) The local geometry description can be determined by using the differential molécules. By putting these differential molecules on the vectors of the panel center the local geometry is described. How this' method works is explained in chapter 6 "Curvilinear panel method". These molecules determine the geometry in the curvilinear coordinate system.

4.4

From curvilinear coordinates 'to reference

plane coordinates.

The functions, which are described at the panel surface, can be transformed into the reference plane. An arbitrary function at the panel surface can be

approximated by a polynomial function in 'the cUrvilinear coordinates ¿, ij:

N Ni

f(,'q) =

+ O(h")

(4.10)

i=O j=O

With

= D7f + O(h"'')

The coefficients of the polynomial can be determined by a finite difference op-erator D' (see chapter 6 "Curvilinear panel method" and appendix H "Diffe-rential molecules.").

Figure 4.1: Indexing of the panels. Relation between the panel index and the coordinate index.

Determination direction reference plane.

The base vectors are determined by the first derivatives of the geometry de-scription:

(N+i)

e = L1

+ o(e"1)

in direction, while i = 0 (4.11)

and

e'2 = Lr±O(r')

in direction, while = 0 (4.12)

Note that these two base vectors are not perpendicular to each other and the length of these base vectors is not equal to one.

The normal vector.

The normal vector of the reference plane can be determined by the cross product

of the base vectors at the origin of the panel (= panel center), instead of the

NW +1.

NPY

NW

ic,

Yjp'

1Î

0 1 2 3 NPX 0 1 2 3

NPX NPX-l-1

- X

ic, jc = index of collocation point. ip, jp = index of panel centre. i,j = index of panel edge.

cross product of the diagonals.

4=

'1®'2

+O(h')

Ile'1 ® e'211

Local base vectors.

Now a set of three orthonormal base vectors can be defined. The first base

vector is defined as the normalized base vector in ¿ direction. The second base vector then is the cross product of the normal vector and the first base vector:

e1 =

IIe'i II

=

(4.14)

Transformation of the Taylor coefficients from curvilinear plane to

reference plane.

The nomenclature of the curvilinear coordinates and the coordinates of reference

plane are both.C, and (. To avoid confusion, the coordinates of the reference

plane are temporarily indicated as n = x',y = 2 _{and z = z3. This}

nomen-clature. is only valid in this section. The curvilinear coordinates are similarly indicated as C =¿l, =¿2 and ( = ¿3.

Now the. derivatives of z = z' and y = z2 to ¿ ¿1 and ?J ¿2 up to the

third order can be determined by:

i

L'Ei =< _-, e ¿> a2 _-;

I.

_-<

X

-Ô.I

X'fjEkçm =< ÔCIOCkOC.m

The Taylor expansion of a function f in the curvilinear coordinate system can be expressed as a Taylor expan8ion in the reference plane. For this calculation the chain rule is used. The Taylor expansion up to the third order is:

(4.13)

linear terms

quadratic terms (4.15)

ftk = >2f,jhz"k

fekem >

(fzhx"em

+

>2fhiXkX'm)

h=1 i=1 2 = >2 frhX"ekEm +

hrl

fxh(xhekxiemen + X

_{mXek +}

xh _{X1k ) +} h=1 i=1

222

>2 >2>2 frh,,ivjx"ekx'emx2e h=li=lj=1

The derivatives f on the right hand side of equation (4.16) are

lin-ear. (Note that f,. =

sa these terms can be added.) All the derivatives of the curvilinear coordinates are known. Aset of equations can be written as a

vector and a matrix. The unknown veëtor f,, contains the terms fh .. . f,,hj. The matrix M contains the coefficients of xhek ....Xhkximxie,. of the right

hand of equation (4.16). The vector f contains the terms of the left hand side

of equation (4.16). The derivatives to the curvilinear coordinates can be

de-termined by the differential operator (see appendix H" Differential molecules."). So equation (4.16) can be transform to:

Mj,=j

(4.17)

The unknown vector L with the coefficient fh ... can be solved by the

inverse relation:

(4.16)

(4.18)

The coefficients' of j, are the Taylor coefficients at the origin of the reference

plane for the same fúnction as defined' in the cur'ilinear plane. So the sin-gularity distribution and geometry can be expressed in the reference plane by substituting the relevant quantities.

4.5

The transformation from global space to

local space.

The global coordinate system is the earth fixed cartesian coordinate system.

All the panel centers, corner points and collocation points are given in this

coordinate system.

The local coordinate system is a cartesian coordinate system, in which the e-axis is along the normal of panel ip, jp. The ¿- and ij axis are the orthogonal

tangents of the paneF surface of panel ip, jp. The origin of this coordinate

To solve the kernel integrals, it is practical to transform the problem into a natural or local coordinate system. In the natural coordinate system the first derivatives of the surface in the environment of the panel center will be small. For the.transformation into the local coordinate system the transformation

matrix L is introduced. This matrix L is composed of the normal vector and

two base vectors, which are perpendicular'to this normal vector. So the normal-vector as calculated in (4.4) is one of the base normal-vectors. The row elements of L are the elements of the base vectors of the local space:

L =(è)T

(4.19)

This is an orthonormal matrix, so the inverse of this matrix is the transpose of the matrix L:

L1 = LT = (è1,i2,;)

(4.20)

This matrix L' can be used for, transformation of the coordinates, defined in the local space to the coordinates defined in the global space.

4..6

The panel surface.

Because the zero and first' order approximation of the reference plane coincide with the panel surface. when there is no curvature, it is only useful to make a Taylor-expansion to the second- or higher-order. The cur.vature expansion of the surface can be described' as a polynomial function of and i in the ¿-space:

NC NCi

a13r) + O(h)

(4.21)

i=O .J=O'

The coefficients a1 of the polynomial are calculated as follows.. First a number of known points at and around the panelin' i-coordinates has to be transformed

to -coordinates:

(k = >13j(Z(ip+dik][jp+djklj - Xip][jp]j) (4.22) i=1

The indices dik and dj represent a relative part of the differential molecule (see appendix H "Differential molecules.").. 'By adding.'these relative indices to the indices of the panel center, the points on the panel surface are given in the global coordinates.

By substituting (k in f of equation (4 16) the local geometry of a panel is determined by using the differential molecule. From' this calculation it follows

that the following coefficients are 0:

000 = a01 = a0 0.;

4.7

The Jacobian.

4.7.1

The Jacobian for the surface integral.

Because the integration domain is transformed to the tangent plane the Jacobian has to be determined The Jacobian only has to be determined if the expansion of the curvature has been approximated by a second or higher order expansion

(For the zero and first order process the.J:acobian will be equivalent to 1, because

the reference plane and the panel surface are coinciding.) For a higher order approximation the Jacobian can be calculated as follows:

and by writingamn as:

Cm_i CJ

+

dm_j .dij (4.28) i=O j=O and: =

Ê Ê a(m_i][n_j]"1

(4.29) i=Oj=O (4.27)

J=

(4.23)

The components of the Taylor-expansion of called c are:

c1 = (i + 1)a[+l] (4.24) and the components of (, called d1 are:

d1 = (i + l)aI[jl]

(4.25)

(a1 are the components of the curvature-expansion) The Jacobian can be developed as a Taylor expansion using:

NC

v'i +

a =

(1 + a) =

k1c +

(4.26) i=O Where:

k,= (

_\ZJ

'(').(.i)....(i+1)

z!

andkò=1

-By substituting in a: