Incompressible Flow Over Airfoils

Prof. SungKi Jung, Ph.D.

Universidade Federal do ABC,
Aerospace Engineering, Centro de Engenharia, Modelagem e Ciências Sociais Aplicadas,
Rua Arcturus 03, 09606-070, São Bernardo do Campo, São Paulo, Brazil.
E-mail: [email protected]

1 Airfoil
1.1 Airfoil Nomenclature

Figure 1: Airfoil nomenclature

Symmetry airfoil generator of NACA 4 Series is

√
yt = 5t(0.2969 x − 0.1260x − 0.3516x2 + 0.2843x3 − 0.1015x4 ), (1)

where, t is the maximum thickness and a ratio based on the chord length. For NACA0012 airfoil,
t = 0.12 or 12 percents of the chord length. For the general NACA 4 series airfoil, we follow the rule,
like as NACAabcd, where a is the maximum camber, b is the location of maximum camber, cd is the
maximum thickness. For example, NACA2412 means that 2 indicate the 2 percents of chord length
for the maximum camber, 4 indicate the 40 percents of chord length for the location of maximum
camber and 12 indicate 12 percents of chord length for the maximum thickness.
A short partial list of airplanes in service in 1982 which use standard NACA airfoils is shown in
table 1. In addition, many of the large aircraft companies today design their own special purpose
airfoils; e.g., the Boeing 727, 737, 747, 757, and 767 all have specialJy designed Boeing airfoils.
Such capability is made possible by modern airfoil design computer programs utilizing either panel
techniques or direct numerical finite-difference solutions of the governing partial differential equations
for the flow field.

1
 Table 1: Commonly used airfoil series
Airplane Airfoil
Beechcraft Sundowner NACA63A415
Beechcraft Bonanza NACA23016 Root
NACA23012 Tip
Cessna 150 NACA2412
Fairchild A-10 NACA6716 Root
NACA6713 Tip
Gates Learjet 24D NACA64A109
General Dynamics F-16 NACA64A204
Lockheed C-5 Galaxy NACA0012

2 Airfoil Characteristics
A movie for flow visualization around the airfoil is Flow behaviors around the airfoil.

2.1 Lift
At the moderate angles of attack, the lift slope is linear. Around the maximum lift coefficient, the lift
curve is non-lnear due to the flow separation induced by the viscous effects. Lift is a function of an
angle of attack, Reynolds number and Mach number, cL = f (Re, M, α).

Figure 2: Schematic of lift-coefficient variation with angle of attack for an airfoil

Definitions:
lif t
ˆ cl = 1
ρu2 S
: Lift coefficient, where S is a reference area
2

dcl
ˆ clα = dα : Lift slope

ˆ αl=0 : An angle of attack at zero lift

ˆ clmax : An angle of attack at the maximum lift coefficient

2
Camber Effects on Lift

Having a camber in the airfoil increases the highest attainable lift coefficient. However, there is a
decrease in the critical angle of attack.

Figure 3: Camber effect

3
Mach Number Effects on Lift

As Mach number increases, the slope of lift coefficient is changed.

Figure 4: Lift slope due to Mach numbers

4
High-Lift Devices on Lift

If using the high-lift devices, the maximum lift coefficient can be increased.

Figure 5: High-lift devices effects on the maximum lift coefficients

5
2.2 Drag
Drag consists of various parameters, like as parasite drag, skin friction drag, induced drag, etc. In
this chapter basic characteristics of drag on airfoil is shown in fig. 6. Also, the drag coefficient is a
function of an angle of attack, Reynolds number and Mach number, Cd = f (α, Re, M ).

Figure 6: Drag coefficient due to cambers

Definitions:
Drag
ˆ cD = 1
ρu2 S
: Drag coefficient
2

