Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

A Three Dimensional Theory for the Design Problem of Propeller Ducts In a Shear Flow I. Falcao de Campos (Maritime Research Institute Netherlands, The Netheriands) ABSTRACT A linearized theory of the three-dimen- sional steady interaction between a ducted propeller system and a radially and circum- ferentially sheared axial onset flow is pre- sented. Following duct lifting surface theory the duct is modelled by a distribution of pressure dipoles and sources on a reference cylinder to represent the effects of loading and thickness. An actuator disk model is used to represent the effects of propeller loading. An integral equation for the pressure distur- bance is derived which may be applied to treat both the effects of loading and thickness. The potential and shear interaction components of the disturbance pressure are treated separate- ly and a computational scheme is applied to solve the integral equation for the interac- tion pressure. The results of sample calcu- lations for the effects of duct loading in axisymmetric and non-axisymmetric wakes are presented and discussed. NOMENCLATURE Amn ~ Ban (cay ~ ~ ) ( ex ~ er ~ e(3) Em F f G _mn = Gmnij,Gmnij Gll, G12, Gal, G22 H(x) hl,h2,h3,h4 I~,Km k L(~) P(lnl ) P,P ~ Ap, Ap Qm+l/2 q R Functions of shear param- eters Shear parameters Duct semichord Unit vectors, Cartesian coordinates Unit vectors, cylindrical coordinates Fourier integral in the downwash calculation External force field Duct camber Kernel function Matrix elements Functions in degenerate kernel function Heaviside unit step function Functions of shear param- eters for duct loading and thickness r S T ),U U U~,U2,U3,U4 u (u,v,w) Modified Bessel functions of order m And, index of radial node i1,i2,ip,kp _ Functions of shear param eters for propeller loading index of radial node Parameter in x-wise Fourier transform Function in downwash calcu lation Downwash kernel function Pressure, respectively its x-wise Fourier transform Strength of pressure dipole, respectively its x-wise Fourier transform Legendre function of second kind and half order Strength of source distribu tion Radius, distance between two joints Transformed radius Right-hand side of type equation Function of source distribu tion Duct thickness Undisturbed axial velocity, respectively its modulus Reference velocity Parameters in analytical defined wake Disturbance velocity Axial, radial and circumfer ential components of distur bance velocity Fluid velocity V(O) TV (x,y, z) (x,r,8),(L,~,4) J.A.C. Falcao de Campos, MARIN, P.O. Box 28, 6700 M Wageningen, The Netherlands 645 Influence function for ra- dial downwash velocity on the duct Average radial velocity in- duced on the duct Radial velocity jump on the duct Cartesian coordinates Cylindrical coordinates

&(x) n A p T 1, V Subscripts d m,n t - Conical angle - Expansion rate - Dirac delta function - Small parameter - Normalised axial distance on the duct - Parameter of x-wise Fourier transform - Fluid density - Integral operator - Integration volumes - Right-hand side of integral equation - Argument of Legendre func- tions - Refers to duct - Index of circumferential harmonics - Refers to propeller - transverse component Superscripts (O) (1) 1. INTRODUCTION - Refers to potential, induc- tion part - Refers to shear interaction part Ducted propellers are a well-established means of ship propulsion. It is well-known, [1] that the use of a ducted propeller with an accelerating type of duct improves the effi- ciency of the propulsor in case of heavy load- ing. Also the use of a decelerating type of duct may be beneficial to reduce the risk of cavitation of the propeller. A large number of conventional duct designs, which have been most successfully applied in practice [2]-[1], are axisymmetric but the application of non- axisymmetric ducts has also been subject to investigation both experimentally [1] and theoretically [3]-[4]. These ducts have been applied to reduce the non-uniformity of the inflow to the propeller in the the ship's wake leading to improved performance from the point of view of efficiency, cavitation and vibra- tions. For the design of ducted propellers a number of analytical tools have become avail- able along the years. Early ducted propeller theories [5]-[6] were based on linearised annular airfoil theory for the singularity representations of the effects of duct loading and thickness, in combination either with an infinite blade number model (actuator disk) or with a finite bladed lifting line model of the propeller. An extensive review of these the- ories was made by Weissinger and Maass in ref. [7]. It is interesting to notice that the truly inverse methods published to date for designing axisymmetric propeller ducts are based on these theories. The methods determine the duct geometry (in the presence of a time averaged propeller induced velocity field) for specified duct pressure distribution [8] or given load and thickness distributions [9]. These methods suffer from the drawback that it is not possible to guarantee a priori that the given pressure or load distribution will lead to an acceptable duct geometric shape. Never- theless, those inverse methods were of great assistance in designing famous ducted pro- peller systematic series, such as the ones ] and [10]. Often to be modified to - published in references [1 the final duct shapes needed meet practical requirements. Following the developments in the numeri- cal methods for the calculation of potential flow on lifting bodies, methods for the hydro- dynamic analysis of ducted propellers evolved to a greater degree of sophistication. More accurate panel representations of the duct geometry have been employed for axisymmetric flow [11] and, more recently, for three-dimen- sional flow [12]-[13]. These last methods for steady three-dimensional analysis have concen- trated on the complex interaction between pro- peller and duct in uniform inflow by incorpo- rating lifting surface or panel representa- tions of the propeller blades. Also complete unsteady potential flow analysis [14] of the ducted propeller system has been attempted. The methods mentioned previously are restricted to potential flows. In reality the ducted propeller operates in the highly non- uniform flow endowed with vorticity in the ship's wake and the interaction with this flow is an important field of research in propulsor design. In dealing with this problem the po- tential flow methods have retained completely their usefulness through the introduction of the concept of the effective onset velocity, which is defined as the total velocity minus the potential velocity induced by the propul- sor. The effective onset velocity has to be computed by some (viscous or inviscid) rota- tional model for the propulsor-hull interac- tion. Examples of the inviscid approach to the computation of the effective velocity for con- ventional propellers in axisymmetric flow can be found in [15]-[17]. There is a considerable amount of ref- erences in the turbomachinery literature deal- ing with the problem of solving approximate forms of the Euler equations for determining the inviscid disturbance flow to parallel shear flows, as can be found in the survey given by Hawthorne [18]. In particular, the large shear - small disturbance approximation, applied along the lines set in the classical works of Karman and Tsien [19] for a lifting line and of Lighthill [20] for a simple source, have been used to fundamentally study the effects of shear in the flow around aero- dynamic shapes [21]-[22], including the an- nular airfoil [23]-[24]. For an infinitely bladed propeller, modelled by an actuator disk, the same approach has been followed to investigate the effects of shear in the incom- ing flow in the axisymmetric case [25]-[26], plane flow [27] and three-dimensional flow 646

