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 Potential Based Panel Method for the Unsteady Flow Around Open and Ducted Propellers S. Kinnas, C.-Y. lIsin, D. Keenan (Massachusetts Institute of Technology, USA) Abstract The unsteady flow around an open marine propeller subject to a spatially nonuniform inflow is analyzed by utilizing a time marching potential based panel method. An efficient algorithm is implemented in order to ensure an explicit I(utta condition (i.e. pressure equality) at the blade trailing edge at each time step. The numeri- ca.l method is shown to be very robust for a broad range of reduced frequencies as well as con.sistent with known analytic solutions and with an existing unsteady lifting surface method. A hybrid panel method is developed for the analy- sis of the unsteady flow around ducted propellers. It combines an unsteady lifting surface method for the propeller with a potential based panel metl~od for the duct. The propeller is essentially treated with an ex- isting time marching vortex lattice scheme which has been developed for open propellers, with the effects of the duct being accounted for via the generalized images of the propeller singularities with respect to the duct. The proposed method is shown to be appropriate for treating unsteady ducted propeller flows in a computa- tionally efficient and robust manner, especially when a given ducted propeller geometry must be analyzed for various inflow conditions. ~ Introcluction In most marine applications the propeller is subject to severe non-axisymmetric wal;es produced from the boundary layer of the vehicle. Therefore unsteady ma- rine propeller flows are a very important aspect of the overall propeller / hull hydrodynamic analysis problem. One of the important steps in determining the co~n- plete propeller/wal;e interaction, is the analysis of the unsteady flow around a propeller in the presence of a. prescribed spatially nonuniform inflow. Accurate pre- dictions of the unsteady pressure distributions on the blade, especially at the blade leading edge and tip, are crucial in determining cavitation inception, unsteady boundary layer separation as well as leading-edge vor 667 tex separation. Numerical lifting surface methods have been exten- sively applied for the analysis of unsteady propeller flows. One of the first investigators to formulate the un- steady propeller lifting surface problem in terms of the acceleration potential was Hanaoka t7,84. The acceler- ation potential formulation was implemented, in terms of finite number of chordwise nodes, by Tsal;onas et al. t32,314. More recently, the unsteady vortex lattice technique was employed by I`erwin and C.S. Lee t174. The unsteady vortex lattice technique was extended by Keenan t14] to include unsteady vortex wake relaxation. A review of the different steady and unsteady lifting sur- face methods as applied to marine propellers has been given by Kerwin t164. One of the major drawbacl;s of the lifting surface methods though, is their inherent failure at the blade leading edge and tip where the blade thickness effects are substantial. On the other hand, panel methods have been applied for the analysis of the steady flow around marine propellers, including the hub 1)Y, among oth- ers, Hess and Valarezo t10], Kerwin, Ixinna,s, J.T. Lee and Shih (who also included the duct) AS] and Hoshino t12~. The first of those methods employs a, source based formulation [9], and the other torso employer a potential based formulation t27~. Panel methods were first extended to treat unsteady flows around two dimensional hydrofoils by Giesing Eli, who employed a source based formulation, and then by Basu and Hancock t1] who employed a. surface vorticity formulation but also allowed for wake relaxation. More recently, panel methods have been extended to treat the unsteady flow around helicopter blades by Maskew t25] and by Morino, I`aprielian and Sipcic t26i, with emphasis given to the evolution of the rotor free wakes. In the present worl;, a time marching potential based panel method is developed for the analysis of the un- stea.dy flow around marine propellers, with emphasis on the accurate predictions of unstea.d~r blade pressures and forces for a large range of blade rate harmonics. A computationally efficient explicit pressure Kutta condi- tion is implemented in the panel method in order to ensure pressure equality at the blade trailing edge at

each time step. The numerics of the panel method are shown to be robust with respect to tl~e size of the time step and the number of panels on the blade for a broad range of reduced frequencies. The nun~erica,1 method is also shown to be consistent with existing analytic solu- tions as well as with an unsteady lifting surface method. Finally, a hybrid lifting surface / potential leased panel method is developed for the analysis of unsteady flows around ducted propellers. The propeller is treated with a time marching vortex lattice method, with the effects of the duct being accounted for via the gener- alized images of the singularities representing the pro- peller and its trailing wake with respect to the duct. The generalized image idea was first applied to the analysis of the flow around a propeller in the presence of a duct by Kinnas and Coney t214. The propeller was modelled in lifting line theory by using a, finite number of semi-infinite helical vortex horseshoes. The effect of the duct on the propeller was accounted for via, the generalized images of the horseshoes with respect to the duct. The flow field of tl~e generalized image of each horseshoe was computed lay sols id for the non- a.xisym~netric potential around the duct in the presence of the horseshoe. Y As Is Figure 1: The propeller in a spatially nonuniform in- flow. 2 The Unsteacly Pane} Methoc! for Open Propellers 2.1 Formulation Consider a propeller subject to ~ spatially nonuniform inflow Us~xs,rs,bs), as shown in Figure 1. The in- flow is taken with respect to the absolute (ship fixed) system of cylindrical coordinates xS,rs,bs. The flow around the propeller will be analyzed with respect to the propeller fixed system (x,y,-), also shown in Fig- ure 1. If the propeller rotates It,] a.r~gula~ velocity-~ (i.e. right-handed), then the inflow relative to the pro- peller, Uin, will be time dependent arid given as _ ~ Uirl (x, y, a, i) = Us(x, r, O-At) + ~ x ~, where r = I/, 0 = arctan(-/y) and .r = (x,y,-). Ale assume at this point, that the interaction be- tween the propeller and the inflow is inviscid and ir- rota.tional i. Then, the time dependent total flow ve- locity relative to the propeller fixed system, Ox, y, z, i), can be written in terns of the perturbation potential, 4(x, y, z, t), as follows: q(x,y,z,t) = t~in(`x,y,z,t) + V¢(~.r,y,z,t). (2) By applying Green's formula. for ¢( .r,y,z,t) at any time t, we will get the following integral equation for the perturbation potential Up at every point p on the propeller blade surface So: 2~>p = Asp [¢>q 0~ -G(p; q) ion ~ dS fsw Gil q, A, t) `' q) dS, (~3) with the subscript q corresponding to the variable point in the integrations; n is the unit vector normal to the propeller surface or to the wake surface, i\: is the po- tential jump across the wake sheet, and G(p; q) is the Green's function. In the case of unbounded three di- mensional fluid domain, G is given as G(<p;q)= Rid >~, (4) with R(<p; q) being the distance between points p and q. Equation (3) expresses the potential on the pro- peller blade as the superposition of the potentials in- duced by a continuous source distribution, G. on the propeller surface, Sp, and a continuous dipole distribu- tion, {3G/0n, on the prop ells surface Sp and its walls Sw. The strength of the source distril:'ution is given, via the kinematic boundary condition, as ion =-tJ~in(Xq' yq, A, t) 'A (5) where xq7yq,zq are the coordinates of point q with re- spect to the propeller fixed system. The strength of the dipole distril~ut.ion is unknown and equal to the perturbation potential on the pro- peller or to the potential jump in the wal;e. The dipole strengths will be determined by inverting integral equa- tion (3~. The position of the trailing wal;e, All,, is assumed to be invariant with time and taller to be the same as that corresponding to the circun~ferentia,ll~ averaged inflow ilk. The dipole strength Afar, O. t) in the walls, is convected along the assumed walls model with angular 1 this is equivalent to assuming that ds is the effective wake. 668

speed as, in order to ensure that the pressure jump in the wake is equal to zero, i.e. i\~(r,b,t) = i\~(r,~T(r),t- co ) = /\; ( ~-bT(ri ) t > ~-bT(r) , ad, ~ '\~(r,Y,t) = /\4S(r); t ~ ~(6) where r,§ are the cylindrical coordinates of the wake surface, Sw, and 0~ (r) is the ~ coordinate of the pro- peller blade trailing edge ati radius r. ~\¢S(r) is the steady flow potential jump in the `val;e when the pro- peller is subject to the circumferentia.lly averaged in- flow. For t < 0 we assume that the propeller is subject to the circumferentially averaged inflow. The unsteady inflow is "turned on" at t = 0. The value of the dipole strength, /~¢T(r,t), at the trailing edge of the blade at time t, will be given by A¢T(,,t) = SbT(~,t) -~T(~',t) = F(r,t) (7) where o+tr,t) and ¢>T(r,t) are the values of the poten- tial at the upper (suction side) and lower (pressure side) blade trailing edge, respectively, at time I. The differ- ence in those potentials is also equal to the circulation It at time t around the Claude section at radius r. The condition (7) is equivalent to requiring the shed vortic- ity from the blade trailing edge to lie proportional to the time rate of change of the circulation around the blade. Me ~ MOW Figure 2: Panel arrangement for a propeller blade and its wake, N = 20, 1~1 = 10, by = 6°. 2.2 Numerical Implementation Equation (3) is a Fredholm integral eclua.tion of the sec- ond kind with respect to op. To solve this equation nu- merically, we discretize the propeller and wake surface into quadrilateral panels. A typical panel arrangement for a propeller and its wal;e is shown in Figure 2. The propeller blade panel arrangement is identical to that used in the steady propeller panel method t244. The time domain is also discretized into equal time intervals /\t. The panels in the wake start from the blade trailing edge and their edges are located on the prescribed wale surface at equal angular spacing AOltr, which is related to the time step At as follows: ~61~,~.' = wAt. (S) On each of the quadrilateral panels, tl~e dipole or source distributions are approximatecl lay constant strength distributions. The discretized aversion of the integral equation (3) is applied at the centroids of each of the propeller panels and at each time step n = t/!\t. This renders the following system of linear equations: NB NP NE] M NW it, it, ail' f ji' (n) + ~ ~ ~ W.~Iim li/\t¢)m l(n) k=1 j=1 K=1 m=1 1=1 NB NP = ~ ~b~;j~K(n); i = 1,NP X NB (9) k =1 j=1 where 1\TB is the number of blades, and for each blade, N is the number of cl~ordu ise panels, .71 I is the number of spanwise panels, ^7\7p-N x -lI is the total number of panels and burl, is the nu~nl~er of chord ise panels in the wake. The influence coefficients ai'J a.ncl b~r\J are defined as the potentials induced at panel i lay unit (constant) strength dipole and source distributions, respectively, located at panel j on blade I`'. The walk influence co- efficients I7V-~]im ~ are defined similarly. The shape of the surface bounded by the edges of each panel (nonplanar in general), is approximated by a hyperboloidal surface and the corresponding influence coefficients are evalu- ated by using analytical expressions thy], [2S;, t134. The use of hyperboloidal surface panels instead of planar panels, has been found crucial for the consistency and convergence of the steady flow propeller panel method t24], especially, when applied to extreme propeller ge- ometries (i.e. high skew and twist) t13~. A similar con- clusion has also been reached by Hoshino t12~. The source strength ~K(n) is defined from equation (5) as it,' (n) = [Jin(~', y, ~ A , 71 /\t) · nK (10) with XK,YA,ZI~ being the coordinates of the centroid of panel j on blade I`' and D`.K being the unit vector normal to that panel. The terms in equation (9) can be regrouped as (the superscript 1, which denotes the lied blade, is omitted): NP M ~'ai,j¢)j~n)+ ~' Jiv~i,m,~(>m,~(l7) = RH,Si(~; i = 1, NP j=1 m=1 (11) NB NP NB HP RHSi(7l) = ~ bI gal (71) - ~ ~ U! j¢'j (Gil K=! j=1 K=2 j=1 ND M NW , . ~l NlV - ~ ~ ~ Tail' if\ ¢)m `(n) - ~ ~ tVi,m,l/~¢m,l(~) k =2 m=1 l=! m=1 l=2 (19) 669

The system of equations (11) must be solaced at each time step n with respect to the potentials on the key blade, K = 1. The potentials on the other blades (It' > 2) and their wal;es, which appear on the right- hand side of equation (11), are taken equal to those on the key blade at an earlier time step, when the blade It' was at the place of the key blade. This scheme, already used in an unsteady propeller lifting surface vortex lat- tice method [19], is essentially an iterative method of solving the system of equations (9) in order to deter- mine the steady state oscillatory solution for a given inflow. The solution, in most cases, was found to con- verge after three propeller revolutions. Solving equation (11) instead of equation 9 requires the inversion of a much smaller matrix and thus reduces substantially the computing time. It also requires less computer storage, since only the solution for the lay blade, rather than for all blades, needs to be stored at each time step. The values of /~¢m,~(n) for I ~ 2 (i.e. for all panels in the wake, except those next to the trailing edge) are determined from the discretized form of equation (6~: i\~)m,'(n) = i\¢m,~(n-1 + 1~; I > 2 ~ n > 1, /\~)m,l~n) = /\¢m; I > 2 ~ 77 < 1. (13) On the other hand, the value of the potential /\¢m,~, at the first wake panel mall be given, by using equation (7), as i&¢ Among + rttl(~7 - 1) Aim with Emend denoting the circulation around the blade strip m at time step n. = tRm(n-3) A¢3,m(n) Figure 3: Expanded view of the dipole distribution in the wake of blade strep m. In deriving equations (13) and (ll), the continu- ous dipole distribution over each wake panel is approx- imated by a constant distribution with strength equal to the mean value of the potentials at the edges of the panel, as shown in Figure 3. Numerically, Em~n) is approximated by laming = Among-amid (15) where ¢+ and Em are the potentials at the upper and lower trailing edge panel, respectively, at blade strip m. Equation (15) is an extension of LIorino's Kutta con dition in steady flow [97], though different from the un- steady Kutta condition employed by NIorino et. al. in [264. The right-hand side of equation (15) approximates the jump in the potential at the trailing edge with the difference of the potentials at the trailing edge control points on the blade, as shown in Figure 3. An i~npro~e- ment on this approximation will be presented in the next section. By substituting equations (15) and (14) into equa- tion (11) we end up with a, linear system of equations with respect to the unknown potentials Ink on the key blade. This system.of equations is inverted at each time n. A modification on the previously presented numeri- cal scheme has been implementecl, in which the dipole sheet on the first wake panel is approximated by a linear instead of a constant distribution. As will be described in section 2.5, that modification is necessary to Natalie the results of the unsteady panel method practically in- sensitive to the size of the time step. In that case, the system of equations (11) will be modified to be Np Al Ad, ai,j f j(n) + Ad, T,;tmEm (n) = RlI.5, (n ~ , i = 1, Np, j=1 m=1 (16) with RHSi(n) = RHSi(n)-Ti,F'mEn'(17-1); i = 1,Np (17) and where TO ~ TR) are the influence coefficients of a linear dipole distribution (in the chord wise direction) with unit value at the left - i.e. trailing edge - (right) edge and zero value at the right (left) edge. Substituting equation (15) in equation (16) we end up, at each time step, with a system of Np linear equations with respect to the Np unknown potentials on the panels of the l;ey blade. Equation (16) can also be written in the following matrix form [A] [I] Jr [T] [I] = [Rl'-I ,s ] ( 18 ) with [A], [<I] and [T] being the matrices of ai,j's, fj(n)'s and T',Lm's, respectively, and [r] = [rl(n),...,r',r(77)]T, [RHS] = [RH3l (n), , RIl ,s'Np (~7 )] ( 1 9) 2.3 Unsteady Pressures and Forces on the Blade The pressure, p, at every point .~. on the blade and for each time step n is computed from the potentials on the blade by applying the unsteady flow Bernoulli's equa- tion with respect to the propeller system: 670

P + ~ q2 + 0; = PI + i q2 (20) where p is the density of the fluid, q is the magnitude of the total surface velocity with respect to the propeller system at point x and time step n, pa and qOO are the pressure and the magnitude of the velocity of the un- perturbed flow, with respect to the propeller system, respectively. The surface velocity q is computed in terms of the derivatives of the perturbation potentials on the blade surface, in the same way as in the steady flow panel method t18~. The b¢/0t term in equation (20) is com- puted numerically, by implementing a, backward finite difference- second order accurate- scheme with respect to time t13~. The inviscid unsteady forces acting on each propeller blade are computed by integrating the unsteady pres- sure distributions on the surface of the blade. 2.4 The Unsteady Pressure Kutta Con dition and the Base Problems In the steady flow panel method foil propellers and ducts t18i, it was found that the nume~ica.l Kutta condition (15) t27] was not sufficient and an iterative scheme was developed in order to modify the solution and ensure pressure equality at the trailing edge. An explicit pressure I`utta condition was also found to be required in the unsteady panel method. For ex- ample, in Figure 4 the unsteady pressure distributions before and after applying the iterative pressure condi- tion, are shown for the N4118 propeller in the wake field described in section 2.5.2. The iterative method to en- sure pressure equality at the tailing edge, is described next. The pressure difference, /\pnt, at the trailing edge of each blade strip m,, is defined as /\pn, (n) = p+ (n)-Pm (n) (21 ) P-P p(up+w ~ ) 8:198 0.250 0.125 O. 000 -O. 1004 , ~ - Before Pressure Kutta Condition . few - ~ i ._ _ ~ r/R = 0. 60 loo ' ' 01.200 ' ' ' 0.450 ' ' ' 01.700> X/C r/R = 0. 30 r/R = 0.80 Figure 4: Pressure distributions vs. choldwise location at different propeller radii, r/R, for the N4118 propeller at blade angle Is = 300°. Jp = 0.S33, ug = 0.2. Before and After applying an iterative pressure I`utta condi- tion. ) = tt*~(~) _ tJ(~)i-~ i,Np*~(k) (24) where k denotes the iteration number and t~p*~(k) are the pressure differences which correspond to the solu- tion t¢*~(k) of the system of linear equations tA]t¢*] (k) = fRHS] _ fTitri (k) (25) The first iteration (k = 1) is taken equal to the solution when the Kutta condition (1.5) is applied as t¢*;(~) = t¢l, Or*; = try (26) In equation (24) the elements of the Jacobian matrix tJ(k)] are defined as with p+ and Pm being the pressures at the upper and J(k) = 0~Pt . (27) lower trailing edge, respectively. We call oj* and F*m the potentials and the circula tions, respectively, which satisfy equation 18 and which, in addition, produce a zero pressure difference (within a desired tolerance) at the blade trailing edge. That is, [A]t¢*] + fT] tr*1 = fRHSi, (22) and /\pm~n) = 0; m = 1,1tI. (2.3) Because of the nonlinear dependence of /\Pm on 4*, an iterative method is used in solving the system of equations (22) and (23) with respect to the unknown 4,1 's and flm Is. The following AI-di~ensional Newton - R.aphson scheme is employed ire ogler to determine the circulations Emus The derivatives in equation (27) are determined numer- ically in a similar way as described in t18~. The present iterative scheme has been found to con- verge rapidly (in 3 to 4 iterations). This would require solving, for each time step, the system of equations (25) as many times as the number of iterations. As a result the already considerable computing time 2 would be increased approximately by a factor equal to the num- ber of iterations, since a substantial part of the total computing time is devoted in inverting the system of equations 22. Instead of inverting the system of equa- tions 25 at each time step, the following technique of base problems is developed and implemented: 2For a discretization with N = 40 and llI = 20 the computing time on an IRIS 4D/25 TO is approximately equal to four hours. 671

