NACA 4 digit series (MPXX) are formulated as follows. The first digit $M$ represents the maximum camber in hundreds of $c$. The second digit $P$ specifies the location of $m$ in tenths of $c$. The third and fourth digit ($XX$) is the maximum thickness of the airfoil in percentage of $c$. The formulation of the upper and lower surface coordinates are as follows
$$x_{u}=x - y_{t}\sin(\theta)$$
$$x_{l}=x + y_{t}\sin(\theta)$$
$$y_{u}=y_{c} + y_{t}\cos(\theta)$$
$$y_{l}=y_{c} - y_{t}\cos(\theta)$$
In here $\theta=\arctan(dy_{c}/dx)$, $y_{c}$ is the coordinates of the mean camber line and calculated as follows
$$y_{c}=\frac{c*m}{p^2}\left[2p\left(\frac{x}{c}\right)-\left(\frac{x}{c}\right)^2\right] \hspace{10pt} \textrm{for}\hspace{5pt} \frac{x}{c}\lt p$$
$$\frac{dy_{c}}{dx}=\frac{2m}{p^2}\left[p-\frac{x}{c}\right] \hspace{10pt} \textrm{for} \hspace{5pt}\frac{x}{c}\lt p$$
$$y_{c}=\frac{c*m}{(1-p)^2}\left[(1-2p)+2p\left(\frac{x}{c}\right)-\left(\frac{x}{c}\right)^2\right] \hspace{10pt} \textrm{for}\hspace{5pt} \frac{x}{c}\ge p$$
$$\frac{dy_{c}}{dx}=\frac{2m}{(1-p)^2}\left[p-\frac{x}{c}\right] \hspace{10pt} \textrm{for} \hspace{5pt}\frac{x}{c}\ge p$$
$y_{t}$ is the thickness distribution and the formula reads
$$y_{t} = (t*c)*(1.4845\sqrt{x}-0.63x-1.785x^2+1.4215x^3-0.5075x^4)$$
Here is an example python code for generating NACA 4 Digit airfoils.
import numpy as np
import matplotlib.pyplot as plot
import math as mt
# NACA 4 Digit Airfoil Generation
m = 2 # maximum camber
p = 4 # maximum camber position
t = 12 # thichkness
c = 2 # chord length
m = 0.01*m # maximum camber in % of chord
p = 0.10*p # maximum camber position in tenths of chord
t = 0.01*t # thickness in % of chord
# Coefficients for 4 digit series
a0 = 1.4845
a1 = -0.6300
a2 = -1.7580
a3 = 1.4215
a4 = -0.5075
n = 100 # number of points along the chord
x = np.linspace(0,c,n) # x coordinate of points along the chord
y = np.zeros(n) # x coordinate of points along the chord
yc = np.zeros(n) # y coordinate of the camber line
dyc = np.zeros(n) # gradient of the camber line
yt = np.zeros(n) # thickness distribution
xu = np.zeros(n) # x coordinate of the upper surface
yu = np.zeros(n) # y coordinate of the upper surface
xl = np.zeros(n) # x coordinate of the lower surface
yl = np.zeros(n) # y coordinate of the lower surface
for i in range(n):
if (x[i]/c < p):
yc[i] = (c*m/p**2)*(2*p*(x[i]/c)-(x[i]/c)**2)
dyc[i] = ((2*m)/p**2)*(p-(x[i]/c))
else:
yc[i] = (c*m/(1-p)**2)*(1-2*p+2*p*(x[i]/c)-(x[i]/c)**2)
dyc[i] = ((2*m)/(1-p)**2)*(p-(x[i]/c))
for i in range(n):
yt[i] = (t*c)*(a0*mt.sqrt(x[i]/c)+a1*(x[i]/c)+a2*(x[i]/c)**2+a3*(x[i]/c)**3+a4*(x[i]/c)**4)
teta = mt.atan(dyc[i])
xu[i] = x[i] - yt[i]*mt.sin(teta)
xl[i] = x[i] + yt[i]*mt.sin(teta)
yu[i] = yc[i] + yt[i]*mt.cos(teta)
yl[i] = yc[i] - yt[i]*mt.cos(teta)
plot.xlim(-0.2,c+0.2)
plot.ylim(-c/3,c/3)
plot.plot(xu,yu,color='deepskyblue')
plot.plot(xl,yl,color='deepskyblue')
plot.plot(x,yc,'g--')
plot.yticks([])
plot.xticks([])
plot.show()
<< 2.4 Conservation of Energy || Contents || 3.1.2 NACA 5 Digit Airfoils >>
[1] Summary of Airfoil Data, I.H. Abbott, A.E. von Doenhoff, L.S. Stevers, Jr, NACA Report No. 824, 1945.
[2] www.airfoiltools.com