Simulation of unsteady motion of a propeller in a fluid including free wake modeling

Engineering Analysis with Boundary Elements 28 (2004) 633–653 www.elsevier.com/locate/enganabound

Simulation of unsteady motion of a propeller in a fluid including free wake modeling Gerasimos K. Politis School of Naval Architecture and Marine Engineering, National Technical University of Athens, 157 10 Zografos, Heroon Politechniou 9, Athens, Greece Received 28 February 2003; revised 20 October 2003; accepted 23 October 2003

Abstract The problem of flow around a marine propeller performing a general 3D unsteady motion in an infinitely extended fluid is formulated and solved using a boundary element method. Hydrodynamic modeling of the freely moving—unsteady—trailing vortex sheet emanating from each blade is achieved, using vortex filaments and a time stepping method. Thus vortex wake – blade interaction can be taken correctly into consideration. The method uses bilinear panels for the representation of blade and free wake geometry and constant panels for the representation of the unknown dipole intensity. The proposed method produces very stable rollup wake patterns for both steady and unsteady problems for a broad range of discretization parameters. Application of the method to a number of specific cases of steady and unsteady propeller motion has shown the very good numerical performance of the method and that good predictions for forces and pressure distributions can be obtained. q 2004 Elsevier Ltd. All rights reserved. Keywords: Boundary element method; Unsteady lifting flow; Free wake modeling; Wake rollup

1. Introduction Boundary element methods (BEMs) have long been applied for the solution of propeller design and analyses flow problems. The first 3D BEM for analyzing the steady flow around a marine propeller can be attributed to Hess and Valarezo [1]. In this original paper the classical Hess and Smith formulation [2,3] has been used for the representation of a steadily translating and rotating propeller. Trailing vortex surfaces emanating from blade trailing edges have been modeled using a prescribed wake shape (PWS) in the form of the classical ‘frozen’ helicoidal model. It is well known [4,5] that boundary layer asymmetry, in case of flows with lift, results in unstable regions of flow behind the body trailing edge, where the tangential component of velocity changes abruptly. These regions are termed trailing vortex sheets. In the context of BEM these regions have to be modeled as free boundaries satisfying certain conditions known as Helmoltz vortex theorems [5 –7]. This introduces an internal geometrical nonlinearity to the BEM since the domain space of E-mail address: [email protected] (G.K. Politis). 0955-7997/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2003.10.004

the integral equation is initially unknown (it has to be determined by solving a free boundary problem). Correct inclusion of the vortex sheet roll-up geometry in propeller flow problems is absolutely necessary since this is the main mechanism of interaction between propeller blades. Historically, the modeling of the free boundary problem for wings started with the PWS method [7], which is based on the assumption of a given wake shape, which can be guessed using either simple linearized models (small perturbation velocities) or experimental information. With the increase in computing power, the wake relaxation method (WRM) has been applied to the solution of steady state flow problems. According to this method some initial wake geometry is assumed for which wake grid panels are defined. Panel geometry is then deformed in a successive iteration scheme until it becomes a flow surface for the steady problem. Finally, for unsteady problems the time-stepping method (TSM) has been applied for the case of flow around wing configurations. In this case the grid representing the free wake surface is dynamically evolving with time [7]. We shall now give a brief review of the state of the art in unsteady calculations of propellers.

634

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

After the publication of the pioneering work of Hess and Valarezo [1] a large number of papers have appeared in the ship hydrodynamics community utilizing different forms of BEM and free wake modeling. Kerwin et al. [8] present a Morino type [9] BEM for the solution of the steady flow problem of a propeller in a duct. They also discuss the applicability of different BEMs for propeller flow problems. Taking into account that blade sectional thickness to chord ratios vary from 20% at blade root to 2% at tip region, they conclude that the Morino type formulation is more efficient due to the lower order of singularity in the kernel of the corresponding integral equation. Trailing vortex surfaces emanating from blade and duct trailing edges have been modeled using PWS in the form of the classical helicoidal model. Belibasakis and Politis [10, 11] present a BEM for the unsteady propeller flow problem based on surface vorticity distributions. The position of the trailing vortex wake in this formulation is assumed known from experimental observations (PWS) [12]. The WRM has been applied to propeller flow problems by a number of authors. Kinnas [13] presents a Morino type (potential based) BEM for the solution of unsteady flows around propellers. The position of the trailing vortex wake is calculated using a WRM with a circumferentially averaged mean flow. Thus resulting vortex geometry is invariant with time. Maitre and Rowe [14] present a similar formulation for the steady state problem. Takinaci [15,16] presents an improved WRM method for the calculation of the steady state vortex wake geometry including roll up. Finally Kinnas [17,18] presents a number of roll up models for the steady case utilizing PWS and WRM. As indicated by Kinnas [18], if his numerical calculation was performed up to a large distance from the blade, the predicted roll-up could eventually become unstable. In the present paper the problem of flow of a propeller performing an unsteady translation with instantaneous ~ velocity vector VðtÞ and unsteady rotation with instantaneous rotational velocity V~ðtÞ is formulated and solved for an observer in an inertia (motionless) coordinate system. A Morino type BEM is used, together with a vortex filament approximation to simulate the trailing vortex system produced by the propeller blades, and a TSM to calculate the time evolution of the propeller unsteady motion. The developed scheme is used to simulate: (a) a propeller moving with steady translational and rotational velocity, (b) a propeller in inclined flow at various angles of inclination and (c) a propeller in heaving motion. Results are presented showing the grid independence of the developed scheme as well as the stability of the free wake patterns.

2. Formulation 2.1. Kinematics Consider a propeller with z blades operating in an unbounded fluid, performing an unsteady translation with ~ and an unsteady rotation with instantaneous velocity VðtÞ

Fig. 1. Geometry of the problem.

