Universal Equations for Continuously Variable

4-Digit Symmetric NACA Airfoils

Tianhou Wang
Department of Engineering, University of Cambridge
August 24, 2018

Contents
1 Background 1

2 Source Definition 1

3 Determination of the Coefficients 2

4 The Final Equations 4

5 Derivation of the Coefficients Formulae 4

5.1 Properties of the Definition . . . . . . . . . . . . . . . . . . . 4
5.2 First Condition of the Equations . . . . . . . . . . . . . . . . 5
5.3 Coefficients d2 and d3 . . . . . . . . . . . . . . . . . . . . . . . 5
5.4 Coefficient d1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
5.5 Coefficient a0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
5.6 Second Condition of the Equations . . . . . . . . . . . . . . . 7

6 Defect of these Formulae 8

1
1 Background
The NACA family of airfoil geometry is usually conveniently generated
through an airfoil plotter available online, which generates a set of data
points through which a smooth curve could be fitted in a CAD programme
or alike. The curve may then be used for further purposes. This approach is
widely used in less demanding modelling scenarios.

However, in parametric modelling where much freedom to modify the

model through the control of parameters is desired, it is required that the
geometry of the airfoil can be varied parametrically within the CAD model
file, which may then be passed onto other software programmes such as FEA
packages or optimisers without further external reference to online resources.
The ability to define the airfoil geometry mathematically within the CAD
programme also provides improved robustness for the model when other pa-
rameters are being varied, without potential pitfalls such as unrestrained
splines assuming reversed curvature, or small gaps originating due to a lack
of fit, giving rise to errors.

Most CAD packages do not allow the complicated algorithm of an air-

foil plotter to be implemented. Instead, they demand explicit equations by
which the curve may be drawn. A good example of this would be SOLID-
WORKS. CREO Parametric, as a more advanced software package, does
allow the implementation of simple mathematical operations, but still, they
lack the ability to solve for the numerical coefficients required to plot an
airfoil. Therefore, it is necessary to have a set of formulae to enable the
determination of airfoil geometry by the substitution of desired properties
(e.g. maximum thickness location) only.

This article aims to present a set of explicit parametric formulae which gov-
erns the shape of a family of NACA symmetric airfoil variants. By explicit
parametric, it is meant that the equation can be determined by a straight-
forward substitution of the control parameters.

2 Source Definition
According to NACA report No.492[1], the NACA 4-digit series airfoil may
be modified in shape by changing the y equation into a piecewise function
dividing at the point where the thickness of the airfoil reaches a maximum.

1
 x = t for 0 ≤ t ≤ 1 (1)
