Plotting the profile of an airfoil generated by random NACA numbers.

When I plot the airfoil profile of the randomly generated NACA parameters, the air foilf seems fine until it sharply cuts off and translates down the axis vertically then continues with the curve. I am not sure what the issue is.

n = 1 % number of air foils generated

% upper and lower limits of camber, position, and thickness
Camberlimits =  [0,9.5];
Positionlimits = [0,90];
Thicknesslimits = [1,40];

% Constants 
a0 = 0.2969;
a1 = -0.1260;
a2 = -0.3516;
a3 = 0.2843;
a4 = -0.1015; % open trailing edge
%a4 = -0.1036; % closed trailing edge

% NACA = MPTT
% M = max camber in percentage of the chord length
% P = max camber position in percentage of the chord length
% TT = thickness in percentage of the cord length

% generate random string to represent an airfoil of given parameters
% M,P,and TT

Var1 = Camberlimits(1) + (Camberlimits(2) - Camberlimits(1))*rand(n,1);
Var2 = Positionlimits(1) + (Positionlimits(2) - Positionlimits(1))*rand(n,1);
Var3 = Thicknesslimits(1) + (Thicknesslimits(2) - Thicknesslimits(1))*rand(n,1);

% Take the decimal values and round to the nearest integer
Var1 = round(Var1,0);
Var2 = round(Var2,0);
Var3 = round(Var3,0);

% Transform the numbers to strings

str1 = num2str(Var1);
str2 = num2str(Var2);
str3 = num2str(Var3);

% Combine the values into one string of element length 5
combinedstr = strcat(str1,str2,str3);

% Since we only want an element length of 4 we will select the value in
% position 1,2,4,and 5. 

var1 = combinedstr(:,1);
var2 = combinedstr(:,2);
var3 = combinedstr(:,4:end);

% Recombine the values into one string that describes the parameters M,P,
% and TT of each randomly generated airfoil.
string = strcat(var1,var2,var3);
stringnum = str2num(string);


%Extract values from NACA string
for m = 1:1:n

M0 = string(:,1)
P0 = string(:,2)
TT0 = string(:,3:4)

M0 = str2double(M0)
P0 = str2double(P0)
TT0 = str2double(TT0)

end 




M = M0/100;
P  = P0/10;
TT = TT0/100;

points = 100; % Number of points used to represent air foil profile

% airfoil grid

x = linspace(0,1,points)';

% Camber and gradient
yc =ones(points,1);
dyc_dx = ones(points,1);
theta = ones(points,1);

for i = 1:1:points
    if (x(i) >= P && x(i) < P)
        yc(i) = (M/P^2)*((2*P*x(i))-x(i)^2);
        dyc_dx(i) = ((2*M)/(P^2))*(P-x(i));
    elseif (x(i) >= P && x(i) <= 1)
        yc(i) = (M/(1-P)^2)*(1-(2*P)+(2*P*x(i))-(x(i)^2));
        dyc_dx(i) = ((2*M)/((1-P)^2))*(P-x(i));
    end 
    theta(i) = atan(dyc_dx(i));
end 

% Distribution of airfoil Thickness 
yt = ones(points,1);
for i = 1:1:points
    term0 = a0*sqrt(x(i));
    term1 = a1*x(i);
    term2 = a2*x(i)^2;
    term3 = a3*x(i)^3;
    term4 = a4*x(i)^4;

    yt(i) = 5*TT*(term0+term1+term2+term3+term4);
end

% Points for upper surface
xu = x(:) - yt(:).*sin(theta);
yu = yc(:) + yt(:).*cos(theta);

% for i = 1:1:points
%     xu(i) = x(i) - yt(i)*sin(theta(i));
%     yu(i) = yc(i) + yt(i)*cos(theta(i));
% end

% Points for lower surface

xl = x(:) + yt(:).*sin(theta);
yl = yc(:) - yt(:).*cos(theta);

