Computation of the Flows over Flapping Airfoil by the - Feng Liu

20 Pages·2004·8.95 MB·English

theAIAA43rdAerospaceSciencesMeeting,Reno,NV,Jan.10–13,2005
Computation of the Flows over Flapping Airfoil by the
Euler Equations
Shuchi Yang∗, Shijun Luo†, Feng Liu‡
Department of Mechanical and Aerospace Engineering
University of California, Irvine, CA 92697-3975
Her-Mann Tsai§
Temasek Laboratories
National University of Singapore
Kent Ridge Crescent, Singapore 119260
To investigate the mechanism of thrust generation by flapping airfoils, the inviscid ver-
sion of a three-dimensional unsteady compressible Euler/Navier Stokes flow solver is used
to simulate the flow field around flapping airfoil NACA 0012 at low speeds. Sinusoidally
plunging or/and pitching oscillations are studied. The compressible Euler code with a
low free-stream Mach number of 0.1 or 0.05 can simulate the incompressible flows without
leading-edge separation. The wake vortex structures are visualized by numerical meth-
ods: vorticity filled-contours, perturbation-velocity vector plots, and streamlines. The
computed wake-vortex structures almost coincide with known low-speed test results. The
time-averaged thrust coefficients, input-power coefficients, and the efficiency over a period
of the oscillation versus the Strouhal number based on the total excursion of the trailing
edge of the airfoil for a number of combined plunging and pitching cases agree well with
the linear and non-linear incompressible potential-flow solutions in the literature. In com-
parison with known test data, the accuracy of the computed forces decreases when the
leading-edge vortices appear and become strong enough to interfere with the trailing-edge
vortices.
I. Introduction
Humanity’s earliest dreams for flight were inspired by birds and bats. Little progress was made on the
understanding of the mechanism of the flapping wings until 1910’s when Knoller1 and Betz2 noted that the
flappingmotionofanairfoilproducesaneffectiveangleofattackresultinganormalforcevectorwhichhasa
component in forward direction, the thrust. This phenomena was observed by Katzmayr3 in his experiment
in 1922. Later Freymuth4 experimentally demonstrated that an airfoil is capable of generating thrust when
itundergoeseitherpurepitchingorpureplunging. Acomprehensivereviewoftheresearchanddevelopment
results of flapping-wing propulsors was given by Rozhdestvensky and Ryzhov.5
Water-tunnel tests of a NACA 0012 airfoil oscillated sinusoidally in plunge were performed by Jone,
Dohring and Platzer6 and Lai and Platzer.7 They provided dye flow visualizations over a broad range of
reduced frequencies and nondimentional plunge amplitudes and compared the experimental data with their
numerical computations by an incompressibe potential flow code.
Neef and Hummel8 computed the two-dimensional airfoils and three-dimensional high aspect ratio wings
byusingEuler equations . The reduced frequency islimited to about 0.1to simulatethe flowof the cruising
of large birds. In that paper, they investigated the asymmetric thrust distribution during each period, and
numerically demonstrated that the additional lift has no influence on thrust output.
∗Post-DoctoralResearcher. AIAAmember
†Researcher.
‡Professor. AssociateFellowAIAA.
§PrincipalResearchScientist,AIAAMember
1of20
AmericanInstituteofAeronauticsandAstronautics
Lewin and Haj-Hariri9 modeled the thrust generation of a symmetric Joukowski airfoil with plunging
motionbysolvingthevorticityequationsforatwo-dimensionalviscousincompressibeflow. Theyinvestigated
how the fate of the leading-edge vortex affects the wake and efficiency.
Anderson10 andAnderson,Streitlien,Barrentt,andTriantafyllou11 obtainedvisualizationandforcedata
for a plunging and pitching airfoil NACA 0012 moving in a water testing tank facility, and found that
agreement between theoretical predictictions of a nonlinear incomressible potential flow computation and
experiment is good when either very weak or no leading-edge vortices form.
Tuncer and Platzer12 used a compressible Navier-Stokes solver to compute the unsteady flow fields and
obtained maximum propulsive efficiency when the flow remains mostly attached over the airfoil oscillated in
a combined pitch and plunge. Recently, Young and Lai13 used a compressible Navier-Stokes solver and an
incompressible potential flow code to simulate the flow over a NACA 0012 airfoil oscillated sinusoidally in
plunge. They found that leading-edge separation appears to dominate the generation of aerodynamic forces
for low reduced frequencies and high dimensionless amplitudes, but becomes secondary for high frequencies
and low amplitudes.
The objective of this study is to explore the application of a compressible Euler solver to simulate the
low-speed flow over an airfoil NACA 0012 oscillated sinusoidally in plunge and pitch. In the following
sections, the Euler solver and the computation model are described. The grid refinement and the time-step
refinementareinvestigatednumerically. TheimplementofthecompressibleEulersolvertosimulatethelow-
speedflowisvalidatedbynumericalexperimentsandincompressiblepotential-flowsolutions. Thecomputed
wake-vortex structures and aerodynamic forces are presented and compared with known experimental data
and numerical potential flow and Navier-Stokes flow solutions. Lastly, conclusions are drawn.
II. Euler Solver
It is known that the Euler solver can capture automatically the shear layer separated from the sharp
trailing edge and its rolling up into a vortex cores while convected downstream. Although the secondary
separations such as those shown in Figure 14 of Young and Lai13 on either or both sides of the airfoil
near the trailing edge from their Navier-Stokes code are absent in the Euler solutions, the gross dominant
characteristicsofthewakeflowfield,i.e. theprimarytrailing-edgevortexconfigurationsandtheirinteractions
with the body surface are reproduced as long as there are no leading-edge separation vortices.
The present Euler solver is based on a multi-block, multigrid, finite-volume method and parallel code for
the three-dimensional, compressible steady and unsteady Euler and Navier-Stokes equations. The method
usescentraldifferencewithablendofsecond-andfourth-orderartificialdissipationandexplicitRunge-Kutta-
type time marching. The coefficients of the artificial dissipation depend on the local pressure gradient. The
orderofmagnitudeoftheaddedartificialdissipationtermsisoftheorderofthetruncationerrorofthebasic
scheme,sothattheaddedtermshavelittleeffectonthesolutioninsmoothpartsoftheflow. Nearthesteep
gradients the artificial dissipation is activated to mimic the physical dissipation effects. The resulting code
preserves symmetry. Unsteady time-accurate computations are achieved by using a second-order accurate
implicitschemewithdual-timestepping. Thesolverhasbeenvalidatedforanumberofsteadyandunsteady
cases.14–17
III. Validation for Low-Speed Flow
The present Euler solver was originally designed for compressible flows. It is known that the numerical
solutionofacompressibleflowsolvermaynotconvergetothephysicalincompressibleflowasthefree-stream
Mach number goes to zero. To remove this problem, a preconditioning techniques to treat the compressible-
flow solver are proposed by Turkel.18 However, in many cases, there exists a small nonzero free stream
Mach number at which the compressible code would yield good approximate incompressible flow. Instead
ofimplementingthe preconditioningtechniques, theappropriete smallfree-stream Machnumber issearched
for by numerical experiments, if it exists.
Asanexample,thesteadyincompressibleinviscidflowovertheairfoilNACA0012atα=10◦iscomputed
by the present Euler code with M =0.2, 0.1, 0.05, and 0.025. The grid of the O-type is 193×33 in the
∞
circumferentialandradialdirections, respectively. Thefarfieldboundaryispositionedawayfromtheairfoil
by 20 times its chord . The initial solution for the Euler computation is the free-stream flow.
Fig. 1 shows the computed pressure distributions over the upper and lower surfaces of the airfoil NACA
2of20
AmericanInstituteofAeronauticsandAstronautics
-7
-6
-5
PanelMethod
-4 M∞=0.20
M∞=0.10
-3 M∞=0.05
Cp M∞=0.01
-2
-1
0
1
2
0 0.2 0.4 0.6 0.8 1
x/c
Figure 1. Pressure distributions cp of steady flow over airfoil NACA 0012 at α=10◦ computed by the Euler
code with different M∞ and by incompressible panel method.
0012 in comparison with those calculated by a panel method for incompressible inviscid flow. It is seen
that the numerical solution of the Euler solver at M =0.1 coincides with the incompressible-flow solution
∞
calculated by the panel method, while the Euler numerical solution at M = 0.01 runs away from the
∞
incompressible-flow result. Hence for this example, the present compressible Euler solver simulates the
incompressible flow very well when M is taken to be 0.1.
∞
For the unsteady low-speed flow considered in this paper, the M -validation of the Euler solver will be
∞
studied after the grid and time-step refinement investigations.
IV. Computational Model
The airfoil NACA 0012 with chord length c performs a harmonic plunging motion h(t) and a harmonic
pitchingmotionθ(t)aroundthepivotpointO locatedbehindtheleadingedgebyadistanceb,inanuniform
flow of low speeds as shown in Fig. 2. The free-stream velocity and Mach number are U and M . The
∞ ∞
pitching movement and the plunging movement for the airfoil are described below.
θ(t)=θ sin(ωt+φ) (1)
0
h(t)=h sin(ωt) (2)
0
θ(t) b c
O
Y h(t)
-dh(t)/dt
U
∞
X
Figure 2. Airfoil in plunging and pitching motion.
where ω is the angular frequency, φ is the phase difference by which the pitching leads the plunging.
The instantaneous angle of attack of the airfoil at the pivot point O consists of two parts: one part due
3of20
AmericanInstituteofAeronauticsandAstronautics
to the pitching motion of the airfoil and the other part due to the plunging motion of the airfoil.
(cid:181) (cid:182)
1 dh(t)
α(t)=θ(t)−tan−1 (3)
U dt
∞
Theangle of attackhas a maximumvalue, α . A nominal angle of attack, α introduced byAnderson
max 0
et al. is defined as.11
(cid:181) (cid:182)
ωh
α =θ −tan−1 0 (4)
0 0 U
∞
The reduce frequency k is defined as
ωc fcπ
k = = (5)
2U U
∞ ∞
where f denotes the frequency of airfoil oscillation in Hz, i.e., f =ω/(2π). For flapping cases, the Strouhal
number is used, St is defined as
fA
St= (6)
U
∞
where A denotes the characteristic width of the wake flow. Usually A=2h .
0
Alternatively, a Strouhal number based on the total excursion of the trailing edge of the airfoil A is
TE
widely used and denoted by St .
TE
fA
St = TE (7)
TE U
∞
The instantaneous aerodynamic force components X(t) and Y(t) are in the x- and y- directions, respec-
tively. M(t)isdefinedastheinstantaneouspitchingmomentaroundthepivotpointO showninFig. 2. The
corresponding non-dimensional coefficients based on the free-stream dynamic pressure and the airfoil chord
are c (t), c (t) and c (t), respectively. The instantaneous pressure coefficient is denoted as c (t).
x y m p
The time-averaged value of −X(t) over a period of the oscillation denoted by F is the thrust generated
by the airfoil.
(cid:90)
1 T
F =− X(t)dt (8)
T
0
The time-averaged power input is
(cid:34) (cid:35)
(cid:90) (cid:90)
1 T dh(t) T dθ(t)
P = Y(t) dt+ M(t) dt (9)
T dt dt
0 0
The thrust coefficient C and the power input coefficient C are based on the free-stream dynamic
T P
pressure and the airfoil chord. The propulsive efficiency η is defined below.
FU C
η = ∞ = T (10)
P C
P
V. Grid and Time-Step Refinement and M -Validation
∞
The plunging oscillation of the airfoil NACA 0012 with reduced frequency, k =4.0 and non-dimensional
plunging amplitude, h /c = 0.0125 is computed. There is no leading-edge separation in this case as shown
0
by Young et al.13 by the Navier-Stokes computation in their Fig. 5.
The C-type grid is used for all the computations in this paper. The far field boundary is located at 20c
away from the airfoil. The coarsest grid is 289×33 with 193 grid points along the airfoil circumference. In
the downstream region of the airfoil, there are 48 vertical grid lines and 65 horizontal grid lines.
Thecoarsestgrid289×33isshowninFig. 3: (left)givestheover-allgridconfigurationwhereonlyevery
fourthgridlinesaredrawnforclarity,and(right)givesaclose-uplookofthegridintheneighborhoodofthe
4of20
AmericanInstituteofAeronauticsandAstronautics
40
0.6
0.4
30
0.2
y/c 20 y/c 0
-0.2
10
-0.4
-0.6
0
-20 0 20 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
x/c x/c
Figure 3. Grid 289×33 configurations: (left) Over-all view with only every 4th grid line shown for clarity,
(right) Close-up look of the neighborhood of the airfoil and its nearby wake.
airfoil and its nearby wake. The grid is refined by inserting a new grid point between every two neighboring
grid points in either or both circumferencial and radial directions and thus the seven grids are formed. The
free-stream Mach number, M = 0.1. The real-time steps in one cycle, nstep = 64. The instantaneous
∞
drag coefficient, c versus h/c is plotted to show the grid dependence of the numerical solutions. The single
x
precision is employed for the Euler computation in this paper.
In the dual-time computation, to ensure a real-time step solution having a sufficient accuray, the steady
modeoftheEulercodeisrequiredtorun120multi-gridcycles. Startingfromafree-streamflowastheinitial
solution, the real-time unsteady solution becomes periodical after 3 cycles of the oscillation. The results of
the fifth cycle is used as the final numerical solution.
Fig. 4 gives c versus h/c for the different grid refinements: (a) 289×33, 577×33 and 1153×33, (b)
x
577×33, 577×65 and 577×129, (c) 289×65, 577×65 and 1153×65. It is seen that the numerical results
of c computed with the grid 289×65 is nearly grid independent.
x
Furthermore, Fig. 5 gives the instantaneous lift coefficient, c and moment coefficient, c versus h/c,
y m
and the pressure coefficient c versus x/c at ωt=0 with three different grids. It is seen that the numerical
p
results of all the aerodynamic forces computed with the grid 289×65 are grid independent.
For the computation of the wake structures, the grid 1153 × 65 is used. There are 769 grid points
along the entire surface of the airfoil. There are 192 vertical grid lines and 129 horizontal grid lines in the
region downstream of the airfoil. Also, a grid-independence experiment for the wake computation is made
by refining the grid downstream of the airfoil while the grid-point number along the airfoil circumference
remains to be 193 as shown in Section VII B.
The time-step refinement is performed with the number of steps in one cycle of the oscillation, nstep=
64,128,256, and 512. The grid, 289×33 and the free stream Mach number, M = 0.1 are used in this
∞
computations. The numerical results show that the numerical solution with nstep = 64 is time-step inde-
pendent.
TousethecompressibleEulersolvertosimulatethelow-speedflowovertheplungingairfoil,theappropri-
ate free-stream Mach number, M is to be determined by numerical experiments with M =0.2,0.1,0.05,
∞ ∞
and 0.025. The grid 289×33 and nstep=64 are used in this computations.
Fig. 6 gives the computed pressure coefficients, c versus x/c when ωt = 0 and π/2 with M =
p ∞
0.200,0.100,0.050, and 0.025. It is seen that the computed pressure distributions have an anomaly behavior
near the trailing edge of the airfoil. The anomaly increases as M is decreased. This anomaly may be
∞
associated with the implementing of the compressible-flow solver to incompressible flow. Except this very
small neighborhood region of the trailing edge, the numerical results show that M = 0.05 would be
∞
an acceptable choice for simulating the low-speed flow by the present compressible code. This conjecture
is further confirmed by the good agreement of the present computation results with the incompressible
potential flow solutions calculated by Anderson et al.11 as shown in Section VII A.
5of20
AmericanInstituteofAeronauticsandAstronautics
0.005
0 258797××3333
1153×33
-0.005
x
c
-0.01
-0.015
-0.02
-0.015 -0.01 -0.005 0 0.005 0.01 0.015
h/c
0.005
0 557777××3635
577×129
-0.005
x
c
-0.01
-0.015
-0.02
-0.015 -0.01 -0.005 0 0.005 0.01 0.015
h/c
0.005
0 258797××6655
1153×65
-0.005
x
c
-0.01
-0.015
-0.02
-0.015 -0.01 -0.005 0 0.005 0.01 0.015
h/c
Figure 4. Instantaneous drag coefficient cx versus h/c computed with different grids: 289×33, 577×33, and
1153×33; 577×33, 577×65, and 577×129 and 289×65, 577×65, and 1153×65.
6of20
AmericanInstituteofAeronauticsandAstronautics
1 289×65
577×65
1153×65
0.5
y 0
c
-0.5
-1
-0.01 -0.005 0 0.005 0.01
h/c
0.4
289×65
0.2 577×65
1153×65
m 0
C
-0.2
-0.4
-0.015 -0.01 -0.005 0 0.005 0.01 0.015
h/c
-1.5
-1
-0.5 527879××6655
1153×65
p 0
c
0.5
1
1.5
0 0.2 0.4 0.6 0.8 1
x/c
Figure 5. Instantaneous lift coefficient cy and moment coefficient cm (about O at b = c/4) versus h/c, and
pressure coefficient cp versus x/c at ωt=π/2 computed with different grids: 289×65, 577×65, and 1153×65.
7of20
AmericanInstituteofAeronauticsandAstronautics
-1.5 -1.5
-1 -1
-0.5 -0.5
p 0 p 0
c c
0.5 0.5
1 MMMM∞∞∞∞====0000....210000520005 1 MMMM∞∞∞∞====0000....210000520005
1.5 1.5
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
ωt=0 x/c ωt=π/2 x/c
Figure 6. Instantaneous pressure coefficient cp versus x/c at ωt = 0 and π/2 computed with M∞ =
0.200,0.100,0.050, and 0.025.
VI. Computational Results for Pure Plunging
PuresinusoidalplungingoftheairfoilNACA0012atlowspeediscomputedbytheEulersolverwiththe
grid 1153×65, nstep = 64, and M = 0.05. The wake structures are shown by the vorticity filled-contour
∞
plots and streamlines when the airfoil is located at typical positions. The angle of attack, lift coefficient and
drag coefficient versus time are plotted to enunciate the mechanism of the aerodynamic-force generation.
Two cases are studied: (1) k = 4.0,h /c = 0.0125, (2) k = 1.5,h /c = 0.2. According to the data in the
0 0
literature, there are no leading-edge vortices in the two cases.
A. k =4.0,h /c=0.0125
0
Fig. 7 gives the vorticity filled-contour plots at ωt = 0,π/2,π, and 3π/2, where + and - denote counter-
clockwise and clockwise directions, respectively. If the wake would move downstream with the free-stream
velocity, the wake wavelength, λ/c should be π/k = 0.785. The above plots match with this value closely.
It is seen that the wake configurations at ωt=0 and π are symmetric with respect to airfoil chord line, but
the vorticities have opposite sign. The wake configurations at ωt = π/2 and 3π/2 are similarly symmetric.
The present results at ωt = π agree qualitatively with the test data and the Navier-Stokes solutions given
in Figure 7 of Young et al.13 The only difference is that in their results two same-sign vortices shed per
half-cycleoftheairfoilandinourresultsoneelongatevortexshedperhalf-cyclevortex. However,thelength
of the latter elongate vortex is equal to the span of the former two same-sign vortices.
Fig. 8givesthecorrespondingstreamlineplotsofthevelocityfieldrelativetotheairfoil. Theyshowthat
the flows are separated smoothly at the trailing edge all the times. The streamlines at ωt = 0 match well
with those of the Navier-Stokes result given by Young et al. in their Fig. 14 except the minor separations
near the sharp trailing edge which are absent in the Euler solution.
Figure 9 gives α,c and c versus ωt. The amplitudes of α,c , and c are 5.71◦,1.418, and 0.007925,
y x y x
respectively. The lift has the same period of the airfoil oscillation. The drag’s period is one half of the lift’s
period. The lift and drag head before the angle of attack. The lead phase angles for c and c are 1.418 and
y x
1.375 rad., respectively. The time-averaged lift over a period is zero. The time-averaged drag is nonzero.
The time-averaged thrust coefficient computed by the Euler solver, C =−0.0081. A skin friction force
T
has to be subtracted from the inviscid-flow solution to simulate the real flow. An estimation of the skin
friction is obtained by the flat-plate skin-friction formulae. At the experimental Reynolds number based on
theairfoilchord,Re=2×104,theskinfrictioncoefficientoftheairfoil,C =0.0188. Andthusforrealflow,
f
C = −0.011, which is a positive drag. Tuncer et al.12 gave the time-average drag coefficient 0.008 from
T
experimental data and Navier-Stokes computations in their Table 1 for the case k = 3.925,h /c = 0.0125.
0
The agreement is good in consideration that the magnitude of C is very small.
T
8of20
AmericanInstituteofAeronauticsandAstronautics
- + - + + - + -
ωt=0 ωt=π/2
+ - + - - + - + -
ωt=π ωt=3π/2
Figure 7. Computed vorticity filled-contour plots for k=4.0,h0/c=0.0125 at ωt=0,π/2,π, and 3π/2, where +
and - denote counter-clockwise and clockwise directions, respectively, with M∞=0.05 and grid 1153×65.
9of20
AmericanInstituteofAeronauticsandAstronautics
0.15 0.15
0.1 0.1
0.05 0.05
y/c 0 y/c 0
-0.05 -0.05
-0.1 -0.1
ωt=0 ωt=1/2π
-0.15 -0.15
0.9 1 1.1 1.2 0.9 1 1.1 1.2
x/c x/c
0.15 0.15
0.1 0.1
0.05 0.05
y/c 0 y/c 0
-0.05 -0.05
-0.1 -0.1
ωt=π ωt=3π/2
-0.15 -0.15
0.9 1 1.1 1.2 0.9 1 1.1 1.2
x/c x/c
Figure8. Computedinstantaneousstreamlinesofthevelocityfieldrelativetotheairfoilfork=4.0,h0/c=0.0125
at ωt=0,π/2,π, and 3π/2, with M∞=0.05 and grid 1153×65.
2 8
Cy
Cx
1.5 α 6
1 4
0.5 2
0
0
cx*1 0 0 αo()
y
c
-0.5 -2
-1 -4
-1.5 -6
-2 -8
0 2 4 6 8 10 12
ωt
Figure9. Instantaneousangleofattackα,liftcoefficient,dragcoefficientforYoungetalcase.13 k=4.0,h0/c=
0.0125.
10of20
AmericanInstituteofAeronauticsandAstronautics