6
2.3 Pitch Moment
In aerodynamics, the pitching moment on an airfoil is the moment produced by the aerodynamic force
on the airfoil if that aerodynamic force is considered to be applied, not at the center of pressure, but
at the aerodynamic center of the airfoil. The pitching moment on the wing of an airplane is part of
the total moment that must be balanced using the lift on the horizontal stabilizer. More generally, a
pitching moment is any moment acting on the pitch axis of a moving body.
The lift on an airfoil is a distributed force that can be said to act at a point called the center of
pressure. However, as angle of attack changes on a cambered airfoil, there is movement of the center
of pressure forward and aft. This makes analysis difficult when attempting to use the concept of the
center of pressure. One of the remarkable properties of a cambered airfoil is that the pitch moment
changes only a little over the operating range of angle of attack of the airfoil. Thus,

dcm
≈ constant, (2)
α

where cm is the coefficient of pitch moment.

Figure 7: A graph showing coefficient of pitching moment with respect to angle of attack. The negative
slope for positive α indicates stability in pitching.

Definitions:

ˆ cm = M oment
1
ρu2 Sc
: Pitch moment coefficient, where c is a reference length
2

ˆ Aerodynamic center (AC): A position that cmα is constant, where cmα = dcm
dα . In general, AC is
at c/4 of airfoil chord for the subsonic flow.

Useful sites:
Airfoil Tools
Airfoil Coordinates Database

7
3 Theoretical Solutions for Low-Speed Flow over Airfoil: the Vortex
Sheet
Vortex filament: The flow induced in any plane perpendicular to the straight vortex filament is
identical to that induced by a point vortex of strength.

Figure 8: Vortex filament.

Vortex sheet: An infinite number of straight vortex filaments side by side, where the strength of
each filament is infinitesimally small. These side-by-sdie vortex filaments is called a vortex sheet.

Figure 9: Vortex sheet.

8
 Now, an infinitesimally small velocity, dV 1 , at point P in fig. 9 is given by

γds
dV = , (3)
2πr

where, γ is the strength of the vortex sheet, per unit length along s, thus γ = γ(s).
In the appendix A-11 , the circulation Γ is
R Z
dΓ = γds −
→Γ= γds. (4)

R
Using the definition of circulation, Γ = u · dl, eq. (4) gives
Z Z
γds = u · dl. (5)

To derive the circulation from eq. (5), consider a vortex sheet as sketched in fig. 10.

Figure 10: Tangential velocity jump across a vortex sheet.

Then, eq. (5) is expressed as,

Z
u · dl = u1 ds + v1 dn − u2 ds − v2 dn = (u1 − u2 )ds + (v1 − v2 )dn. (6)

Finally, if dn ≈ 0, eq. (6) gives

Z Z
γds = u · dl = (u1 − u2 )ds (7)

Since γ is constant, a relationship is derived using eq. (6) and a property of γ,

Z s2
γ ds = γ(s2 − s1 ) = γds = (u1 − u2 )ds → γ = (u1 − u2 ). (8)
s1

A combination of vortex sheet on the airfoil and uniform flow gives us a lift force as mentioned in
the Kutta-Joukowski theorem, L = ρ∞ V∞ Γ. Figure 11 shows a schmatic in regard to the combination.

Figure 11: Thin airfoil approximation.

1
Appendix A-1

9
4 The Kutta Condition
The Kutta condition is imposed to the trailing edge of the airfoil. In turn, the Kutta condition is
insufficient to explain the flow behavior, but it provides good approximations at low angles of attack.
In principal, the Kutta condition is
γ(T E) = 0. (9)

In the real physics, the Kutta condition is invalid at high angles of attack. Figure 12 shows the
flow behaviors to explain why the Kutta condition has a drawback at high angles of attack.

Figure 12: Effect of different values of circulation on the potential flow over a given airfoil at a given
angle of attack. Point 1 and 2 are stagnation points.

5 Kelvin’s Circulation Theorem and the Starting Vortex

The mathematical definition of Kelvin’s circulation theorem2 is

DΓ
= 0. (10)
Dt

Equation (10) gives a following property as

Γ1 = Γ2 , (11)