[28]. The latter reference includes the ef- fects of shear of a radially and circumfer- entially varying axial inflow to the actuator disk and incorporates some non-linear terms in the equations of motion adequate for extending the method for heavier loadings. For ducted propellers the effects of shear have been investigated in [11] using a numerical vortex panel method in axisymmetric flow. Recently, Lee [29]-[30] presented a lin- earised analysis of the ducted propeller sys- tem in axisymmetric shear flow suitable for solving the duct design problem in the pres- ence of the propeller. He showed that the shear significantly affects not only the duct camber and ideal angle of attack but also the duct induced velocity at the propeller. The objective of the present paper is to present a three-dimensional theoretical anal- ysis of the steady interaction between a ducted propeller and a radially and circum- ferentially sheared axial inflow. Consistent with the steady flow assumption, an infinitely bladed propeller (actuator disk) is consi- dered. Although not strictly necessary to the analysis, an axisymmetric loading over the propeller disk will be assumed in this paper for simplicity. The duct loading and thickness may vary in the circumferential direction. The theory may be used in the analysis of a ducted propeller with a given duct shape but it will be most readily applied (non-iteratively) in the inverse mode i.e. to determine the duct section camber and angle of attack for given thickness and loading distributions. As in the method followed by Lee [29] in the axisymmetric case, the analysis developed herein approaches the solution of the problem from the theory of the Poisson equation to derive a system of coupled integral equations for the disturbance pressure harmonics. How- ever, some essential differences with the formulation used by Lee are noteworthy: First, a single type of integral equation in the dis- turbance pressure is used throughout, encom- passing both the effects of loading and thick- ness. Second, by separating shear interaction effects from induction effects, the method recovers the potential formulation of the duct lifting surface theory, see for instance [6], [31]-[33], and actuator disk theory [34], while, at the same time, the presence of any singularities of the pressure field are re- moved from the problem for the interaction with shear. The paper is organized as follows: In section 2 the theoretical analysis is pre- sented. Section 3 deals with the numerical procedures employed so far to solve the integral equation and compute the velocity field. In section 4 the results of sample calculations illustrating the effects of shear in the velocity field due to a non-symmetric duct in a wake field are presented and dis- cussed. The paper closes with some remarks regarding the basic limitations of the method and its further development. 2. THEORETICAL ANALYSIS 2.1. Equations of Motion for General Steady Disturbances to a Shear Flow We start by deriving the linearised Euler equations for the steady flow of an inviscid and incompressible fluid in the presence of an external force field. Anticipating the use of sources and sinks to represent thickness ef- fects, we will assume the rate of expansion to be zero except at the points where such singu- larities will be present. The continuity equation for an incompres- sible fluid reads, see for instance ref. [35]: V.~= A, (1) where ~ is the fluid velocity and ~ is the lo- cal rate of expansion. The Ruler equations for the steady flow of an incompressible ideal fluid are (it V)4 + 9(~) = pp , (2) where p is the pressure, p the fluid density and F the external force field per unit vol- ume. We introduce a Cartesian coordinate sys tem (x,y,z), with unit vectors (e ,e ,e ) and a cylindrical coordinate system with Yunii vec tors (e ,er,e~), Fig. 1. We consider a radi ally And circumferentially non-uniform axial flow D(r,8), independent of the axial coordi nate, to be disturbed baby the presence of the ducted propeller. Let u denote the disturbance velocity so that the fluid velocity may be written as ](x,r,0) = U(r,8)ex + u(x,r,6) . (3) ~\~ Y / \ Fig. 1. Coordinate system and axial inflow. 647 x

The equations of motion (2) may be lin- earised by assuming the disturbance velocity to be small, say of 0(~), in comparison with the undisturbed velocity 0=0(1). Substituting eq. (3)2 into (2) and neglecting squared terms of 0(s ) in the disturbance velocities we ob- tain u aU (A v)(u~ ) v(~) F (4) since (Uex.V)(Uex) = U a(Uex)/ax = 0 We further decompose the disturbance ve- locity into its axial and transverse compo- nents u ueX + Ut ~ ith ( ) ~ ~ With (5) we fin that (u. 9) (Uex) = U ax (Uex) + (Ut. v) (Uex) = = eX(ut. MU) and (4) may be written in the form u aU + e~x(u~t. W) + V(~) = Fp (6) Taking the divergence of (6) and using (1) we obtain 92(~) + 2 ax (U-t.vu) = 9.(Fp) - U ax · (7) We may Eliminate the transverse velocity com- ponent ut between (6) and (7) to obtain a sin- gle equation for the pressure. To do so we first take the transverse vector component of (6) U ax + ~t(~) = p ' where Vt - -ex x (e x V) is the gradient op- erator In the transverse plane, and substitute (8) into (7) to obtain V2(~) + _ :9tU.Vt(P)] = ( pt. VtU) - U ax Equation (9) is a linear partial differ- ential equation for the pressure disturbance which needs to be solved for given U. F and distributions. 2.2. Ducted Propeller Model Let us consider now the specific form of the right-hand side of eq. (9) for the external force fields F representing the duct and pro- peller loadings and for the rate of expansion field ~ representing the duct thickness in the linearized theory. We note that, since ]=0(1), from eq. (9) we will require F=O(~) and A=0(~) Denoting by S(x,r,O) the right-hand side of eq. (9) (multiplied by p), we obtain in cylin- drical coordinates: aFX 1 a~rFr) 1 aF S(x,r,6) = aX + r ar + r ae 5) - U (Fr ar + rat ax) - Pu ax . (10) Consistent with the linearization applied to the Euler equations, we will place the force field singularities and the rate of ex- pansion singularities (sources) producing the disturbance velocities, on a reference surface aligned with the undisturbed flow d. In the present ducted propeller application we choose the reference surface to be a cylinder of con- stant radius and chord. We represent the duct loading by a radi- ally directed force field F=(O,F ,0) distri- buted on a cylinder of radius R and chord 2c, extending from x=-c to x=+c, Fig. 2, in the form: Fr = Apd(x,8)[H(x+c)-H(x-c)] &(r-Rd) , (11) where H(x) is the Heaviside unit step func- tion, H(x)=0 if x<O, H(x)=1 if x>O, and 8(x) is the Dirac delta function. Substituting eq. (11) into eq. (10) we obtain S(x,r,8) = _ Apd(x,0)[H(x+c)-H(X-C)] Fir- - U ar) &(r-Rd) + &'(r-Rd)] r l -C 648 . (12) +c Fig. 2. Duct geometry conventions x

It is seen from equations (9) and (12) that the duct loading is represented by a distribu- tion of pressure dipoles of strength Apd(X,0) on the reference cylinder. The duct thickness is represented by a source distribution q(x,6) on the reference cylinder. The rate of expansion becomes A(x,r,0) = = q(x,0)[H(x+c)-H(x-c)] &(r-Rd) (13) Substitution of eq. (13) into (10) yields S(x,r,8) = - pU aX (q(x,8)[H(x+c) - H(x-c)]) 6(r-Rd) . (14) The propeller loading is assumed axisym- metric over the propeller disk. Neglecting the radial and circumferential components, the force field is assumed to be of the form {=(Fx,O,O), with Fx(x,r) = App(r) 6(x-xp) , r<R , (15) where Ap (r) is the radial distribution of propellers loading, x the axial coordinate of the propeller plane find R the propeller ra- dius. In this case the function S takes the form S(x,r) = App(r) 8'(x-xp) . (16) The expression (16) provides the usual repre- sentation of the force field associated with the actuator disk as a distribution of pres- sure dipoles on the disk. 2.3. Formulation of the Integral Equation 2.3.1. Derivation of Kernel Function By application of Green's third identity we may transform the Poisson-type eq. (9) into an integral equation for the pressure distur- bance. We first note that the function S is only different from zero in a bounded region of space. Assuming the pressure disturbance to vanish at infinity to sufficient order, eq. (9) can be recast in the form (using cylindri- cal coordinates): p(x,r,6) + 2~ ~ U ma ap + (a a¢)(a aP)]R do = = _ 1 ~ S((,a,¢) do (17) where R is the distance of the field point (x,r,0) to the integration point (L,a,4) and do = ad~dad¢. In eq. (17) the integrations are to be carried out over the whole region of space v extending to infinity. The integral eq. (17) can be reduced by the use of Fourier transforms to a system of coupled one-dimensional integral equations for the pressure disturbance harmonics. To perform such reduction we represent the pressure field by a double Fourier expansion in the axial and circumferential coordinates as follows: p(x,r,6) = no a' = 2~ £ Fine ~ pn(r;~)e dX , (19) n=-- _= with inverses) pn(r;A) = 1 ~ ~ p(x r ~)e-i(\X+nO) dude . (20) The "source" potential 1/R given by eq. (18) admits the expansion 1 _ 1 ~ eim( Ott) R ~ m=- a) ~ Ik(m)(r'a;k)eik(X-~) dk (21) _co with _(Iklr) K~(lkla) , r<a r) lm~lKl~) , r>~ ~ (22) Ik(m)(r.a;k) = ~ m ; ; m;; ; where I and K are the modified Bessel func . m m lions of order m. We introduce the notation a(r,8) = U ar ~ b(r,8) = U r at (23) 1. The following convention was adopted in de- fining Fourier transforms: so f(~) = ~ f(x)e ibex dx _a) with inverse f(x) = 2~ ~ f(~)ei\X dA _co for a continuous F.T.; f(~) = Fin with inverse f = 1- ~ fee ins dO n z~ O for a discrete Fourier series. 649