% for i = 1:1:points
%     xl(i) = x(i) + yt(i)*sin(theta(i));
%     yl(i) = yc(i) - yt(i)*cos(theta(i));
% end

f1 = figure(1);
hold on;
grid on;
axis equal;
plot(xu,yu,'r-');
plot(xl,yl,'k-

This is the graph that I am given from the code.

Any help would be great. Thanks!

Each time the code is run it will give a random 4 digit NACA profile or should at least. There are still a couple glitches.

[Keaton Looper]

How would I generate multiple random airfoils at once and then graph each one? I am trying to use a CFD solver to find the drag and lift coefficients to then evaluate the fitness of each airfoil and use a genetic algorithm to do crossover and mutation of each selected parent airfoil and optimize the airfoil design.

Answer 1

[David Goodmanson]

Hello Keaton,

I changed the line

TT0 = string(:,3:4)     to     TT0 = string(:,3)

since the program does not run otherwise. The real change was in the Camber and gradient section going from

if (x(i) >= P && x(i) < P)    to     if (x(i) >= 0 && x(i) < P)

in which case you get airfoils without the jogs. The problem is that the 'if' test in the upper half of the for loop is never true. Fixing the lower limit allows that part of the for loop to run.

[Keaton Looper]

Is there a way that you can see to generate random NACA variables more efficiently?

[David Goodmanson]

The main thing is that you do not have to use a for loop and step through the x values one by one. Rather than suggest specific changes to your code, here is a function that produces the airfoil. I have not tested it vs. known results but I believe it is correct. The function is fairly quick so you should be able to invoke MPTT values one by one and save the results. For x>p, the expression for the camber yca looks different from the usual expression but is equivalent. I just wanted to show the symmetry between the pair x1,p and the pair x3,q

function [xup yup xdn ydn yca yth nacastr] = naca(M,P,TT,x,opt)
% naca four digit airfoil, MPTT
%
% input:
% M is max camber as percentage of the chord
% P is distance aft of leading edge for max camber, in tenths of chord
% TT is max thickness as percentage of the chord
% the chord runs from 0 to 1 in x
% opt is an optional input variable. The trailing edge
% is open if opt='o' and closed if opt = 'c'.  The default is 'c'.   
%
% output:
% xup,yup describe upper surface. xup has different spacing than x 
% xdn,ydn describe lower surface. xdn has different spacing than x
% yca is the camber mean line, function of x
% yth is the thickness, function of x, and does not depend on MPTT
% nacastr is the 4 digit string 'MPTT'
%
% [xup yup xdn ydn yca yth nacastr] = naca(M,P,TT,x)
% 
% David Goodmanson
if TT<=9
    nacastr = [num2str(M) num2str(P) num2str(0) num2str(TT)];
else
    nacastr = [num2str(M) num2str(P) num2str(TT)];
end
m = M/100;
p = P/10;
q = 1-p;
t = TT/100;
a0 = 0.2969;
a1 = -0.1260;
a2 = -0.3516;
a3 = 0.2843;
if nargin == 5 & opt == 'o'
    a4 = -0.1015; % open
else
    a4 = -0.1036; % closed
end
% thickness
yth = 5*t*(polyval([a4 a3 a2 a1 0],x) + a0*sqrt(x));
% camber
x1 = x(x<=p);
x2 = x(x>p);
x3 = 1-x2;
yca = [(m/p^2)*(x1.*(2*p-x1)) (m/q^2)*(x3.*(2*q-x3))];
dycadx = [(2*m/p^2)*(p -x1) -(2*m/q^2)*(q -x3)];
% rotation perpendicular to camber line
th = atan(dycadx);
xup = x - yth.*sin(th);
xdn = x + yth.*sin(th);
yup = yca + yth.*cos(th);
ydn = yca - yth.*cos(th);

Answer 2

[Aya]

vrinstall-check