User:Eas4200c.f08.aero.E/Matlab 1: NACA Airfoil

Matlab: NACA 4 Digit Airfoil Series
This results in a total area of 0.1027 m2 with the centroid at (0.4197, 0.0160).

% filename: airfoil % % PURPOSE: %      To plot an airfoil, calculate it average area and %      to calculate and plot the centroid. % % AUTHOR: %      Stephen Roberts function airfoil % Estabolish chord length and NACA parameters: % c = chord length % m = maximum camber in percentage of chord length c % p = position of maximum camber in tenths of the chord % ta = maximum airfoil thickness in percentage of the chord c = 1; m = 0.02*c; p = 0.4*c; ta = 0.15*c; % ns = number of segments ns = 100; % y = axis along the chord % zc_before = camber mean line before maximum camber % zc_after = camber mean line after maximum camber % ddy_zc_before = derivative of zc_before with respect to y % ddy_zc_after = derivative of zc_after with respect to y % zt_positive = thickness distribution above the camber mean line % zt_negative = thickness distribution below the camber mean line % theta_before = arc tangent of ddy_zc_before % theta_after = arc tangent of ddy_zc_after y = linspace(0,c,ns+1)'; zc_before = (m/(p^2))*(2*p*y - y.^2); zc_after = (m/(1-p)^2)*((1-2*p) + 2*p*y - y.^2); ddy_zc_before = -2*m*(y - p).*1/(p.^2); ddy_zc_after = -2*m*(y-p).*1/((p - 1).^2); zt_positive = (5*ta)*(0.2969*(y.^(1/2)) - 0.1260*(y) - 0.3516*(y.^2) + 0.2843*(y.^3) - 0.1015*(y.^4)); zt_negative = -(5*ta)*(0.2969*(y.^(1/2)) - 0.1260*(y) - 0.3516*(y.^2) + 0.2843*(y.^3) - 0.1015*(y.^4)); theta_before = atan(ddy_zc_before); theta_after = atan(ddy_zc_after); % determine location of maximum camber in terms of i i = 1; while y(i) < (0.4*c) %y(i) %i i = i + 1; end i=i-1; % input the nodal coordinates for j = 1:1:i position(:,j) = [y(j) - zt_positive(j)*sin(theta_before(j)); zc_before(j) + zt_positive(j)*cos(theta_before(j)); 0]; end for k = (i+1):1:(ns+1) position(:,k) = [y(k) - zt_positive(k)*sin(theta_after(k)); zc_after(k) + zt_positive(k)*cos(theta_after(k)); 0]; end for r = (ns+2):1:(ns+2+ns-i-1) position(:,r) = [y(2*ns-r+2) - zt_negative(2*ns-r+2)*sin(theta_after(2*ns-r+2)); zc_after(2*ns-r+2) + zt_negative(2*ns-r+2)*cos(theta_after(2*ns-r+2)); 0]; end for q = (ns+2+ns-i):1:(2*ns) position(:,q) = [y(2*ns-q+2) - zt_negative(2*ns-q+2)*sin(theta_before(2*ns-q+2)); zc_before(2*ns-q+2) + zt_negative(2*ns-q+2)*cos(theta_before(2*ns-q+2)); 0]; end % set up the nodal coordinate arrays for s = 1:1:(2*ns) x(s) = position(1,s); v(s) = position(2,s); z(s) = position(3,s); end % set up the element connectivity array node_connect for t = 1:1:(2*ns-1) node_connect(1,t) = t;    node_connect(2,t) = t+1; end node_connect(1,2*ns) = 2*ns; node_connect(2,2*ns) = 1; % Calculate total area for u = 1:1:(2*ns-1) dL(:,u) = (position(:,u+1) - position(:,u)); dA(:,u) = cross(position(:,u),dL(:,u)); end dL(:,2*ns) = (position(:,1) - position(:,2*ns)); dA(:,2*ns) = cross(position(:,2*ns),dL(:,2*ns)); dummy = sum(dA,2); Total_Area = abs((1/2)*dummy(3)) % Calculate centroid x_centroid(1) = 0; y_centroid(1) = 0; for u = 2:1:(2*ns) x_centroid(u) = (2/3)*(norm(position(:,u)))*position(1,u)/norm(position(:,u)); y_centroid(u) = (2/3)*(norm(position(:,u)))*position(2,u)/norm(position(:,u)); end X = sum(x_centroid.*abs((1/2)*dA(3,:)))/Total_Area Y = sum(y_centroid.*abs((1/2)*dA(3,:)))/Total_Area % plot the whole system for u = 1:1:(2*ns) node_1 = node_connect(1,u); node_2 = node_connect(2,u); yy = [x(node_1),x(node_2)]; zz = [v(node_1),v(node_2)]; xx = [z(node_1),z(node_2)]; plot3(yy,zz,xx,'-',X,Y,0,'+') hold on end % display the figure the way I want to axis([0,1,-.5,.5]) view([0 0 1]) grid title('NACA 2415')
 * }