and we expand a(r,8) and b(r,8) in Fourier with series a(r,8) = ~ anteing n=-= a, b(r,8) = £ bn(r)ein~ , n=-- with an(r) = 2~ ~ a(r,0)e ins dO , bnfr) = 2~ ~ b(r,8)e ins dO . Substituting eq. (19) and eqs. (21)-(23) into the first integral of (17) and carrying out the integration in a we obtain (a ap + by an) R do = a' co co = ~ £ e ~ dk e ~ da a Ik( )(r,a;k) _= O co co x ~ dX £ pa (a) P'(a;A) + L m-n- - -n _00 n__m ~ n=-= 00 lion bm-n(a) pn(a)] ~ dL ei(\ k)t (26) _co where Pn = dpn/da. By noting that co 8(~-k) = 2 ~ ei(A-k)( dL (27) _- the integration in ~ can be carried out to give the integral, the value ao co Go a' £ imO ~ dk eikx 2 £ ~ do Ik(m)(r,a;k) m=-= -= n=-= O a am_n(a) pn(a;k) + in bm_n(a) pn(a;k) . (28) The integral involving p' can be further re- duced by partial integration. Evaluating the derivatives of the terms involving the Bessel functions we finally obtain ~ (a at + a at) R dt ~ co = - £ Him ~ do eiAx m=-= O x £ ~ Gmn(r~a;A) pn(a;A) da ~(29) n=-0 0 650 (24) Gmn(r,a;A) = and A (a) = (1-1ml) a (a) + 2 Km(lAlr) [Amn( a) Im( I Al a) + Bmn(a)l\I Im_1(1\Ia)] , r>a 2 Im(lAlr) [Amn(a) Km(lAla) Bmn(~)l\I Km_1(1\Ia)] , r<a (30) (25) a a' (a) - in b (a) , (31) B (a) = a a (a) an men In eq. (31) the prime denotes differentiation with respect to the function argument. We note that the kernel function Gmn(r,a;A) has a dis- continuity at r=a given by - 2Bm_n(r)~A~ tIm(~A~r) Km_1(~A~r) + Im-l(~A~r) Km(~A~r)] = 2am-n( ) (32) 2.3.2. Duct Loading The right-hand side of eq. (17) for the effect of duct loading may be reduced in a similar way. We assume the duct loading to be represented by the double Fourier expansion co a, APd(X,8) = 2~ £ e ~ APd (A)e dA , (33) n=-= -= n with inverse APd (A) = 2~ ~ ~ APd(x'8)e-i(n8+\x) dxdO n O -c (34) Substituting eqs. (33), (12) and (21) in (17) we obtain, after carrying out the integrations in a, ~ and A, R do a, co = - 2 £ elms ~ dA eiAx x [Rd IK( )(r,Rd;A) APd (A) + m a) 2Rd IK( )(r,Rd;A) £ am_n(Rd) APd (A)] n=-= n (35)