First, the observation is made that for each time step the matrix f~RIIS] iS identical in equations (29) and (253. Thus, substracting equations 9.5 from 92 we get fA] [my ]( )-t44] =-ALTO [~.~*~(~)-God] (~>8) By defining t[~] = thy*] (k) _ t¢] t[~] = tf *](k) _ tri (29) equation (9S) can then be written as [A] tb44 =-tT] t{~] . (30) At this point, we define the base potentials, [q>]m, which correspond to the system of equations (30), and which are the solutions to the base problems tA]i4~]m = _iT]iBim ~ m = 1,.1~1, (31) where tB]m = IBM = 0, B2 = 0, ~ Bm = 1, ..., BM = COT. (32) Physically, the base potentials Brim correspond to the potentials on the propeller blade, when there is no inflow and the potential jumps are equal to zero in all the wake panels except the first panel in the wale at blade strip m, in Chicle there is a linear dipole distri- bution with potential equal to unity at the left (trailing edge) end, and equal to zero at the right end. The base problem is also shown schematically in Figure 5. Notice that the base solutions depend only on the pro- peller discretization and that they are independent of the propeller inflow and the time step n. Due to the linear character of the system of equa- tions (30), the solution, Ah], can then be expressed as a linear superposition of the base potentials Al [it] = ~ firm [] (33) m=1 and by using the definitions (29), the solution to the system of equations (25) is expressed as: Al [~*](k) = [I] + ~ (Em( ) -EnI) [I?] · (34) m=1 In conclusion, by employing the previously described technique of the base problems, we avoid solving the sys- tem of equations (25) for each Kutta, condition iteration and at each time step. Instead, we express the solution to equations (25) in terms of the base potentials. The base potentials are determined before the unsteady so- lution process starts, by solving the so stem of equations (31) as many times as the number of the spanwise blade sections M, which is far less than the number of times the system of equations (22) is solved. As a result the Figure 5: Definition of the base problem at blade strip m. increase in the total computing time, when the present iterative pressure I(utta condition is applied, is insignif- icant. 2.5 Numerical Validation of the Unsteady Pane] Methoc] In this section? the convergence of the results from ~lif- ferent applications of the present. unsteady flow panel method is investigated by wearying the main discretiza- tion parameters. In addition, the consistency of the results from the presented panel method versus ana- lytical or numerical results from lifting surface theory is investigated. Due to the facts tl-~a.t potential based panel methods cannot treat zero thickness wings (equa- tion (3) degenerates to an identity for zero thickness) and that lifting surface theory is exact only for zero thickness, a direct comparison of the presented method with lifting surface theory is not possible. Instead, the following consistency test has been proposed by Kinnas t204. The unsteady panel method is applied on a series of wings (or propellers)' with identical mean camber surfaces but with thickness distributions scaled from an original distribution by a (non-zero) uniform factor. The resulting spanwise circulation distributions are ex- tra.polated to zero thicl;ness and Glen compared with the results from applying an unsteady lifting surface method to the zero thickness wing (or propeller) from the above series. When applying the consistency test for several wings and propellers in steady flow it has been found that the effect of thiel~ness on the spanwise circulation distribution is almost linear with thickness for a large range of thicknesses [~)0], [13], and thus lin- ear extrapolation with thickness of the results from the panel method is employed. In the next two sections the convergence and consis- tency of the unsteady panel method is investigated first when applied to two dimensional hydrofoils and then to propellers. 672

