User:Eas4200c.f08.gator.edwards/Matlab Program

The following MATLAB function was created to plot the shape of a NACA 4-digit series airfoil. It also calculates the area and centroid of the cross section. Following the function is the implementation and output, copied out of the command window.

function airfoil(NACA, varargin)
 
% Plot NACA 4-digit airfoil 
% Inputs:
%   NACA - 4-digit NACA airfoil description
% Optional Inputs:
%   c - chord length 
%     default:  c = 1
%   ns - number of segments used to calculate area (must be multiple of 10)
%     default:  ns = 1000
%   ycen, zcen - coordinates for center point of area calculation
%     default:  ycen = zcen = 0
% Outputs:
%   A - area of airfoil cross section
%   Plot of airfoil shape with camber line and centroid marked
 
% Process optional inputs
if nargin == 1
    c = 1;
    ns = 1000;
    ycen = 0;
    zcen = 0;
elseif nargin == 2
    c = varargin{1};
    ns = 1000;
    ycen = 0;
    zcen = 0;
elseif nargin == 3
    c = varargin{1};
    ns = varargin{2};
    ycen = 0;
    zcen = 0;
else
    c = varargin{1};
    ns = varargin{2};
    ycen = varargin{3};
    zcen = varargin{4};
    if floor(ns / 10) - ns / 10 ~= 0
        error('ns must be a multiple of 10')
    end
end
 
% Convert inputs to percentages
NACA = num2str(NACA);
m = str2num(NACA(1)) / 100;
p = str2num(NACA(2)) / 10;
ta = str2num(NACA(3:4)) / 100;
 
% Divide y-axis into ns segments
y = 0:1/ns:1;
 
% Calculation of mean camber line (zc) and dz/dy
for i = 1:length(y)
    if i <= p * length(y)
        zc(i) = (m / p^2) * (2 * p * y(i) - y(i)^2);
        dzdy(i) = (m / p^2) * (2 * p - 2 * y(i));
    else
        zc(i) = (m / (1 - p)^2) * ((1 - 2 * p) + 2 * p * y(i) - y(i)^2);
        dzdy(i) = (m / (1 - p)^2) * (2 * p - 2 * y(i));
    end
end
 
% Thickness distribution calculation
zt = 5 * ta .* (.2969 .* sqrt(y) - .1260 .* y - .3516 .* y.^2 + .2843 .* y.^3 - .1015 .* y.^4);
 
% Calculation of angle from y-axis to tangent to mean camber line
theta = atan(dzdy);
 
% Calculate coordinates of points on airfoil surface
yu = y - zt .* sin(theta);
zu = zc + zt .* cos(theta);
yl = y + zt .* sin(theta);
zl = zc - zt .* cos(theta);
 
% Scale to cord length
zc = zc * c;
y = y * c;
yu = yu * c;
zu = zu * c;
yl = yl * c;
zl = zl * c;
 
% Calculate area of cross section
A = 0;
for i = 1:ns
    r = [yu(i) - ycen, zu(i) - zcen, 0];
    dl = [yu(i+1) - yu(i), zu(i+1) - zu(i), 0];
    dA = cross(dl, r) / 2;
    A = A + dA;
end
for i = 1:ns
    r = [yl(i) - ycen, zl(i) - zcen, 0];
    dl = [yl(i+1) - yl(i), zl(i+1) - zl(i), 0];
    dA = cross(r, dl) / 2;
    A = A + dA;
end
 
% Calculate centroid location
ybarA = 0;
zbarA = 0;
sumdA = 0;
for i = 1:ns
    dA = (yu(i+1) - yu(i)) * (zu(i+1) + zu(i)) / 2;
    dcenty = (yu(i+1) + yu(i)) / 2;
    dcentz = (zu(i+1) + zu(i)) / 4;
    ybarA = ybarA + dA * dcenty;
    zbarA = zbarA + dA * dcentz;
    sumdA = sumdA + dA;
end
for i = 1:ns
    dA = abs((yl(i+1) - yl(i)) * (zl(i+1) + zl(i)) / 2);
    dcenty = (yl(i+1) + yl(i)) / 2;
    dcentz = (zl(i+1) + zl(i)) / 4;
    ybarA = ybarA + dA * dcenty;
    zbarA = zbarA + dA * dcentz;
    sumdA = sumdA + dA;
end
ybar = ybarA / sumdA;
zbar = zbarA / sumdA;
 
% Create crosshair on Centroid 
a = ybar + .025 * c;
b = ybar - .025 * c;
e = zbar + .025 * c;
d = zbar - .025 * c;
 
a = [a b];
b = [zbar zbar];
d = [e d];
e = [ybar ybar];
 
% Combine upper and lower airfoil outlines
ytot = [yu fliplr(yl)];
ztot = [zu fliplr(zl)];
 
% Plot mean camber line, airfoil outline, and centroid crosshair
hold off
plot(y, zc, 'b-')
hold on
plot(ytot, ztot, 'k-')
plot(a, b, 'g-', e, d, 'g-')
axis([0 c -c/2 c/2])
legend('Mean Camber Line', 'Airfoil Outline', 'Centroid')
hold off
 
% Output area and centroid location
fprintf('\nCross Sectional Area A = %5.4f\n', A(1,3))
fprintf('Centroid Location C = (%5.4f, %5.4f)\n\n', ybar, zbar)

---

EDU>> airfoil(2415);

Cross Sectional Area A = 0.1013 Centroid Location C = (0.4203, 0.0157)

File:Airfoil Plot.jpg