Clear Filters
Clear Filters

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

18 views (last 30 days)
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
% 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)
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));
theta(i) = atan(dyc_dx(i));
% 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);
% 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;
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.
  1 Comment
Keaton Looper
Keaton Looper on 26 Feb 2024
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.

Sign in to comment.

Answers (2)

David Goodmanson
David Goodmanson on 26 Feb 2024
Edited: David Goodmanson on 26 Feb 2024
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
Keaton Looper on 26 Feb 2024
Is there a way that you can see to generate random NACA variables more efficiently?
David Goodmanson
David Goodmanson on 28 Feb 2024
Edited: David Goodmanson on 28 Feb 2024
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)];
nacastr = [num2str(M) num2str(P) num2str(TT)];
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
a4 = -0.1036; % closed
% 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);

Sign in to comment.

Aya on 28 Feb 2024

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!