1.00 ~ c=1 ~ tt~8 pow O. 60 ~O.ZO /~=T r 1. · ~ Unit -0.20 Figure 6: Two dimensional symmetric foil subject to a -o 60 sinusoidal gust -1.00 ~ 2.5.1 Sinusoidal Gust in Two Dimensions The two dimensional version of the presented method is applied on several symmetric modified NACA66 forms t3] with chord c = 1 and various maximum thicknesses Tmaa: The foils are subject to a sinusoidal gust normal to their chord, of amplitude By = 0.2 and of frequency, w. The gust is carried downstream lay a uniform flow lr = 1 and reaches the leading edge of the foil at time t = 0, as shown in Figure 6. The vorticity shed from the foil is positioned along the trailing edge bisector line. The reduced frequency of the gust is defined as usual as 2U (35) First, the sensitivity of the method to the size of the time step, i\t, is investigated, when a, constant strength dipole distribution is utilized in the first wake panel. The results, shown in Figure 7, appear to be very de pendent on the ratio of the length of the first wake panel UL\t to the length of the trailing edge panel ~XT. De pendency on the time step is a very undesirable charac teristic of any unsteady panel method, especially if the method is expected to be allele to handle low as well as high frequency components of the incoming flow. However, when the first panel in the walls is treated with a linear rather than constant dipole distribution, as described in section 2.9, the results, shown in Figure 8, become independent of the time step. Replacing the rest of the constant strength panels in the wake with linear was found to have no significant effect on the results. The poor performance of the constant wake dipole panel is attributed to the fact that the poten tials induced on the trailing edge panels by the saw tooth dipole distribution on the first wake panel, shown shaded in Figure 6, are not negligible. Those saw-tooth induced potentials become larger, the larger the ratio of the first wake panel to the trailing edge panel, and the larger the slope of the dipole distribution at the trailing edge (i.e. frequency) are. Due to the previous investigation, a. linear dipole distribution is always employed in the first wale panel in the unsteady panel method for hydrofoils as well as propellers. - ~ = 2.0 ~ ~u,zr = 8.u ~ of = 4.0 ~ 'I ' = 16.0 - T 8. 0 2. 0 16. 0 . 0 Ut C Figure 7: Convergence with time step size. Constant dipole distribution on the first wake panel. Modified NACA66 form with loYo thickness to chord ratio. Si- nusoidal transverse gust with amplitude vg = 0.2 and reduced frequency k = 0.5. Number of panels on the foil N = 40. 1'00L ' ' ' r 0.20 -O. 20 -O. 60 -1.00 - - ~ = 2.0 Ago Uz t = 8 0\(; Act - . gut` = 4,0 ~ ~ = 16.0 ~ ; 0 4. 0 B. 0 12. 0 16. 0 2t . 0 Ut Figure 8: Convergence with time step size. Linear dipole distribution on the first wake panel. Modified NACA66 form with 10% thickness to chord ratio. Si- nusoidal transverse gust with amplitude vg = 0.2 and reduced frequency k = 0.5. Number.of panel on the foil N = 40. 1.00 0.60 r 0.20 -0. 20 -0. 60 -1.00 j\ ~ ~ _ - 200 panels ~ 60 panels ~ . 120 panels + 40 panels l l l l l l l l l to C .0 4.0 8.0 12.0 16.0 2C .0 Ut c Figure 9: Convergence with number of panels on the hydrofoil. Lou' frequency case. Modified NACA66 form with logo thickness to chord ratio. Sinusoidal transverse gust with amplitude vg = 0.2 and reduced frequency k = 0.5. 673

0.12 0.08 0.04 r BOO -O. 04 -O. Oa -O. 121t ~ I ~ ~ ~ ~ ~ ~ 1 n 1:, . - 200 panels ~ 60 panels . 120 panels + 40 panels . 0 18. 4 18. 8 19. 2 19. 6 2( . o -0. 12 Ut c Figure 10: Convergence with number of panels on the hydrofoil. High frequency case. Modified NACA 66 form with 10% thickness to chord ratio. Sinusoidal transverse gust with amplitude vg = 0.2 and reduced frequency 1c = 5. The convergence with number of panels on the foil, N. is then shown in Figures 9 and 10 for a low and a high reduced frequency, respectively. The convergence appears to be very good. Finally, in order to compare the results from the panel method with the analytic results for a flat plate in a sinusoidal gust - the well known Sears problem t29] - the consistency test, as described in the beginning of this section, is applied. The unsteady panel method is applied to three foils with man = 0.05, 0.1 and 0.2. The results (in this case the time history of the circu- lation around the foil), appear to be linear in thickness and thus, are linearly extrapolated to zero thic.l;ness and compared against the steady state analytic circu- lation of a flat plate in a gust. The comparisons are shown in Figures 11 and 12 for a low and a high re- duced frequency. In both those cases, the consistency test appears to be valid, within acceptable accuracy 3. The analytic complex Circulation, r, with respect to 1.00 F 0.60 0.20 r -0. 20 -0.60 _ -1.00 ~ .0 -my A ~ : - Analytical - Panel Method (Extrapolated) . . . . 4.0 8.0 12.0 16.0 20.0 Ut Figure 11: Consistency test of unsteady panel method with analytic result - Low frequency. Sinusoidal trans- verse gust with amplitude vg = 0.2 and reduced fre- quency k = 0.5. 30nly the steady state numerical circulation should be com- pared to the analytical. . . 0.08 0.04 ~ 0.00 -0. 04 -o. 08 'if - Analytical - Panel Method (Extrapolated) _ .0 18.4 18.8 ~.2 ~.6 2 Ut C Figure 12: Consistency test of unsteady panel method with analytic result - Hiy1' frequency. Sinusoidal trans- verse gust with amplitude vg = 0.2 and reduced fre- quency k = 5. the leading edge, for the flat plate in a sinusoidal gust is r e 2i Jo(~)-iJl(~) (36) 2,rcv9 irk Ho(k)-iHl(~) where Jo, J1 and Ho, Hi are the Bessel and Hankel func- tions, respectively, and i = Am. 2.5.2 The Propeller in Nonuniform Inflow The unsteady panel method is applied to the DTRC N4118 propeller, of which the geometry is given in All. The incoming wake is assumed to be axial, with a once per revolution circumferential variation: Use = Up(1 + ng COS dS) (37) with Up being the circumferentially averaged inflow and ug the amplitude of the wake variation, taken equal to 0.2. The advance coefficient, defined as Jp = 2,rUp/~wD'), with D being the diameter of the propeller, is taken equal to 0.833. 7.00 6.00 5 00 4.00 rayon 27rRUp 3.00 2.00 lOO t o Go , . . - ^8W = 3-0° OR = 0. 60 0~8w=6.0° /~ ~ E] ^8W = 12.:r/1l = 0 Ul:~ S`, of' Wf r/R = 0. 30 I\ i. 0 60. 0 120. 0 180. 0 240. 0 300. 0 wt Figure 13: Convergence of the unsteady propeller panel method with time step size. Circulation vs. blade an- gle at different propeller radii, r/R, for the N4118 pro- peller. Jo = 0.833, ug = 0.2. 674

7 Do - - - 6.00 _ 5.00 _ r*loo 4-00 21rRUp 3.00 2.00 1 nrl 0.00 ).0 r/R = 0.,:~ - 40c.30s panels 40c.25s panels 40c.20s panels 40c~15s panels 60. 0 t - . 0 leO. 0 240. 0 ~0. 0 36 ). 0 but .... . . Figure 14: Convergence of the unsteady propeller panel method with chordwise number of panels. Circulation vs. blade angle at different propeller radii, r/R, for the N4118 propeller. Jp = 0.833, ug = 0.2. . . ~ nr, . BOO r*loo 4 00 DRUB 3.00 2.00 1.00 l r/R = 0.6/,~ - 60c"20s panels ~ 50c.20s panels e' 40c.20s panels + 30ca20s panels O. 00 ). O 60. 0 120. 0 18O. 0 240. 0 300. 0 36 ). 0 cut Figure 15: Convergence of the unsteady propeller panel method with spanwise number of panels. Circulation vs. blade angle at different propeller radii, r/R, for the N4118 propeller. Jp = 0.833, ug = 0.2. 1.50 r*loo Wrap o.so - .50 I. ~ -2. SC l r/R= 0.80/// /~r/R = 0. 30 I- Lifting Surface - Panel Method (Extrapolated) ). O 60. 0 1aD. 0 18O. 0 240. 0 at. 0 ~ ,. It Figure 16: Consistency test of the unsteady propeller panel method with the unsteady lifting surface method. Circulation vs. blade angle at different propeller radii, r/R, for the N4118 propeller. Jp = 0.S33, ug = 0.2. The convergence of the results for this case is inves- tigated by varying the time step size, the chordwise and the spanwise number of panels. The results are shown in Figures 13, 14 and 15, respectively. The convergence with all those parameters is shown to be very good. Finally, the consistency test for Else unsteady panel method versus the unsteady lifting surface method [11], is applied for the Ar4118 propeller. The linearly extrap- olated results from the panel method and those from the lifting surface method are in very good agreement, as shown in Figure 16. costs, rs, Bs) ~ ED ~ ~_~: Figure 17: Ducted propeller in spatially nonuniform in- flow. 3 The Unsteady Panel Method for Ducted Propellers 3.1 Formulation Consider now, a ducted propeller, shown in Figure 17, subject to the spatially nonuniform inflow Us(`xs, rs, bs) with respect to the ship fixed system as described in Section 2.1. The duct is defined~by its chord length, ID, the thickness and camber distributions of its merid- ional section, the angle of attack, °rD, and, the axial distance, ED, of its leading edge from the midpoint of the propeller at the hub section. The formulation of the unsteady ducted propeller flow problem is identical to that in the case of open propellers, which was described in Section 2.1. The perturbation potential, A, is again determined bar us- ing Green's formula, equation (3~. Both the propeller and duct surfaces must be panelled and unsteady wal;es must be shed from the trailing edges of the propeller blades and duct. In the present work, we model the duct with a potential based panel method and the pro- peller with a lifting surface vortex lattice method which will be described in Section 4. A typical arrangement of panels on the duct, the propeller and their trailing wakes is shown in Figure 18. 675

