Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 667
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
OCR for page 668
OCR for page 670
OCR for page 672
OCR for page 673
OCR for page 674
OCR for page 675
OCR for page 676
OCR for page 677
OCR for page 678
OCR for page 679
OCR for page 680
OCR for page 681
OCR for page 682
OCR for page 683
OCR for page 684
OCR for page 685
OCR for page 686
Representative terms from entire chapter:
lifting surface
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(~, (4)
with R(
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], [
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
