NacaAirfoil.m
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
endNACAAirfoilMesh.m
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);