In eq. (35) we have used IK(a )(r,a;~) to de- note m) X Im( X r) Km( X a) , r<a IK(a (r,a;~) = {X K (jar) I,(~X~a) , r>a (36) 2.3.3. Duct Thickness In the case of duct thickness, inserting (14) and (21) into the right-hand side of eq. (17), we obtain S(L,a,¢) do = v _ P ~ eimb J do J da a ~ m=-- -= O J dk IK(m)(r,a;k)eik(X~~) ~(a-R~ _= a:{[H(~+c)-H(~-c)] J U(a,¢) q(L,¢)e ¢> ~ . (37) Introducing the function Tm(L,a) = 2~ J U(a,¢) q(t,¢)e i ~ do , (38) and carrying out the integration in a we have for the integral (37) -2p ~ e Rd J dX IK( )(r,Rd;~) Go x ~ dL ei\(X a) dad {Tm(L,Rd)[H(~+c)-H(~-c)]} . -a (39) In the previous expression the integration in ~ may be performed by introducing the Fourier expansion Tm(A) = J Tm(L,Rd)e ~ dL (40) -c to yield S(L,a,¢) do = -2 ~ Tim ~ dX eii v m=-- -= x ipX Rd IK( )(r,Rd;~) Tm(~) (41) 2.3.4. Propeller Loading With the propeller loading given by eq. (15), the right-hand side of eq. (17) is re- duced in the same way as the previous cases. The integral is ~ R v = - J dL JP da a J do £ eim(~~~) -= O O m=- I dk IK(m)(r,a;k)eik(X~~) Ap (a) (-x ) _= (42) Since the load distribution is assumed inde- pendent of +, the integration in ~ gives In . J e lmt df = an 8mO , (43) where ~ is the Kronecker delta, and, hence, the integral becomes R v ~ aim J dk JP da IK( )(r,a;k) m=-= -= O alp (a) J eik(x-~) (-x ) d: = 2 £ aim J do Minx m=-= _ -i~x R ,^` x iX e P JP IK`U'(r,a;~) App(a) a da. (44) 2.3.5. Integral equation Gathering the previous results expressed by eqs. (29), (35), (41) and (44), by equating the corresponding harmonics of the left and right-hand sides of eq. (17) we may write (symbolically) the integral equation in the form ~ ~ ~ Pm T Pn = Em ' with the operator T defined by (45) ao co Pn ~ J Gmn(r~a;~) pn(a;~) da (46) n=-0 0 The function ~ takes different forms for the effects of dumct loading, thickness and propeller loading. In the case of duct loading it is given by 651

~m(r;A) = = Rd IK( )(r,Rd;~) APd (a) + m a, 2Rd IK )(r,Rd;~) ~ am_n(Rd) APd (I) n=-= n (47) It is easily seen from the properties of the modified Bessel functions that ~m(r;~) is discontinuous at the duct reference cyl- inder r=R . The discontinuity amounts to the (Fourier transformed) pressure jump on the duct, i.e. the strength of the pressure dipole Id ( A) m In the case of duct thickness, eq. (41), the function takes the form (r;~) = ipA Rd IK( )(r,Rd;~) Tm(~) · (48) It is also seen that this function is continuous at r=R but has discontinuous first derivatives at that radius. Finally, in the case of propeller load- ing, the function Am is (r;A) = -ink Ret Ink = -id e ~ [~ IK`V'(r,~;~) App(~) ~ do, (49) which is a continuous function of the radius with continuous first derivatives. The eq. (45) constitutes an infinite sys- tem of one-dimensional integral equations, with a discontinuous kernel, for the distur- bance pressure harmonics. At this stage it should be remarked that, as shown in the classical work in shear flow problems [19], [20], the advantage of having introduced the x-wise Fourier Transform is to reduce the three-dimensional problem to a set of de- coupled two-dimensional problems, each one for a single value of the parameter ~ in eq. (45). However, it should be noted that, in the case of a circumferentially sheared inflow, the circumferential pressure harmonics do not decouple. In fact, they are coupled through the harmonics of the shear parameters, as shown by the form of the integral operator in eq. (46). We would like also to stress that the integral operator T is a function of the shear parameters a and b only, being independent of duct loading, thickness or propeller loading, which affect only the right-hand side of eq. (45). 2.4. Separation of Potential and Interaction Effects The numerical approach to the solution of the integral equation (45) has to be done with care. First, the kernel G is discontinuous at ran, as given in (32). This may not consti- tute a major problem since a number of numer- ical techniques are available to handle this type of kernel, as it may be found, for in- stance, in ref. [36]. Second, the right-hand side of eq. (45) is either discontinuous at r=Rd, as in the case of duct loading, or has the discontinuity in its first derivative, as in the case of duct thickness. The integral equation (45) is of the second kind and this means that the solution will inherit the dis- continuous behaviour at r=Rd. This fact has major consequences for the computation of the velocity field because, as shown in the next section, differentiation of the pressure in the radial direction is required to derive the radial velocity component from the radial mo- mentum equation. In the pursuit of an accurate numerical solution of the integral equation (45) and of an accurate computation of the associated ve- locity field, it is quite natural to distin- guish two contributions to the pressure field, respectively the potential pressure p(O), which satisfies the Poisson equation (9) in the absence of shear, ~tU=O, and the interac tion pressure p(1) which satisfies the remain- ing part of the equation (9). For the trans- formed equation (45) in Fourier space this means that p = p(O) + p(1) with p(O) = ¢(0) m m and (50) (51) Pm ) - T pml) = ¢(1) + T ¢(0) (52) In eqs. (51) and (52) we have decomposed the right-hand side into a shear-independent term 9(0) and a shear-dependent term 9(1). In the case of duct loading we have m Em ) = Rd IK( )(r'Rd;~) Ad (I) (53) m and a) .(1) = 2Rd IK(m)(r,Rd;A) ~ am_n( d) Pd n=-= n (54) In the cases of duct thickness and pro- peller loading the shear-dependent term ~ml) is inexistent and Amy) is given by eqs. (48) and (49), respectively. 652

The potential flow problem eq. (51) for the duct loading and thickness leads to the formulation of the duct lifting surface the- ory. Solution of the loading (vortex distri- bution) and thickness (source distribution) problems has been given in refs. [6], [31], [32], where expressions can be found for the evaluation of the radial velocity (downwash) at the duct reference cylinder. Expressions for the duct induced axial and tangential ve- locities are also given in ref. [6]. We will discuss briefly the solution of this problem in the section devoted to the computation of the velocity field. The solution of the shear interaction problem eq. (52) requires the computation in the right-hand side of eq. (52) of the term T M(°). A straightforward computation allows this term to be expressed in terms of func- tions which can be easily evaluated numeri- cally. The results of these computations for all cases are given in the Appendix A. It may be readily verified using the expressions of the Appendix A that the right-hand side of eq. (52) in the cases of duct loading and thickness is continuous and has continuous first derivatives at r=Rd, a property which is shared by the solution pml) of the shear in- teraction problem. 2.5. Velocity Field In accordance with the decomposition of the pressure field into its potential and in- teraction parts, we write for the disturbance velocity ~ -(0) -(1) u = u + u (55) where u(O) is the velocity associated with the potential pressure p(O), which we call the in duction velocity, and u(1) is the shear inter- action velocity -(0) By definition u satisfies eq. (1) and the momentum equation U aaX + V(P(O)) = Fp, (56) while u(1) is a solenoidal velocity field V U-(1) o satisfying (57) U aaX + e-XP(ut°) + utl)).~7U] + V(~)) = 0 (58) If the pressure fields p(O) and p(1) are known, respectively from eq. (51) and the so- lution of eq. (52), the velocity fields u(O) and u(1) can be evaluated by integrating eqs. (56) and (58). We note that, in using (58) for the determination of the interaction velocity, the potential part of the solution affects the radial and tangential components of the inter- action velocity only through the interaction pressure p(1). For the axial velocity com ponent, however, there is an additional term u(O). VU coupling the two problems and this requires the calculation of the induction ve locity utO) first. Integrating eq. (56) from x=-= to x, we obtain for the induction velocity components u(O) = _ p(O) 1 ~ x dE, (59) v( ) = - U ~ ar ( p ) d: + U ~ p dE, (60) W(°) 1 ~X 1 a (P(°)) d: (61) In writing eq. (59) to (61) we have made u(O) and the use of the fact that the velocity disturbance pressure p(O) vanish at infinity upstream, x=-~. In the cases of duct loading and thick- ness, F =0 and the evaluation of the axial and tangential velocity components u(O) and u(1) do not pose particular problems. They can be computed from the inverse Fourier transform of ~m°), attention being paid to the calculation on both sides of the duct reference cylinder r=R _0. However, in the evaluation of the ra- dia] velocity component from eq. (60) special care has to be taken in carrying out the re- quired differentiation and integration at r=Rd. Such evaluation constitutes the subject of duct lifting surface theory dealing with the evaluation of the Cauchy singularity in the kernel function for the radial downwash. The solution of this problem can be found in refs. [6], [31]-[33] and will not be treated in detail here. In Appendix B the relevant expressions are collected for the case of duct loading, together with an outline of their derivation from the present formulation. In the case of propeller loading F =0 and with F given by eq. (15), we obtain [or the axial coXmponent u(O) = - Pu(O) + ~U H(x-xp), r<Rp (O) pU , r>R (62) the well-known result from actuator disk For the interaction velocity u(1), inte- gration of eq. (58) from x=-= to x yields 653

u(1) = _ p(1) _ a ~x(v(O)+v(l)) do - b ~ (w(O)+w(1)) do (63) _co ', ( 1 ) = ~ u ~ a (p(l)) dE, (64) _- w = - U ~ r ae ( p ) d: ' (65) _a) with the definition (23) of a and b. Again, in deriving (63) to (65) we have made use of the fact that the interaction pressure and veloc- ities vanish at infinity upstream. In eqs. (63)-(65), which hold for the cases of loading and thickness we have been considering, v(1) and w(1) have to be evaluated first from the interaction pressure; u(1) is then computed from eq. (63) with known values of v(O), v(1), w(O) and w(1) 2.6. Duct Boundary Conditions We describe the duct surface by speci- fying the deviations of the outer and inner surfaces from the reference cylinder due to the conical angle, camber and thickness dis- tributions of the duct sections, Fig. 2. The outer and inner surfaces may be given by the expression r(x,0) = Rd ~ a(~)(x-c) + f(x,0) + 2 t(x,8) , -c<x<c, (66) where a(~) is the conical angle, f(x,8) is the camber distribution, satisfying f(-c,8) = f(c,0) = 0, (67) and t(x,8) is the thickness distribution. In eq. (66) the plus and minus signs apply to the duct outer and inner surfaces, respectively. The conical angle a(~)is defined positive as shown in Fig. 2. The duct surface may be de- scribed by the equation Fd(x,r,0) = r - Rd + a(~)(x-c) - f(x,0) + 2 t(x,8) = 0 . (68) The unit vector normal to the surface may be expressed as ~Fd/~9F ~. The kinematic boundary condition on the ducdt surface reads (~.Wd)/~9Fd~ = 0, which can be written, using eq. (68), (69) (U+u) ta(~) _ af + _ at] + v +Wr ta'(~)(x-c) _ af + 2 a~] = 0 . (70) The first order approximation to (70) is V(X'Rd+°'6) = = u(Rd,6)~-a(~) + aX + 2 a-x] (71) Writing v(x,R '+O,8) = v(x,Rd,8) + 2 Av(x,Rd,0) , (72) where v(x,Rd,8) is the average radial velocity at r=R between the outer and the inner sur- faces and Av(x,Rd,8) is the radial velocity jump, we conclude from (71) and (72) that Av(x,Rd, O) = U(Rd, O) ax ' which gives to the first order the strength of the source distribution q(x,0) = U(Rd,8) at · (74) From eq. (71) we derive also to first order v(x,Rd,6) ax- a(~) = U(Rd,~) (75) Integrating (75) and using (67), we obtain =(~) = ~ 2 J. v dx (76) -c and f(x,0) = oc(x+c) + ~U dE . (77) -c It is well-known that the use of the first order boundary condition (71) can be rather crude for design purposes at increasing load- ings, since the perturbation velocities will no longer be small compared with the undis- turbed velocity U(Rd,8). A simple refinement to eq. (75) can be used if we neglect the contribution of the circumferential~ velocity w to the component of the velocity ~ normal to the surface in eq. (70). If, in addition, we retain the first order approximation to the source distribution (74), we obtain af v ax ~ ~= U+u (78) which includes the axial perturbation velocity in the calculation of the conical angle and camber distributions. 654

3. NUMERICAL SOLUTION PROCEDURE For the purposes of numerical analysis a dimensionless form of the equations is used. The duct radius Rd is taken as reference length, velocities are made dimensionless by a reference velocity DO, taken here as the (fi- nite) asymptotic value of the velocity U at large radii, and pressures are made dimension less by pU 2. To solve numerically the integral equa- tion (52) we first apply a transformation to the radial coordinate r=r(r) which maps the integration interval (0,-) onto the interval (0,1). Various possibilities exist to choose such transformation, but we have applied a simple exponential transformation defined by r = 1 - exp(-ar) , (79) with the constant a=ln(1/2) chosen as to map the duct radius r=1 to r=0.5. The integral equation is solved by a quadrature method us ing the trapezoidal rule. Introducing the equidistant nodes r; = (i-1)/(NR-1) , i=l,,NR , (80) the integral in eq. (46) evaluated at the node i is approximated by (omitting the parameter for simplicity) ~ mn(ri'~) Pn(~) do = = 2 Gill ~ wiJ G12j Pn J (2) 2 G21i ~ wij G22j Pnj ' (81) with ill Km( I HI ri ) ~ G21i = Im( I HI ri ) G12j = [A n(ri ) Im( 1 Al r; ) + Bmn(rj ) I \1 Im_l( I A| rj ) ] (dr/dr) G22j = [Amn(rj ) Km( | )~| rj ) Bmn(rj ) 1 x1 Km_l( l Ad rj ) ] (dr/dr)j (82) w(.J) and wi2) being the weights of trape zoidal rule. For the first N+1 harmonics of the disturbance pressure, the discretized form of the left-hand side of equation (52) reads with N NR = ( Pn )mi = ~ N jI1 Gmnij Pn; , (83) G .. = G .. for man or id mn~J man - G mii = 1 - Gmmii ' 655 mnlJ 2 Glli G12; WiJ , j<i , 2(Glli G12i wii G21i G22i Wii ) G21i G22j WiJ I, j=i, , j>i . (84) With the mapping (79) and the equidistant node distribution an equal number of nodes are placed inside and outside the duct. The map- ping ensures a larger concentration of nodes near the axis r=O, where the Bessel functions K are rapidly varying. To achieve an accept- a~le accuracy near r=0 and at large radii ex- ponential scaling of the Bessel functions is applied. The integrals appearing on the right- hand side of equation (52) (see Appendix A) are evaluated by trapezoidal integration using the same set of nodes (80). The system of equations is solved in sequence for the set of values of the param- eter ~ coinciding with the nodes of Laguerre integration in the interval (0,=) which is used for inverting the pressure x-wise Fourier Transform. A LU factorization of the matrix is carried out for each value of X, and then used to obtain the solution vectors for the differ- ent right-hand sides (47), (48) and (49). With the interaction pressure obtained from the inversion of the double Fourier Transform (20) at a specified number of axial planes, the interaction velocity field is com- puted with the eqs. (63)-(65). The integrals are computed by trapezoidal integration. Cen- tral differencing is applied to evaluate the radial and circumferential derivatives. 4. RESULTS OF SAMPLE CALCULATIONS To illustrate the effects of shear a num- ber of sample calculations were carried out for an analytically defined wake field. Here we will only discuss the results concerning the effects of duct loading on the velocity field. The wake field chosen is a sinusoidal perturbation superposed to the axisymmetric wake field used by Lee in ref. [30]. It is defined by the expression U(r,8) = 1 - U1 e 2 (1 - U. r2 sing U O) (85)

With U1=0.72, U2=1.70 and U3=0 we recover the axisymmetric wake of Lee. For the wake field (85) the functions Amn and Bmn approach zero at the axis sufficiently fast to ensure a zero value of the element G22j at the axis, j=1. In all the calculations a chord/diameter ratio c/R =0.5 has been chosen. In general the duct loading and thickness distributions are specified by a sum of chordwise modes which are used in the solution of the duct lifting surface problem. In the present application we have chosen the NACA a=0.8 chordwise loading, which is constant from the leading edge to 0.8 of the chord and decreases linearly to zero from 0.8 of the chord to the trailing edge. The Kutta condition of zero loading at the trailing edge is automatically satisfied by this load distribution. The duct loading is expressed in terms of the duct section lift coefficient defined by ZZUO ~Zc' _1 d (86) Since the theory is linear, the distur- bance pressure and velocities are proportional to the loading coefficient CL, and the results hold for an arbitrary loading. The results are presented for a considerably high duct loading and a strongly sheared inflow which may be considered as representatives of a typical ducted propeller application. In such cases the disturbance velocity is no longer small in comparison with the velocity of the incoming flow. Of course the assumption of small disturbance velocities does not hold in such cases and non-linear effects will certainly affect (in an unknown manner) the accuracy of the predictions. In any event, the results are easily scaled to smaller loadings. To examine the convergence of the numeri- cal solution of the integral equation a compu- tation for axisymmetric flow was carried out first for a section lift coefficient CL=0.9. Fig. 3 shows the convergence of the interac- tion pressure at the duct inlet plane x=-0.5 with the number of radial nodes. The conver- gence of the Fourier inversion procedure was checked by performing the calculations with 7 and 15 nodes of the Laguerre integration, the results remaining unchanged. It can be seen that with a number of 65 nodes the solution coincides with the solution with 129 nodes except near the axis r=0. It must be remarked that none of these solutions accurately sat- isfy the condition of a zero radial derivative at the axis, as shown in the enlarged view of the region near the axis in Fig. 4, although with increasing number of nodes the numerical error decreases. This was to be expected since, with the present quadrature solution method, the condition is not enforced explic- itly and, therefore, the accuracy of the solution will depend on the discretization error of the integrals appearing both on the °°5°~""1""1""1""1""~""1""1""1""1""~ 0.045 0.040 OD35 Cat to =) 0.030 , Q - Q 0.025 a, 0.020 Lo 0.015 m OD11 OD05 ~1 ~ ~ ~ !, , 1 I 1 ., I I 1 1 ODOR ODD O. 20 0.40 0.60 0.80 1.00 170 1.40 DO 1.BO 2.00 -,,~L~ . \+ x + NR=32 x KR=65 NO x x x x AX i; r/Rd Fig. 3. Convergence of interaction pressure at x=-0.5 with the number of radial nodes. Axisymmetric flow. c/Rd=0.5, a=0.8, CL=O.9. + - 32 X N1--65 029 OD40 _ <Jo 0.047 _ : Q. 0.046' i _ Q 0.045 a) : ~ OD44 _ on : a) OD43 [L : OD42 _ 0.041 r/Rd Fig. 4. Detail of computed interaction pres- sure near the axis at x=-0.5 for dif- ferent number of radial nodes. Axisymmetric flow. c/Rd=0.5, a=0~8, CL=O.9. 656

left and right-hand sides of the integral equation. It should also be noted that the numerical solution exhibits the expected de- gree of smoothness at the duct radius r=1. We insist in the importance of obtaining a smooth pressure distribution since the velocity field is derived from it by numerical differentia- tion. The radial component of the interaction velocity on different planes is shown in Fig. 5 for the solution with 65 nodes. The inter- action radial velocity is rather small. In Fig. 6 the axial velocity distribution due to the effect of the duct loading inside and out- side the duct is shown for this case. The ve- locity distribution of the oncoming flow is also shown in this figure. The induction ve- locity is simply the potential flow velocity in uniform flow divided by the local inflow velocity. This of course produces the large velocities close to the symmetry axis. The interaction velocity is negative and corrects to a certain extent this extreme behaviour. Nevertheless, for this strongly sheared inflow and, in contrast with the uniform flow case, the present method predicts in shear flow an increase to the axis of the disturbance veloc- ity due to the duct. Again near the axis the interaction velocity and thus, the disturbance velocity due to the duct is influenced by the local error in the pressure distribution. As a second application we considered a non-axisymmetric wake field defined by (85) with U1=0.72, U2=1.70, U3=2.0 and U4=0.5. The 0.030 n non O 0.010 - . _ ° 0.000 - ._ -0.010 -n non \/ \ ~` ,, . 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Radius r/Rd Fig. 5. Interaction radial velocity due to duct loading. Axisymmetric flow. c/Rd=0.5, a=0.8, CL=O.9. wake field is shown in Fig. 7. At all radii the velocity is lowest at O=0. (say the upper part of the propeller plane) and highest at O=180 deg. (the lower part of the propeller plane). To examine the capability of the duct to change this incoming velocity, two types of duct loading were considered: an Axisymmetric loading with the section lift coefficient CL=0.9, as in the previous Axisymmetric case, and an asymmetric loading defined by CL = CLo + ACL cos 0, with ACL = 0.15CLo; this cor responds to a 30% variation of duct section loading around the circumference, the highest loading being placed where the incoming veloc- ities are lowest and vice versa. The same chord-diameter ratio c/R =0.5 and the same chordwise load distributidon are assumed. The integral equation was solved using a number of 65 radial nodes, for 7 nodes of the Laguerre integration. A number of 3 circumferential harmonics was sufficient to carry out the so- lution in the present case. The induction, interaction and the total disturbance axial velocities at x=0 inside the duct are shown, respectively in Figs. 8, 9 and 10 for the duct with Axisymmetric loading. Similar results are shown in Figs 11, 12 and 13 for the case of asymmetric duct loading. In the case of Axisymmetric loading oppose site tendencies are found for the inner radii up to about 0.5 and the outer radii. For the outer radii the total disturbance velocities are lower at O=180 deg., where the incoming -- xlRd = +0.3 x/Rd= 0.0 ,' x/Rd = -0.3 / \ '' ',,, . it\ -0.1 -n ~ " / " . . . / ~ 1.0 0.9 0.8 0.7 0.6 - O 05 - 0.4 c' 0.3- o > 0.2 ct 0.1 0.0 - "' - UIUO/ .. / U/UO ~ -- Total=lnd.+lnt. Induction - Axisymmetric Wake 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Radius r/Rd Fig. 6. Axial velocity due to duct loading at x=O. Axisymmetric flow. c/Rd=0.5, a=0~8, CL=O.9. 657

velocity is large while for the inner radii large velocities are found at that angular po- sition. This is caused by the large positive interaction velocity calculated for the small- er radii, see Fig. 9 for r/Rd=0.30. The domi- nant term in eq. (63) is the second term which couples the negative radial velocities with a large positive value of the shear parameter a (note that for this wake field this parameter is larger than on the base axisymmetric wake). For the outer radii the interaction velocities are very small, and the total disturbance ve- locities are dominated by the induction part. In the case of asymmetric loading the induction part shows a more pronounced vari- ation, Fig. 11, as expected from the duct load variation. The level of interaction velocities decreases considerably at the inner radii, as shown in Fig. 12, especially at the angular position O=180 deg. where the local duct load- ing is lowest. The interaction velocities at the outer radii remain small. Still an in- crease of the total disturbance velocity to the axis is found in this case, Fig. 13. Clearly, an important result of practical in- terest is the better capability of the second duct with asymmetric loading to make the in- coming flow more uniform. out 1 0.7 l 0.6 o 1 = 0.5 ~ 1 1 ~ 1.0 0.9 Go 0.8 - = - - 0.7 o 0.6 a) > 0.5 Cal ._ X 0.4 0.3 0.2 n 1 1 ~r/Rd = 0.83 ~r/Rd = 1.00 )( r/Rd =0.61 ~r/Rd =0.71 r/Rd=0.42 )( r/Rd= 0.51 ~r/Rd = 0.00 (9 r/Rd = 0.30 o.o 30 60 90 120 150 180 210 240 270 300 330360 Angular Position (deg.) Fig. 7. Axial velocity of non-axisymmetric wake field. Eq. (85) with U1=0.72, U2=1.7, U3=2, U4=0.5. 0.45: 1 ~ r/Rd = 0.61 () r/Rd = 0.51 [S r/Rd = 0.30 0.40 ~ =~0.35 ~ 1 -_ ~ ) O0.4 0.2 0.1 0.0 - _ . _ o a) ct 0.30 . _ x in: n25 ~TT~ I I ___ 0 30 60 90 120 150 180 210 240 270 300 330 360 0 30 Angular Position (deg.) r/Rd = 0.96 r/Rd = 0.91 r/Rd = 0.83 r/Rd = 0.71 \\ ,; . . . . . 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) Fig. 8. Induction part of the axial velocity due to duct loading at x=O. Non-axisymmetric wake. Axisymmetric duct loading. c/Rd=0.5, a=0.8, CL=0.9. 658

5. SUMMARY AND CONCLUSION A linearized analysis of the three-dimen- sional steady interaction between a ducted propeller system, modelled by a suitably cho- sen external force field and source distribu- tion, and a non-axisymmetric sheared axial onset flow has been given. The basic assump- tion of the theory is that the disturbance velocities should remain small compared with the velocity of the onset flow. Consistent with this assumption, linearised boundary con- ditions are applied on the duct surface. As a consequence of the linearization the effects of duct loading and thickness and propeller loading can be treated separately and super- posed. For the three-dimensional problem a for- mulation in terms of the pressure seems an obvious choice and a linear integral equation has been derived for the disturbance pressure, which can be applied to describe both loading and thickness effects. By separating the po- tential pressure from the interaction pres- sure, the integral equation governing the in- teraction part may be in principle solved with great accuracy using suitable numerical proce- dures. The solution of the potential part and, more specifically, the computation of the cor- responding induced velocity field can be ob- tained from the results of duct lifting sur- face theory. As it stands, the analysis may be applied to solve the design problem of a duct in the ship's nominal axial wake field in the pres- ence of the propeller. This requires the spec- ification of the duct loading and thickness distributions, both circumferentially and 0.34 on o2l ~ 0.1 1 . _ 0 0.0 > Cal -0. 1 in: 1 -0.2 -0.3 -n4 / ~ / \ / ~ r/Rd = 0.61 \ An/ (9 r/Rd = 0.51 \ [S r/Rd = 0.30 , . . . . . 0 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) chordwise, and the propeller radial load dis- tribution. The computation of the disturbance velocities on the duct reference cylinder en- ables the determination of the duct camber and conical angle distributions. Due to the non- linear character of the interaction thrust forces between duct and propeller some itera- tion will be required to meet given duct and propeller thrust values. A computational scheme has been developed to solve the integral equation and evaluate the velocity field. The scheme has been ap- plied to investigate the influence of shear on the velocity field due to the effects of duct loading for axisymmetric and non-axisymmetric wake fields. Relatively low computational costs are associated with the utilization of the scheme when compared with a more direct numerical approach to the solution of the problem. Although the present results have not yet been validated by a proper comparison with other calculation methods or experimental data, the following conclusions may be drawn from the results of the sample calculations presented in this paper: - The effects of shear on the disturbance ve- locity inside the duct are found to be rather significant for the strongly sheared wake fields used in the calculations. Both for the axisymmetric wake and the non-axi- symmetric wake fields, the axial disturbance velocity distribution considerably differs from the distribution in uniform flow. - For the axisymmetric case a disturbance ve- locity distribution increasing to the sym- metry axis is found as a result of the ef- fect of shear. Also for the non-axisymmetric it< r/Rd = 0.96 r/Rd = 0.91 r/Rd = 0.83 r/Rd = 0.71 O o.oo .~ o > ._ -0.05 0 30 . . . . . 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) Fig. 9. Shear interaction part of the axial velocity due to duct loading at x=O. Non-axisymmetric wake. Axisymmetric duct loading. c/Rd=0.5, a=0.8, CL=0.9. 659

wake the effects of shear considerably af- fect the disturbance velocity distribution both radially and circumferentially at the inner radii up to r/Rd=0.5. For the outer radii the effects of shear become much smaller. - A large effect on the disturbance velocities on the inner radii arising from the effects of shear, is found if the duct circumferen- tial load distribution is changed. The foregoing results certainly indicate the need for pursuing the development of a model for the interaction between the propel- ler, duct and wake field, if the design of wake-adapted ducts using design criteria for flow rectification is aimed at. At present, work on this linearised model will concentrate on the application of the method to duct de- sign including the effect of duct thickness and propeller loading. In addition, steps to validate the numerical results of the model will be undertaken. REFERENCES 1. Oosterveld, M.W.C., "Wake Adapted Ducted Propellers," Doctor's Thesis, Netherlands Ship Model Basin Publ. No. 345, 1970, Wageningen, The Netherlands. 2. Van Manen, J.D., and Oosterveld, M.W.C., "Analysis of Ducted Propeller Design," Trans- actions of the Society of Naval Architects and Marine Engineers, Vol. 74, 1966, pp. 522-562. 3. Turbal, V.K., "Theoretical Solution of the Problem on the Action of a Non-Axisymmet- rical Ducted Propeller System in a Non-Uniform Flow," Proceedings of the Symposium on Ducted Propellers, The Royal Institution of Naval Architects, 1373, pp. 11-19. °9 owl 0.7 l o 0.6. - ~ 0.5 o > 0.4 Ct x 6 0.3 0.2 0.1 ~' / \ ~ r/Rd = 0.61 C) r/Rd = 0.51 r/Rd = 0.30 . o 4. George, M.F., "Linearised Theory Applied to Annular Aerofoils of Non-Axisymmetrical Shape in a Non-Uniform Flow," Transactions of the Royal Institution of Naval Architects, Vol. 120, 1977, pp. 201-210. 5. Dickmann, H.E., and Weissinger, J., "Beitrag zur Theorie optimaler Dusenschrauben (Kortdusen)," Jahrbuch Schiffbautechnischen Gesellschaft, Band 49, 1955, pp. 253-305. 6. Morgan, W.B., ''Theory of the Annular Air- foil and Ducted Propeller,' Fourth Symposium on Naval Hydrodynamics, Office of Naval Re- search, ACR-92, 1962, pp. 151-197. 7. Weissinger, J., and Maass, D., "Theory of the Ducted Propeller - A Review," Seventh Sym- posium on Naval Hydrodynamics Office of Naval , Research, DR-148, 1968, pp. 1209-1264. 8. Morgan, W.B., "Some Results From the In- verse Problem of the Annular Airfoil and Ducted Propeller," Journal of Ship Research, Vol. 13, No. 1, March 1969, pp. 40-52. 9. Dyne, G., "A Method for the Design of Ducted Propellers in a Uniform Flow," Publ. No. 62 of the Swedish State Shipbuilding Experimental Tank, 1967, Goteborg, Sweden. 10. Dyne, G., "An Experimental Verification of a Design Method for Ducted Propellers," Publ. No. 63 of the Swedish State Shipbuilding Experimental Tank, 1968. 11. Falcao de Campos, J.A.C., "On the Calcu- lation of Ducted Propeller Performance in Axi- symmetric Flows," Doctor's Thesis, Netherlands Ship Model Basin Publ. No. 696, 1983, Wagen- ingen, The Netherlands. 12. Kerwin, J.E., Kinnas, S.A., Lee, J.-T. and Shih, W.-Z., "A Surface Panel Method for the Hydrodynamic Analysis of Ducted Propel- lers," Transactions of the Society of Naval Architects and Marine Engineers, Vol. 95, 1987, pp. 93-122. 0.50: 1 0.451 0.40 >, 0.35 .~ 0 1 o - 0.30 Ct . _ ~: 0.25 n20 r/Rd = 0.96 ~ r/Rd = 0.91 C) r/Rd = 0.83 r/Rd = 0.71 . . . . . . . . 0 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) . . . . . . . . . . . . 0 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) Fig. 10. Total disturbance axial velocity due to duct loading at x=O. Non-axisymmetric wake. Axisymmetric duct loading. c/Rd=0.5, a=0.8, CL=0.9. 660