instantaneous rotational velocity V~ðtÞ: Let Oxyz denote an inertia (motionless) coordinate system and let O0 x0 y0 z0 denote a system moving with the propeller, Fig. 1. Let also ~rs ðtÞ denote the relative position of the second system with respect to the S first. Let L ¼ i¼1;z Li denote the propeller volume as union of volumes of the z blades, L , R3 : Let ›Li ; i ¼ S 1; …; z denote the bounding surface of each blade and ›L ¼ i¼1;z ›Li : Let li [S›Li denote the trailing edge line on each blade and l ¼ i¼1;z li : Consider also z surfaces Wi ; i ¼ 1; …; z embedded in R3 with boundaries (lines) ›Wi ; where ›Wi > ›Li ¼ li ; i.e. part of the boundary of ›Wi lays on Li and coincides with the trailing edge. We shall call Wi the ‘trailing vortex sheet’Sor ‘free wake’ emanating from the ith blade. Let also W ¼ i¼1;z Wi : Assuming inviscid, incompressible and irrotational flow at points P [ V þ ; R3 2 L 2 W and introducing a virtual inviscid, incompressible and irrotational flow for the region V 2 ; L (inside blades), a potential function exists for the perturbation velocity ~v created by the propeller movement in both regions V þ and V 2 (~v ¼ 7w; where w is the perturbation potential). Let wþ and w2 be the potential at V þ and V 2 ; respectively. Let also wþ and w2 denote the potential in the two sides of trailing vortex surface W (i.e. W is considered a two sided surface). Then the following representation theorem holds for the potential function, Katz and Plotkin [7] 1 ð s 1 ð n~ ·~r dS þ wðPÞ ¼ 2 m 3 dS ð1Þ 4p ›V r 4p ›V r where s and m denote the source and dipole intensities, respectively, defined on points of the surface ›V ; ›L þ W and given by the following relations

s ¼ n~·7ðwþ 2 w2 Þ

ð2aÞ

m ¼ w2 2 wþ

ð2bÞ

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

where n~ ; n~þ ; Fig. 1 denotes a unit normal vector to the ! bounding surface and r ¼ PQ ; P denotes the (constant) control point and Q the variable (integration) point. Note that Eq. (1) holds uniformly for all points P [ R3 2 ›L 2 W: Applying Morino ideas [9], an integral equation for the perturbation potential due to propeller movement can be obtained by requiring a homogeneous internal condition

w2 ðPÞ ¼ 0

ð3Þ

i.e. the perturbation potential inside blades should be zero. In this case due to a well known theorem of potential theory, Gunter [19], the perturbation velocity everywhere inside the moving body is zero. Thus

w ðP [ V Þ ¼ 0 ) 7w ðP [ V Þ ¼ 0 ) ~v ¼ 0 2

2

2

2

2

ð4Þ

The no penetration boundary condition on solid body surfaces requires that the relative fluid-body velocity should be zero. A point P [ L moves relative to the inertia Oxyz system, with velocity ! ! ! ð5Þ qðtÞ ¼ VðtÞ þ VðtÞ £ ð~r0 2 ~rs Þ Assuming that no flow at infinity exists, the no penetration boundary condition on body surface ›L becomes n~·7wþ ¼ n~ ·~q

ð6Þ

Substituting relations Eqs. (4) and (6) in relation (2a), we get þ

s ¼ n~·ð7w 2 7w Þ ¼ n~·~q 2

ð7Þ

635

Note that time derivative of the potential function is evaluated in the inertial coordinate system. Introducing the perturbation velocity jump by the relation v2 k~vm l ; ~vþ m2~ m

ð10Þ

and the mean vortex sheet perturbation velocity by the relation ~vþ þ ~v2 m l~vm k ; m ð11Þ 2 and substituting relations (2b), (10) and (11) to relation (9), we get

›m D m þ l~vm k·7m ¼ 0 , m ¼ 0 ›t Dt

ð12Þ

where by Dm =Dt ; ›=›t þ l~vm k·7 we denote a substantial derivative based on the mean perturbation velocity at points of the trailing vortex sheet. Relation (12) indicates that trailing vortex sheets are moving with the mean perturbation velocity while dipole intensity m is a material property when mean velocities are considered. Let ~xðu; v; tÞ denote the equations describing the geometry of a time evolving vortex sheet parameterized by u; v (do not confuse curvilinear coordinate v with perturbation velocity ~v). The time evolution of the vortex sheet can then be calculated by integrating in time the relation d~x ¼ l~vm k dt

ð13Þ

Substituting relation (7) in representation theorem (1) and taking the limit as P2 ! P0 with P [ V 2 ; P0 [ L; we get (we assume smooth boundaries) 1 1 ð 1 ð n~·~q n~ ·~r dS ð8Þ mþ m 3 dS ¼ 2 4p ›V r 4p ›V r

To find the mean perturbation velocity l~vm k on the vortex sheet we take the gradient of relation (1) and then take the limits as Pþ ! P0 [ W and P2 ! P0 [ W; i.e. from both sides (þ ) and (2 ) of W: We then take the mean value of the calculated limits. This procedure results in the following relation

which is a Fredholm integral equation of the second kind with a weakly singular kernel. Existence and uniqueness of solutions for smooth domains for such type of equations can be found in Kress [20] and Mikhlin [21].

l~vm k ¼ ~vðP [ WÞ ~r ~r 1 ð 1 ð s 3 dS 2 g~ £ 3 dS ¼2 4p ›V r 4p ›V r ! 1 ð dl £ ~r 2 m£ 4p ›W2l r3

2.2. Dynamics The domain space for S m in Eq. (8) consists of two parts. The first part is ›L ¼ i¼1;z ›Li (blades), which is a solid boundary with known motion, Eq. (5). TheSsecond part is consisting of the trailing vortex sheets W ¼ i¼1;z Wi which carry the generated vorticity at the blade trailing edges to the downstream region satisfying Helmholtz’s equations, Kochin et al. [5]. An equivalent form of the Helmholtz dynamic equations for the case of free vortex sheets and in terms of m can be found as follows: Applying the unsteady Bernoulli equation for two neighboring points Pþ and P2 taken on the two sides W þ and W 2 of a trailing vortex sheet and in the vicinity of a point P0 [ W; we get

› wþ 1 ›w 2 1 þ rl~vþ l2 ¼ þ rl~v2 l2 2 2 ›t ›t

ð9Þ

ð14Þ

where the last integral is a line integral over the free part of the line bounding W (i.e. ›W 2 l) and g~ is the surface vorticity intensity vector defined by the equation

g~ ; n~ £ ð~vþ 2 ~v2 Þ ¼ n~ £ ð7wþ 2 7w2 Þ

ð15Þ

with n~ ; n~ ; Fig. 1. Derivation of relation (14) is somewhat cumbersome in Hess [3]. A more elegant proof of Eq. (14) is given in Appendix 1. Surface integrations in Eq. (14) are taken in the Cauchy principal value sense. General existence theory for such types of integrals can be found in Mikhlin [21]. Eq. (12) can also be used to calculate the rate of generation of vorticity at points on the trailing edges. The problem of the introduction of a proper ‘generating vorticity relation’ in BEM formulations of lifting flows has been þ

636

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

extensively discussed, e.g. by Sarpkaya [22], Cottet and Koumoutsakos [6] and Katz and Plotkin [7]. The corresponding relation is usually termed ‘Kutta condition’. By applying relation (12) at points on the trailing edge for two successive time steps and taking into account that a material point P [ l at time t would travel with the mean flow velocity to become a point P0 [ W at time t þ dt; we get

wu ðtÞ 2 wl ðtÞlP ¼ mðt þ dtÞlP0

ð16Þ

where wu ðtÞ; wl ðtÞ denote perturbation potentials at points Pu ! P; Pl ! P; i.e. Pu and Pl approaching P [ l from different sides of l; and ~rP0 ¼ ~rP þ l~vm k·dt

ð17Þ

Relation (16) is the type of Kutta condition, which shall be used in our formulation. 2.3. Pressures, forces and moments Bernoulli equation for an observer in the moving coordinate system O0 x0 y0 z0 can be found in Katz and Plotkin [7] p 2 p1 1 ›w ¼ 2 ð7wÞ2 2 ~vref ·7w 2 2 r ›t

ð18Þ

where: ~vref ¼ 2~q given by relation (5), and time derivative of potential refers to the moving system. Forces and moments can then be calculated using ð ! ð ! ~¼ ~ ~r1 £ðp·~n þ D Þ·dS F ðp·~n þ DÞ·dS; M¼ ð19Þ ›L ›L ! where D denotes a viscous correction drag force given by ~ ¼ 1 CD rl~vtot l~vtot D ð20Þ 2 and ~vtot is total velocity relative to an observer in the moving system, given from ~vtot ¼ ~vref þ ~v

ð21Þ

and CD is the skin friction coefficient of a smooth plate in turbulent flow given by, Blevins [23] 0:0986 CD ¼  ð22Þ 2 Ux 2 1:22 log10 n where Ux=n is the local Reynolds number at the referenced radial position of the blade.

Fig. 2. Grid points for blades and free wake. The ‘Kutta-strip’ can be identified by the chord-wise size of its panels which is different from that of the blades.

bilinear elements have been selected for this purpose, Fig. 2. The same type of elements have been selected for discretizing the Sdynamically evolving trailing vortex sheet W : Wi ¼ m¼1;M EWi;m ; where EWi;m denotes the mth element on the ith trailing vortex sheet (each blade has its own trailing vortex sheet) and M is the total number of elements at time t used to describe the trailing vortex sheet of each blade. Note that number of elements describing W increases with time while their position should be calculated as part of the solution, through relation (13). In this respect, at each time step, a new ‘strip’ with K elements (termed also ‘kutta strip’) directly adjacent to each trailing edge li is created, Fig. 2. Let EWi;k ; k ¼ 1; K denote the elements of the kutta strip of the ith blade. For reasons of grid continuity K coincides with the number of elements used to describe each blade in the radial direction. 2. With the assumptions of the previous paragraph, integral Eq. (8) at each time step t takes the following descritized form 1 1 ð n~·~r mþ m dS 2 4p ›LðtÞþkutta_stripðtÞ r 3 1 ð 1 ð n~ ·~q n~ ·~r dS 2 ¼ m dS 4p ›LðtÞ r 4p WðtÞ2kutta_stripðtÞ r 3 ð 1 1 X n~·~r ) mi0 ;n0 þ mi;n dS 3 2 4p ELi;n r i¼1;z n¼1;N

3. Description of the time stepping algorithm The following natural steps have been used in materializing the solution algorithm of a propeller in general 3D unsteady motion 1. Discretization of the boundary Li ; i ¼ S 1; z of each blade into N boundary elements: ›Li ¼ n¼1;N ELi;n ; where ELi;n denotes the nth element of the ith blade. Four node

þ

ð 1 X n~ ·~r mi;k dS 3 4p EWi;k r i¼1;z m¼1;K

¼

1 X ð 1 n~·~q dS 2 4p 4p r ELi;n i¼1;z n¼1;N

X

ð

m EWi;m

n~ ·~r dS r3

i¼1;z m¼Kþ1;M

ð23Þ

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

where i0 ; n0 denote control point position, i0 ¼ 1; z; n0 ¼ 1; N: Constant dipole intensity mi;n for both blade and trailing vortex sheet region has been used for the discretization of Eq. (8). Regarding source distribution (first term in the right hand side), its intensity is analytically known at each time step, relation (7), and thus the surface integral can be calculated numerically. Considering the second term on the right hand side, both m and the geometry of its support W (trailing vortex sheets) are assumed to be known from previous time steps (see fourth step). As shown in Eq. (23) a kutta-strip is used as an additional unknown at each time step. Dipole intensities are then calculated by satisfying Eq. (23) at the z £ N centroids of blade elements simultaneously with Eq. (16) relating dipole intensity of kuttastrip elements with dipole intensities of the directly adjacent body elements (z £ K more equations). Thus body and kutta-strip dipole strengths can be calculated at each time step. 3. Using a constant dipole intensity, formula (14) for the calculation of the mean perturbation velocities on free vortex sheet, degenerates to ~vðP [ WÞ 1 X ¼2 4p

ð

i¼1;z n¼1;N

1 2 4p

X ð i¼1;z m¼1;M

! ! ð ~r dl £ ~r s 3 dS þ mi;n £ dS r3 r ›ELi;n ›ELi;n

! dl £ ~r mi;m £ r3 ›EWi;m

ð24Þ

where source and dipole intensities and integration surface W are known from previous time steps (see fourth step). More specifically the first term is a surface integral with regular kernel for points P [ W while the second and third terms in the right hand side are line integrals of known dipole strengths (second step) over the boundary line of each element. Note that these integrals are equivalent to the well known Biot –Savart formula for the calculation of induced velocities from singular vortex lines, Katz and Plotkin [7]. Thus the integrals on the right hand side of Eq. (24) can be calculated and mean induced velocities on the trailing vortex sheet can be found. In our implementation we use relation (24) to calculate velocities directly at grid points of the grid representing the free wake of each blade. Since Biot – Savart formula gives infinite self-induced velocities a smoothing technique is necessary to desingularize Eq. (24), Sarpkaya [22]. Details of the scheme used in the present paper are given in Section 4 (see numerical implementation of the method). 4. With the induced perturbation velocities known at each grid point ~xi;l ; i ¼ 1; z l ¼ 1; LðtÞ; where LðtÞ denotes the number of nodes representing the free wake of each

637

blade at the current time step, the new position of the nodes can be calculated using relation (13), which in discretized form becomes1 ~xnew xi;l þ ~vðP ; ~xi;l Þdt; i;l ¼ ~

i ¼ 1; z ð25Þ

l ¼ 1; LðtÞ where x~new i;l denotes the new position of the i; l grid point and ~vðP ; ~xi;l Þ denotes the mean perturbation velocity calculated at this point using Eq. (24) (third step). Note again that dipole strengths are invariant quantities during the motion with the mean perturbation velocity, relation (12). 5. Knowing the dipole strengths on the blade surface at each time step (second step), perturbation velocities can be calculated by taking the surface gradient of the dipole strength ~v ¼ 7wþ ¼ 27s m

ð26Þ

where 7s is used to denote the surface gradient. Relation (26) is numerically implemented using a central difference scheme. Details are given in Section 4. With perturbation velocities known blade surface pressures, and blade and propeller forces and moments can be calculated at each time step, using Eqs. (18) and (19). 6. As a final step to our algorithm we move the propeller as a rigid body to its new position and return to step 2. Calculation of the new propeller position can be achieved easily as a composition of the known ~ translation with instantaneous velocity VðtÞ and an instantaneous rotation with rotational velocity V~ðtÞ:

4. Numerical implementation 4.1. Grid generation Grid generation is the necessary starting point for every BEM. Since we are using a TSM method the trailing vortex sheet grid is produced as part of the solution. So grid generation is necessary only for the propeller solid boundary. This is a serious advantage of the method in comparison to PWS or WRM methods, where discretization of the trailing vortex sheets is also necessary. The coordinates of the discretized surface should be sufficiently accurate since any inaccuracies can lead the hydrodynamic pressure to become noisy. According to our experience accuracy in six significant digits seems to be necessary in order to have satisfactory results in grid description. A detailed discussion of propeller geometry is given by Klein, [24]. There are two physical directions in blade surface along which nodes have 1

This is a first-order TSM. Although higher-order Runge-Kutta integration schemes are expected to improve accuracy, we make this choice to simplify computer implementation of the method.

638

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 5. Equal-s; equal-c surface grid 20 £ 10. Fig. 3. Equal-r; equal-c surface grid 20 £ 10.

to be distributed, a radial (or spanwise) direction and a chordwise direction. The most natural (and also common) selection in discretizing blade surface is to use equally spaced radial intervals and equally spaced or cosine spaced chordwise intervals. Figs. 3 and 4 present the grid representation of a highly skewed blade subdivided in 20spanwise £ 10-chordwise boundary elements using equally spaced radial intervals (equal-r) and either equally spaced or cosine spaced chordwise intervals (equal-c or cos-c). By carefully looking at those figures we observe that cosine chordwise spacing is better since more nodes are distributed at the leading edge where blade geometry changes abruptly. Regarding radial spacing, equally spaced radial intervals suffer at the blade tip region. Facing this problem we have developed an alternative grid generation scheme in which the spanwise and chordwise curvilinear distances at the surface of the blade are evenly subdivided to produce the grid points (equal-s scheme). This can be achieved using the following steps: (a) Produce a usual equally spaced radial-chordwise grid as that shown in Fig. 3. (b) Calculate the curvilinear lengths along chordwise direction and redistribute evenly or according to a given pattern (e.g. cosine pattern) the grid nodes. (c) Calculate the spanwise curvilinear lengths for the grid points found in bth step and redistribute evenly or according to a given pattern the nodes at this direction.

Fig. 4. Equal-r; cos-c surface grid 20 £ 10.

(d) Return to ath step and repeat the process. The proposed algorithm usually converges in two iterations. The result of applying this idea is shown in Figs. 5 and 6. Comparing Figs. 3 and 4 with Figs. 5 and 6 we conclude that for the same number of elements much better description of the propeller surface is possible. Finally note that all necessary surface interpolations are performed using parametric cubic splines. For a propeller operating at off design conditions a sudden pressure change is expected near the leading edge region. In such cases the use of cosine spacing along chordwise direction with evenly spaced nodes in the spanwise direction (equal-s) is preferable such as in Fig. 6. Our grid generation routine allows someone to decide regarding the type of subdivision independently in each direction. It can also produce grids for propeller with a hub, Fig. 7, or a propeller in a duct. Grid generation for the hub and its intersection with the propeller blades is a more complicated problem. Since the present work is mainly directed to the simulation of unsteady propeller blade motion and the formation of the trailing vortex wake we shall leave the discussion of hub grid generation for a future paper. 4.2. Treatment of integrals Almost in every BEM the researcher is faced with the calculation of some type of singular integrals.

Fig. 6. Equal-s; cos-c surface grid 20 £ 10.

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

639

error bounds rather than relative, as in Ref. [24]. Furthermore for the evaluation of the improper integrals for dipoles and sources the cutoff radius used is cutoff_radius ¼ 0:01·minðside_u; side_vÞ

Fig. 7. Hub grid.

The appearance of the singularity in the kernel is always connected with the relative position of the control point P and the integration point Q; Fig. 1. In our formulation we have two types of singular integrals † Weakly singular integrals (improper integrals), appearing in the calculation of element self induction factors (or self influence coefficients) due to dipoles (left hand side of Eq. (23)) and in the calculation of the self part of the source distribution term (first term in the right hand side of Eq. (23)). † Non-integrable (at P ; Q) in the usual sense line integrals of the Biot – Savart type, appearing in the right hand side of relation (24) (second and third terms). With the selected bilinear boundary elements the integration range of those integrals is a straight line (side of the bilinear element). Thus integrations as well as treatment of singularity can be performed analytically as will be explained in the sequel. The first term in the right hand side of Eq. (24) is regular since this integral is always evaluated on the trailing vortex wake whereP – Q: (otherwise this term presents a strong Cauchy type surface singularity as P ! Q; Mikhlin [21]). The numerical calculation of surface improper integrals is a simple procedure and can be achieved using common numerical integration rules and a cutoff function symmetric around the singularity (cutoff radius). In our implementation all improper integrals are calculated using adaptive Simpson quadrature. According to this the integrals are first transformed to iterated and then the integrations (in each direction) are performed using adaptive Simpson quadrature, Press et al. [25]. The only difference between our quadrature method with that presented in Ref. [24] is in the convergence criterion of the integrals where we use absolute

ð27Þ

where, side_u and side_v are the lengths of the two sides of the element for which the self integration is performed. Propeller blade is a ‘difficult’ boundary with length scales changing abruptly in the spanwise direction. Thus at hub region blade thickness is comparable to chord and thus boundary elements at opposite sides of the blade are relatively ‘far’. At the tip region thickness/chord tends to zero. Thus boundary elements at opposite sides of the blade are very near. Adaptive integration has been proved to be a fast and reliable method which behaves well in all these regions. Regarding numerical handling of Biot – Savart type integrations (second and third terms in relation (24)), consider a straight line vortex of unit strength starting at ~ Let also the position of the point a~ and ending at point b: ~ control point be denoted by p: Then the induced velocity ~v at p~ from the line vortex defined by a~; b~ can be calculated analytically and is given by the following relation (for details see Appendix 2) ! ~v ¼ s2·ð~l £ a1 Þ=ð4pÞ ð28Þ ! where s2; ~l and a1 can be calculated as follows ! ! ! a1 ¼ a~ 2 p~; b1 ¼ b 2 p~ ð29Þ ! ! ! ~l ¼ ð ! b1 2 a1 Þ=d; d ¼ l b1 2 a1 l ð30Þ ! ! s0 ¼ l a1 l2 ; s1 ¼ ~l· a1 ð31Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f ¼ d2 þ 2·d·s1 þ s0; g ¼ s0 2 s12 ð32Þ pffiffiffi pffiffiffi s2 ¼ ðs1·f 2 s0·ðd þ s1ÞÞ=ð s0·f·gÞ ð33Þ ~ induced For a non-degenerate vortex line (i.e. a~ – b), velocities calculated from relation (28) can become infinite only in the case when g!0

ð34Þ

This observation indicates that the scalar variable g is the proper cutoff parameter when relation (28) is used. It is very interesting that using g as the cutoff parameter very stable patterns of unsteady wake rollup are predicted with large variations in the cutoff selection. 4.3. Calculation of potential surface gradient Let fn21 ; fn ; fnþ1 denote the values of the potential at three consecutive nodal points P1 ; P2 ; P3 (centroids of elements) on blade surface along either u (chordwise) or v (spanwise) directions. Let also l1 ¼ lP1 P2 l; l2 ¼ lP2 P3 l be the curvilinear physical distances of P1 ; P2 ; P3 along either u or v direction. These distances can be easily calculated using the element shape functions and the surface metric tensor.

640

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 8. Rollup pattern for propeller DTMB P4119 at angle ¼ 1808, J ¼ 0:833; angular step angle ¼ 108.

Derivative at P2 is calculated using the following formula (weighted central difference)

ev are unit vectors along the element local u where ! eu and ! and v directions. ! eu and ! ev can be calculated using

›f f 1·l2 þ f 2·l1 ¼ ›lu or v l1 þ l2

~eu ¼

ð35Þ

fn 2 fn21 ; l1

f2 ¼

fnþ1 2 fn l2

and lu or v denotes the physical length along direction u or v: At boundary points of ›Li a usual forward or backward difference scheme has been used. By applying relation (35) along u and v directions, the corresponding physical potential derivatives can be found. Tangential perturbation velocity can then be calculated using

›f ›f ~v ¼ ! eu þ! ev ›lu ›l v

~ev ¼

d~r=dv ld~r=dvl

ð37Þ

where ~r is the position vector at point u; v of the element. Some details regarding computer implementation of the method, using object oriented features of FORTRAN 95, can be found in Appendix 3.

where f1 ¼

d~r=du ; ld~r=dul

ð36Þ

5. Numerical examples 5.1. Impulsively started propeller As a first example we shall apply our method to an impulsively started propeller, that is a propeller which at t ¼ 0 is starting to move with translational steady velocity V~

Fig. 9. Rollup pattern for propeller DTMB P4119 at angle ¼ 3608, J ¼ 0:833; angular step angle ¼ 108.

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

641

Fig. 10. Rollup pattern for propeller DTMB P4119, 20 £ 15 blade grid, J ¼ 0:833; angular step angle ¼ 88.

and rotate with steady rotational velocity V~: For our example we shall use the model propeller DTMB 4119, for which detailed measurements of blade pressure distributions and total axis forces exist, Jessup [27], and which has been used in the past, Gindroz et al. [28], as a benchmark for evaluating boundary element codes. Basic geometric characteristics of this propeller are given by Jessup [27]. Results shall be presented at two different propeller advance coefficients J ¼ 0:833 and J ¼ 0:5 (J ; V=nD; n denotes the propeller rotational speed and D the propeller diameter). Advance coefficient J ¼ 0:833 is the propeller design point where blade sections are at the ideal (shock free) angle of attack. When reducing the advance coefficient below this ideal value, the loading of the propeller increases and the rollup of the free wake is expected to increase as well. To evaluate the grid convergence of our algorithm, results shall be given for different grid densities and different values of the angular rotational step. Fig. 8 shows the rollup patterns of the trailing vortex sheet for a propeller at 1808 angular position (measured with regard to its position at t ¼ 0) and for two different grid arrangements: (a) a 15 £ 10 grid (i.e. 15 elements spanwise and 10 elements chordwise at each side of the blade) and (b)

a 20 £ 15 grid. Fig. 9 shows similar results for a propeller at 3608 angular position. Both figures refer to the design advance coefficient J ¼ 0:833: We see a strong starting vortex shed from each blade and travelling downstream with the perturbation velocity creating strong spiral type regions at the back of the wake for each blade. Furthermore strong hub and tip vortices are continuously shed from the respective blade regions and ‘absorb’ the weaker vorticity regions at inner blade radii producing also spiral type patterns. Comparing the rollup patterns in the right and left part of Figs. 8 and 9 we observe that they are practically identical. Figs. 8 and 9 have been produced with a 108 angular step. To examine the effect of the angular step in the rollup pattern, calculations have been performed for a smaller step of 88 and for a 20 £ 15 grid and for a greater step of 158 with a coarser grid of 8 £ 5. Fig. 10 shows the corresponding rollup patterns for the first case while Fig. 11 for the second case. Comparing Figs. 10 and 11 with Figs. 8 and 9 we observe that rollup patterns are practically the same, independently of the angular step. Furthermore systematic calculations with various grid densities and angular steps, not given here, have

Fig. 11. Rollup pattern for propeller DTMB P4119, 8 £ 5 blade grid, J ¼ 0:833; angular step angle ¼ 158.

642

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 12. Propeller thrust coefficient as a function of propeller angular position, Propeller DTMB P4119, J ¼ 0:833:

shown that the predicted rollup pattern is very stable independent of the selection of the discretization parameters. Fig. 12 presents the thrust coefficient Kt of the propeller as a function of its angular position for the four different combinations of grid densities and angular steps (Kt ; T=ðrn2 D4 Þ; T denotes the propeller thrust and r is the fluid density). From this figure we observe that stable numerical results are obtained with convergence achieved from a grid density of 15 £ 10 and up. In the same figure the steady state experimental thrust coefficient is shown, Jessup [27]. Comparing the experimental thrust coefficient with the calculated one, we observe that the present method predicts

the experimental thrust with error less than 1%. From the same figure we observe that the impulsively started propeller reaches a steady state thrust at less than half a revolution. Figs. 13– 15 present results of the chordwise pressure distributions for the DTMB 4119 propeller at the design advance coefficient and at three radial positions: r=R ¼ 0:3; 0.7 and 0.9. Our results refer to the propeller at angle 4008. Pressure coefficients in those figures are nondimensionalized with the total undisturbed velocity: CP ;

p 2 p1 : 1 rl~ql2 2

Fig. 13. Chordwise pressure distribution for propeller DTMB P4119, J ¼ 0:833; Propeller at angle 4008, r=R ¼ 0:3:

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

643

Fig. 14. Chordwise pressure distribution for propeller DTMB P4119, J ¼ 0:833; Propeller at angle 4008, r=R ¼ 0:7:

More specifically Figs. 13 –15 compare pressure distribution obtained with the present method with experiments, Jessup [27], and calculations made by Hoshino [28]. The present method seems to predict very satisfactory the pressure distribution at r=R ¼ 0:3 and r=R ¼ 0:7: At r=R ¼ 0:9 good predictions are obtained at blade face (pressure side) while fair results are obtained for the suction side at leading and trailing edge regions. This can be attributed to complex 3D viscous flow at the blade tip region which

cannot be modeled effectively using BEM. Note that differences in Cp at the tip region do not affect seriously propeller thrust since blade loading tends to zero there. Fig. 16 shows the rollup patterns of the trailing vortex sheet for the P4119 propeller operating at advance coefficient J ¼ 0:5: Comparing the rollup patterns with the corresponding patterns at J ¼ 0:833; we observe that a stronger spiral type deformation is produced from the starting vortex as well as from the hub and tip vortices

Fig. 15. Chordwise pressure distribution for propeller DTMB P4119, J ¼ 0:833; Propeller at angle 4008, r=R ¼ 0:9:

644

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 16. Rollup pattern for propeller DTMB P4119, 20 £ 15 blade grid, J ¼ 0:5; angular step angle ¼ 108.

(as expected). From the same figure we also observe that starting vortices affect the structure of the vortex wake since they approach its upstream portion using their own induced velocities. They also seem to create something like a vortex ring which accelerates the flow in the propeller region. We believe that this phenomenon, that is the creation of an axial jet due to the redistribution of wake vorticity, is of crucial importance for correctly predicting the propeller developed thrust. Fig. 17 shows the calculated propeller thrust coefficient as a function of propeller angular position for three

different grid densities. From this figure we again confirm the grid independence of the numerical results for relatively coarse grids. By a more careful observation at this figure we detect a small instability in the calculated Kt for the 15 £ 10 grid. Such instabilities can occur due to a ‘resonance’ between the shed vorticity mechanism at the trailing edge and the developed circulation around blade sections as a result of satisfying a Kutta condition. Usually by changing the grid density and/or the angular step this instability disappears. In the same figure it is shown that the experimental (steady state) thrust

Fig. 17. Propeller thrust coefficient as a function of propeller angular position, Propeller DTMB P4119, J ¼ 0:5:

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

645

Fig. 18. Pressure distribution ðCp Þ for propeller DTMB P4119, J ¼ 0:5; Propeller at angle 4008.

coefficient is reached in less than half a revolution. Observe also that the present calculation method predicts the experimental steady thrust with error less than 1%. Finally Fig. 18 shows the calculated pressure distributions on blade face and back at J ¼ 0:5: 5.2. Propeller in inclined flow As a second example we shall apply our method to a flow around an inclined axis propeller. We shall use the model propeller DTMB 4718 for which detailed measurements of blade pressure distributions exist at 08

and 7.58 inclination angles, Jessup [29]. We shall perform calculations at 0 and 7.58, for comparison with experimental results, and for 108, 208, 308 and 408 to investigate the effect of angle inclination on the performance of the propeller. For the case of 08 inclination, calculations should be performed for three different advance coefficient values J ¼ 0:751 (design), J ¼ 0:6 and J ¼ 0:4: Fig. 19 shows the thrust coefficient Kt of the propeller as a function of its angular position at the three selected advance coefficients and at propeller inclination 08. The experimental Kt for the design J is also shown on this figure. We observe that the present

Fig. 19. Propeller thrust coefficient as a function of propeller angular position with the advance coefficient J as parameter, Propeller DTMB P4718.

646

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 20. Chordwise pressure distribution for propeller DTMB P4718, J ¼ 0:715; Propeller at angular position 4008, r=R ¼ 0:7:

method predicts the experimental thrust with error less than 1%. Figs. 20 – 22 present comparisons of the calculated chordwise pressure distributions with experimental results at the design advance coefficient and at radial position r=R ¼ 0:7 at 08 inclination. The present method again seems to predict satisfactorily the pressure distributions. Figs. 23 and 24 present comparisons of first harmonic amplitude of pressure coefficients with the experimental ones for 7.58 inclined propeller and at r=R ¼ 0:7 and 0.9. Some differences exist in this case but we have to note that unsteady pressure measurements on model propellers are a challenging task with errors not

entirely under control. Fig. 25 shows the propeller thrust coefficient for different propeller inclinations. We observe a non-linear variation of propeller thrust with inclination. Also at inclinations greater than 208 a small oscillation in the calculated trust is starting to occur. This oscillation can be attributed to the increasing asymmetry of the free wake with regard to propeller axis (Fig. 27). Furthermore to investigate the non-linear thrust increment with respect to the axis inclination we made some additional calculations at 08 angle inclination but with reduced advance coefficient according to the relation J ¼ cosðincl: angleÞ·Jdesign ðJdesign ¼ 0:751Þ: This relation is usually

Fig. 21. Chordwise pressure distribution for propeller DTMB P4718, J ¼ 0:6; Propeller at angular position 4008, r=R ¼ 0:7:

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 22. Chordwise pressure distribution for propeller DTMB P4718, J ¼ 0:4; Propeller at angular position 4008, r=R ¼ 0:7:

Fig. 23. Chordwise first harmonic amplitude pressure coefficient with propeller DTMB P4718, inclined 7.58, r=R ¼ 0:7; J ¼ 0:715:

Fig. 24. Chordwise first harmonic amplitude pressure coefficient with propeller DTMB P4718, inclined 7.58, r=R ¼ 0:9; J ¼ 0:715:

647

648

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 25. Propeller thrust coefficient as a function of propeller angular position with the propeller axis inclination angle as parameter, Propeller DTMB P4718, J ¼ 0:751:

used in practical propeller calculations to represent the effect of inclination. It can be derived simply by taking the velocity vector triangle at the propeller at inclined condition. The results of calculations are shown in Fig. 26. We observe that at smaller angles the use of the reduced J in calculations can correctly predict the thrust coefficient. Greater differences occur at inclinations over 208 which are not of practical interest. Fig. 27 shows the rollup patterns of the trailing vortex sheet for the P4718 propeller operating at advance coefficient J ¼ 0:751: More specifically Fig. 27e and f show the dipole strength distribution at 408 inclination angle. Observe the characteristic asymmetry which is equivalent to an increase in blade loading at 10 o’clock position. Fig. 28 shows the single blade thrust force for the propeller at various inclinations where for the reading of the horizontal axis of this figure we have use the convention that angular position 08 denotes blade at 12 o’clock. Observe that although axis forces are steady, Fig. 23, single blade forces oscillate at axis frequency with amplitude increasing with axis inclination. Note also that maximum and minimum of single blade force occur at the positions of about 908 and 2708 (as expected) with a small shift to the right with increasing inclination angles due to the free wake asymmetry, Fig. 27.

coefficients: (a) J ¼ 0:751 corresponding to a propeller translating with a velocity v0 ¼ 3:756 m/s and rotating with frequency n ¼ 8:2 s21 and (b) J ¼ 0:4 corresponding to v0 ¼ 2 m/s and n ¼ 8:2 s21. The propeller performs a vertical oscillation (heaving) according to the following relation   D vt yðtÞ ¼ sin 4 2 where, D ¼ 0:61 m is the propeller diameter and v the angular frequency of propeller rotation. Thus in two propeller turns one full heaving oscillation occurs. Fig. 29 contains rollup patterns for the heaving propeller at the two advance coefficients. Figs. 30 and 31 show the single blade forces and Fig. 32 the total propeller thrust at the two advance coefficients. In Fig. 32 we have also plotted

5.3. Propeller in heaving motion As a final example of the use of the developed method we shall present unsteady calculations for the propeller DTMB 4718 in heaving motion. This operating condition could happen when a ship is travelling in waves. Calculations have been performed at two advance

Fig. 26. Comparison of thrust coefficient in inclined axis condition with the corresponding thrust coefficient calculated with properly reduced J; Propeller DTMB P4718.

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

649

Fig. 27. Rollup patterns for propeller DTMB P4718, 10 £ 10 blade grid, J ¼ 0:751; angular step angle ¼ 108, at various propeller axis inclination angles.

the corresponding steady propeller thrust at the same advance coefficients for comparison. From Figs. 30 and 31 we observe that due to heaving motion a serious asymmetry of the instantaneous thrust produced by each blade occurs. Finally Fig. 32 shows

that the propeller thrust is also unsteady with a minimum value equal to the corresponding steady (w/o heaving) and extremum values occurring at angles nearly equal to 1808 and 3608 with a mean value greater than the corresponding w/o heaving. This can be explained by

650

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Fig. 28. Single blade thrust coefficient as a function of blade angular position with the propeller axis inclination angle as parameter, Propeller DTMB P4718, J ¼ 0:751:

Fig. 29. Rollup patterns for propeller DTMB P4718 in heaving motion, 10 £ 10 blade grid, angular step angle ¼ 108.

Fig. 30. Single blade thrust coefficient as a function of blade angular position for the propeller at heaving motion, Propeller DTMB P4718, J ¼ 0:751: (angular position 08: blade no 1 at 12 o’c lock).

the fact that the existence of heaving tends to increase the instantaneous rotational velocity of the propeller blade sections while extremum values for thrust coincide with corresponding extremum values of the heaving velocity as expected.

6. Conclusions An efficient time stepping method is presented for the simulation of the general unsteady motion of a propeller in a 3D infinite fluid. Unsteady free wake rollup and contraction

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

651

Fig. 31. Single blade thrust coefficient as a function of blade angular position for the propeller at heaving motion, Propeller DTMB P4718, J ¼ 0:4:

Fig. 32. Propeller (total) thrust coefficient as a function of blade angular position for the propeller at heaving motion, Propeller DTMB P4718.

has been included in the modeling allowing thus blade – wake interaction at higher propeller loadings to be taken correctly into consideration. The method is based on the use of four-node bilinear boundary element representation of the blade surface and the free wake. A Morino type formulation is used for the derivation of the corresponding integral equation with constant dipoles for each element. Free wake motion is simulated by tracking in time the grid nodes with the total fluid velocity for an inertial observer. Application of the method to a number of cases shows that very stable rollup wake patterns are produced for both steady and unsteady cases accompanied by good predictions of propeller blade pressures and forces.

Acknowledgements The work described in this paper has been carried out during a sabbatical leave of absence. I would like to thank the National Technical University of Athens, the School of Naval Architecture and Marine Engineering and the

department of Ship and Marine hydrodynamics for the opportunity they gave me.

Appendix A. Derivation of Eq. (14) A derivation of Eq. (14) is given by Hess, [3]. We shall present here a more elegant (and fast) proof based on cartesian tensor formulation of the representation theorem (1). The usual Einstein convention for repeated and free indices shall be used. Lemma 1. Let S be a Lyapunov surface, Gunter [19], embedded in R3 with regular boundary ›S: Let m be a sufficiently smooth function defined on S: Then with the terminology given previously in this paper, the following identity holds ð ð ~r ~r d~l £ ~r ð m 3 ¼ ð~n £ 7mÞ £ 3 dS þ ðmn~ 7Þ 3 dS ðA1Þ r r r ›S S S where, direction of line integration in the left hand side of Eq. (A1) has to conform with the direction of the unit normal vector n~ to the surface S:

652

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

Proof. We start from Stokes formula for ðS; ›SÞ in tensor form ð ð Fi dli ¼ e ijk dSj Fi;k ðA2Þ

representation theorem (1), apply Eq. (A1), take the limits as Pþ ! P0 [ W and P2 ! P0 [ W and finally take the mean value of the calculated limits.

where e ijk denotes the permutation symbol, Fi a sufficiently smooth vector function defined on S and Fi;k the covariant derivative of Fi : Select Fi according to the following relation

Appendix B. Derivation of Eq. (28)

›S

def

S

Fi ¼ e ikj m

rk Aj r3

ðA3Þ

where Aj denotes a parallel vector field in R : Then relation (A2) becomes   ð ð rk rm 1ikj m dli 3 Aj ¼ 1ijk dSj 1imn m 3 An ;k ðA4Þ r r ›S S 3

But 1ikj and An can be considered as constants in covariant differentiation, McConnell [30], thus     r r 1imn m m3 An ¼ 1imn m m3 An r r ;k ;k   rm r ¼ 1imn m;k 3 An þ 1imn m m3 An ðA5Þ r r ;k Using (A5) relation (A4) becomes ð ð r r 1ikn m dli k3 An ¼ 1ijk dSj 1imn m;k m3 An r r ›S S   ð r þ 1ijk dSj 1imn m m3 An ðA6Þ r ;k S and eliminating the parallel vector field from both sides, we get   ð ð r ð r r m1ikn dli k3 ¼ 1ijk dSj 1imn m;k m3 þ 1ijk dSj 1imn m m3 r r r ;k ›S S S ðA7Þ

Consider a straight line vortex of unit strength, starting at ~ Then the induced velocity point a~ and ending at point b: from the line vortex, ~v; at a control point p~; is given by 21 ðd ~l £ ~r ds; ðB1Þ 4p 0 r 3 where ~l is given by Eq. (30), s is the physical length along the line vortex (0 # s # d; where d is defined as in Eq. (30)), and ~r ¼ ~rðsÞ ¼ q~ðsÞ 2 p~; where q~ðsÞ is a position vector that scans the line vortex. Now q~ ¼ a~ þ s~l; so ~r is given by ! ~r ¼ a~ þ s~l 2 p~ ¼ a1 þ s~l ðB2Þ ~v ¼

and the numerator of the integral in Eq. (B1), becomes ! ~l £ ~r ¼ ~l £ ð ! a1 þ s~lÞ ¼ ~l £ a1 ; ðB3Þ which is independent of s: So the integral to be calculated is ðd 1 ds ¼ s2 ðB4Þ 3 0 r Using the definitions for s0 and s1 given by Eq. (31), we have r 2 ¼ s2 þ 2s s1 þ s0

ðB5Þ

and so the integral (B4) is given by the analytical relation (33). Using relations (33) and (B3), Eq. (B1) takes the form of Eq. (28).

But

m1ikn

 rk d~l £ ~r  r dli 3 ¼ m 3  ; 1ijk dSj 1imn m;k m3  r r r n   ~r ¼ ð~n £ 7mÞ £ 3 dS r n

Appendix C. Computer programming hints

ðA8Þ

and

       djm djn    rm dSj m rm 1ijk dSj 1imn m 3 ¼  r ;k  d d  r 3 ;k km kn    ~r  r ¼ m m3 dSm ¼ ðmn~ 7Þ 3 dS r ;n r n

ðA9Þ

where djm denotes the Kronecker delta and the following identity has been used: ðrm =r 3 Þ;m ¼ 0: Substituting relations (A8) and (A9) to relation (A7) we get identity (A1) A Using Lemma 1 the proof of Eq. (14) is a simple task. We only have to take the gradient in both sides of potential

A necessary final step in the application of a BEM is the computer implementation of the method. Our code has been written in Fortran 95 using the Lahey LF95 V5.5 Fortran compiler running on a PC (Pentium 4 –1800 MHz). The object oriented features of Fortran 95 have been fully implemented in our code, Chapman [26]. More specifically the heart of the code is a module which extends the Fortran language in a way more suitable for handling BEM calculations. By inspecting BEM calculations we observe that handling of vectors and vector operations is a very common task. Thus a ‘type vector’ object is defined and the operations between vectors are defined (i.e. addition, subtraction, multiplication by a number, inner and outer product, equality) by ‘overloading’ the corresponding fortran operators for reals. By using the Fortran 95 attributes ‘PURE’ and ‘ELEMENTAL’ in the corresponding routines, all operations can be performed

G.K. Politis / Engineering Analysis with Boundary Elements 28 (2004) 633–653

automatically on whole arrays or array sections. This results in enormous simplification in programming by avoiding doing loops and reducing, as a result, the size of the code. In the special case of an unsteadily moving propeller, ‘ELEMENTAL’ routines can be defined performing the rigid body translation and rotation of the propeller, task which has to be performed at each time step. Thus the hole propeller can be moved with a single command, e.g.: rotate(ths,cb), where ths denotes the vectorial instantaneous rotational angle and cb denotes the array containing the nodal points of the propeller. A quadrilateral bilinear element is an object character~ ~c; d~ and the dipole di and/or source s ized by its nodes a~; b; strength which it caries. We have found useful to introduce a ‘type boundary_element’ object with the following structure Type boundary_element Type(vector),pointer::a,b,c,d Real::di,s End type boundary_element The attribute pointer to the nodes of the boundary element has to do with the observation that the nodes of a boundary element are also nodes of the grid describing the body. That is there is an injective mapping of element nodes to the grid nodes which can be introduced in the code by making the element nodes to ‘point’ to the grid nodes. This results in better utilization of computer memory since only grid nodes for propeller and free wake are stored. On the other hand a boundary element appears continuously in our calculations as an individual object and thus it is easier to be handled as such in computer programming. The same observation holds for the sides of each element which appear as individual objects in the Biot – Savart formulas. Introduction of another object of pointer type to describe nodes and vorticity of element sides simplifies the computer implementation of the method. Note also that with the use of pointers any deformation of the nodes describing our free boundary or of the nodes of the moving propeller is automatically transferred to the element geometry. Closing this paragraph note that dipole influence coefficients are independent from the motion of the propeller. Note also that they are symmetric for all propeller blades. Thus one can calculate the dipole influence coefficients for the ‘key’ blade and then shift them cyclically to obtain the corresponding dipole coefficient matrix for the propeller. This can be done before starting the time stepping algorithm. References [1] Hess JL, Valarezo WO. Calculation of steady flow about propellers by means of a surface panel method. 23rd Aerospace Sciences Meeting, AIAA, Reno, Nev.; January 1985. [2] Hess JL, Smith AMO. Calculation of Nonlifting potential flow about arbitrary 3-D bodies. Journal of Ship Research 1964;8(2):xx.

653

[3] Hess JL. Calculation of potential flow about arbitrary three dimensional lifting bodies. Report no. MDC J5679-01, McDonnell Douglas Corporation; 1972. [4] Glauert H. Elements of aerfoil and airscrew theory. Cambridge: Cambridge University Press; 1926. first pub. [5] Kochin NE, Kibel IA, Roze NV. Theoretical Hydrodynamics. New York: Interscience Publishers; 1964. [6] Cottet GH, Koumoutsakos PD. Vortex methods—Theory and practice. Cambridge: Cambridge University Press; 2000. [7] Katz J, Plotkin A. Low speed aerodynamics. New York: McGrawHill; 1991. [8] Kerwin EK, Kinnas SA, Lee JT, Shih WZ. A surface panel method for the hydrodynamic analysis of ducted propellers. Transactions SNAME 1987;. [9] Morino L, Kuo CC. Subsonic potential aerodynamics for complex configurations: a general theory. AIAA Journal 1974;12(2). [10] Belibasakis KA, Politis GK. A boundary integral equation formulation of the Neuman problem for a vector field in R3 with application to potential lifting flows. Engineering analysis with Boundary elements 1995;16:5–17. [11] Belibasakis KA, Politis GK. A nonlinear velocity based boundary element method for the analysis of marine propellers in unsteady flow. International Shipbuilding Progress 1998;45. [12] Politis GK, et al. Selection of blade section characteristics in computer aided propeller design, SNAME. Marine computers 86 Symposium. [13] Kinnas SA, Hsin CY, Keenan DP. A potential based panel method for the unsteady flow around open and ducted propellers. 18th Symposium on Naval Hydrodynamics, Ann Arbor, Michigan, USA; 1990. [14] Maitre TA, Rowe AR. Modeling of flow around a marine propeller using a potential-based method. Journal of Ship Research 1991;35(2). [15] Takinaci AC. A wake rollup model for heavily loaded marine propellers. International Shipbuilding Progress 1996;43. [16] Takinaci AC, Atlar M. On the importance of boundary layer calculations instead of viscous correction in heavily loaded marine propellers while using a surface panel method. Ocean Engineering 2001;28. [17] Pyo S, Kinnas SA. Propeller wake sheet roll-up modeling in three dimensions. Journal of Ship Research 1997;41(2). [18] Kinnas SA, Pyo S. Cavitating propeller analysis including effects of wake alingment. Journal of Ship Research 1999;43(1). [19] Gunter NM. Potential theory and its applications to basic problems of mathematical physics. New York: Frederick Ungar Pub. Comp; 1967. [20] Kress R. Linear integral equations. Applied mathematical series 82, Berlin: Springer; 1989. [21] Mikhlin SG. Multidimensional singular integrals and integral equations. New York: Pergamon Press; 1965. [22] Sarpkaya T. Computational methods with vortices—the freeman scholar lecture. Journal of Fluids Engineering 1988;111(5):1989. [23] Blevins RD. Applied fluid dynamics handbook. Princeton, NJ: Van Nostrand Reinhold Company; 1984. [24] Klein JL. A rational approach to propeller geometry. Propellers 75, SNAME, TandR Symposium; 1975. [25] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran. Cambridge: Cambridge University Press; 1986. [26] Chapman SJ. Fortran 90/95 for scientists and engineers. New York: McGraw Hill; 1998. [27] Jessup SD. An experimental investigation of viscous aspects of propeller blade flow. PhD Dissertation. The catholic University of America, Washington, DC; 1989. [28] Gindroz B, Hoshino T, Pylkkanen JV, editors. 22nd ITTC propulsion committee, propeller RANS/Panel method workshop, Proceedings, Grenoble, France; 1998. [29] Jessup SD. Measurement of the pressure distribution on two model propellers. David Taylor Naval Ship Research and Development Center, Report no DTNSRDC-82/035; 1982. [30] McConnell AJ. Applications of Tensor Analysis. New York: Dover Publications; 1957.