and if the airfoil is initially at rest, velocity u around the airfoil is zero, u = 0. It results that the
circulation Γ is also zero, Γ = 0. When the airfoil is set into motion, it will develop a nonzero
circulation. Figure 13 shows the mechanism of circulations around the airfoil.

Figure 13: The creation of the starting vortex and the resulting generation of circulation around the
airfoil.
2
Appendix A-2

10
 By construction of the circuits, we have

Γ2 = Γ3 + Γ4 , (12)

and by Kelvin’s theorem, Γ2 must be the same as the outer circuit’s Γ1 befor the airfoil started to
move,
Γ2 = Γ1 = 0. (13)

Therefore, the airfoil and the starting vortex must have equal and opposite circulations,

Γ3 = −Γ4 . (14)

6 Classical Thin Airfoil Theory: The Symmetric Airfoil

If an airfoil is very thin, we can assume that the aerodynamic characteristics of the thin airfoil is
similar with that of the camber line of the thin airfoil as shown in fig. 14.

Figure 14: The thin airfoil and its camber line.

Then, let us employ the vortex sheet on the camber line. On the camber line to be a streamline,
the component of velocity normal to the camber line must be zero at all points along the camber line.
Thus,
V∞,n + w0 (s) = 0, (15)

where, V∞,n and w0 (s) are the component of velocity normal to the camber line of freestream and
induced (by the vortex) velocities, respectively. Figure 15 shows the schematics of the vortex and
freestream. The geometry of fig. 16 for the component of velocity normal to the camber line of
freestream velocity in eq. (15) yields
  
−1 dz
V∞,n = V∞ sin α + tan − . (16)
dx

For the thin airfoil at small angle of attack, both α and tan−1 (−dz/dx) are small values. Then, eq. (16)
gives  
dz
V∞,n = V∞ α − . (17)
dx
In eq. (15), w0 (s) can be then defined by the infinitesimally small velocity as mentioned in eq. (3).

11
 Figure 15: Placement of the vortex sheet for thin airfoil analysis.

Figure 16: Determination of the component of freestream velocity normal to the camber line.

At first, keep in mind that if the airfoil is thin, the camber line is close to the chord line, which is
shown in (b) of fig. 15. Then, it makes the approximation that

w0 (s) ≈ w(x). (18)

Equation (18) gives a relationship between the induced velocity and a distance ξ, which are shown
in fig. 17. Here, the velocity dw induced by the vortex is derived by eq. (3) as the infinitesimally small
velocity,
γds γ(ξ)dξ
dV = → dw = − . (19)
2πr 2π(x − ξ)

12
 Equation (19) provides the velocity w(x) induced at point x by all the elemental vortices along
the chord line by integrating eq. (19),
Z c
γ(ξ)dξ
w(x) = − , (20)
0 2π(x − ξ)

where, c is a chord length.

Figure 17: Calculation of the induced velocity at the chord line.

Substituting eq. (20) and eq. (17) into eq. (15), we obtain
  Z c
0 dz γ(ξ)dξ
V∞,n + w (s) = V∞ α − − = 0. (21)
dx 0 2π(x − ξ)

Considering a symmetric airfoil, eq. (21) becomes

Z c
1 γ(ξ)dξ
= V∞ α, (22)
2π 0 x−ξ

dz
since, for the symmetric airfoil, dx = 0 in eq. (21). To deal with the integral in eq. (22), let us
transform ξ into θ via the following transformation3 ,

c
ξ = (1 − cos θ), (23)
2

then, the differential form for dξ in eq. (23) is obtained as

c
dξ = sin θdθ. (24)
2

As a same manner with eq. (23), x, where is a fixed point, in eq. (22) can be defined by a particular
value of θ, namely, θ0 , such that
c
x = (1 − cos θ0 ). (25)
2
Substituting eqs. (23) to (25) into eq. (22), we obtain
c π γ(θ) 2c sin θdθ π
Z Z Z
1 γ(ξ)dξ 1 1 γ(θ) sin θdθ
= V∞ α → c c = = V∞ α. (26)
2π 0 x−ξ 2π 0 2 (1 − cos θ0 ) − 2 (1 − cos θ) 2π 0 cos θ − cos θ0