13. Kinnas, S.A. and Coney, W.B., "A Syste- matic Method for the Design of Ducted Propel- lers " Fourth International Symposium on Prac tical Design of Ships and Mobile Units, 1989, Varna, Bulgaria. 14. Feng, J. and Dong, S., "A Panel Method for the Prediction of Unsteady Hydrodynamic Performance of the Ducted Propeller with a Finite Number of Blades," Proceedings of the International Symposium on Propeller and Cavi- tation, 1986, Wuxi, China. 15. Huang, T.T., Wang, H.T., Santelli, N. and Groves, N.C., "Propeller/Stern Boundary Layer Interaction on Axisymmetric Bodies: Theory and Experiments," DTNSRDC Report 76-0113, Dec. 1976. 16. Huang, T.T. and Groves, N.C., "Effective wake: Theory and Experiment," Thirteenth Sym- posium on Naval Hydrodynamics, The Shipbuild ing Research Association of Japan, 1980, pp. 651-673. 17. Dyne, G., "A Note on the Design of Wake Adapted Propellers," Journal of Ship Research, Vol. 24, No. 4, 1980, pp. 227-231. 18. Hawthorne, W.R., "On the Theory of Shear Flow," Gas Turbine Laboratory Report No. 88, Oct. 1966, Massachusetts Institute of Tech- nology, Cambridge, Massachusetts. 19. Von Karman, T. and Tsien, H.S., "Lifting- line Theory for a Wing in Nonuniform Flow," Quarterly of Applied Mathematics, Vol. 3, 1945, pp. 1-11. 20. Lighthill, M.J., "The Fundamental Solu- tion for Small Steady Three-dimensional Dis- turbances to a Two-dimensional Parallel Shear Flow," Journal of Fluid Mechanics, Vol. 3, 1957, pp. 113-144. 21. Weissinger, J., "Linearisierte Profil- theorie bet ungleichformiger Anstromung," Acta Mechanica 10, 1970, pp. 207-228. 0.7 l 0.6 t o 0.5 .~ O 0.4 X 0-3 0.2 0.1 0.0 22. Weissinger, J., "Linearisierte Profil- theorie bei ungleichformiger Anstrdmung. Teil II: Schlanke Profile," Acta Mechanica 13, 1972, pp. 133-154. 23. Weissinger, J. and Overlach, B., "Grund- lagen zu einer Theorie des Ringflugels in axialsymmetrischer Scherstromung," ZAMM 55, 1975, pp. 413-421. 24. Overlach, B., "Linearisierte Theorie der axialsymmetrischen Stromung um Ringflugel bei ungleichformiger Anstromung," Dissertation, 1974, Karlsruhe. 25. Goodman, T.R., "Momentum Theory of a Pro- peller in a Shear Flow," Journal of Ship Re- search, Vol. 23, No . 4, Dec. 1979, pp . 242- 252. 26. Falcao de Campos, J.A.C. and Van Gent, W., "Effective Wake of an Open Propeller in Axisymmetric Shear Flow," Netherlands Ship Model Basin Report No. 50030-3-SR, May 1981. 27. Van der Vegt, J.J.W., "Actuator Disk in a Two-Dimensional Non-Uniform Flow," Interna- tional Shipbuilding Progress, Vol. 30, No. 348, Aug. 1983, pp. 158-178. 28. Van Gent, W., "A Model of Propeller Ship Wake Interaction," Proceedings of the International Symposium on Propeller and Cavitation, 1986, Wuxi, China. 29. Lee, H., "Ducted Ship Propellers in Ra- dially Sheared Flows," Ph.D. Thesis, Stevens Institute of Technology, 1985. 30. Lee, H., "Effects of Radially Sheared Inflow on the Design of Propeller Ducts," Third International Symposium on Practical Design of Ships and Mobile Units, 1987, Trondheim, Norway. 31. Weissinger, J., "Zur Aerodynamik des Ringflugels I. Die Drukverteilung dunner, fast drehsymmetrischer Ringflugel in Unterschall- stromung," D.V.L. Bericht Nr. 2, 1955. 0.45 . r/Rd = 0.61 r/Rd = 0.51 - Cl r/Rd = 0.30 ~.< 0.40 o , 0.35 .~ o > ct 0.30- x 0.25 r/Rd = 0.96 r/Rd = 0.91 r/Rd = 0.83 [1 r/Rd=0.71 _ ~ 1 1 1 1 0 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) , . . . . . 0 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) Fig. 11. Induction part of the axial velocity due to duct loading at x=O. Non-axisymmetric wake. Asymmetric duct loading. c/Rd=0.5, a=0.8, CLo=0.9. 661