em Figure 18: Panel arrangement of the duct, propeller and their trailing wakes in unsteady flow. Only half of the duct panels are shown. One way of analyzing the unsteady flow around the ducted propeller is by treating the duct and propeller as one body and by employing a time marching scheme similar to that described in Section 2.2. Instead, in the present work, we only treat the propeller with a time marching lifting surface scheme with the effects of the duct on the unsteady propeller flowfield being a,c- counted for via the generalized images of the propeller singularities with respect to the duct. Some of the ad- vantages of the latter Versus the former scheme will be outlined later in this section. It is well l;nown that the potential flow around a, body in the presence of an infinite `~-all, can be modelied with singularities distributed on the surface of the body, as if the flow domain were unbouncled, with the effects of the wall being accounted for via the images of the singularities with respect to the wall. The generalized image idea, as introduced by Kinnas and Coney [21], is an extension of the classical image idea for the case that the infinite wall is a generally shaped body, for example a duct. The flow field of the generalized image of a sin- gularity (i.e. source, vortex or dipole) `vith respect to a generally shaped body A, is defined as the modification to the flow field of the singularity due to the presence of body A. In mathematical terns, as described in the next section, the combined flow field due to the singu- larity and its generalized image with respect to a body A, consists the modified Green's function which satis- fies the boundary conditions on A. A similar idea has been applied in analyzing the flow around a body in the presence of a free surface, where the modified Green's function is required to satisfy the free surface condition [224. The implementation of the generalized image idea in the analysis of the unstea(ly flow around ducted pro ~_ ~.~ q (a) q (b) ~N ~ qi (c) Figure 19: Decomposition of the flow field around a duct and the propeller vortex loops, subject to a, spa- tially nonuniform inflow. pelters consists of the following steps, which are also shown schematically in 'Figure 19: Step 1 Solve for the flow around the duct (without the propeller) subject to the nonuniform inflow Us~xs~rs,bs), by using a potential based panel method t184. The resulting flow field is in general non-axisymmetric, but steady with respect to the ship fixed system. The total velocity, qD, is com- puted at the propeller control points. Step 2a For each unit strength vortex loop on the pro- peller key blade and its trailing wale, the gen- eralized image flow field is computed. This is achieved by solving for the non-axisymmetric po- tential on the duct in the presence of that loop and in zero inflow, bar apply ing a potential based panel method [914. The same panel arrangement on the duct is used as in Step 1, thus avoiding re- computing the duct to duct influence coefficients 676

for each vortex loop. Step 2b The velocities, qi, induced by tile generalized image of each of the unit strength vortex loops at the propeller control points are computed and stored. Step 3 The unsteady propeller lifting surface method, described in Section 4, is then applied as if the propeller was open, but with the following modi- fications in order to take into account the duct: · The propeller inflow is taken equal to qD, as computed in Step 1. The propeller to propeller influence coeffi- cients, velocities A, are modified by adding to them the corresponding generalized im- a.ge influence coefficients, qi, which were com- puted in Step 21~. Steps 2a and 2b consume the largest and Step 3 the smallest part of the total computing time. If the same ducted propeller must be analyzed for different inflows, Steps l and 3 must be repeated, but Steps 2a and 2b, which depend on the duct and propeller geometry but not on the inflow, must not, as long as the generalized image influence coefficients have been stored. Thus, the same configuration can be analyzed for different inflow conditions, rather quickly, after the first condition has been analyzed. The matrices that have to be inverted for each vor- tex loop in Step 2a are identical, and the corresponding right-hand sides are independent from each other. This makes Steps 2a and 2b amenable to parallel processing. This would not be feasible, if the duct and propeller were treated with a tinge marching scheme, since the right-hand sides at each time step are not independent from each other. U . nA A (ha, =~ _ By \ \ / ~ \ \ // \~/ 6 G WB Figure 20: Potential flow around bodies A and B. The flow field due to the modified source G = G + Gi is shown with dashed lines. 3.2 The Generalized Image Mode] A mathematical justification l~elli~:ld the generalized im- age model is given in this section. bite consider two bod- ies A and B subject to a genera.! inflow Uin as shorten in Figure 20. Assume that the floral around A and B is lifting and that the geometry of their trailing wakes WA and WB, respectively, is known. The flow can be analyzed by considering bodies A and B as one and proceed in the usual way by apply- ing Green's formula, equation (3), with respect to the perturbation potential ~ on A or B: 27r~ = | ~ dS + | ¢) dS SA On B SO 071. A + | (~¢,)A {3G dS~ | (~¢,)B dS WA art WD on, - | G(-US nA)dS-~ G(-[Nn t7B)dS (38) where G is the infinite fluid domain Green's function, defined by equation (4), SA and SIB are the surfaces of bodies A and B. nA and nB are the normals on the bodies A and B. respectively, nit and nWB are the normals on the wake surfaces LISA and I4{B, respectively, (~)A and (~\~)B are the jumps in the potential across the wake sheets 1174 and T1rg, respectively. An alternative way of formulating Green's formula on B. in the presence of A, can be found by introducing the Green's function G which satisfies V2G(p; q) = -4~(xp - xq)~(yp - yq)~(zp - zq) outside A (39) with ~ being the delta function and, (xp,yp,zp) and (Xq~yq~Zq) being the Cartesian coordinates of points p and q, respectively. Also ~ FIG = a, on A, (40) nA G(p;q)~O as p loo, (41) VG= finite at tic trailing eclge of A. (42) ~ Physically, G(p; q) corresponds to the potential due to a point source of strength-4~, placed at q, in the presence of body A, as shown in Figure 20. Notice that ~ G must satisfy the I`utta condition (42) at the trailing edge of A. It can easily be shown that G satisfies 677

G( ) | G( ) [JG d5 I (~G)A DIG d5 + 41,G(P; q) WA IOWA (43) for p ~ A. G is the infinite domain Green's function. The jump of the potential (~\G)A in the wake of A is set equal to the difference of the potentials at the trailing edge of A in order for the I`utta condition (42) to be satisfied. ~ The function G will be given in the fluid domain as ~ G(P; q) = G(P; q) + + -[| G-dS + J (I) dS] (44) 41r SA Ulna H/A IOWA for points p outside body A. By defining i(P;q) 47E [,ISA Lana t ~ G can then be expressed as 6,,| (~G)A dS] IVA On WA (45) ~ G(P; q) = G(P; q) ~ Gi(P; q) (46) Gi is the modification to the infinite fluid domain Green's function due to the presence of body A. We call Gi the generalized image of the source G with respect to body A. In the case the body A is an infinite wall then Gi is identical to the actual image of C; `srith respect to the wall. It can be shown that the pertu~ba.tion potential, ¢, in the presence of both A and B can then be expressed as ~ = (A + LAB (47) where HA iS the perturbation potential outside A in the presence of the flow field Uin but in the absence of body ~ B. and SIB iS the perturbation potential on B in the presence of the flow field Uin + V¢4. The potentials jA and SIB satisfy the following equations: 270A = ,/ ¢)A-dS + .| (COMA) dS- SA 071 A IV 4 (fin II{A - | Gt-Uin 77.A)dS oil A; (4S) 2~)B = | OB-dS + | (~(~)B)-dS- SB SUB WB ~nWB - | G(\-[Jill 7IB- (9 HIS on B. (49~) The physical meaning behind equation (47) is that the flow in the presence of bodies A and B can be de- composed into two parts: · The flow a.rouncl body A in the presence of the incoming floor [din · The flow around body B in the presence of the modified flow field resulting from the previous step. ~ The potential RIB can be expressed as the superpo sition of distributions of modified sources G and modified dipoles DING on body B and its wake. The modified dipole Rig can easily be shown, by using equations (40) and (42), that satisfies the kine- matic boundary condition and the I`utta condition on A. Physically, the modified dipole consists of the infi- nite fluid domain dipole and its generalized image with respect to body A. It is obvious that the role of bodies A and B in the decomposition, equation (47), is interchangeable. The special case of A being the dock and B being the propeller, has been utilized in analyzing a ducted propeller in a spatially nonuniform flow, as described in Section 3.1. 4 Unsteady [Lifting Surface Theory for Ducked Propellers The panel method for ducts described earlier has been combined with a lifting surface representation of the en- closed propeller to permit time-domain solution for un- steady blade forces. We describe here the formulation of the propeller boundary value problem, the method by which forces are obtained and review the conver- gence behavior of the model. Some findings based on the application of the code are given in the next section. 4.1 Formulation Keenan [14] has shown that the vortex lattice lifting surface representation of the propeller blades described here follows naturally from the potential-based method described in Section 2. Differentiating equation (3) to obtain the perturbation velocity at the field point x and then passing to the limit of vanishing blade (and wake sheet) thickness one finds is ,u(~)n~xj tV,~(n(~) Vex _' mud -Is ~w(~)n(X) tv=(n(~) V`: AX _ ':~ )] do -2~rn(x) vl(x). (50) This is an expression for the velocity normal to the blade at point x due to a distribution of dipoles of un- t;nown strength ill on the blade surfaces Sp and known low on the wake surfaces Sw. The dipole strength in the wake, ,`~, is known from the history of blade cir- culation. We require that 678

v1 n = (U-v2) n - on Sp, (51) vat ~ 0 at infinity (52) where vat is the perturbation due to the propeller and wake sheets, v2 is the background velocity and U is the velocity of the body itself. The total fluid velocity is V = vet + V2 There are two integrals in equation (50); one, Ip, over Sp and one, Ill,, over Sit,. The discretization is the same for both so we consider only Ip, namely Ip = ~ P(~)n(x) tV=(n(~) ~ \7~ `,~)1~da. (53) The surface Sp is divided into quadrilateral panels whose vertices are placed on the camber surface of the blade as described later. Within each panel, the dipole strength ~ is assumed constant. Thus Ip is approxi- mated as Ip=~/li/s n(x) tV=(n(~) <~x-'id)] The surface Sp is replaced by the collection of panels Si. Two problems arise at this point. The required integrals are hyper-singular and must be performed on arbitrarily curved surfaces Si. By recalling the equiv- alence between a constant strength dipole patch and a vortex loop around its perimeter (see, e.g., t94~) one can escape these difficulties. One then can write the discrete form of Ip as Ip =~,llin~x)~{ ':~ xdli. The integral on Si is replaced with a simple application of Biot-Savart's law for the straight vortex segments bounding each panel. The arrangement of the panelling on the blades is shown in Figure 21. The spacing follows the so-called quasi-continuous method (QCM) described by Lan t23] and greatly enhances the method's ability to capture the square root like behavior near the blade edges as compared to classical vortex lattice methods (VLM). This arrangement of panels is often referred to as "c.o- sine spacing". QChI has the added benefit over VLLI of placing panels closer to the erlges of the blade, thins increasing resolution. For a. ducted propeller, the spa.n- wise panelling is modified. The typical clucted propeller will not exhibit squa.re-root-lil;e behavior of the span- wise loading near the tip (unless there is a large tip gap) so QCM spacing is inappropriate in the spanwise direction. Figure 21 shows the uniform spacing used in the spanwise direction for duct-fitted propellers. The tangency boundary condition on the blades is satisfied by collocation at control points given by QCM. The geometry of the wale sheets is prescribed. The angular extent of the panels is equal to that swept out by the blade in one time step except that the wake panel adjoining the trailing edge is one quarter this size. Figure 21: Typical blade panellings. On the left, fol- lowing Lan's Quasi-Continuous Method. On the right, using uniform spacing spanwise as for a ducted pro- peller. On the right, sections near the root have been trimmed to fit a tapered hub. The discretization of equation (50) is described more fully in [14~. Eventually, it leads to the linear system K M N ~ ~ ~ ai,n,m,krin,m,k k=! m=~1 n=! K Al Nw ~ ~ ~ I,, - 2- 2~ 2-ai,n,m,kitT] n,m,k k=~1 m=! n=2 -ni v~(xi). (~56) Each ai,n,m,k is the Hence at control lit i due to the (n., m)th unit strength dipole panel/`ortex loop plus it's exact image in the duct. Having computed these matrix elements, it is left to solve for the Es's ``rhich are simply (55) the jumps in dipole strength going from one panel to its neighbor along a chordwise strip. With h Al N.1, be = - ~ ~ ~ al n m i. lawn m k - ni v~(xi). (57) k=! m=} n=2 one can recast equation (56) into the form ~ (All ) (Alla ) ~ · . 1\ (Awl) .. banjo;) ~ ~ fir,, ~ ~ (b,,) ~ (58) Each block in the matrix corresponds to the influence of one blade on another. Suppose that the matrix is separated so that only the diagonal blocks appear on the left side. Then equation (5S) breal;s down into the following series of I; sub-equations: K (Akk) (rk) = (bk)-At, (ski) (ri ), ~ = 1, 2' , It t=1 (59) The form of equation (59) hints at a substantial compu- tational savings. If the [S's appearing on the right side were known the matrix solution could be reduced from an (N x M x It' )2 effort to It x (N x AI )2. Fortunately, for a normal propeller, the blade-to-blade influence rep 679

resented by the summation on the right is small. This circumstance permits use of a block Ga,uss-Seidel itera- tive solution, thus: k-~ (`Akk~) (l~k`)( A) = (be; - ~ (A~i) (rS)(n+~) i=! - ~ (Audi) (rs)(n), ~ = 1,~,...,It'. (60) i=k+1 The superscripts in parentheses indicate the iteration level. The iteration continues until the desired tolerance for Us is achieved. The speed-up over direct solution of the full matrix equation is, in fact, even greater than suggested above. Each (Aij) is constant in time (by virtue of the place- ment of the 1st shed vortex, cf.~144) so an LU decompo- sition may be performed once at the beginning of the calculation. Only back-substitution is required at each time step thereafter. For normal propellers, a further advantage accrues from the fact that each blade is the same shape. That means that each (A',,k) is the same only one block needs to be decomposed. The resulting scheme is very fast. Computation of the right-hand-side vector b tales far longer than does extracting Us once b is known. 4.2 Forces on the propeller Once the discrete vortex strengths ale l;nown, the blade forces may be calculated. Previous implementations of vortex lattice lifting surface models have obtained blade forces from the "rotating bedspring" analogy (a term coined by J. E. Kerwin) 5,11,17. In that ap- proach, one imagines the lifting surface to be replaced by the array of singularities. The forces computed are those arising from the application of JouLowski's rule to each vortex segment using velocities evaluated at the vortex segment (usually on its midpoint). Guermond [6] has correctly criticized this approach. Its fault lies in the fact that velocities calculated on the body any- where away from the control points do not satisfy the kinematic boundary condition (51~. Kerwin [15] has observed that while the leading edge suction acts prin- cipally at the first vortex, its effect is distributed over the entire chord. Thus, moments will lee incorrect even if total forces are accurate. Two alternative methods of calculating force distri- butions may be suggested as being correct. One may in- terpolate the velocities, correctly computed at the con- trol points, to the vortex elements or one may evalu- ate vorticity at the control points. The latter approach is taken here. Also, rather than apply the Joul;owski rule, local pressures will be evaluated. This is seen as a convenience since one often requires surface pressures anyway for additional studies. 4.2.1 Blade surface pressures To obtain surface pressures we employ Bernoulli's equa- tion in the form P-Pa = 2p(VcO VOO -V V -2 .9 ) (61) where Pa is the ambient pressure far enough upstream so that do/0t may be assumed small, VOO is the total fluid velocity at the same location and V is the total fluid velocity at the point of interest. This differs from the normal Bernoulli equation derived for unsteady, purely potential flow in that the velocities involved are not merely To but include the contribution from the rota- tional background flow. Equation (61) also differs from Bernoulli's equation for steady rotational flow because of the do/0t term. A more expansive discussion of this point can be found in [144. Having solved equation (50), the velocities at the control points on the blade are immediately available as ~ M N V (Xi ~ = V2- 27r ~ ~ [I V~,n,m,krn,m,k Nw Jr 2_ v~ n make wn,m.~ (62) n=1 where V2,n,m,k, V',n,m,k are the vecto' il~fltlences at xi of the unit strength vortex loops. The velocity on the left is written as VM to emphasize that it is the mean velocity at the control point. There is also a local jump in velocity caused by the local vortex density Rim k which is 2~Vn,m,k = An m,k X nn,m,~ = Vspiln~m,k (63) so that Vi = V3l ~ i\v on the suction/pressure sides. Equation (61) also includes the term b<~/0t,which must be evaluated for unsteady flows. Unsteady in this sense refers to the ship-fixed, inertial coordinates so, in this model, propellers operating in axisymmetric inflows are unsteady. This term can be written as ,~ = rQut + ,~ | (64) The quantity at is the tangential velocity induced by the singularities on the propeller and wale sheets and Q is the shaft angular velocity. The second term on the right is the potential change at the blade-fixed point of interest. It turns out that this last term can be related sim- ply to the circulation around that portion of the blade forward of the control pointt14~. In discretized form, it is: fit ~ ~ ,'~j ~ (653 the sum being from the leading edge to the 77th panel along the chord. Thus, the pre,ssllre jump across the blade Sp is p+ -p rip (V TV+ b! (O O )) assuming pi+ = pOO anal that the vorticity in the bacl< 680

ground is weak. The pressure jumps are combined with the corresponding directed panel area. to obtain an in- cremental force OF which acts at the control point. 4. 2.2 Leading Edge Suction As a consequence of idealizing the propeller blades as zero thickness surfaces we introduce a singularity in pressure at the leading edge. This gives rise to the well known leading edge suction force. Derivations of this force may be found in many classical references. One such is by Milne-Thomson t304. The leading edge suc- tion force Fs is where Fs = - 47rpCs i (67) Cs = 1i~ ~ gyps ~ 1 (68) and l is a unit vector directed towarcls the tip along the leading edge, i is directed strean~wise normal to the leading edge and lying in the surface So. The quantity s is the arc length from the leading edge measured along a curve on Sp in the i direction; fly is the vortex density on Sp. Lan showed in his analysis of the QCM that the value of Cs could be obtained directly by computing the downwash at the leading edge. In this worl;, a direct ap- plication of the limit process expressed in equation (6S) is employed instead. To evaluate the limit, the vorticity fly is calculated at the control points along a stream`~ise strip. At each control point, we set An = An 1. Now, apart from a factor related to the leading edge sweep, we have a. sequence in n, Cs" = \f~'Yn, (69) which must lee evaluated for n ~ 0. This requires extrapolation since the quantities called for are only known inboard of the leading edge. The extrapolation can be made reliable if it is assumed that fly may be expressed as gypsy = ~;P(s) (70) where P(s) is a polynomial to any degree in s. If this is the case, the slope of the curve /y vanishes at the leading edge. The assumption is not very restrictive. It is true for a flat plate and for parabolic camber. It strictly fails for NACA a-series loadings but these are mathematical idealizations, requiring a logarithmic sin- gularity in camber at the leading edge. The loading actually achieved is probably closer to something ex- pressible by (70~. This approach is computationally ef- fective in that it obtains Cs directly front already known quantities: no additional influence coefficients need be calculated. 4.3 Results After the duct and propeller pa.nelling is specified the duct problem is solved for each panel on the propeller and wake sheets in turn. This gives the duct contri- bution to the propeller panel Green's functions. The remainder of the code is virtually the same as for an open propeller except that the influence coefficients in the matrix include the image effects. The time-domain solver then proceeds just as fast for the ducted problem as for an open propeller. meanline: a = 0.8 thickness form: NACA66 r/R 0.182 0.300 o.ioo o.soo 0.600 0.700 0.800 0.900 .950 .000 P/D rake skew c/D [/c . 1.423 o.ooo o.ooo 0.179 0.018 _ 1.402 o.ooo o.ooO 0.208 0.020 __ .389 o.ooo o.ooO 0.232 _ 0.022 .380 o.oOo o.oOo 0.254 0.023 .379 O.ooo o.ooo 0.273 0.023 .386 0.000 ~ o.OOO 0.288 0.023 .408 0.OoO O.oOO-0.299 0.020 .446 0.000 o.OOO 0.306 0.014 .472 ~O.ooo O.oOo 0.308 o.oos .502 ~0.000 0.000 0.309 0.00s Test Duct ID = R ED = 0° t/lD = 0.10 NACA66 Table 1 t/D 0.039 0.031 0.023 0.016 0.015 0.0ls 0.0ls 0.015 . 0.015 0.015 As long as the panelling of the duct, propeller or wake is not changed, one may run the time-domain solver repeatedly for various inflow configurations with- out returning to the duct presolver. Thus one may amor- tize the presolver overhead according to the number of inflow cases one considers. VVe have tested the code for a simple propeller/duct combination described in Table 1. In l~igure 22 we show the spanwise distribution of circulation for uniform flow at a nominal advance coefficient of J = 1.0. Three curves are presented for various panellings of the duct and propeller combination. The 6 x 5 curve has six pan- els over the span of the propeller, five over the chord. The duct has 24 panels circumferentia.lly and 40 around the chord. The curves labelled 10 x 5 and 10 x 10 have 3.00 . ~1~ so so O.WO L :~" r ~- ~/ "' .~. 6 x 5 ducted ·--10 x 5 ducted -10 x 10 ducted ·~10 x 10 open 1.00 Figure 22: Convergence of span`` ise distribution of cir- culation with propeller and duct panelling. 681

l --- r/R=O. 6368 6 x 5 case as.. r/R=0.6959 10 x 5 case _r/R=0.6959 10 x 10 case 00 of ~ 4. 0(, 4.00 , Kazoo ~ W~~/~ ~I5z-oo~ ~=:J O O """t,, ""'''''"'''- '""""" ,, 360. 0 oo~ooo 180. 0 360. 0 180. 0 360. 0 n° n°° Figure 23: Convergence of unsteady circulation for the ducted propeller. 40 panels circumferentially on the duct. For compa.r- ison, we include a fourth curve which is for the same propeller but without a duct. All have twelve stream- wise along the propeller Wolfe corresponding to 30 time steps per revolution. Predictions of steady I' agree well with corresponding results from the code based on the work of Kerwin el. al. t18~. To study the unsteady behavior we ran the same three cases in the artificial inflow given by Vx(~) = 1-c ~ ~ y ~-~ + ~-b = 0 Vt = 0 where c = as(1+~) b = as-asi/2 (1 +£)a,rctan=-~/2. (~72) ~Ox1Op~pdbr~d _ oxso~m ~J'1.O 1 2 3 4 6 Hamlonk d Watt rats did _ open Figure 24: Unsteady torque computed directly for the ducted propeller and inferred from an open propeller case. , . -10 x 10 ducted, r/R=0.9609 -10 x 10 Acted, r/R=O. 6959 ---10 x 10 open. r/R=O. 6959 · ~ ~ 10 x 10 own, r/R=O. 9609 Figure 25: Unsteady circulation for ducted and open propeller. The quantities a and ~ permit adjustment of Else depth of the wake defect and its width, respectively. Here they are set to a = 0.5 and ~-0.05. The result is a wake with a 50% defect at zero degrees and has a volumetric mean velocity of unity. Figure 23 shows unsteady circulation near 0.7R as a function of shaft rotation angle for the three ducted cases presented alcove. Notice convergence appears to be very rapid. The lifting surface code used in this model generally gives very good convergence behavior for circulation both spatially and temporally. We find here that the ductecl case is nearly as well behaved once the duct is sufficiently resolved circumferentially. Finally, to examine the eRect tl~e duct has on the un- stea.dy forces we look at Figure 91 which compares the (71) 10 x 10 ducted propeller case with a. run using the same propeller without the duct. For the unducted propeller we include the augment to the steady flow which the duct produces; only the panel image effects are turned off. Given only an open propeller code one might ima.g- ine taking this latter approach to get unsteady forces on a ducted propeller. As the figure shows, one would underpredict unsteady torque by a, significant fraction. Thrust behaves similarly. Figure 25 compares circula- tion for these same runs at two radii. Near 0.7R, the curves for the open and ducted cases are, apart from the expected mean offset, essentially the same. How- ever, the curves near 0.96R show sigllific.antly different behavior. Thus, the sectional forces can be expected to show even greater divergence in their spectra than do the global forces shown alcove. 5 Conclusions A time marching potential based panel method was pre- sented for the analysis of the unsteady flow around open marine propellers subject to spatially nonuniform in- flows. An efficient algorithm is implemented in order to ensure an explicit Rutty condition (i.e. pressure equa,l- ity) at the blade trailing edge at each time step. The 682

numerics of the method are shown to be very robust for a broad range of reduced frequencies. The method is also shown to be consistent with known analytic solu- tions as well as with an existing unsteady lifting surface method. The method provides the user with accurate unsteady pressure distributions which may be coupled with a viscous / inviscid interaction scheme to predict unsteady boundary layer separation and / or leading edge separation. A hybrid panel method is also developed for the analysis of the unsteady flow around ducted propellers. The method combines an unsteady lifting surface method for the propeller with a potential based panel method for the duct. The propeller is treated with a time march- ing vortex lattice scheme as if the propeller were open. The effects of the duct on the propeller are included via the modified inflow which is due to the presence of the duct (in the absence of the propeller) and via gen- eralized images for the propeller and its trailing wake. The proposed method is shown to be appropriate for treating unsteady ducted propeller flows in a computa- tionally efficient and robust manner, especially when a given ducted propeller geometry must lee analyzed for various inflow conditions. 6 Acknowledgments Support for this research was provided by the MIT Sea Grant College Program and the David Taylor Re- search Center, Department of the Navy, Grant Num- ber NA86AA-D-SG089, and by the Office of Naval Re- search Contracts N00014-S9-J-3194 and N00014-87-K- 0-192. The authors would like to thank Professor Justin E. Kerwin of the Department of Ocean Engineering at MIT for his valuable comments and suggestions during the course of this work, and also Lisa Shields, a former graduate student in the Department of Ocean Engineer- ing, for her help during the initial development of the unsteady ducted propeller computer codes. References [1] B.C. Basu and G.J. Hancocl;. The unsteady motion of a two-dimensional aerofoil in incompressible in- viscid flow. Journal of Fluid Mechanics, Vol. 87:pp 159-178, 1978. t2] R.J. Boswell and M.L. Miller. Unsteady Pro- peller Loading - Aleasurement, Correlation, with Theory and Paramet'ic Steady. Technical Re- port DTNSRDC 2695, DTNSRDC, October 1968. t3] T. Brocl;ett. Minimum Pressure Envelopes for AIodiJied NA CA-~6 Sections with ETA CA a=0.8 Camber and Buships Type I and Type II Sections. Report 1780, DTNSRDC, Teddington, England, Feb 1966. Joseph P. Giesing. Nonlinear two-dimensional un- steady potential flow with lift. J. Aircraft, 5~2), Mar-Apr. 1968. D.S. Greeley and J.E. Kerwin. Numerical methods for propeller design and analysis in steady flow. Trans. SNrAAlE, vol 90, 1982. t6] J. L. Guermond. About collocation methods for marine propeller design. In Propellers '88, 1988. i7] T. Hanaol;a. Hydrodynamics of an oscillating screw propeller. In Proc. 4th Symp. Slav. Hydrodyn., Natl. Acad. Press, 1962. tS] T. Hanaoka. Numerical Lifting Surface Theory of a Screw Propeller in Non-uniform Flow (Part 1: Fundamental Theory). Technical Report 6~5~:1-14, Ship Res. Inst., Tokyo, 1969. t9] J. L. Hess. Calculation of Potentialf7OW About Ar- bitrary Three-Dimen.sional Liftiny Bodies. Tech- nical Report MDC J5679-01, McDonnell Douglas, October 1972. t10] J. L. Hess and Exalter O. Valarezo. Calculation of steady flow about propellers by means of a surface panel method. In 23rd Aerospace Sciences Fleeting, AIAA, Reno, Nevada, January 1985. t11] T. Hoshino. Application of quasi - continuous method to unsteady propeller lifting surface prob- lems. J. Soc. Nav. Arch. Japan, INS, 1985. T. Hoshino. Hydrodynamic analysis of propellers in steady flow using a surface panel method. In Proceedings of the Spring Meeting, The Society of Naval Architects of Japan, May 19S9. t13] Ching-Yeh Hsin. Development and Analysis of Un- steady Propeller Panel A~Iethod. PhD thesis, De- partment of Ocean Engineering,~MIT, 1990. Under preparation. t14] D. P. Keenan. AIarirle Propellers in TJn.steady pow. PhD thesis, Massachusetts Institute of Technology, May 1989. , J. E:. Kerwin. Private communication. t16; J. E. Kerwin. Marine propellers. Annual Review of Fluid Mechanics, 18:367-403, 1986. 683 J. E. Kerwin and C. S. Lee. Prediction of steady and unsteady marine propeller performance by nu- merical lifting-surface theory. SHAME Transac- tions, 86, 1978. J.E. Kerwin, S.A. Kinnas, J-T Lee, and W-Z Shih. A surface panel method for the hydrodynamic analysis of ducted propellers. T'ans. SNA ME, 95, 1987.

[19] I`erwin.J.E. and C-S Lee. Prediction of steady and unsteady marine propeller performance by numer- ical lifting-surface theory. Trans. $ PIE, vol S6, 1978. t20] S.A. I`innas. A General Theory for the Coupling Between Thickness and Loading for NVings and Propeller. Submitted for Publication, June 1990. [21] Spyros A. Kinnas and William B. Coney. On the optimum ducted propeller loading. In Proceedings of the Propellers '88 Symposium, SNAME, Virginia Beach, VA, September l9S8. t22] F.T. Korsmeyer, C.H. Lee, J.N. Newman, and P.D. Scalvounos. The analysis of wave effects on tention- leg platforms. In OAIAE '88 Conference, 1988. t23] C. E. Lan. A quasi-vortex-lattice method in thin wing theory. Journal of A ircr`/ft, 11~9), September 1974. t24] J. T. Lee. A potential based panel method for the analysis of marine propellers in steady flow. Tech- nical Report 87-13, Dept. of Ocean Engineering, Massachusetts Institute of Tecl~nology, July 19S7. (25] Brian Masl~ew. Influence of rotor blade tip shape on tip vortex shedding- an unsteady inviscid anal- ysis. In Proceedings of 36th Annual AHS Forum, 1980. t26] Luigi Morino, Zaven Jr. Kaprielian, and Slobo- dan R. Sipcic. Free lVake Aerodynamic Analy- sis of Helicopter Rotors. Technical Report CN. DAAG29-S0-C-001C, U.S. Army Research Office, May 1983. [27] Luigi Morino and Ching-Cl~iang Kuo. Subsonic potential aerodynamic for complex configurations : a general theory. AIAA Journal, vol 12(no 2~:pp 191-197, February 1974. [2S] J.N. Newman. Distributions of sources and normal dipoles over a quacl~ilateral panel. Journal of En- gineering Alalhematics, vol 90:pp 113-126, l9S6. t29] William R. Sears. Some aspects of non-stationary airfoil theory and its practical application. Journal of the Aeronautical Sciences, vol. S(No. 2), 1941. t30] L. M. Milne - Thomson. Theoretical Aerodynamics. Dover, fourth edition, 1966. [31] S. Tsakonas, W.R. Jacobs, and M.R. All. An exact linear lifting-surface theory for a marine propeller in a nonuniform flow field. Journal of Ship Re- search, vol. 17(No. 4~:pp 196-207, Dec. 1973. t32] S. Tsakonas, NV.R. Jacobs, and P.H. Rank. Un- steady propeller lifting-surface theory with finite number of chordwise modes. Journal of Ship Re- search, vol. 12(No. l):pp 14-45, 196S. 684

