User:Eas4200c.f08.aero.E/Matlab 2: Convergence and Multi Cell

{| class="toccolours collapsible collapsed" style="width:100%" !Single Cell Airfoil Convergence

% filename: airfoil_convergence_study %    % PURPOSE: %      To find the value of J such that the difference between the %      current value and the previous value is less than 1% %      (Also plot graph showing J vs ns) %   % AUTHOR: %      Stephen Roberts function airfoil_convergence_study % NS = final value of number of segments % ns = intermediate value of number of segments End_num_Segments = 100; Start_num_Segments = 5; for ns = Start_num_Segments:1:End_num_Segments % 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 = .5; m = 0.02*c; p = 0.4*c; ta = 0.15*c; % 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_original = linspace(0,c,ns+1)'; y = y_original/c; 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 for i = 1:1:(ns+1) position(1,i) = y_original(i); end for i = (ns+2):1:(2*ns) position(1,i) = y_original(2*ns+2-i); 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 total length of exterior for i = 1:1:(2*ns) Length_Vector(i) = norm(dL(:,i)); end Total_Length = sum(Length_Vector); J(ns) = 4*((Total_Area)^2)*0.002/Total_Length; end % Fix J vector to start at Start_num_Segments for i = 1:1:(End_num_Segments-Start_num_Segments) j(i) = J(i+Start_num_Segments-1); end J=j; % Create vector "Number_of_segments" to represent the number of segments for i = 1:1:(End_num_Segments-Start_num_Segments) Number_of_segments(i) = Start_num_Segments+i-1; end % Calculate the percent difference between each J value for i = 1:1:(End_num_Segments-Start_num_Segments-1) Percent_Difference(i) = (((J(i+1))-(J(i)))/(((J(i+1))+(J(i)))/2))*100; end i = 1; while abs(Percent_Difference(i)) > 1 i = i+1; end % Plot J versus the number of segments plot(Number_of_segments,J,Number_of_segments(i+1),J(i+1),'+') xlabel('Number of Segments') ylabel('J') title('Single Cell Airfoil: Convergence Study') grid




 * }