32. Weissinger, J., "Zur Aerodynamik des Ringflugels III. Der Einfluss der Profil- dicke," D.V.L. Bericht Nr. 42, 1957. 33. Ordway, D.E., Sluyter, M.,M. and Sonnerup, B.O.U., "Three-Dimensional Theory of Ducted-Propellers," TAR-TR-602, Aug. 1960, Therm Advanced Research, Ithaca 1, New York. 34. Hough, G,.R., and Ordway, D.E., "The Generalized Actuator Disk," Developments in Theoretical and Applied Mechanics, Vol. II, Pergamon Press, Oxford [etc.], 1965, pp. 317-336. 0.5 0.4 0.3 o , 0.1 .~ o 0 Ct FAX 0.2 0.0 -0. 1 -no ._ -0.3 -0.4 rn m / \ r/Rd = 0.61 () r/Rd = 0.51 r/Rd = 0.30 0 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) 35. Batchelor, G.K., An Introduction to Fluid Dynamics, Cambridge University Press, - Cambridge, 1967. 36. Baker, C.T.H., The Numerical Treatment of Integral Equations, Clarendon Press, Oxford, 1977, pp. 375. 37. Erdelyi, A., ea., Tables of Integral Transforms, Vol. 1, McGraw-Hill, New York, 1954, pp. 106. O o.oo .~ o ct <: -0.05 r/Rd = 0.96 r/Rd = 0.91 C) r/Rd = 0.83 r/Rd = 0.71 . . . . . . . . . . . O 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) Fig. 12. Shear interaction part of the axial velocity due to duct loading at x=O. Non-axisymmetric wake. Asymmetric duct loading. c/Rd=0.5, a=0.8, CLo=O.9. 0. o.8 0.7 o 0.6 - ~ 0.5 o > 0.4 Ct ._ x 0.3 0.2 0.1 \ ~ / ~ r/Rd = 0.61 () r/Rd = 0.51 r/Rd = 0.30 o.so n45 n~n o - >~ 0.35 o ~D 0.30 0.25 0.20 ) ~ ' 5 ~ _] . . . . . . . . . . . O 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) Fig. 13. Total disturbance axial velocity due to duct loading at x=O. 0 30 60 90 120 150 180 210 240 270 300 330 360 Angular Position (deg.) Non-axisymmetric wake. Asymmetric duct loading. c/Rd=O.S, a=0.8, CLo=O.9. 662