- DISCUSSION Jin7hang Feng Pennsylvania State University, USA (China) 1. Regarding the implementation of Kutta C.D., an analytical expression of Eq. 27 can easily be formulated. In fact, such an expression has been successfully used in a 3-D panel code developed at Penn. State. Why the authors chose to use a numerical approach to determine the Jacobin, which is inevitably more expensive, at least in terms of CPU time? 2. Theoretically, when using the generalized image model, the immense numerical effort needed to evaluate the modified Green function ~ may well surpass~any likely advantage in solving the potential at a later stage after G is known. More important I think, any simplified model derived as a result of employ such scan be equivalently established with the original Green Function. Would the authors explain their consideration of using such a procedure? 3. A non-uniform inflow for propeller is usually rotational. The inflow vortex will redistribute when propeller disturbance is introduced. Further, wake vortex shed from the blade trailing edge might deform, including rolling up as one often sees in a similar two dimensional lifting flow. Why is it important to use an unsteady approach to tackle the problem if the vortex transportation feature cannot be accounted for properly? Would it not be easier, for example, just to use a quasi-steady approach to attempt the problem? AUTHORS' REPLY On the implementation of the iterative pressure Kutta condition: to determine the elements of the Jacobin matrix (Eq. 27), we evaluate the derivatives involved numerically by perturbing the circulation distribution at each strip by a small amount and then computing its effect on the pressure difference at the trailing edge of all strips (the details are described in [18]). This operation involves only differentiations of potentials to find pressures, and thus requires very minimal CPU time, especially as compared to the CPU time required to solve for the propeller potentials at each time step. The discusser mentions an alternative method to determine the elements of the Jacobin matrix, but he does not give a reference and thus we cannot assess the validity of his approach for our problem. On the use of the generalized image model: as described in Section 3.1 of the paper, the generalized image of each unit strength vortex loop on the propeller or on its wake corresponds to the solution on the duct in the presence of that vortex loop and in zero inflow. Thus, the generalized image coefficients depend only on the geometry of the duct and propeller and do not need to be recomputed for each time step in the propeller solution or for a different inflow. The primary reasons for selecting the generalized image technique as opposed to a direct duct and propeller time marching solver are mentioned in Section 3.1. They are as follows: (1) the same geometry can be run with considerably less CPU time for different inflows and, (2) the computations of the generalized images of each propeller vortex loop are independent of each other and can be performed very quickly on a parallel processing computer. On his suggestion to use a quasi-steady technique instead of our fully unsteady method: it is well established in the hydrodynamic community that the quasi-steady theory fails in predicting unsteady forces on propellers, especially at blade rate (for example look at Fig. 15 of [16]). We prefer to be systematic in accounting for each of the effects (e.g., potential flow effect, effective wake effect). At present, we have developed a computationally reliable method to account for the potential flow effect. We plan to couple our method with a Euler solver in determining the effective wake contribution resulting from the vorticity transport Mr. Feng mentions. DISCUSSION Ali H. Nayfeh Virginia Polytechnic Institute and State University, USA It would be very useful to compare your results with the detailed experimental results of D. Telionis and co-workers and the panel results D.T. Mook and co-works for the case of sinusoidal gust in two dimensions. D.T. Mook, A.H. Nayfeh, and co-workers developed steady and unsteady discrete and continuous vortex methods for lifting surface and rotor blades at high angles of attack, including interference effects. They used a Kutta condition that does not require iteration. The results were published in the AIAA Journal of Aircraft. AUTHORS' REPLY On comparing the results of our method in two dimensions against existing experiments: this is something that we plan to do in the near future, though not before we include the tunnel wall effects as well as the unsteady boundary layer effects. The primary goal of the research reported in the paper was to produce a boundary element method (BEM) for the unsteady potential flow around hydrofoils or propellers, which method would be accurate, robust, and consistent for a broad range of reduced frequencies. We will perform an experiment at the MIT water tunnel on a two dimensional hydrofoil subject to a sinusoidal gust in reduced frequencies (k= 10) which are closer to propeller applications than those used in previous experiments. We also plan to compare our results to this experiment. On the application of an iterative pressure Kutta condition: it has been found to be necessary when applying velocity based BEM formulations (Hess[9]) and more recently when applying potential based BEM formulations (Kerwin, Kinnas, Lee & Shih [18]). However, when a lifting surface (zero thickness) vortex lattice technique, such as that of Mook & Nayfeh, is applied, then the Kutta condition is equivalent to requiring the bound vorticity at the trailing edge to be equal to zero (which can be imposed either explicitly or implicitly) and thus an iterative condition is not needed (for example look in Kerwin & Lee[18]). 685