NACA5桁翼型のプロット

en.wikipedia.org

function naca5(n1,n2,n3,n4,n5)
L = 0.15*n1
P = 0.05*n2
S = n3
TT = 0.1*n4+0.01*n5
t = TT;
n = 100;
xs = linspace(0,1,n);
xus = [];
xls = [];
yus = [];
yls = [];
if S==0 %non-reflex
    if n2==1
        p = 0.05;
        r = 0.0580;
        k1 = 361.40;
    elseif n2==2
        p = 0.10;
        r = 0.126;
        k1 = 51.640;
    elseif n2==3
        p = 0.15;
        r = 0.2025;
        k1 = 15.957;
    elseif n2==4
        p = 0.2;
        r = 0.29;
        k1 = 6.643;
    elseif n2==5
        p = 0.25;
        r = 0.391;
        k1 = 3.230;
    end
    for idx=1:n
        x = xs(idx);
        yt = 5*t*(0.2969*sqrt(x)-0.1260*x-0.3516*x^2+0.2843*x^3-0.1015*x^4);
        if x <= r
            yc = k1/6*(x^3-3*r*x^2+r^2*(3-r)*x);
            dyc_dx = k1/6*(3*x^2-6*r*x+r^2*(3-r));
        else
            yc =k1*r^3/6*(1-x);
            dyc_dx = -k1*r^3/6;
        end
        theta = atan2(dyc_dx,1);
        xu = x - yt*sin(theta);
        xl = x + yt*sin(theta);
        yu = yc + yt*cos(theta);
        yl = yc - yt*cos(theta);
        xus(end+1) = xu;
        xls(end+1) = xl;
        yus(end+1) = yu;
        yls(end+1) = yl;
    end
else %reflex
    if n2==2
        p = 0.10;
        r = 0.130;
        k1 = 51.990;
        k2_k1 = 0.000764;
    elseif n2==3
        p = 0.15;
        r = 0.217;
        k1 = 15.793;
        k2_k1 = 0.00677;
    elseif n2==4
        p = 0.20;
        r = 0.318;
        k1 = 6.520;
        k2_k1 = 0.0303;
    elseif n2==5
        p = 0.25;
        r = 0.441;
        k1 = 3.191;
        k2_k1 = 0.1355;
    end
    for idx=1:n
        x = xs(idx);
        yt = 5*t*(0.2969*sqrt(x)-0.1260*x-0.3516*x^2+0.2843*x^3-0.1015*x^4);
        if x <= r
            yc = k1/6*((x-r)^3-k2_k1*(1-r)^3-r^3*x+r^3);
            dyc_dx = k1/6*(3*(x-r)^2-k2_k1*(1-r)^3-r^3);
        else
            yc = k1/6*(k2_k1*(x-r)^3-k2_k1*(1-r)^3*x-r^3*x+r^3);
            dyc_dx = k1/6*(k2_k1*3*(x-r)^2-k2_k1*(1-r)^3-r^3);
        end
        theta = atan2(dyc_dx,1);
        xu = x - yt*sin(theta);
        xl = x + yt*sin(theta);
        yu = yc + yt*cos(theta);
        yl = yc - yt*cos(theta);
        xus(end+1) = xu;
        xls(end+1) = xl;
        yus(end+1) = yu;
        yls(end+1) = yl;
    end
end

figure();
plot(xus(xus<=p),yus(xus<=p),'r');
hold on;
plot(xus(xus>p),yus(xus>p),'b');

plot(xls(xus<=p),yls(xus<=p),'r');
plot(xls(xus>p),yls(xus>p),'b');
axis equal;
title(sprintf('NACA%d%d%d%d%d',n1,n2,n3,n4,n5));

5桁の意味
（NACA-LPSTT）
0.15*L　最大の揚力係数
0.05*P　最大キャンバー位置
キャンバーがSimple（0）かReflex（1）かどうか
TT　厚み

NACA23015の例
最大揚力係数0.3
最大キャンバー位置0.15
キャンバーはSimple
厚み15%

キャンバー位置の影響

キャンバーラインの影響（SimpleとReflex）
Simple

Reflex

厚みの影響