APPENDIX A Shear-Dependent Pressure Disturbance Term T ~n°) Duct Loading The term is w T ~n° ) = Rd ~ G( )(r;~) APd ( \) n=-- mn n with Gd )(r;~) = 21XlKn(l~lRd) [Km(l~lr) h1(r) + mn Im(lAlr) h2(r)] 21Al In(lklRd) + {m(l~lr) h4(Rd) ' r<Rd ~ (88) (L) dmn(r;~) = 21\l Kn(l~lRd) Km(lklr) h1(Rd) + 21Xl In(lklRd) [Km(lAlr) h3(r) + Im(lklr) h4(r)] with hi(r) = ~ G12(a) In( 1 Nl a) da, (90) R h2(r) = ~ G22( a) In( 1 Nl a) da , (91) r h'(r) = ~ G12(~) Kn(lkl a) da ~(92) Rd co (r) = r G22(a) Kn( l ~l a) da r r , r>Rd ~ (89) and (93) G12(r) Amn(r) Im(l~lr) + Bmn(r) 1\l Im 1(l~lr) , (94) G22(r) = Amn(r) Km(l~lr) Bmn(r) |~| Km 1(l~lr) Duct Thickness T ~n°) = ipX Rd where Gd )(r;~) = 2 Kn(~Rd) [Km(~r) h1(r) + Im(~r) h2(r)] + 2 In(~Rd) Im(~r) h4(Rd) ' r<Rd ' (97) (T) (87) dmn(r;~) = 2 Kn(~Rd) Km(~r) h1(Rd) -2 In(~A~Rd) [Km(~r) h3(r) Im(~\lr) h4(r)] , r>Rd . (98) Propeller Loading (O) ~ p (L) m with ~ + (99) G( )(r;~) = 2 Km(|~|r) i1(r) + m 2 Im(lAlr) i2(r) + 2 Im(l~lr) x ip(Rp) ~ G22(a) Ko( 1 ~l a) da, p r<Rp , (100) G( )(r;~) = 2 Km(|~|r) [i1(Rp) + ip(Rp) ~ G12( a) Ko( | \| a) da] + 2 Im(l~lr) where co x ip(Rp) ~ G22( a) Ko( I \1 a) da ~ r>Rp , (101) r i1(r) = ~ G12( a) [Ko( | \| o) ip( a) + · (95) a) ~ Gd )(r;A) Tn(~) (96) n=-= mn l2( and Io( 1~1 a) kp(a)] da ~ r<Rp , (102) Rp (r) = r G22( a) [Ko( 1 \1 a) ip( a) + Io( I \1 a) kp( a) ] do , r<Rp , (103) r m(r) = ~ Io(|~|~) lipp(a) a da, (104) R kp(r) = [P KO(|~|a) lipp(a) a do (105) 663

APPENDIX B Radial Downwash v(O) Due To Duct Loading The radial momentum equation is U TV ~ = 1 ~ _ + p ar lip p [H(x+c)-H(x-c)] 6(r-Rd) Taking Fourier transforms we obtain iApU v(O) = ~; mn r a9( ) ~ . (106) L ar + ~Pdm(~) 6(r~Rd)] . (107) m=-= Evaluating the derivative in eq. (107) with (53) we obtain with iNpU v(O) = a, = - ~ aim 2 Rd \2 APd (A) m=-= m x [Im(|)~|r) Km(IAIRd) + Km(l~lr) Im(I I d)] (108) At the duct r=Rd we obtain, after inverting the Fourier transforms v(O) and APd ' m v(O) (x,Rd, O) = = - 4i U ~ mime [2 Rd J. Ad ( A) do m=-= -c m co x 1~ ~ IK( ?( I AIR ~ e \( () d\] , (109) ~(m-l)'' ~'-~d' with IK(m 1) = (Im+1 + Im-1)(Km + Km-1) (110) In terms of non-dimensional velocities and pressures v(O)(x,R,,8) = - ~ VmO)(x,Rd) e me, (111) m=-= with Vm0) = 4~ R- r And ( i) Fm(n) do, (112) d -1 m where 0o m(b) = - 2- ~ ~ IK(m+l)(|)~|) elan dA (113) and = -(x-~) Rd (114) The Fourier integral (113) can be written in the form co Fm(n) = J. ~ IK(m_l)(~) sin (\n) dX (115) and may be evaluated in terms of the Legendre functions of second kind and its derivatives. Using the results for the Fourier transforms of the modified Bessel functions [37] it may be shown that P (ink) Fm(n) = 2 m rat + mu, (116) Pm( ~ n ~ ) = 2 [(m-1)Qm+l/2(°) (m+l)Qm_3/2(~) + 2mb Lm(~)] , (117) where Lm(cc) is given by the recurrence rela Lm(~) = Lm_1(~) - 2 [0m-1/2( ) em 3/2(°)] ~m22 , (118) Lo(~) = 0, L1(~) = ~ 2- [Q 1/2(~) + Qi/2( )] ' (119) and a= 1 + 2~ . With eq. (116), eq. (112) becomes (120) m d 27i Rd ~_~ ~ do m '21 ~ Ad (a) d(] (121) -1 m Qm+1/2(~) = -A + O(ln~ nit ) for small values of A, Pm(0) = 1 and the first integral in eqs. (3-36) has a Cauchy singularity. The result (121) agrees with the result given by Weissinger in ref. [31]. Further reduction of eq. (121) including the treatment of the Cauchy singularity can be found in the orig- inal reference. 664

