叶片（airfoil）横截面方程及案例

本人一直潜水，现贡献一把。希望有分加。

付上一程序用MatLab 编写的，计算叶片的横截面方程及叶片形状。希望各位笑纳。有用的请用，没用的别骂。

运行后的点坐标可导入CAD。


clear;
clc;
syms V c alpha
alpha=2*pi/45; % the attack angle is 8 degree and change to radia;
c=1; % the nondimensional chord length;
n=200; % the half number of numerical panel quantity;
t=0.12*c; % the max thickness of the airfoil
X(1)=c;
for i=1:n/2+1;
        X(i)=c*(n/2+1-i)/(n/2);
        X(n+2-i)=X(i);
end
% the following is to form camber line for NACA 2412;
% m = maximum camber as fraction of chord,p = chord wise position of maximum
% camber. e.g. NACA 2412 means m = 0.02, p = 0.4 and tmax = 0.12.
m=0.02;
p=0.4;
for i=1:n+1
    if X(i)<=p
        YC(i)=(m/p^2)*(2*p*X(i)-X(i)^2);
    else
        YC(i)=(m/(1-p)^2)*((1-2*p)+2*p*X(i)-X(i)^2);
    end
end
for i=1:n+1
    if i<=n/2+1
        Y(i)=-5*t*(0.29690*X(i)^0.5-0.12600*X(i)-0.35160*X(i)^2+0.28430*X(i)^3-0.10150*X(i)^4)+YC(i);
    else
        Y(i)=+5*t*(0.29690*X(i)^0.5-0.12600*X(i)-0.35160*X(i)^2+0.28430*X(i)^3-0.10150*X(i)^4)+YC(i);
    end
end
for i=1:n
    x(i)=(X(i+1)+X(i))/2;
    y(i)=(Y(i+1)+Y(i))/2;
end
for i=1:n
    phi(i)=atan2((Y(i+1)-Y(i)),(X(i+1)-X(i)));
    if phi(i)<0
        phi(i)=phi(i)+2*pi;
    end
end
for i=1:n
    for j=1:n
        A(i,j)=-(x(i)-X(j))*cos(phi(j))-(y(i)-Y(j))*sin(phi(j));
        B(i,j)=(x(i)-X(j))^2+(y(i)-Y(j))^2;
        C1(i,j)=sin(phi(i)-phi(j));
        C2(i,j)=cos(phi(i)-phi(j));
        D(i,j)=(x(i)-X(j))*cos(phi(i))+(y(i)-Y(j))*sin(phi(i));
        S(j)=((X(j+1)-X(j))^2+(Y(j+1)-Y(j))^2)^0.5;
        E(i,j)=(x(i)-X(j))*sin(phi(j))-(y(i)-Y(j))*cos(phi(j)); % (B(i,j)-A(i,j)^2)^0.5;
    end
end
for i=1:n
    for j=1:n
        F(i,j)=log((S(j)^2+2*A(i,j)*S(j)+B(i,j))/B(i,j));
        G(i,j)=atan((S(j)+A(i,j))/E(i,j))-atan(A(i,j)/E(i,j));
        J(i,j)=-(C2(i,j)/2)*F(i,j)-C1(i,j)*G(i,j);
        K(i,j)=-C2(i,j)+(1/S(j))*(F(i,j)*(A(i,j)*C2(i,j)+D(i,j)/2)+G(i,j)*(A(i,j)*C1(i,j)+E(i,j)*C2(i,j)));
    end
end
for i=1:n
    for j=1:n
        if j==1
            JK(i,j)=J(i,j)-K(i,j);
        else
            JK(i,j)=K(i,(j-1))+J(i,j)-K(i,j);
        end
    end
    JK(i,n+1)=K(i,n);
    JK(n+1,1)=1;
    JK(n+1,n+1)=1;
end
for i=1:n
    RHS(i,1)=2*pi*sin(alpha-phi(i)); % 2*pi*V*sin(alpha-phi(i));
    RHS(n+1,1)=0;
end
gamma=(inv(JK)*RHS)';
for i=1:n+1
    Cp(i)=1-(gamma(i))^2; % 1-(gamma(i)/V)^2;
    vpa(Cp,6);
end
hold on
plot(X,Y);axis equal
xlabel('X');
ylabel('Y');
title('Figure of the NACA 2412 airfoil');
figure
hold off
hold on
plot(X,Cp,'k-+');
xlabel('X');
ylabel('Cp');
title('x~Cp');
hold off
% the following is to calculate the whole circulation Gamma;
for i=1:n
    Gamma=sum(((gamma(i)+gamma(i+1))/2)*S(i));
end
% % Validation on an elliptic airfoil
t=10;
k=(t-1)/(t+1);
for i=1:n/2+1
    thetal(i)=pi*(n/2+1-i)/(n/2);
    thetal(n+2-i)=thetal(i);
     v(i)=-2*sin(thetal(i))/((1-k)^2+4*k*sin(thetal(i))^2)^0.5;
     Cp1(i)=1-v(i)^2;
end
for i=1:n+1
    v(i)=-2*sin(thetal(i))/((1-k)^2+4*k*sin(thetal(i))^2)^0.5;
    Cp1(i)=1-v(i)^2;
end
for i=1:n+1
     thetau(i)=pi*(i-1)/n;
     v(i)=-2*sin(thetau(i))/((1-k)^2+4*k*sin(thetau(i))^2)^0.5;
     Cpu(i)=1-v(i)^2;
end
figure
hold on
plot(thetal,Cp1,'R-*');
plot(thetal,v,'k:.');
legend('thetal~Cpl','thetal~velocity');
% plot(thetau,Cpu);
xlabel('theta');
ylabel('Cp');
hold off

不错，谢谢分享了

嗯，很好！支持一下！

看不懂，但是支持一下楼主！

支持一下