4位NACA翼型：
主要参考：
1. 知乎评论区：https://zhuanlan.zhihu.com/p/621483821
2. 网友个人博客：https://blog.hantaotao.top/archives/32/
% 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)
% 尾缘有厚度，最后一个系数是-0.1015；尾缘是尖的，要把它改成-0.1036
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);
运行后，直接输入0012、2412等即可
NACA 5-digit 五位翼型
主要参考：
1.MATLAB版本：https://ww2.mathworks.cn/matlabcentral/fileexchange/23241-naca-5-digit-airfoil-generator/files/naca5gen.m
2.MATLAB版本：https://github.com/bvb-06/NACA-Airfoil-generator/blob/master/NACA.m 实现了reflex的功能
3.Python版本：https://github.com/dgorissen/naca/blob/master/naca.py
4.文字资料 https://web.stanford.edu/~cantwell/AA200_Course_Material/The%20NACA%20airfoil%20series.pdf
5. http://airfoiltools.com/airfoil/naca5digit
我把reflex==0和reflex==1的版本分开了，下面是reflex==0:
% NACA 5 digit airfoil generation code.
% Input 5 digit from command line, eg. "23012"
clear
s=input('\nInput four digit mumber of the NACA airfoil:\n?','s');
Cl=str2num(s(1))*3/20; % 升力系数
p=str2num(s(2))/20; % 4digit是max camber position最大弯度的位置，5digit就只是单纯分隔点
reflex=str2num(s(3)) % 这里只能是0，reflex为1的话，见下文
t=str2num(s(length(s)-1))/10+str2num(s(length(s)))/100; % thickness 最后两位是百分比
x=0:0.01:1;
% 如果需要两头密集、中间稀疏的
% x=0.5*(1-cos(x*pi));
a0= 0.2969;
a1=-0.1260;
a2=-0.3516;
a3= 0.2843;
a4=-0.1015; % For finite thick TE 后缘有厚度
% a4=-0.1036; % For zero thick TE 后缘无厚度
yt=(t/0.2)*(a0*sqrt(x)+a1*x+a2*x.^2+a3*x.^3+a4*x.^4);
P=[ 0.05 0.1 0.15 0.2 0.25 ]; %http://airfoiltools.com/airfoil/naca5digit 表格数据
M=[ 0.0580 0.1260 0.2025 0.2900 0.3910]; % 网址里面对应r
K=[361.4 51.64 15.957 6.643 3.230 ]; % 网址里面对应k1
m=spline(P,M,p); %三次样条插值
k1=spline(M,K,m);
xc1=x(find(x<=p));
xc2=x(find(x>p));
xc=[xc1, xc2];
if p==0 % 对称的
xu=x;
yu=yt;
xl=x;
yl=-yt;
zc=zeros(size(xc));
else
yc1=(1/6)*k1*( xc1.^3-3*m*xc1.^2+m^2*(3-m)*xc1 );
yc2=(1/6)*k1*m^3*(1-xc2);
zc=(Cl/0.3)*[yc1,yc2];
dyc1_dx=(1/6)*k1*( 3*xc1.^2-6*m*xc1+m^2*(3-m) );
dyc2_dx=repmat((1/6)*k1*m^3,size(xc2));
dyc_dx=[dyc1_dx,dyc2_dx];
theta=atan(dyc_dx);
xu=x-yt.*sin(theta);
zu=zc+yt.*cos(theta);
xl=x+yt.*sin(theta);
zl=zc-yt.*cos(theta);
end
x=[fliplr(xu), xl(2:end)];
z=[fliplr(zu), zl(2:end)];
upper_lower_separator=min( find(x==min(x)) );
xU=x(1:upper_lower_separator);
xL=x(upper_lower_separator:end);
zU=z(1:upper_lower_separator);
zL=z(upper_lower_separator:end);
plot(xU,zU,'b-',xL,zL,'r-',xc,zc,'k-.');
grid on;
axis equal;
运行后，直接输入23012即可
reflex==1的版本：
% NACA 5 digit airfoil generation code.
% Input 5 digit from command line, eg. "23112"
clear
s=input('\nInput four digit mumber of the NACA airfoil:\n?','s');
Cl=str2num(s(1))*3/20; % 升力系数
p=str2num(s(2))/20; % 4digit是max camber position最大弯度的位置，5digit就只是单纯分隔点
reflex=str2num(s(3)) % 1请用这个函数，0请用上面一个
t=str2num(s(length(s)-1))/10+str2num(s(length(s)))/100; % thickness 最后两位是百分比
x=0:0.01:1;
% 如果需要两头密集、中间稀疏的
%x=0.5*(1-cos(x*pi));
a0= 0.2969;
a1=-0.1260;
a2=-0.3516;
a3= 0.2843;
a4=-0.1015; % For finite thick TE 后缘有厚度
% a4=-0.1036; % For zero thick TE 后缘无厚度
yt=(t/0.2)*(a0*sqrt(x)+a1*x+a2*x.^2+a3*x.^3+a4*x.^4);
P=[ 0.1 0.15 0.2 0.25 ]; %http://airfoiltools.com/airfoil/naca5digit 表格数据 camber position
M=[ 0.1300 0.2170 0.3180 0.4410]; % 网址里面对应r
K1=[51.990 15.793 6.520 3.191]; % 网址里面对应k1
K2_K1 = [0.000764 0.00677 0.0303 0.1355];
m=spline(P,M,p); %三次样条插值 airfoiltools里面的r
k1=spline(M,K1,m);
k2_k1=spline(M,K2_K1,m);
xc1=x(find(x<=p));
xc2=x(find(x>p));
xc=[xc1, xc2];
if p==0 % 对称的
xu=x;
yu=yt;
xl=x;
yl=-yt;
zc=zeros(size(xc));
else
yc1=(1/6)*k1*( (xc1-m).^3-k2_k1 * (1-m).^3 * xc1 - m .^3 * xc1 +m.^3);
yc2=(1/6)*k1*( k2_k1 *(xc2-m).^3-k2_k1 * (1-m).^3 * xc2 - m .^3 * xc2 +m.^3) ;
zc=(Cl/0.3)*[yc1,yc2];
dyc1_dx=(1/6)*k1*( 3*(xc1-m).^2-k2_k1 * (1-m).^3 - m.^3 );
dyc2_dx=(1/6)*k1*( 3*k2_k1*(xc2-m).^2-k2_k1 * (1-m).^3 - m.^3 );
dyc_dx=[dyc1_dx,dyc2_dx];
theta=atan(dyc_dx);
xu=x-yt.*sin(theta);
zu=zc+yt.*cos(theta);
xl=x+yt.*sin(theta);
zl=zc-yt.*cos(theta);
end
x=[fliplr(xu), xl(2:end)];
z=[fliplr(zu), zl(2:end)];
upper_lower_separator=min( find(x==min(x)) );
xU=x(1:upper_lower_separator);
xL=x(upper_lower_separator:end);
zU=z(1:upper_lower_separator);
zL=z(upper_lower_separator:end);
plot(xU,zU,'b-',xL,zL,'r-',xc,zc,'k-.');
grid on;
axis equal;
例如，输入23112，即可生成。
浙公网安备 33010602011771号