DISCUSSION All H. Nayfeh Virginia Polytechnic Institute and State University, USA The convergence problem near the axis may be alleviated if one uses an analytical rather than a numerical solution there. Such a procedure was used by Nayfeh, Kaiser, and co-workers in the seventies to treat acoustic waves propagating in circular ducts. The results were published in the AIAA Journal, Journal of Sound and Vibration, and the Journal of the Acoustical Society of America. AUTHORS' REPLY The author thanks Prof. Nayfeh for pointing out the possibility of using an analytical solution in the vicinity of the axis. Indeed, see for instance references [1] and [2]; the form of the solution near the axis is determined by the potential part of the solution, even in the absence of a uniform flow core and certainly in the case of bounded shear parameters a and b at the axis. In this case, the analytical solution for each circumferential pressure harmonic behaves like the modified Bessel function Im multiplied by a constant. Its asymptotic behavior for small arguments can then be used to determine the constant by collocation of the integral equation (45) at a point off the center line but sufficiently close to it. [1] Eversman, W., Effect of Boundary Layer on the Transmission and Attenuation of Sound in an Acoustically Treated Circular Duct," Journal of the Acoustical Society of America, Vol. 49, No. 5, 1971, pp. 1372-1380. [a] Nayfeh, A. H., Kaiser, J. E. and Telionis, D. P., Acoustics of Aircraft Engine-Duct Systems," AIAA Journal, Boll 13, No. 2, 1975, pp. 130-153. DISCUSSION William B. Morgan David Taylor Research Center, USA This is a very interesting paper on the design problem of ducted propellers. I have a question concerning the optimum dueled propeller in a shear flow. Can you say anything about how to calculate the optimum ducted propeller considering both the optimum propeller and the optimum duct shape in a shear flow? AUTHORS' REPLY Dr. Morgan raises the question of how to optimize both the propeller and the duct in a shear flow. Although the interaction with shear covers only a particular aspect of the propulsor-hull interaction, the question may be addressed without considering in detail the shear producing mechanism which is the presence of the ship's hull with its boundary layer and wake. From an untheoretical point of view, it would be interesting to know for a given sheared inflow velocity field what are the load distributions on the duct reference cylinder and on the propeller disk which minimize the kinetic energy of the fluid left far behind the ducted propeller system and which satisfy a given thrust constrains". Such optimization would lead to the d,ucted propeller loading for maximum recover of the kinetic energy of the fluid present in the sheared incoming flow. The formulation of such optimization problem has not been attempted. From a more practical point of view, the present model may be conceivably used to determine the duct shape for maximum attainable rectification of the inflow to the propeller, for given propeller radial load distribution in a given wake field. This is a desirable feature from the standpoint of avoiding cavitation on the propeller blades which may also lead to efficiency improvement. 665