( √
a0 t + a1 t + a2 t2 + a3 t3 for 0 ≤ t ≤ m
y= (2)
d0 + d1 (1 − t) + d2 (1 − t)2 + d3 (1 − t)3 for m < t ≤ 1

where m is a parameter that controls the fractional position of the max-

imum thickness point along the chord, and t which varies from 0 to 1, is a
parameter for both the x and y equations.

In the appendix of the NACA report, there is a table of the eight coef-
ficients to be used for several discrete shapes of the airfoil. This, however,
is not suitable for parametric modelling and optimising purposes, as it is
required that the shape can vary continuously. Besides, it is not possible
to implement a switch function in most CAD packages to allow for the use
of discrete coefficients. Therefore, the following section presents the formu-
lae to determine these coefficients from control parameters, with the rest of
this article presenting the mathematical basis on which these formulae are
derived.

3 Determination of the Coefficients

The following definitions are used for the rest of this article:

C chord length of the airfoil section

m fractional position of the thickest point along the chord
T fractional thickness of the airfoil
D fractional thickness of the trailing edge

It should be noted that C is the master scaling factor of the geometry,

with all the other parameters being fractional.

Coefficient d0
1
d0 = D (3)
2

Coefficient d1

d1 = −2.5m4 + 7.1667m3 − 2.725m2 + 0.5033m + 0.155 (4)

2
Coefficient d3
−0.2 + (1 − m)d1 + 2d0
d3 = (5)
(1 − m)3

Coefficient d2
−d1 − 3d3 (1 − m)2
d2 = (6)
2(1 − m)

Intermediate coefficient R
(1 − m)2
R= (7)
2d1 (1 − m) − 0.6 + 6d0

Coefficient a0
a0 = 0.2969 for normal nose radius (8)

Intermediate coefficient β
1 a0
β= + √ 3 (9)
R 4( m)

Coefficient a3
√
0.1 − 21 a0 m + 12 βm2
a3 = (10)
m3

Coefficient a2
1
a2 = β − 3ma3 (11)
2

Coefficient a1
a0
a1 = − √ + 3a3 m2 − βm (12)
2 m

By substituting the above coefficients into Equation (1) and 2, it is possible

to obtain a curve for the upper half of the airfoil of unity chord length. In
order to obtain the final geometry:
• scale both expressions with the factor C
T
• scale the y expression with the factor 0.2

• connect the trailing edge with a straight line

3
4 The Final Equations
With the coefficients determined as described above, and applying the
appropriate scaling, the form equations of the modified airfoil may be written
in the following form:

x = Ct for 0 ≤ t ≤ 1 (13)
(
CT
√
0.2
[a0 t + a1 t + a2 t2 + a3 t3 ] for 0 ≤ t ≤ m
±y = CT (14)
[d + d1 (1 − t) + d2 (1 − t)2 + d3 (1 − t)3 ]
0.2 0
for m < t ≤ 1

Alternatively, in the explicit form:

(
CT
[a Cx + a1 Cx + a2 ( Cx )2 + a3 ( Cx )3 ]
p
0.2 0
for 0 ≤ x ≤ mC
±y = CT
[d + d1 (1 − Cx ) + d2 (1 − Cx )2 + d3 (1 − Cx )3 ] for mC < x ≤ C
0.2 0
(15)

5 Derivation of the Coefficients Formulae

This section presents the derivation of the coefficients formulae used above.

Please take note that x = t (from Equation (1)) in the following discus-
sions.

5.1 Properties of the Definition

NACA report No.492 has demanded that the family of the modified sym-
metric 4-digit airfoils be described in the piecewise form by two third-order
polynomials (Equation (2)). It can be immediately seen that this elegant
choice of equation form always satisfies the following properties:

1. the curve starts from the origin.

2. the curve terminates at x = 1, y = d0 , so the trailing edge thickness
can be controlled by the coefficient d0 alone.
3. the curve has an infinite gradient at x = 0, i.e.
dy
lim =∞
x→0 dx
this guarantees that the curve has a round nose unless the coefficient
a0 is set to zero, which would result in a sharp leading edge.

4
5.2 First Condition of the Equations
Because the form equations are piecewise, in order for the two y equations
to match at the point t = m, it is necessary to enforce the following:

1. y = 0.1 at t = m, i.e. the curve reaches a specified value of thickness

at the point of x = m. Given that the thickness, when scaled, is to
be done globally, i.e. as a multiplier in front of the y equations, it
is convenient to choose a constant 0.1 here while the scaling can be
applied afterwards. This is the reason why the scaling factor as seen in
Equation (14) and (15) has a denominator of 0.2.

2. dy


dt x=m
= 0, i.e. the curve reaches a maximum at the point of x = m.

Both y equations must satisfy these two conditions.

5.3 Coefficients d2 and d3

Consider the second half of the y equation:

y = d0 + d1 (1 − t) + d2 (1 − t)2 + d3 (1 − t)3
dy
= −d1 − 2d2 (1 − t) − 3d3 (1 − t)2
dt

By enforcing the condition as described in Section 5.2, the following simul-

taneous equations are arrived at:
(
d0 + d1 (1 − m) + d2 (1 − m)2 + d3 (1 − m)3 = 0.1
−d1 − 2d2 (1 − m) − 3d3 (1 − t)2 =0

where d0 = D2 is specified. It is then possible to arrive at Equation (5) and

(6), where d2 and d3 are calculated in terms of d0 and d1 .

5.4 Coefficient d1
NACA report No.492 has chosen 5 discrete values for coefficient d1 for
reference. The value of d1 represent the angle of the curve at the trailing
edge, as the following equation makes clear:

dy 
= −d1
dx x=1

5
 Value of d1 must be chosen so that a reversal of curvature is avoided in the
after-portion of the curve. By performing a polynomial fit to the data chosen
in NACA report No.492, Equation (4) is arrived at. The following table
presents the values of d1 chosen in NACA report No.492, and the respective
values as calculated from Equation (4).

Table 1: Value of Coefficient d1

m d1 from NACA d1 from Eqn.(4)
0.2 0.200 0.1999936
0.3 0.234 0.2339909
0.4 0.315 0.3149888
0.5 0.465 0.4649875
0.6 0.700 0.6999872

It can be seen that Equation (4) is very precise for the range 0.2 ≤ m ≤ 0.6.
However, this method is not robust, and the use of this polynomial fitting
method is a known defect of the equations presented in this article. Please
refer to the next section for further discussion on this issue.

5.5 Coefficient a0
Please take note that, in Equation (8), the value of coefficient a0 has been
specified. This is no random real number. The coefficient a0 is related to the
radius of the nose of the airfoil.

The formula for radius of curvature of a curve in Cartesian coordinates is:

3
(1 + y 0 2 ) 2
Radius of Curvature =
y 00
2
dy 00 d y
where y 0 = ,y = 2
dx dx

In this case, it can be determined that:

a0
y 0 = √ + a1 + 2a2 x + 3a3 x2
2 x
a0
y 00 = − √ 3 + 2a2 + 6a3 x
4( x)

6
 Substituting into the radius of curvature formula gives:

a0 2
lim R =
x→0 2
Therefore, it is further possible to define a fifth control parameter to control
the nose radius, from which a0 may be determined, and the formulae pre-
sented above will work for a range of values for a0 . However, this is out of
the scope of this article. It must be noted that a0 shall not take a negative
value.

5.6 Second Condition of the Equations

With the coefficients of the second half of the curve determined, it is now
necessary to examine the four coefficients of the first half of the curve. The
same approach as discussed above may be used, where the First Condition
would produce a set of two simultaneous equations.
( √
y = a0 m + a1 m + a2 m3 + a3 m3 = 0.1
dy 
 a 2
dt t=m
= 2√m + a1 + 2a2 m + 3a3 m = 0
0

With a0 determined as described above, there must still be one more con-
dition before these four coefficients can be solved for.

The Second Condition must be enforced: radius of curvature must match

exactly at t = m for the two y equations.

Radius of curvature of the curve at x = m, calculated from the second

half, is given by Equation (7), which is derived from the definition, taking
note of the following:

dy 
= −d1 − 2d2 (1 − m) − 3d3 (1 − m)2
dt t=m
d2 y 

= 2d2 + 6d3 (1 − m)
dt2 t=m

where d2 and d3 should be expressed in terms of d1 and d0 by Equation (5)

and (6). Note that, since all four d coefficients are now determined, Equation
(7) evaluates into a constant R0 .

7
 Radius of curvature of the curve at x = m, calculated from the first half,
is given by:
1
R= = R0
− 4(√a0x)3 + 2a2 + 6a3 x

and this will give the third simultaneous equation required:

a0 1
− √ 3 + 2a2 + 6a3 x =
4( x) R0

It is now possible to solve for a1 , a2 , and a3 , and the result is given by

Equation (10), (11), and (12).

6 Defect of these Formulae

The method and results presented above contain a defect, which is related
to the determination of the coefficient d1 . It is unclear how was the values
presented in NACA report No.492 chosen, and the only condition that d1
must satisfy is that a reversal of curvature in the second half of the curve
must not be allowed.

If the curvature is not to reverse, the sign of the radius of curvature,

evaluated at any point between x = m and x = 1, must be the same as the
sign of R0 used in the previous section, which is negative. Therefore, this
condition can be written as follows:
3
(1 + y 0 2 ) 2
< 0 for m ≤ x ≤ 1
y 00
where y 0 = −d1 − 2d2 (1 − x) − 3d3 (1 − x)2 , y 00 = 2d2 + 6d3 (1 − x)

The condition above is very hard to be solved for analytically, but it should
yield an inequality, which the value of d1 must satisfy. It should be noted that
d2 and d3 is present in the inequality, which is related to d0 . In NACA report
No.492, the value of d1 is only a function of m, and for all airfoils reported,
d0 is kept at a constant value. Therefore, it is unclear if the variation of d0
will make d1 , which, under the current framework, being only a function of
m, violate this condition at some point. Hence, the variation of d1 must be
applied only around the value of 0.002 (which is the value reported in NACA
report No.492) in a very cautious manner.

8
 Further work should be done to numerically study the extent to which d0
can be allowed to vary.

References
[1] John Stack and Albert E. von Doenhoff. Tests of 16 related airfoils at high
speed. Available at https://ntrs.nasa.gov/search.jsp?R=19930091566