NACA四位数翼型公式及代码

NACA四位数翼型的型面曲线本身并不高深，只是一些多项式曲线而已。它是依据理想流动理论和风洞吹风数据，通过对一系列性能好的翼型拟合得到的。用户可以自定义最大弯度（m）、最大弯度位置（p）和最大厚度（t）3个参数，生成相应的翼型。

四个数字的意义举例如下：

从上面的定义可以看到NACA四位数翼型的两个局限：

1. 最大弯度只能是0%，1%，2%......9%之一；

2. 最大弯度位置只能是0%，10%，20%......90%之一；

比如，想要生成最大弯度在35%处的翼型就做不到。当然这影响不大，因为一般情况下气流对翼型的最大弯度位置并没有那么敏感。

只有这三个参数显然还无法生成翼型，还需要知道厚度和弯度的分布规律才行。这两个规律是给定的，都是多项式曲线，对于所有NACA四位数翼型都相同。分别如下：

弯度分布用中心线的y坐标表示：

厚度分布用半厚度y_t表示：

有了弯度和厚度分布，就可以给出翼型上下表面的曲线了：

最后多说一句，还有一种常用的低速翼型厚度分布规律，即C4翼型，和NACA翼型不一样，但低速时的性能相差无几，C4翼型的厚度分布如下：

大家如果感兴趣可以自己把NACA和C4的翼型厚度分布画在一张图上，看看它们的区别。总体来说C4的前部要比NACA胖一点。

% NACA four digit airfoil generation code.

% Input four digit from command line, eg. "2412"

clear

s=input('\nInput four digit mumber of the NACA airfoil:\n?','s');

t=str2num(s(length(s)-1))/10 str2num(s(length(s)))/100;    % thickness

p=str2num(s(length(s)-2))/10;     % max camber position

m=str2num(s(length(s)-3))/100;    % max camber

if p==0                     % To avoid p=0

    p=0.0001;

end

x=0:0.001:1;

for i=1:length(x)

    yt(i)=(t/0.2)*(0.2969*x(i)^0.5-0.1260*x(i)-0.3516*x(i)^2 0.2843*x(i)^3-0.1015*x(i)^4);

        % thickness distribution

    if x(i)<=p

        yc(i)=m/p^2*(2*p*x(i)-x(i)^2);  % center line (front half)

    else

        yc(i)=m/(1-p)^2*((1-2*p) 2*p*x(i)-x(i)^2);  % center line (rear half)

    end

    if i==1

        theta=pi/2;

    else

        theta=atan((yc(i)-yc(i-1))/(x(i)-x(i-1)));

    end

    xu(i)=x(i)-yt(i)*sin(theta);        % upper surface

    yu(i)=yc(i) yt(i)*cos(theta);

    xl(i)=x(i) yt(i)*sin(theta);        % lower surface

    yl(i)=yc(i)-yt(i)*cos(theta);

end

% plot airfoil

 plot(xu,yu,'b-',xl,yl,'b-',x,yc,'k-.');

 grid on;

 axis equal;

% output airfoil data

 filename=['NACA',s,'.dat'];

 fid=fopen(filename,'w');

 fprintf(fid,'   x_u      y_u      x_l      y_l\n');

 for i=1:length(x)

     fprintf(fid,'%8.4f %8.4f %8.4f %8.4f\n',xu(i),yu(i),xl(i),yl(i));

 end

 fclose(fid);

运行后，输入四位数字：

还可以试试其它的，比如

0015 ▼

4412 ▼

9401 ▼

0000 ▼

6666 ▼

9999 ▼

9110 ▼

哈哈，后面这些怪异的翼型当然不能用。