NACA翼型发生器对于导入icem建模有非常大的帮助。

function [y,x] = NacaAirfoil(varargin)

switch nargin
    case 0
        AF= '0012';   % designation
        M = AF(1);   % The maximum value of the mean line as chord/100.
        P = AF(2);   % Is the chordwise position of the maximum chamber as chord/10.
        XX= AF(3:4); % Is the maximum thickness, t/c, in chord precent (x 100%).
        TE= 'zero';  % Trailing edge (TE)
        n = 30;      % x-nodes to describe airfoil
        xs='halfcos';% x-nodes spacing method
        wf= 0;       % writte solution file
    case 1
        AF= varargin{1}; 
        M = AF(1);
        P = AF(2);
        XX= AF(3:4);
        TE= 'zero';
        n = 30;      
        xs= 'halfcos';
        wf= 0;
  case 2
        AF= varargin{1}; 
        M = AF(1);
        P = AF(2);
        XX= AF(3:4);
        TE= varargin{2};
        n = 30;      
        xs= 'halfcos';
        wf= 0;
  case 4
        AF= varargin{1}; 
        M = AF(1);
        P = AF(2);
        XX= AF(3:4);
        TE= varargin{2};
        n = varargin{3};      
        xs= varargin{4};
        wf= 0;
  case 5
        AF= varargin{1}; 
        M = AF(1);
        P = AF(2);
        XX= AF(3:4);
        TE= varargin{2};
        n = varargin{3};      
        xs= varargin{4};
        wf= 1;
        pth=varargin{5};
    otherwise
        error('wrong number of input arguments, check again dummy.')
end

% set airfoil main characteristics
m=str2double(M)/100;   % mean line coord
p=str2double(P)/10;    % maximun chamber position
t=str2double(XX)/100;  % relative thickness (t/c)

% Fixed parameters
a0= 0.2969;
a1=-0.1260;
a2=-0.3516;
a3= 0.2843;
switch TE 
    case 'finite'   % For finite thick TE
        a4=-0.1015;  
    case 'zero'     % For zero thick TE
        a4=-0.1036;
    otherwise
        error('not a valid trailing edge')
end

%% Build Airfoil cross section
% Set spacing for discrete values in x-direction
switch xs 
    case 'halfcos'  % Half cosine based spacing
        beta=linspace(0,pi,n)'; x=(0.5*(1-cos(beta))); 
        header=['NACA',AF,':[',num2str(2*n),'panels,Half cosine x-spacing]'];  
    case 'uniform'  % uniform based spacing
        x=linspace(0,1,n)';
        header=['NACA',AF,':[',num2str(2*n),'panels,Uniform x-spacing]'];
    otherwise
        error('not a valid spacing method')
end

% find points before and after the maximum chamber position 'p'
xc1=x(x<=p);  % for t <= p
xc2=x(x> p);  % for t > p
xc=[xc1;xc2]; % concatenate vectors

% compute thickness function yt(x) 
yt=(t/0.2)*(a0*sqrt(x)+a1*x+a2*x.^2+a3*x.^3+a4*x.^4);

% 
if p==0  % symmetric airfoil
    % thickness embelope with camber or mean line.
    xu=x;    yu= yt;    % upper surface
    xl=x;    yl=-yt;    % lower surface
    yc=zeros(size(xc)); % camber line function
else  % non-symmetric airfoil
    % camber lines
    yc1=(m/p^2)*(2*p*xc1-xc1.^2);               % for t <= p
    yc2=(m/(1-p)^2)*((1-2*p)+2*p*xc2-xc2.^2);   % for t > p
    yc=[yc1;yc2];       % camber line function yc(x)

    dyc1_dx=(m/p^2)*(2*p-2*xc1);                % for t <= p
    dyc2_dx=(m/(1-p)^2)*(2*p-2*xc2);            % for t > p
    dyc_dx=[dyc1_dx;dyc2_dx];   % camber line slope dyc/dx
    theta=atan(dyc_dx);  % camber line slope angle 

    xu= x-yt.*sin(theta);  yu=yc+yt.*cos(theta); % upper surface
    xl= x+yt.*sin(theta);  yl=yc-yt.*cos(theta); % lower surface
end

% Plot airfoid embelope
plot(xu,yu,'bo-'); hold on  % upper surface
plot(xl,yl,'ro-'); hold off % lower surface
axis equal

% Output points
x=[flipud(xu) ; xl(2:end)];
y=[flipud(yu) ; yl(2:end)];


%% Write Points to file
if wf == 1
    F1=header;
    F2=num2str([x,y]);
    F=strvcat(F1,F2);
    fileName=[pth,'NACA',AF,'.dat'];
    dlmwrite(fileName,F,'delimiter','')
end

end

close all;
 
%-------------------------------------------------------------------------% 
% Define all the parameters 
%-------------------------------------------------------------------------% 
 
Na = 50;
M = 30;
A = 10; %Clustering points near leading edge, inf = no clustering
B = 2; %Clustering points near the surface, 1 = no clustering
t = .12;
c = 1;
L = 1;
N = floor(2*c/(2*c+L)*Na);
x1 = c*ones(1,N+1);
y1 = zeros(1,N+1);
x2 = x1;
y2 = y1;
y2(1) = c;
s1 = y1;
s2 = s1;
N2 = 250;