The solution of eq. (26) is purely obtained in terms of the mathematical concerns (mathematical
3
Appendix A-3

13
theory of integral equation). Simply, the solution is

1 + cos θ
γ(θ) = 2αV∞ . (27)
sin θ

In order to derive the lift L = ρ∞ V∞ Γ (Kutta-Joukowski theorem), the total circulation Γ around
the airfoil should be given by the integral of the elemental vortex of strength γdξ (it is a same manner
with eq. (4), Z c
Γ= γ(ξ)dξ. (28)
0

Substituting eq. (25) and eq. (27) into eq. (28), we obtain
Z c Z π Z π Z π
c 1 + cos θ c
Γ= γ(ξ)dξ = γ(θ) sin θdθ = 2αV∞ sin θdθ = αcV∞ (1 + cos θ)dθ
0 0 2 0 sin θ 2 0 (29)
= παcV∞ .

The lift is
2
L = ρ∞ V∞ Γ = ρ∞ V∞ παcV∞ = παcρ∞ V∞ , (30)

then the lift coefficient is

L παcρ∞ V∞ 2
cl = = 1 2
= 2πα, (31)
q∞ c 2 ρ∞ V∞ c

and the lift slope is

dcl d(2πα)
clα = = = 2π. (32)
dα dα
Based on the same principal with the lift of eq. (30) and the lift coefficient of eq. (31), we can
obtain the moment and moment coefficient. At first, the moment at the leading edge of the airfoil is
derived by the increment of lift dL and the distance ξ,
Z c Z c Z c Z c
MLE = − ξ(dL) = − ξ(ρ∞ V∞ dΓ) = − ξ[ρ∞ V∞ γ(ξ)dξ] = −ρ∞ V∞ ξγ(ξ)dξ. (33)
0 0 0 0

Substituting eq. (23) and eqs. (28) and (29) into eq. (33), the moment is given as

1 2 2 πα
MLE = − ρ∞ V∞ c , (34)
2 2

then, the moment coefficient, as mentioned in subsection 2.3, with the lift coefficient eq. (31) is

MLE πα cl
cm,LE = 1 2
=− =− , (35)
2 ∞ V∞ Sc
ρ 2 4

where, S = c for the two dimensional airfoil. The moment coefficient about the quarter-chord point
is defined as
cl
cm,c/4 = cm,LE + , (36)
4
substituting eq. (35) into eq. (36), then we obtain

cl cl cl
cm,c/4 = cm,LE + = − + = 0. (37)
4 4 4

Figure 18 shows comparisons between the theory and experiment for the lift and moment coeffi-
cients for an NACA0012 airfoil.

14
Figure 18: Comparison between theory and experiment for the lift and moment coefficients for an
NACA0012 airfoil. (Experimental data are from Abbott and von Doenhoff, “Theory of Wing Section”)

7 The Cambered Airfoil

The cambered airfoil is started from the general solutions of the thin airfoil as mentioned in eq. (21),
Z c  
γ(ξ)dξ dz
= V∞ α − . (38)
0 2π(x − ξ) dx

dz dz
In case of a symmetric airfoil, dx = 0, but for the cambered airfoil dx 6= 0. For the left side of
eq. (38), we can employ eq. (26), thus eq. (38) is replaced to
Z π  
1 γ(θ) sin θdθ dz
= V∞ α − . (39)
2π 0 cos θ − cos θ0 dx

The general solutions of eq. (39) is purely placed on the mathematics area. Just we employ the
solution of eq. (39) as
∞
!
1 + cos θ X
γ(θ) = 2V∞ A0 + An sin nθ . (40)
sin θ
n=1

Let us now obtain expressions for the aerodynamic coefficients for a cambered airfoil. Using eq. (24)
and eq. (28), the total circulation due to the entire vortex sheet from the leading edge to the trailing
edge is Z c Z π
c
Γ= γ(ξ)dξ = γ(θ) sin θdθ. (41)
0 2 0

15
Substituting eq. (40) into eq. (41), we obtain

∞
" #
Z π X Z π
Γ = cV∞ A0 (1 + cos θ)dθ + An sin nθ sin θdθ . (42)
0 n=1 0

From any standard table of integrals,


Z π Z π π/2 for n = 1
(1 + cos θ)dθ = π, and sin nθ sin θdθ = (43)
0 0 0 otherwise.

Hence, eq. (42) becomes  π 

Γ = cV∞ πA0 + A1 . (44)
2
Substituting eq. (44) into the Kutta-Joukowski theorem (L = ρ∞ V∞ Γ), the lift per unit span is

2
 π 
L = ρ∞ V∞ Γ = ρ∞ V∞ c πA0 + A1 , (45)
2

and the lift coefficient is

2 c πA + π A


L ρ∞ V∞ 0 2 1
cl = 1 2
= 1 2
= π(2A0 + A1). (46)
2 ρ∞ V∞ S 2 ρ∞ V∞ S

For the lift coefficient in eq. (46), we must define the coefficients of A0 and A1 . The coefficients
are derived in section 7.1, and the solutions of A0 and An are just used in the following equations.
Z π Z π
1 dz 2 dz
A0 = α − dθ0 , and An = cos nθ0 dθ0 . (47)
π 0 dx π 0 dx
Substituting A0 and A1 for n = 1 of eq. (47) into eq. (46), we obtain
Z π
2 π dz
  Z 
1 dz
cl = π(2A0 + A1) = π 2 α − dθ0 + cos θ0 dθ0
π 0 dx π 0 dx
(48)
1 π dz
 Z 
= 2π α + (cos θ0 − 1)dθ0 ,
π 0 dx

and the lift slope is

dcl
clα = = 2π. (49)
dα
Equation (49) is very interesting because the lift slope for cambered airfoils and symmetrical airfoils
are exactly same. Also, it is a general result from thin airfoil theory that the lift slope is constant
for any shape airfoil. However, the lift coefficients between the cambered and symmetrical airfoils are
different due to the fact that eq. (48) has an integral formula for the geometry of the camber. Next,
we will derive the angle of zero lift, αL=0 , as only difference between the cambered and symmetrical
airfoils, because of αL=0 = 0 for the symmetrical airfoil. From eq. (49),

dcl = clα dα → cl − cl0 = 2π(α − αL=0 ) → cl = 2π(α − αL=0 ), (50)

where, cl0 = cl (αL=0 ) = 0 as the zero lift. Then, substituting eq. (50) into eq. (48), we obtain
Z π
1 dz
αL=0 = − (cos θ0 − 1)dθ0 . (51)
π 0 dx

16
 For the moment of the cambered airfoil, eq. (33) can be used, hence the resulting coefficient of
moment of the cambered airfoil is
 
π A2
cm,LE =− A0 + A1 − , (52)
2 2

then, substituting eq. (46) into eq. (52), we have

hc π i
l
cm,LE = − + (A1 − A2 ) . (53)
4 4

At the quard point of the chord, the coefficient of the pitch moment is derived using eq. (36),

cl hc
l π i c
l π
cm,c/4 = cm,LE + =− + (A1 − A2 ) + = (A2 − A1 ) . (54)
4 4 4 4 4

Unlike the symmetric airfoil, where cm,c/4 = 0, eq. (54) demostrates that cm,c/4 is constant (inde-
pendent of angle of attack) because A1 and A2 are a function of the geometry at the given airfoil shape.
Thus, the quarter-chord point is the theoretical location of the aerodynamic center for a cambered
airfoil. Next, the center of pressure can be obtained. Prior to obtain the center of pressure, note that
the definition of the center of pressure is the average location of the pressure acting on the surface, as
shown in fig. 19.

Figure 19: Definition of the center of pressure [NASA web page]

Mathematically, it is defined as
R
xp(x)dx MLE cm,LE · c
xcp = R =− =− . (55)
p(x)dx L cl

17
Substituting eq. (53) into eq. (55), we obtain
 
c π
xcp = 1 + (A1 − A2 ) . (56)
4 cl

7.1 Coefficients, A0 and An

In eq. (40), the coefficients, A0 and An must be defined to close the system of eq. (40). To find these
specific values of coefficients, eq. (40) is substituted into eq. (39), then we have
P∞
2V∞ A0 1+cos θ
π π

 
sin θ + n=1 An sin nθ sin θdθ
Z Z
1 γ(θ) sin θdθ 1 dz
= = V∞ α − . (57)
2π 0 cos θ − cos θ0 2π 0 cos θ − cos θ0 dx

Rearranging eq. (57), eq. (57) gives

π ∞ Z
1 X π An sin nθ sin θdθ
Z
1 A0 (1 + cos θ)dθ dz
+ =α− . (58)
π 0 cos θ − cos θ0 π 0 cos θ − cos θ0 dx
n=1

The left and right sides of eq. (58) is derived in appendix A-44 . The integral formula is purely
included in the mathematic area. Here, we just employ the solutions of the integral formula. The left
side of eq. (58) is
π  Z π
1 π
Z Z 
1 A0 (1 + cos θ)dθ 1 dθ cos θdθ 1
= A0 + = A0 (0 + π) = A0 , (59)
π 0 cos θ − cos θ0 π 0 cos θ − cos θ0 π 0 cos θ − cos θ0 π

and the right side of eq. (58) is

Z π
sin nθ sin θdθ
= −π cos nθ0 . (60)
0 cos θ − cos θ0

Substituting eqs. (59) and (60) into eq. (58), eq. (58) can be reduced to
∞
X dz
A0 − An cos nθ0 = α −
dx
n=1
∞ (61)
dz X
or = (α − A0 ) + An cos nθ0
dx
n=1

Just from eq. (61), the coefficients of A0 and An may not be obtained, hence we employ the Fourier
cosine series expansion for the function of dz/dx. The general form of the Fourier cosine series is given
by
∞
X
f (θ) = B0 + Bn cos nθ, (62)
n=1

where, the coefficients B0 and Bn are given by

Z π Z π
1 2
B0 = f (θ)dθ and Bn = f (θ) cos nθdθ. (63)
π 0 π 0

dz
Once we compare eq. (61) and eq. (62), assuming f (θ) = dx , (α − A0 ) and the cosine function of
eq. (61) can be treated as B0 and the cosine function of eq. (62), respectively. Then, the coefficients
4
Appendix A-4

18
of A0 and An are given using the eq. (63) as

1 π 1 π dz 1 π dz
Z Z Z
α − A0 = B 0 = f (θ)dθ = dθ0 → A0 = α − dθ0 ,
π 0 π 0 dx π 0 dx
(64)
2 π 2 π dz
Z Z
An = B n = f (θ) cos nθdθ = cos nθ0 dθ0 .
π 0 π 0 dx

7.2 Example
Consider an NACA23012 airfoil. The mean camber line for this airfoil is given by
   x 
z x 3  x 2 x
= 2.6595 − 0.6075 + 0.1147 for 0 ≤ ≤ 0.2025
c c c c c (65)
z  x x
= 0.02208 1 − for 0.2025 ≤ ≤ 1.0.
c c c

Calculate (a) the angle of attack at zero lift, (b) the lift coefficient when α = 4◦ , (c) the moment
coefficient about the quarter chord, and (d) the location of the center of pressure in terms of xcp /c,
when α = 4◦ . Compare the results with experimental data.

Solution We need dz/dx for the solutions of all questions. From the given shape of the mean
camber line, this is
   
dz x 2 x x
= 2.6595 3 − 1.215 + 0.1147 for 0 ≤ ≤ 0.2025
dx c c c
(66)
dz x
= −0.02208 for 0.2025 ≤ ≤ 1.0.
dx c

Transforming from x to θ, where x = 2c (1 − cos θ) in eq. (25), we have

 
dz 3 2
= 2.6595 (1 − 2 cos θ + cos θ) − 0.6075(1 − cos θ) + 0.1147
dx 4
(67)
= 0.6840 − 2.3736 cos θ + 1.995 cos2 θ for 0 ≤ θ ≤ 0.9335 rad
= −0.02208 for 0.9335 ≤ θ ≤ π.

Note: For simplicity, we have dropped the subscript zero from θ. θ0 is the variable of integration
- it can just as well be symbolized as θ for the variable of integration.

(a) From eq. (51), Z π

1 dz
αL=0 = − (cos θ − 1)dθ. (68)
π 0 dx
Substituting eq. (67) into eq. (68), we have

1 0.9335
Z
αL=0 = − (−0.6840 + 0.30576 cos θ − 4.3686 cos2 θ + 1.995 cos3 θ)dθ
π 0
(69)
1 π
Z
− (0.02208 − 0.02208 cos θ)dθ.
π 0.9335

From a table of integrals, eq. (69) becomes

19
 1
αL=0 = − [−2.8683θ + 3.0576 sin θ − 2.1843 sin θ cos θ + 0.665 sin θ(cos2 θ + 2)]0.9335
0
π
1
− [0.02208θ − 0.02208 sin θ]π0.9335
π (70)
1
= − (−0.0065 + 0.0665) = −0.0191rad
π
= − 1.09◦ .

(b) α = 4◦ = 0.0698rad. From eq. (50),

cl = 2π(α − αL=0 ). (71)

Substituting eq. (70) into eq. (71), we obtain

cl = 2π(α − αL=0 ) = 2π(0.0698 + 0.0191) = 0.559. (72)

(c) The value of cm,c/4 is obtained from eq. (54). For this, we need the two Fourier coefficients,
A1 and A2 . From eq. (64) Z π
2 dz
An = cos nθdθ. (73)
π 0 dx
For A1 , substituting eq. (67) into eq. (73) gives

π
2 0.9335
Z Z
2 dz
A1 = cos θdθ = (0.6840 cos θ − 2.3736 cos2 θ + 1.995 cos3 θ)dθ
π 0 dx π 0
2 π
Z
+ (−0.02208 cos θ)dθ (74)
π 0.9335
= 0.0954

For A2 , substituting eq. (67) into eq. (73) gives

2 π dz
Z
2 dz
A2 = cos 2θdθ = (2 cos2 θ − 1)dθ
π 0 dx π dx
2 0.9335
Z
= (−0.6840 − 2.3736 cos θ − 0.627 cos2 θ − 4.747 cos3 θ + 3.99 cos4 θ)dθ
π 0 (75)
2 π
Z
+ (0.02208 − 0.0446 cos2 θ)dθ
π 0.9335
= 0.0792

Finally, The value of cm,c/4 in eq. (54) is obtained using eqs. (74) and (75),

π π
cm,c/4 = (A2 − A1 ) = (0.0792 − 0.0954) = −0.0127. (76)
4 4

(d) From eq. (56) with eqs. (72), (74) and (75),
 
c π ch π i
xcp = 1 + (A1 − A2 ) = 1+ (0.0954 − 0.0792) , (77)
4 cl 4 0.559

20
then, xcp /c is
xcp
= 0.273. (78)
c

Comparisons with experimental data

Table 2: Comparisons of analytical solutions and experimental data for the NACA23012 airfoil
Calculated Experiment
αL=0 -1.09 -1.1
cl (at α = 4◦ ) 0.559 0.55
cm,c/4 -0.0127 -0.01

8 Lifting Flows over Arbitrary Bodies: The Vortex Panel Numerical

Method
Our talks for the thin airfoil theory show meaningful results in cases that the thickness of the airfoil is
less than 12 percent and the angles of attack are small. Hence, thin airfoil theory is quite restrictive
when we consider the whole spectrum of the aerodynamic application. To overcome the limitation of
the theory mentioned in previous sections, we treat the vortex panel method, which is a numerical
technique that has come into widespread use since the early 1970s. In turn, there exists no closed-form
solution for γ. Hence, the solution of γ must be obtained numerically. This is the purpose of the vortex
panel method.
For the vortex panel method, we start from the velocity potential. The velocity potential for the
vortex flow at a specific point P is
Γ
φ=− θ. (79)
2π
For the increment in velocity potential, dφ, induced at point P by the elemental vortex γds is

γds
dφ = − θ, (80)
2π

where, the details of eq. (80) can be found in the appendix A-55 . Then the velocity potential at P
due to the entire vortex sheet is Z
1
φ(x, z) = − θγds. (81)
2π
Let us approximate the vortex sheet by a series of straight panels, as shown in fig. 20. Let the
vortex strength γ(s) per unit length be constant over a given panel, but allow it to vary from one
panel to the next. These panel strengths are unknowns; main thrust of the panel technique is to solve
for γj j = 1ton, such that the body surface becomes a streamline of the flow and such that the Kutta
condition γT E = 0 is satisfied. The midpoint of each panel is a control point at which the boundary
condition is applied. Also, at each control point, the normal component of the flow velocity is zero.
Let P be a point located at (x, y) in the flow, and let rpj be the distance from any point on the
jth panel to P , as shown in fig. 20. The distance rpj makes the angle θpj with respect to the x axis.
5
Appendix A-5

21
 Figure 20: Panel distributions over the surface of a body of arbitrary shape.

From eq. (81), the velocity potential induced at P due to the jth panel, ∆φj , is
Z
1
φj = − θpj γj dsj . (82)
2π j

In eq. (82), γj is constant over the jth panel, and the integral is taken over the jth panel only.
The angle θpj is given by
y − yj
θpj = tan−1 − , (83)
x − xj
then, the potential at P due to all the panels is summed over all the panels using eq. (81)
n n Z
X X γj
φ(p) = φj = − θpj dsj . (84)
2π j
j=1 j=1

Since point P is just an artitrary point in the flow, let us put P at the control point of the ith panel
shown in fig. 20. The coordinates of this control point are (xi , yi ). Then eqs. (83) and (84) become

yi − yj
θij = tan−1 , (85)
xi − xj

then, θpj in eqs. (83) and (84)is replaced to θij of eq. (85). Thus, eq. (84) becomes

n n Z
X X γj
φ(xi , yi ) = φj = − θij dsj . (86)
2π j
j=1 j=1

At the control points, the normal component of the velocity is zero (V∞,n +w0 (s) = 0), as mentioned
in eq. (15) of section 6. V∞,n of the ith panel is given by

V∞,n = V∞ cos βi , (87)

and the details is described in the appendix A-66 . The normal component of velocity induced at
6
Appendix A-6

22
(xi , yi ) by the vortex panels is
∂
w0 (s) = Vn = [φ(xi , yi )]. (88)
∂ni
∂φ ∂φ
In turn, eq. (88) was mentioned in our previous lecture, u = ∂x and v = ∂y . If the coordinates are
∂φ ∂φ
transformed to the normal and parallel, Vn = ∂n and Vp = ∂p . Combining eqs. (86) and (88), we have

n Z
X γj ∂θij
Vn = − dsj . (89)
2π j ∂ni
j=1

Substituting eqs. (87) and (89) into eq. (15), we obtain

n Z
X γj ∂θij
V∞ cos βi − dsj = 0. (90)
2π j ∂ni
j=1

Equation (90) is allowable to obtain the vortex strength, γ, since eq. (90) is purely a function
of geometry and inflow. Once we try to solve eq. (90) numerically, the Kutta condition should be
imposed on the trailing edge, γ(T E) = 0, to meet the number of equations and unknowns, which is
happened by the inconsistency between the grid points and control point on the surface.
Once we obtain vortex element γj sj , the total circulation is

n
X
Γ= γj s j . (91)
j=1

Lastly, the lift using eq. (91) is obtained from

n
X
L = ρ∞ V∞ γj sj . (92)
j=1

23