An improved panel method for unsteady aerofoil flows

AN IMPROVED PANEL METHOD

FOR UNSTEADY AEROFOIL FLOWS

N. Vlachos1 and S.P. Fiddes2 ,

Department of Aerospace Engineering, University of Bristol Bristol, UK

Abstract

This paper presents some improvements to a low-order panel method for computing unsteady, inviscid, in-compressible potential flows. An explicit Kutta con-dition for equalization of the trailing edge pressures is defined and implemented iteratively in a linearised form yielding improved prediction of pressures at the trailing edge. The wake shedding model, which pos-sesses an arbitrary parameter for this type of low-order method, is also examined and a method for removing the indeterminacy of the wake-shedding parameter on a rational basis is put forward. These improvements are incorporated into a two-dimensional panel method and results obtained for a NACA0012 section in tran-sient and oscillatory motions. Results are compared with those from a moving-mesh finite-volume Euler solver for incompressible flow, and the two methods are shown to agree closely for the lift and pitching mo-ment predictions, while comparison of the computed wakes shows major differences. Additional results are also shown for the panel method for the case of the convection of a point vortex past a moving aerofoil and its wake.

Nomenclature Symbols

c aerofoil chord (solution characteristic length)

section lift coefficient

section pitching moment coefficient (about

"jc

= 0.25)

N, Nw number of panels around the aerofoil sec-tion and number of wake panels

npc solution time-steps per cycle of motion

p pressure

q surface-tangent perturbation velocity r0 intrinsic position vector