%-------------------------------------------------------------------------% 
% Solve Algebraic Method
%-------------------------------------------------------------------------%
for i = 2:N2+1
    x1(i) = c*(1-(i-1)/N2);
    y1(i) = max([t/.2*(.2969*x1(i)^.5-.126*x1(i)-.3516*x1(i)^2+.2843*x1(i)^3-.1015*x1(i)^4) 0]);
    s1(i) = s1(i-1)+sqrt((x1(i)-x1(i-1))^2+(y1(i)-y1(i-1))^2);
    x2(i) = c*(1-2*(i-1)/N2);
    if x2(i) >= 0;
        y2(i) = c;
    else
        y2(i) = sqrt(c^2-(x2(i))^2);
    end
    s2(i) = s2(i-1)+sqrt((x2(i)-x2(i-1))^2+(y2(i)-y2(i-1))^2);
end
S1 = s1(end)*(1-exp(-(0:N)/A))/(1-exp(-(N)/A));
S1(end)/s1(end)
S2 = (0:N)*s2(end)/N;
for i = 1:N+1
    X1(i) = interp1(s1,x1,S1(i));
    Y1(i) = interp1(s1,y1,S1(i));
    X2(i) = interp1(s2,x2,S2(i));
    Y2(i) = interp1(s2,y2,S2(i));
end
r = exp(((1:(M+1))-M+1)/((M+1)/B));
r = r-r(1);
r = r/max(r);
for i = 1:N+1
    for j = 1:M+1
        x(i,j) = X1(i) + r(j)*(X2(i)-X1(i));
        y(i,j) = Y1(i) + r(j)*(Y2(i)-Y1(i));
    end
end
N = Na-N;
x = [((L+c):-L/N:(c+L/N))'*ones(1,M+1); x];
y = [(c*r'*ones(1,N))'; y];
x = [x; x(end-1:-1:1,:)];
y = [y ;-y(end-1:-1:1,:)];
[n,m] = size(x);

%-------------------------------------------------------------------------% 
% Plot Algebraic Method
%-------------------------------------------------------------------------%
mesh(x,y,zeros(n,m))
view(0,90)
colormap([0 0 0])
axis('equal','tight');
xlabel('Length [unitless]','FontSize',14);
ylabel('Length [unitless]','FontSize',14);
title('Algebraic Method ','FontSize',18);

%-------------------------------------------------------------------------% 
% Predetermined Constants
%-------------------------------------------------------------------------%
de = 1/M;
de2 = de^2;
dn = 1/Na;
dn2 = dn^2;
dedn = 2/(M*Na);
P = 100; %number of iterations for the elliptic solver.

%-------------------------------------------------------------------------% 
% Solve Elliptic Method
%-------------------------------------------------------------------------%
 
for k = 1:P
    for i = 2:2*Na-1
        for j = 2:M-1
            xn = (x(i+1,j)-x(i-1,j))/(2*dn);
            yn = (y(i+1,j)-y(i-1,j))/(2*dn);
            xe = (x(i,j+1)-x(i,j-1))/(2*de);
            ye = (y(i,j+1)-y(i,j-1))/(2*de);
            a(i,j) = xn^2+yn^2;
            b(i,j) = xe*xn+ye*yn;
            c(i,j) = xe^2+ye^2;
        end
    end
    for i = 2:2*Na-1
        for j = 2:M-1
            x(i,j) = (a(i,j)/de2*(x(i+1,j)+x(i-1,j))+c(i,j)/dn2*(x(i,j+1)+x(i,j-1))...
                -b(i,j)/dedn*(x(i+1,j+1)-x(i+1,j-1)+x(i-1,j-1)-x(i-1,j+1)))...
                /(2*(a(i,j)/de2+c(i,j)/dn2));
            y(i,j) = (a(i,j)/de2*(y(i+1,j)+y(i-1,j))+c(i,j)/dn2*(y(i,j+1)+y(i,j-1))...
                -b(i,j)/dedn*(y(i+1,j+1)-y(i+1,j-1)+y(i-1,j-1)-y(i-1,j+1)))...
                /(2*(a(i,j)/de2+c(i,j)/dn2));
        end
    end
end
xn = ones(Na,M)/0;
yn = xn;
xe = xn;
ye = xn;
for i = 2:Na-1
    for j = 2:M-1
        xn(i,j) = (x(i+1,j)-x(i-1,j))/(2*dn);
        yn(i,j) = (y(i+1,j)-y(i-1,j))/(2*dn);
        xe(i,j) = (x(i,j+1)-x(i,j-1))/(2*de);
        ye(i,j) = (y(i,j+1)-y(i,j-1))/(2*de);
        J(i,j) = yn(i,j)*xe(i,j)-xn(i,j)*ye(i,j); % Jacobian
    end
end
figure(2)
mesh(x,y,zeros(n,m))
xlabel('Length [unitless]','FontSize',14);
ylabel('Length [unitless]','FontSize',14);
title('Elliptic Method ','FontSize',18);
axis('equal','tight');

%-------------------------------------------------------------------------% 
% Plot Elliptic Method
%-------------------------------------------------------------------------%
 
view(0,90)
colormap([0 0 0])

%-------------------------------------------------------------------------% 
% Plot Jacobian
%-------------------------------------------------------------------------%
 
figure(3)
x=2:30; y=2:50; [X,Y]=meshgrid(x,y);
surf(X,Y,J);
view(-45,45)
title('Jacobian ','FontSize',18);