(x',

y') frorn center

of rotation

r vector from singularity to point of influence S, i1 surface-panel tangent and normal vectors

1 _{graduate student} 2_Profcssor

Sn, Sw body and wake surfaces

u, v flow velocity components in (::v,

y)

V0 inertial velocity of intrinsic frame origin

x 1 y inertial frame coordinates

t1 r time and non-dimensional time r :::: tU j c

a aerodynamic incidence

r

circulation

1 vorticity

!:'J.CPTE pressure coefficient difference between

up-per and lower surfaces at the Os panel chord-length trailing edge (w wake shedding parameter

p surface doublicity strength (equivalent to ifJ - ¢;)

¢ perturbation velocity potential cr surface source density

( eqmva ent to . 1 7);; -

oq,, oq,

Bn

l

no

rate of rotation of intrinsic frame

w angular frequency Superscri ptsjSu bscripts

intrinsic (body-fixed) frame quantity previous time level quantity

quantity norrnaliscd against steady-state value

quantities internal to Sn or body panel index

u,D upper and lower surface at the trailing edge n,w body and wake quantities respectively ss steady-state

Introd11etion

In many panel rnethods1 the I<utt.a condition of classical <:1.crofoil theory has been applied in a man-ner that. implicitly accounts for the physical effects of viscosity at points of high surface curva.t.ure, simulta-neously fixing the level of bound circulation of lifting surfaces. In its analytic application for steady flow about hvo-dimcnsional acrofoils, the Kutta condition

specifies that the -flow separates from the aerofoil sur· face at the aerofoil trailing edge such that the flow speeds on both surfaces of the trailing edge remain finite[l]. From this follow restrictions on the flow ve-locity at the edge different for finite-angle and cusped trailing edge geometries

[6],

both cases having in com-mon equal pressures on upper and lower surfaces at the trailing edge.

This equalization of pressure on either side of a trailing edge is also used to define the I<utta condition in unsteady aerofoil problems where the conditions on the velocity at the trailing edge are then modified by the sense of vorticity being shed from the edge[4,

6].

It has been suggested that there is an upper limit to the frequency range for which this interpretation of the I<utta condition can be applied [10, 11, 12, 13, 14], but such extreme conditions will not be considered in this work.

A simple \implicit~ I\.utta condition1 is used in some well known low-order panel methods for steady and unsteady flows [11 71 81 151. This paper will present an implementation of the \equal-pressure' I<utta con-dition and treatment of the aerofoil trailing edge in the framework of a low-order panel method. This im-proved form of the Kutta condition is needed fol' un-steady flow to improve the prediction of surface pres-sure especially in the vicinity of trailing edge. The goal of this work is an implementation of this condi-tion that docs not require higher order singularities or geometry representation on any part of the solid or wake surfaces modelled and must not incur heavy computational cost.s. Furthermore, it is neccs:sary to tackle the problem of the free wakc>sheclding parame-ter which is a feature of this type of low-order panel method [1, 5].

Low-order panel methods arc very versatik1 and developrnents of the approach described here can ul-timately he applied to the solution of unsteady in-compressible potential flows about multiple thrce-dirncnsional deforming geometries. For sirnulations of such complexity, simplicity and computational ef-ficiency arc paramount and are the prirnary incentives for pursuing improvcrnent.s for tl!is panel method.

Solution Method

Doublet and source singularities are used to repre-sent the surface of two-dimensional profiles in incorn-pressihle potential flow where the governing cq11at.ion (for continuity) is Laplace's equation

( 1) and the potential ¢ is related to the now velocity hy

'V¢

=

(u, v) (2)

Since Eq. (1) is a linear differential equation the principle of superposition applies and the complete so-lution can be assembled from a combination of elemen-tary solutions. The elemenelemen-tary solutions used here in-clude the two-dimensional vortex, doublet and source singularities which satisfy Laplace's equation in there-gion of the flow around them[1]. Uniform distributions of sources and doublets on flat panels are used for the body, while the wake is represented by doublet panels only (equivalent to point vortices at the doublet panel junctions.)

The total potential at a point in the flow is the summation the global ¢ influence of all the (J and p

singularities distributed over Sn and Sw

¢(x, y, z)

=

_!__

r

J1:

(~)

ds 47r}sa.w n r

- _!__

r

(J

(~)

ds (3)

47r}so r

The surfaces Sn and Sw arc modelled by N and Nw discrct.e panels respectively, so the discrete form of Eq. (3) is needed

N N Nw

¢;(x,y)

=

I>JA;,J

+

L(JJB;,J

+

LPkC;,k (4)

j::::l j::::l k=l

The influence coefficients Ai,j,Bi,J and Ck₁J which re-place the integrals in Eq. {3) for constant unit strength doublet and source panels are given by

A · - -1

)1[}(1)

- - ds '1 - 4rr ₅

,8n

ri,j ( 5) B; ·

= _

_!__

1-

1-ds ,.J 47r _s, 7'· . l,J ( 6) C; k

=

_!__ {

.!!._(-

1- ) ds ' 1rr Js~<Dn l'i,k (7) where the integrals arc taken over the each of the body and wake panels denoted by

si

and

sk.

'The solution aims to find the distributions of u and fl that produce a flow-field which satisfies the physically correct boundary condition, i.e. zero flow through Sn. In the numerical solution, this is strictly true only at the discrete points where the boundary condition is enforced) in this cctsc points just inside Su at panel centres (the collocation points.)

The boundary condition is tackled in two steps in this method. Firstly, to nullify the local onset flow clue to the kincrnat.ic velocity norrnal to Sn at each collocat.ion point., rr set so t.hat

(To= -(V,

+

n

X r,) ·ii (8)

at each panel. This only satisfies the boundary con-dition

locally

and the (J distribution induces a

pertur-bation potential to the global flow which does not yet

satisfy the tangency condition on Sn. Therefore, in the second step, the 11 distribution must be adjusted to nullify the global perturbation due to ff. Since the above also means that

8¢,

=

0

8n

(9)

and there are no singularities internal to

Sn

the Dirich-let boundary condition can be enforced[!]

r/!;

=

0 (10)

where r/!; is given by Eq. (4). This results in the fol-lowing matrbc equation for finding the needed combi-nation of surface doublet strengths (including a wake doublet panel) where JlN flW [ lUISI

l

RJIS,

=

RIISN N W-1 (11)

lUIS;=

'f:;u;B;,;

+

L

JtkC'i,k (12)

j=l k=l

The above matrix equation has one more unknown than equations and requires an additional equation for a unique solution provided by applying a trailing edge condition.

Unsteady solution

The unsteady solution proceeds by time-stepping through a series of instantaneous solutions of Eq. {11) and a trailing edge condition. The solution is aware of the flow history through the strengths and positions of the wake singularities.

Changes in the level of aerofoil bound circulation result in the shedding of vorticity at the trailing edge, the strength of which is determined with the applica-tion of the Kutta condiapplica-tion. The is wake shed at points where the Kutta condition is enforced, one new free vortex being created for each solution time-step rep-resenting the vorticity shed during that tirnc interval. The wake points arc, in the time-steps subsequent t.o their creation, convected with the local flow velocity, the new positions being calculated by a simple Eu\e:r integration scheme which assumes a constant veloc-ity history of each point during each tirne step. This scheme is known to exhibit cumulative error (worsened

with large time-step size) and instability with small time-step size[16]. However, this scheme is considered adequate for the current computations. To limit the instabilities arising from close vortex spacing, a solid-body rotation vortex core of radius typically 10-5 _of the aerofoil chord is applied in the calculation of ve-locity influence between free vortices.

The condition of zero net circulation production in the solution (the Kelvin condition) is automatically satisfied by this method since each wake doublet panel created is equivalent to two equal and opposite point vortices at the panel ends which together contribute zero net circulation.

Unsteady pressure

The form of the Bernoulli equation which relates the flow velocity and acceleration (in the inertial frame) to the pressure, in the absence of any force fields (e.g. ncglect.ing gravity), is (from Ref. [1])

1' (V¢)' 8,P 1'=

+ + =

-p 2 8t p (13)

This defines the pressure 111 terms of the local flow conditions and the conditions at infinity where the flow perturbation is zero.

Since pressure on the surface of the body is of

in-terest, the analysis is restricted to the external velocity potential and its derivatives on Sn. The r/!(x,y) field used in Eq. {13) is an instantaneous sampling of the perturbation field induced by the singularities bound to and moving with S'n and Sw. In addition to the

8¢/8l

due to the fluctuation in singularity strengths, <;\ point fixed in (a:, y) w\H, in general, also perceive a change in ¢ with time clue to the movement of the sing11laritics. To obtain t.he total

8¢/8l

for a point fixed in inertial space the following transformation is

needed

8¢1

8t J:,y

iJr/!1

D¢1

Osl

Dt

s,n

+

EJs n,t

Dl

.t·,y

+ -

Oq\1

-

8nl

8n s ,t Dt x ,!J ( 14)

The potential carried by the body on the outside of Su with cloublicity ft is sirnply

q\

=

r/!;

+

Jl

=

Jl

frorn Eq. (10), which means lhat

f)£ I

()l 3,1!

8p

8t ( 15)

As S'n moves through inertia! space the c:ornpo-nent of flow velocity in the surface normal direction is

y'

x'

Sw

Figure 1: Inertial and aerofoi1 intrinsic coordinates

x'

Figure 2: Detail of surface quantities

known directly from the local velocity of the surface since the surface will carry the fluid in the ii direction

8

¢

=

(V,

+

D,

x r,) ·

t1

=

-<r

iJn

( 16)

In the § direction the fluid moves relative to the (x, y) \vith the local tangential perturbation velocity (irrespective of (V,

+

0, X r,) ·

s)

so in the

s

direction we have

. 8¢ 8jl. \ l ¢ · s = - = - = q

as

as

( 17)

Combining the expressions of Eq. (16-17) and Eq. (15) yields the total inertial 8¢/81

8¢

I

8n

I

as

I

8''1

- = -<r -

+

q -

+ -,

Dt x,y Dl ;1:,y 8l :c,y 8l s,n

( 18)

where (s, n) arc surface coordinates local to and fixed to each point on Sn and known in terms of (x', y'). Given the motion state of the body (Vo,Oo) through

(x,yL

it can be shown that

-

as

I

=-

(V o

+

" Ho X ro · S . ) .

8t

x,y

( 1 9)

~~~

=

-(V,+D,

x r,)

·n

="

x,y (20) The full expression relating the local surface con-ditions with far-fteld concon-ditions now becomes (from Eq. 13)

Poo-p_!(' _{- - q -(;} ')+&p

-p 2 8t

- q(V,

+

D,

x r,) ·

s

(21) To obtain an expression for Cp from the above, it is necessary to choose a velocity which defines a ref-erence stagnation pressure. For the calculation of the force/moment coefftcients to be possible via integra-tion of the Cp around the section the reference velocity must be constant over Sn. The most usual choice is the translational velocity of the whole body through (x,

y)

although other reference speeds must be chosen for cases of pure rotation.

Eq. (21) in terms of the reference flow speed now yields the expression for the pressure coefficient on the body surface C'r

=

+("'

-q2

+

2q(V,

+

12, xr,)-8- 2:) V:·cj where C' _{p -} P-1 v2 Poo 2P ref

Unsteady Kutta Condition

(22)

(23)

The commonly used implicit. form of the Kutta con-clition, implemented here for comparison, specifies for steady flow [1]

lu - ''/,.

= ()

In the case of this low order panel method, this condi-tion is replaced with

(see Fig. 3) and is implerncntcd by

(24) In conjunction with the Kutt.a-Joukowski theorem this equation irnplies that the loading at the trailing edge is zero. Hence, the pressures at the trailing edge should he equal and the implicit condition of Eq. (24) should affect the nearby fl solution such that this re-sults. The solution has a direct effect on the trailing

Figure 3: Detail of trailing edge panels

edge pressures through

OJ.l/&s.

This is approximated by first order finite differencing of the step-wise contin-uous

J.l·

However, it has been found that Eq. (24) does not provide sufficient control over the distribution of the discrete values of fJ. to result in equal trailing edge pressures as calculated in the post-processing.

In unsteady flow, the inadequacy of Eq. (24) is worsened by the omission of the 8¢ j &t component of the unsteady pressure. The condition which mnst now be applied, derived from Eq. (22) for /!,.Cp"' = 0 at the trailing edge, is

q' -

0''

+

2

OJ.ll -

q'

+ "'' -

2

OJ.ll

U U OtU L L O t £

- 2qu (V, + 0 X

r,) ·

su

+

2qL

(V, + 0

X

r,) ·

ii,,

=

0 (25) This can not be implemented directly in the solu-tion of the linear boundary condisolu-tion equasolu-tions since there are unknown quadratic terms of Jli in the veloc-ity q2• In order to avoid the computational overhead of solving with a non-linear equation an approxirnation is needed which willlinearise Eq. {25). The technique employed by this panel method, introduced by Bose[2], is to replace the quadratic perturbation velocity terms with products of the current (unknown) velocity q and the perturbation velocity at the previous tirne level ij,

Le. to assume that

(

()'')'""

&s

(&'')

&

1

_,1

_{_}

- - -qq

&s

8s

t-t>t - (26)

The linearised trailing edge condition implemented ls now

qu (iju -2(V, +0 x ro) ·su)

- qL (ij,,- 2(V,

+

0 x r,). iii-) - 0'2

+

0'2

+

2 OJI

I -

2 OJI

I

= 0

U L ()t U Q/ [_ (27)

The unknown q terms arc replaced by their

nurncr-ical analogues (second order)

_ ( 2 2Js, )

qu-- Ss, +Ss, + Ss,(Ss, +Ss,) J.l,

-t ( 2

+

2Ss,

+

2Js, ) ·

- Ss, + Ss,

Ss,(os, + Ss,)

os,(os, + os,) J.l,

( 2<ls, )

- Jl (28)

Ss,(Ss, +Ss,)

'

(29) The application of the explicit trailing edge con-dition involves only unknO\vn quantities (in this case, pi) on the aerofoil surface in an equation of the form

AN-2PN-2

+

AN-tflN-1

+

AN-2PN

- A,p,- A,p,- A,/1,

=

1U!S'N+! (30)

Note that the strength of the shed wake panel Jlw is not directly coupled to the solution through Eq. (27) in this formulat.ion. Tnstead is determined indirectly by its influence on the aerofoil bounda1·y condition, i.e. through t.he coeflicient.s C;,w in Eq. (11). This makes it possible to usc any type of singularity to model the wake without aiTccting the implementation of the trail-ing edge condit.ion, although it might then be neces-sary to explicitly enforce the I<elvin condition.

The approxirnation in Eq. (2G} relics on the as-sumption that the perturbation velocities at the trail-ing edge do not change rapidly with time, and so Eq. (27) is less accurate with increasing rates of change and tirne-st.ep size. In such cases, sub-iterations arc performed within the same titnc level, the new values of q bci11g used to update the approximation (j, un-til the pressure difference is below some set tolerance, typically 0.005 of the reference dynamic head. The iterative process involves solving the matrix equation Eq. (11) changed only by the current set of coefficients from Eq. (~W), and can be stopped once !::l.CP'n-: is be-low some preset magnitude.

It has been fonnd that this schcrne works for most unsteady motions where t.he arnplitudes and frequen-cies of the motions arc not too great. The scheme is not convergent for all situations. At points in the solu-tion where there arc rapid accelerasolu-tions at the trailing edge point, more iterations arc required to attain the required difference in pressure, and in some cases the solution locks before the required pressure difference

is reached. It is sometimes useful to reset the coef-ficients of the iterative process with a sub-iteration using Eq. (24)), although after this the solution of-ten locks in the same way. Doubling the precision of the computations gives only a slight improvement. Other more successful techniques for enhancing the sub-iteration process are being sought. One candidate which has not been sufficiently explored is adaptive time-stepping. Reducing the time step size will always make the approximation of Eq. (26) more accurate and convergence may be improved if the time-step size is reduced in proportion to rates of change of flow quan-tities at the trailing edge.

Wake Shedding Parameter

The current method does not provide the velocity exactly at the trailing edge or just off the trailing edge which would be needed to determine the position of each new shed wake vortex. Velocity calculation on or infinitesimally close to the trailing edge is as accurate

as the interpolation/extrapolation scheme used to in-terpret the singularities in the vicinity. To bypass the lack of solution detail the average of the upper and ]ower surface trailing edge velocities is taken as the most representative of the time-average of the velocity with which the newest wake point is convected during the time-step it is generated. The wake doublet panel created docs not constitute part of a streamline. It is assumed that the new vortex has followed a streak line to its current position (as computed by above approx-imation) but that the direction of the panel does not dictate the exact direction of flow off the trailing edge. Trying to prescribe or solve for such detail is inconsis-tent with the order of the current method

+v

ot

new wake point

previous wake po1nt -~

Figure 4: Detail of wake shedding at aerofoil trailing edge

The lack of detail in the wake shedding process introduci'S the pararncter (w to the method which is defined as the fraction of the local velocity-convected distance from the trailing edge at which the new wake point is placed. The same type of indeterminacy was encountered for the the PM ARC code by Martin et al.[5]. In the latter case, the optimum value of (w

was sought by comparison with analytic results (Wag-ner and Theodorscn functions) for a thin flat plate in transient and oscillatory motions. Given that the lat-ter are themselves approximations) this type of com-parison is not entirely satisfactory as a calibration of

(w and another route is taken here.

More continuous representation of the shed vor-ticity may remove the (w indeterminacy and it was

attempted to replace the doublet panel at the trail-ing edge by a flat linear -y panel in order to verify the optimum (w. This introduced an additional unknown d-yjds which it was hoped would allow the solution to detcrm:ne automatically the center of vorticity within the new wake panel. The result was an oscillatory feed-back instability and further work in this direction was abandoned. Another candidate approach was based on the work of Graham[17] using the complex-plane solution for the flow off an infinite wedge corner as a local solution at the trailing edge to determine the strength and positioning of the shed vorticity. This too was not taken to completion due to difficulty in inter-facing the local solution with the global panel method solution and due to the prospect of poor extensibility to three dimensions. In general) it is not desirable to use higher order singularity distribution or geometry in the trailing edge region in order to preserve the sim-plicity of the method and to allow direct. extension to a three-dimensional vortex lattice/panel method.

Fortunately) it. becomes apparent that the method is) to a good extent) self-calibrating from the following observation: as !:ir is decreased the range of position tha.t the new vortex point. can take and the strength of the vortex both diminish. Therefore) as the time-step is decreased (w becomes less significant. This seems consistent within the fnm1ework of the time discrcti-sation since as !:ir -> 0 a continuous vorticity sheet would result in making (w practically ineffective in

the limit (the shed vortex point. would always coin-cide with the trailing edge). Obviously) due to nu-merical considerations (rnachi11e precision) vortex wake stability) rnemory) CPU tirne) thi::; limit can not be ap-proached to the extent that would fully calibrate the method, i.e. to the extent \vhcrc (w variation would have no effect at all on the results. However) calibra-tion results can be obtained with small enough !:ir to give the clcsirccl rnaximum cleviat.ion in the results for a given change in (w so that larger tirne-step runs can be used for longer ~irnulations.

To show that. (w bccornes less significant with de-creasing D.r a series of runs were performed with the NACA0012 aerofoilundcrgoing one cycle of a= 5° sin{2kr) pitch oscillation in steady forward flight. This was done for

Cv

values of 0.30 to 0.60 in equal increments of 0.1 ()) repeating tl1is series of calculations four times, initially with C,r"' 0.097 and halving C,r each time (sec Fig.

G).

The difference between the

solutions between successive values of (w at the differ-ent /!,.r clearly diminishes with decreasing !!,.r show-ing that the effect of (w becomes less significant. The solution

eM

is used for the sensitivity analysis since it is most strongly affected by changes in (w as seen in Fig. 5. In this case, which exhibits a maximum

eM

"'! 0.030, reducing /!,.r from 0.097 to approxi-mately 0.012 reduces the maximum change in solution

for /!,.(w = 0.1 from 10% to approximately 1% of

max-imum

eM.

Assuming that the method converges to the correct solution with decreasing /!,.r and that the discrete vortex shed must lie somewhere in the rniddlc of the curve of continuous vorticity it is intended to represent, using any (w in the range between O.a and 0.7 will yield a solution within a few percent of the datum solution, even for

eM

which is sensitive to (w

variation.

The decreasing effect of (w can only be exploited if for decreasing ~r the method converges to one solu-tion. The same set of runs used for the above analysis are used to show the convergence characteristics of the method. The

eM

results of the smallest /!,.r solution is used as a datum and the deviation of the other solu-tions frorn this is calculated within each (w series. Dif-ferent convergence characteristics are seen at difDif-ferent values of (w· In all cases, as expected, the smaller .6.r runs show smaller deviation from the dat.nm, although not with uniformity within the cycle (see Fig. 7). The inability of the iterative Kutta condition to converge at the ends of the cycle is believed to be the reason for this non-uniformity and, on the whole, the method tends to a single solution as /!,.r

->

0.

To find the optimum (w for use with larger 6-r for each type of motion modelled a reference ;:;olut.ion should be computed with /!,.r ;;) 0.01 and (, "' 0.5. Then, results at the desired 6-r can be computed for a range of (w values so that the optimum can be selected by comparison with the reference run which is known, from the above previous discussion, to produce results sufHciently independent of (w. This calibration proce-dure would also be applicable to the three-dimensional extension of the current method. Alternatively, the calibration may be bypassed by obtaining the solu-tions with very fmc time-step size and (w ~ 0.5.

Discussion of Results Performance of the unsteady I<:utta condition

Figure 9 compares the Cp error at. the {.railing edge with that of the linearised explicit I<utta condition of Eq. (27). The

eP

error is plotted against the strength of the wake vortex shed during the time-step that the 6-CPTr:: was exhibited. This essentially plots the er-ror in terms of a measure of the degree to which the flow is unsteady since the strength of the wake vor-tex being shed is inversely proportional to the rate of

change of bound circulation (by the Kelvin theorem). Therefore, it is apparent that Eq. (24) produces error increasing in proportion to the degree to which the flow is unsteady and that the explicit Kutta condition of Eq. (27) performs significantly better.

Impulsive start at incidence

The transient response of the NACAOOI2 section started impulsively from rest into constant. speed for-ward flight at a

=

5° was computed. The results for some values of (w arc shown in Fig. 8. It is appar-ent that (w has a stronger effect on the transient (;M

response. The C~., response approaches Wagner's an-alytic indicia! response for a flat plate [1] which pre-dicts a starting (:L of 0.50. Perfect agreement is not expected since the panel method includes the effects of the time discrctisat.ion and of aerofoil thickness. Oscillatory pitch motion

The NACA0012 section was used for a set of com-putations spanning a range of k from 0.02 to 0.30 shown in Fig. 5. The upper limit is considered low enough t.o avoid the regions of k where the validity of the unst.eady Kutta condition is questionable. The NAC:A0012 section exhibits the collapse and reversal of the

e/,

loop at. approximately k

=

0.2 as predicted by Geisslcr[9] among others. The time~step is kept low in order to minimise (w effects with npc

=

400 used in the cornput.ations. Since .6.r is related to the number of time-steps used by 6-r

=

1rjknpc, it not constant for t.he runs at. different. frequencies. This is most no-ticeable fork

=

0.02 from the magnified effect of (w in the C M response.

Wake convection sensitivity

In addition t.o recalculating the location of the wake as the acrofoil moves (solid-body translation and rota-tion), keeping the wake completely force-free requires the calculat.ion of velocity influence ( u and v compo-nents) at each of the current wake points from all other wake vortices and acrofoil singularities. This results in distortion of the wake and amounts to ap-proximately 2Nw (Nii

+

Nw) velocity influence calcu-lations per t.irne-st.ep. Each solution tirne-st.ep a new wake-point is created so Nw is increrncntecl and the computation is slowed nearly in proportion to

Na,.

Therefore it. is desirable to reduce the portion of the wake for wl!ich t.hc distorting components of flow ve-locity are applied, say for t.hc Kma:.r,· vortices nearest

the trailing edge. Using the case of the NACA0012 section in ±5° sinusoid pitching motion it. is evident frorn Fig. 15 that. the solution is relatively insensitive to the calculation of t.he dist.ort.ion of t.he far wake. H.cdncing !Cnrcr results in significantly different wake

shapes (Fig. 16) but the effect on the solution is not noticed until I<max/npc 1'0

1%.

Multiple body prob-lems and simulations where the nature of the motion keeps the wake in proximity of the body may not be eligible for this computational saving.

Euler method comparison

The steady-flow aspects of the method can be veri-fied by comparison with analytic results. However 1 for the unsteady a..<5pects1 such as the wake treatment, it is useful to compare with another method which solves for similar flow but which does not deal ex-plicitly with the creation and convection of the shed vorticity. The method of Gaitonde[3] is used for com-parison, as it solves the Euler equations (inviscid, in-compressible but rotational flow) so fluid rotation is implicitly included. Unsteady solutions are obtained using artificial compressibility and a moving mesh on which a cell-centered finite volume scheme is imple-mented.

Initially, the resulting vortex wake frorn the two methods is compared to assess their ability to capture the ·flow history. The vorticity is extracted from the so-lutions ofGaitonde[3] for the k = 0.30 pitching motion case. This is done only at the solution instants corre-sponding to the quarter-phase positions of the pitching motion cycle for the NACA0012 section at the high-est frequency compared (k = 0.3). The panel method wake points arc superimposed on the resulting iso-1 contours in Fig. 14 with a graphical indication of the sign and strength of the panel method wake circula-tion.

The comparison indicates that the Euler solver is less capable of capturing and retaining the vorticity generated by the aerofoil1 especially beyond five aero-foil chord lengths downstream, Decreasing grid den-sity and increasing cell aspect ratio further from the acrofoil make it impossible for the Euler solver to re-solve the wake as well as the panel method and it ex-hibits appreciable diffusion of the vorticity. Further-more) since this is a deforming mesh method) the fluid behind the aerofoil is sarnp1ed by a grid which rnoves cyclicly through it as the trailing edge moves normal to the direction of the onset flow. This means that the wake is sampled successively hy grid regions of varying density. toss of How dct<-Li\ because of thls grid varia-tion and the numerical dissipavaria-tion incorporated in the method may be the main contributing factors to the

clifrusion and loss of vorticity.

The difference in wake-capturing capability docs not seem to significantly affect the agreement. for C'L and CAJ seen in Fig. 12 where results are presented for the oscillatory pitch motion cases. The C'L agreement is excellent with both methods predicting the drop in hysteresis and the reversal of the

c/,

loop direction at k

=

0.20. The agreement. in

eM

is good but.

de-teriorates with decreasing frequency. This appears to be because the two methods predict slightly different steady-state solutions for dCM / dcx of the NACA0012 profile used. In order to gauge the agreement in a way which eliminates some of the effects of the post-processing (in this case, pressure integration around the profile), Cp distributions were also compared at the quarter-phase positions of the motion cycle. Fig. 13 shows the almost identical pressure distributions, t.he Euler solver capturing more detail of the Cp variation near the leading edge since over 90 points are used to sample the profile as opposed to the 60 panels used by

the panel method. Vortex interaction studies

The panel method was used to simulate the case of the NACA0012 profile at a= 0' flying past a free vortex. The interaction is simulated for the two cases (shown in Fig. I 0 and Fig. 11) of a vortex of positive and neg-ative circulation of thf, same magnitude starting five chord lengths in front and 0. JOe below relative to the aerofoil centre. The oscillatory characteristics exhib-ited by the CL and eM results during the interactions indicate that the presence of the aerofoil wake damps the influence of the vortex and that modeling the wake

is important in the correct prediction of the pressures during the interaction.

Condnsion

The lincarisecl unsteady Kutta condition described here irnprovcs significantly the pressure prediction at the trailing edge for the simple transient and pitch-oscillation motions without excessive computational cost. However, in1provement of the numerical imple-rncntation is needed in cases where the sub-iteration process fails to converge.

Comparison of panel method results using the new l<utta condition with the results of the cell-centered moving-mesh Euler solver of Gaitoncle(3] shows good agreement between the two very different methods. This establishes the ability of the panel method to compute the unsteady incompressible flow about an acrofoil section with accuracy comparable to a finite volume method solving the Euler equations but. in a

fraction (typically 5%) of the computing time. Fur-thermore, the panel method, by 'fitting, the wake rather than capturing it (as the finite-volume method docs}, is capable of producing superior solutions for the wake shape. This superiority is more apparent at a distance from the trailing edge and is therefore rnore significant in multiple-body problems. It does not appear to he a factor in the current comparison which involves the computation of flow about a single aerofoil.

The wake convection process, which is one of the more computationally demanding aspects of the cur-rent rnethod, is more significant for the wake points in the vicinity of the trailing edge and can be disabled for the wake beyond a certain distance (typically 5% of the motion cycle length) to reduce the computatiomtl effort, without noticeably affecting the solution on the aerofoil.

The effect of variations in the essentially 'free' wake-shedding parameter ((w) which afflicts low-order methods ofthis type, can be controlled simply through reduction of the solution time-step size. It has been shown that decreasing D..r nearly eliminates the ef-fect of (w (allowed to vary within a reasonable range) on the solution. This is true even for the aerofoil CM which was shown to be especially sensitive to (w. Therefore, selection of the optimum (w for larger time-step calculations (for computational efficiency) can be selected by comparison (of a sample of the motion modelled) with a very fine time-step calibration so-lution.

With these improvements in the most critical ar-eas of the Kutta condition and wake shedding treat-ment, the type of low-order singularity method de-scribed here promises to provide useful solutions for the unsteady potential flow about arbitrary three-dimensional geomet.ries with relatively small compu-tational cost. For unsteady flows in particular, tl1c fact that the method only requires quantities on the body and wake surfaces makes it very attractive for more complex motions where surface deformations arc involved.

References

[1] Katz, J. and Plotkin, A., Low Speed

Aerody-namics: From Wing Theory to Panel kfclhods, McGraw-Hill, 1991.

[2] Bose, N., "Explicit Kut.ta Condition for an Unsteady Two-Dimensional Constant Potential Panel Method," AIAA Journal, Vol. a2, No. 5, May 1994, pp. 1078-1080.

[a] Gaitondc, A.L., "A Dual-Tirnc Moving Mesh Method for 2D Unsteady Tncomprcssihlc Flows Past Aerofoi!s," Bristol University, Aerospace En-gineering Dept. Heport No. 735, December 1 99f). [4] Mangler, I<.W. and Smith, .J.II.B., "Behaviour of

the Vortex Sheet at the Trailing Edge of a Lifting \Ving,l) Royal Aero. Soc. Journal, Vol. 71, Novem-ber 1970, pp. 906-908.

[5] Martin, and I< roo, "Comparison of Pi\! ARC and Analytic Ttcsults for Two-Dimensional Unsteady Airfoils," AIAA Paper 93-06:!6, a 1st Aerospace Sciences Meeting and Exhibit, Reno, NV, .Jan. 1993.

51.9

[6] Maskell, E.C., "On the I<utta-Joukowski Condi-tion in Two-Dimensional Unsteady Flow/' Royal Aircraft Establishment Tech. Memo. Aero 1451, Sept. 1972.

[7] Morino, L. and Kuo, C.C., "Subsonic Potential Aerodynamics for Complex Configurations: A General Theory,'' AIAA Journal, vol. 121 no. 2,

pp. 191-197, 1974.

[8] Richason, T.F., Katz, J. and Ashby, D.L., "Un-steady Panel Method for Flows with Multiple Bodies Moving Along Various Paths," AIAA Journal, vol. 32, no. I, pp. 62-68, Jan. 1994.

[9]

Geissler, W., "Berechnung der Druckverteilung

an oszilliercnden Tragfliichen mit endlicher Dicke in inkompressibler Stromung," DFVLR-AVA-Internal Hept. 76-48, 1976.

[10] Archibald, F.S. Unst.eady Kutta Condition at High Values of the Reduced Frequency Pararn-cter", Journal of Aircraft, vol. 12, pp. 545-550, 1975.

[11] Commerford, G.L. and Carta, F.O., "Unsteady Aerodynamic Response of a Two~Dimensional

Airfoil at High Values of Ttcduccd Frequency," AIAA Journal, vol. 12, no. 1, pp. 43-48, 1974.

[12]

Satyanarayana, 13., and Davis, S., "Experimental

Studies of lJnsteady Trailing Edge Conditions," A!AA Journal, vol. 16, no. 2, pp. 125-·129, 1978. [13] Flcct.cr, S., "Trailing Edge Condition for

Un-steady Flows at High Reduced Frequency," AIAA Paper 79-0152, Jan. 1979.

(11] Poling, D.Ft., and Telionis, D.P., "The Response of Airfoils to Periodic Disturbances ~·· The Un-. steady J{uttUn-.a Condition,'' AIAA Journal, volUn-. 24,

no. 2, pp. 193-199, 1986.

(15] Katz, .J., and Maskcw, B., "Unsteady Low-Speed Aerodynamic i'vfodel for Complete Aircraft Con-figurations," Journal of Aircraft, vol. 25, no. 4, pp. 302-310, 1988.

[16] Moore, D.W., "A Numerical Study of the Roll-up of a F\nit.e Vortex Sheet," Journol of Fluid

Mechanics, vol. G:l, pp. 225, 1974.

[17]

Graham, ,J.ivl.R.,"Application of Discrete Vor-tex i\'fcthoc!s to the Computation of Separated Flov·:s", in Nmncrical fvlcthocls for Fluid Dynam-ics H, Oxford University Pre-ss, Oxford, pp. 273-302, 1986.

0 ~ 0 II -"' 0

"'

0 II .<

s

0 II -"'

8

0 II -"' CL 0.6 0.4 0.2 0.0 ··--·---- - r- . -0.2 -0.4 -0.6 -5 0.6 0.4 0.2 0.0 -0.2 -0.4 ----0.6 -5 0.6 0.4 0.2 0.0 -0.2 --0.4 --0.6 -5 0.6 0.4 0.2 0.0 -0.2 --2.5 0 --2.5 0

---.--- -~--25 0 0.600 . 0.500 0.400 0.300 ---, ~2.5 0

n ( dcg.)

2.5 2 ,.

_"

2.5 2.5 -0.040 '--~--1--~--' 5 -5 -2.5 0 2.5 5 0.040 --~--,.--~-0.020 0.000 -0.020 -0.040 _L....-._,,______1~--....L...-5 --5 0 2 5 5 0.020 0.010 0.000 -0.020 '--~-~-~-5 -5 --2.5 0 2.5 5 0.008 0.004 0.000 -0.004 ---2.5 0 2.5 5

n (dcg.)

!; lg"llrC 5: CL and eM of the NACAOO! '2 scc!.i(l!l in n ::::: sin '2./.:r

pitching motion at 'Zero lllC<Hl incidence for value~ of Cw fr01n

(UOO to O.GOO (N

=

90, n1"

=

100). 54.10

6CM

0.004 0.002 0.000 -0.002 -0.004 0.004 0.002 0.000 -0.002 -0.004 0.004 0.002

··-'''::;;;:;:it::::.::>

/',r = 0.097' -5.0 -2.5 -5.0 -2.5 0.0 0.0 I I 2.5 2.5 5.0 5.0 - - c = - -0.000

c

---"---0.002 -0.004 0.004 0.002 0.000 -0.002

·rr

!

/',r = 0.021! -5.0 -2.5 0.0 2.5 5.0

.----.----,--,---r·

i

/',r = 0.012! -0.004 ...l...L...___...l. -5.0 -2.5 0.0 2.5 5.0 (Y (w 0.60 t0 0 50 (,, 0.50 to 0.40 (w 0.40 to 0.30

Figure 6: Effect. of timc-st.ep size on CAt difference for incrcntents in (w of 0.10 (NACAOOI2 o

=

sin 2kr, k

=

0.:300).

-0.0!0 .... L~~-~~-~ -5.0 -2.5 0.0 2.5 5.0 0. 0 1 0 1·---....,....--- r~-~~-0.005 0 000 -0.005 -Q.QJQ L-~----· _.L.~.___l _____ j _ _ _ _ -5.0 --2.5 0.0 2.5 5.0 0.010 ---- r ---.. ·r ·-~--~r----0.005 0.000 -0.005 -0.0!0 -5.0 -2.5 0.0 2.5 5.0 0.010 0.005 0.000 -0.005 -S.O -2.5 0.0 2.5 5.0 { \ t:ir -0.097 L\r .. O.CHD L\T -- 0.021₁

Fig me 7: Efrcct of (w 011

eM

CO!l-vcrgct\('C for decreasing tinH>step.

Tl1c CM solution for ~r

=

0.012

IS used as datutn (NACi\00!'2 <> CC'/ sin 2h, k I I o.:lOO). 0.8 0.6 0.4 0 2 4 0.8 0.4 0.0 -0.4 -0.8 -1.2 0 2 4 0 0 (w = 0.300 (w

=

0.500 (w=0.700 Wagner function 6 8 T 6 8 T 10 !0

Figure 8: Transieut. C,, and

Cu

for NIIC:AOOJ2

section started impulsively into constant. speed

for-ward motion at o· :::: 5°. The normalising

stcady-strltc values arc C1.,. = 0.585 and CM,. = 0.{){)!):l2 with c':.r

=

0.02. f'iC,_,TB 0.2 · -0.1 0.0 -0.1 0.2 ---O.O?S explicit k=0.3 implicit k:::::0.3 implicit k::::0.2 implicit k:::::O.l -0.000

l'w

0.025

Figure D: Pnf'onnancc of the 1\utta condition::: con1par('d for

the Nt\Ci\0012 S(~ction in harmonic pit.ch tno!.ion (n::::: sill Lkr. n.11c:::::::: tj()(), DC1,n,· tol(:rancc for the P.Xplicit 1\utt.a cottdit iott set

to O.OOo).

Cr

eM

0.20 -.,---...,.--- 0.050 0.10

1

0.025

)\

I ...

----·

\

0.00

\

0.000 !

'

'

2 4 G 8 10 ? 4 6 8 10 T T

Figure !0: Vortex trajectory Cl!ld

c,_

and(.'"' f!tJct.uat.iOII during the intnactioll of a free vortex with a NACA0012 section flying at zero incidence rclat.ivc to the nndisturhccl flow. Tlw vortex

passeS t!Jc <lCfOfoi) leading edge at. T ~ iJ.7;. ( J' = --·0. J)

CIII 0.025 ---,---.--- ··-...---,--0.000 - __ --0.025 2 4 6 B 10 T T

Figure 11: Vortex I rajcct.ory :1JHI (,'1, ;utd (',\t flt~ctuat.ion dming tiH' inkradion of a free vor!.r'.\ wi Lh a N ACA 001 '2 ~eel. ion flying aL '/,('fn i llcid('!!Ci' rch 1.ive t.o the l! 11d islurhcd flow. Tl!i~ vorl n.:

passes !.he acrofoil !c:1ding, ul~c at r ~-I ;J. ( ]' 0.1)

0

"'

0 II

"'

0 N 0 II

"'

0 II

"'

N 0 0 II

"'

CL 0.6 0.4 PANEL 0.2 ULER 0.0 --0.2 -0.4 -0.6 ~--·5 -2.5 0 2.5 0.6 0.4 PANEL EULER 0.2 0.0 --0.2 .()_4 ··5 -2.5 0 2.5 0.6~-~-0.4 0.2 0.0 --0.2 ---0.4 -0.6 .,{/

~9)

)f···

,.:(j/:.-~

iPANEL .,., EULER

;/}{{~-~--PANEL

S-S ---'---l _______ [_ _______ _ . -5 --2.5 0 2.5 n· (dcg.) 5 5 5 Figure l2:

(;a i Loll •.I •fl)

Con1parison \\'i!.ll

1'11r lhc Ni\C·\111112 111o1.i(lll :d. zero 111e;ln incidr·nn•.

eM

0.040 0.020 0.000 -0.020 --0.040 ---l....----__ L...__ _ _ _J. ______ ··5 -2.5 0 2.5 0.030 ---...---0.020 0.010 0.000 --0.010 -0.020

>

--0.030 _________ t _________ L__ _ _ _ l____ --5 -2.5 0 2.S 0.015 0.010 0.005 I 0.000

\

--0.005 ----0.010 -0.015 ______ j - - · - - - - ._L__~---5 -2.5 0 2.5 0.008 0.004 0.000 0.004

--o.ooa ________

j_ ___ _ --5 --2.5 ·-.. -, __ 0 2.5 5

)

,.

·'

·~

)

5

EULI~H IP.cl.li<ld !'l'~idL:-;

or

:-;ec1.io11 in n :--·· ~i11 "2/,·r pitching

51.13 0. u I 0. u I 0. u I 0. u I 0'

=

oo

dcxjdl

=

Omax 0.6 0.4

~...__.._

0.2 0.0

... 0 ..

~~~ -0.2 PANEL -0.4 _EULER -0.6 -0.8 -1.0 ·---' 0.0 0.2 0.4 0.6 0.8 1.0 xfc Q' =5° do:jdt = 0 1.0 0.5 0.0 PANEL -0.5 EULER -1.0

L_-~---'~-~---~-J

0.0 0.2 0.4 0.6 0.8 1.0 x/c n = 0° dnfdt = .. -n'""" 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 --0.8 -1.0 0.0 1.5 1.0 0.5 00 PANEL EULER 0.2 0.4 0.6 0.8 1.0 xfc Q' =-5° do:jrlt. = () i I

'

I -·O.S _I I .. 1.0 0.0 0.2 0.'1 0 6 O.B 1.0 x/c

Figure 1:>: Co!llparisotl of pn•:-::--un·

di:---tribtJLions ;d. t.hc quartcr-ph;:1:--e 1h1~itiot1:-;

of the N:\Ci\0012 section ill pit.cl1 n

" ' " ' 0 dOf/dl=O

dOf/dt = 0

dOf/dl = -n

dOfjdt ::; 0

....

<>

"f > 0 (E:ULE:R CONTOURS) < ::> "r < 0 (EULER CONTOURS) • PAN8L VORTEX POINT r > 0 0 PANE:L VORTEX POINT r > o

Figure 14: Euler-method vorticity contours with panel method wake geometry (wake strength indicated by point size) for the quarter-phase positions of the pitching motion (The vertical scale is exaggerated by x2.)

-/\',,,., ;:c: ]()()

CL

0.50 ~-~--~--~-~ 0.25 /(max

=

•100: /\'max= 1001 f(m<u: = •!1 · I 0.00 --0.25 -0.50 '---'---~~ -5 -2.5 0 (Y 2.5 5 ()M 0.050 ~-·---~---.---·

~~=~-0025 0.000

i

-0025

_--1----//

!

_./

j -0.050 t _ _ _ L,_, _ _____j___ _____ L.,. -5 -2.5 0 n 2.5 5

Figure 15. !(mar sensitivity cornp<-~r­

isons for NACJ\0012 inn-= sin2kr motion at k

=

o.:wo

with 11pc = t100 .

...

.-.;·:---~

/ _/',.

"~---~--Fig1₁rc [(): /\'

111a.r sCJlsit.ivii.J' COill[lariso!l of \\"akr• sllii[H's for lllc f\1,\(',.\001'2 ;-1fln t.I11Tr'

cycles of n- -:-.:: sinwf. I!IOI.iOH ;d, ~· ·:_:()_;\()()with Ill''":::;[()()_