%Ronald %all alfas in DEGREES! %by ascorrea, with help from Jacob Huffman %"filler" variable just makes it easier to switch between Ronald and Grimace function Ronald(paramfile, pranformfile, filler) testrange={'20'}; for testnumber=1:length(testrange) prefix='FDR0'; name=strcat(prefix, testrange{testnumber}); paramfile=strcat(name, '.param'); pranformfile=strcat(name, '.plan'); %DIRECTORY LOCATIONS--------------------------------------------------------------] rootLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Matlab Code'; airfoilLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Airfoils'; pranformLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Planforms'; pranformConstraintLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Planform Constraints'; paramLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Parameter Files'; outputsLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Ronaldplus Outputs'; %-----------------------------------------------------------------------] n=50; cd(paramLOCATION); [variablenames, variablevalues]=textread(paramfile, '%s %f'); inputvalues=num2str(variablevalues); input=strcat(variablenames, inputvalues); [constantnames, constantvalues]=textread('constants.param', '%s %f'); cd(pranformLOCATION) [spans, chords, offsets, airfoil]=textread(pranformfile, '%f %f %f %s'); cd(airfoilLOCATION) [sref swet]=WhettedArea(spans, chords, offsets, airfoil, n); for k = 1:length(variablenames) eval([variablenames{k} '= variablevalues(k);']); end for k = 1:length(constantnames) eval([constantnames{k} '= constantvalues(k);']); end Dragredux_BLI=(1-Dragredux_BLI/100); wingspan=2*spans(end); AR=wingspan^2/sref; % wakeplot=findwake(airfoil, chords, spans, offsets, n, alt_cruise, M_cruise, laminarflow_cruise, fuse_width); %TAKEOFF------------------------------------------------------------------] % cd0_to=findcd0(c_f, M_to, sref, swet); [cd0_to D_wing_to D_fuse_to]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_to, alt_to, M_to, sref, c_f, D_engine, cd_landinggear, sfront_landinggear, fuse_width); % finddragpolar(cd0_to, cLalpha_to, cL0, efficiency_to, AR, cd_correction, 'Takeoff') % cLreq_to=findcLrequired(alt_to, sref, M_to, weight_to); cLreq_to=.48; alphareq_to=findalpharequired(cLalpha, cL0, cLreq_to); cD_to=findcD(AR, efficiency_to, cd0_to, cLreq_to, cd_correction); eloverdee_to=cLreq_to/cD_to; cM_to=cMalpha*alphareq_to+cM0; q_to=finddynpress(alt_to, M_to); cD_to=cD_to*Dragredux_BLI; cd0_to=cd0_to*Dragredux_BLI; D_to=cD_to*q_to*sref; %CRUISE------------------------------------------------------------------] [cd0_cruise D_wing_cruise D_fuse_cruise]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_cruise, alt_cruise, M_cruise, sref, c_f, D_engine, 0, 0, fuse_width); % finddragpolar(cd0_cruise, cLalpha, cL0, efficiency_cruise, AR, cd_correction, 'Cruise') cLreq_cruise=findcLrequired(alt_cruise, sref, M_cruise, weight_cruise); alphareq_cruise=findalpharequired(cLalpha, cL0, cLreq_cruise); cD_cruise=findcD(AR, efficiency_cruise, cd0_cruise, cLreq_cruise, cd_correction); cM_cruise=cMalpha*alphareq_cruise+cM0; q_cruise=finddynpress(alt_cruise, M_cruise); cD_cruise=cD_cruise*Dragredux_BLI; cd0_cruise=cd0_cruise*Dragredux_BLI; D_cruise=cD_cruise*q_cruise*sref; eloverdee_cruise=cLreq_cruise/cD_cruise; %LOITER------------------------------------------------------------------] [cd0_loiter D_wing_loiter D_fuse_loiter]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_loiter, alt_loiter, M_loiter, sref, c_f, D_engine, 0, 0, fuse_width); % finddragpolar(cd0_loiter, cLalpha_loiter, cL0, efficiency_loiter, AR, cd_correction, 'Loiter') cLreq_loiter=findcLrequired(alt_loiter, sref, M_loiter, weight_loiter); alphareq_loiter=findalpharequired(cLalpha, cL0, cLreq_loiter); cD_loiter=findcD(AR, efficiency_loiter, cd0_loiter, cLreq_loiter, cd_correction); cM_loiter=cMalpha*alphareq_loiter+cM0; q_loiter=finddynpress(alt_loiter, M_loiter); cD_loiter=cD_loiter*Dragredux_BLI; cd0_loiter=cd0_loiter*Dragredux_BLI; D_loiter=cD_loiter*q_loiter*sref; eloverdee_loiter=cLreq_loiter/cD_loiter; %LANDING------------------------------------------------------------------] [cd0_landing D_wing_landing D_fuse_landing]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_landing, alt_landing, M_landing, sref, c_f, D_engine, cd_landinggear, sfront_landinggear, fuse_width); % finddragpolar(cd0_landing, cLalpha_landing, cL0, efficiency_landing, AR, cd_correction, 'Landing') % cLreq_landing=findcLrequired(alt_landing, sref, M_landing, weight_landing); cLreq_landing=.401; alphareq_landing=findalpharequired(cLalpha, cL0, cLreq_landing); cD_landing=findcD(AR, efficiency_landing, cd0_landing, cLreq_landing, cd_correction); cM_landing=cMalpha*alphareq_landing+cM0; q_landing=finddynpress(alt_landing, M_landing); cD_landing=cD_landing*Dragredux_BLI; cd0_landing=cd0_landing*Dragredux_BLI; D_landing=cD_landing*q_landing*sref; eloverdee_landing=cLreq_landing/cD_landing; cd(outputsLOCATION) filenameplan=regexprep(pranformfile, '.plan', '_planform.tiff'); filenameoutput=regexprep(pranformfile, '.plan', '.txt'); plottitle=regexprep(pranformfile, '.plan', ''); fid = fopen(filenameoutput, 'w'); plottitle fprintf(fid, plottitle); fprintf(fid,'\nOutputs--------------------------------------------------\n'); fprintf(fid, 'Wingspan: %f, AR: %f, Planform Area: %f, Wetted Area: %f \n', wingspan, AR, sref, swet); fprintf(fid, '\n----------Takeoff--------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_to, cd0_to, cLreq_to, cD_to, D_to, cM_to, eloverdee_to); fprintf(fid, '\n----------Cruise---------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n D_Wing:%f, D_fuselage:%f', alphareq_cruise, cd0_cruise, cLreq_cruise, cD_cruise, D_cruise, cM_cruise, eloverdee_cruise, D_wing_cruise, D_fuse_cruise); fprintf(fid, '\n----------Loiter---------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_loiter, cd0_loiter, cLreq_loiter, cD_loiter, D_loiter, cM_loiter, eloverdee_loiter); fprintf(fid, '\n----------Landing--------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_landing, cd0_landing, cLreq_landing, cD_landing, D_landing, cM_landing, eloverdee_landing); fprintf(fid, '\nInputs--------------------------------------------------\n'); for k=1:length(input)-1 fprintf(fid, '%s %s \n',input{k}, input{k+1}); end fclose(fid); cd(rootLOCATION) end end % function cd0=findcd0(c_f, M, sref, swet) % if M>.7 % c_f=c_f*(1-.09*M^2); % end % cd0=c_f*(swet/sref); % end function [cd0 D_wing D_fuse]=findcd0(spans, chords, offsets, airfoil, n, laminarflow, alt, v_mach, sref, c_f, D_engine, cd_lg, Sfront_lg, width) wingSpan = 1; % n = 50; width=width/2; [q v_fps dynvisc]=findconstants(alt, v_mach); drag=0; drag_fuse=0; drag_wing=0; Chords = findChords(n, chords, spans); for i = 1:(size(chords)-1) for j = 1:n dist = wingSpan*(spans(i+1)-spans(i))./n; chordloc=linspace(0, Chords(i,j)); delta_x = Chords(i,j)/100; if spans(i)>width for k = 1:100 x(k)= k * Chords(i,j)/100; Re_x(k)= v_fps * x(k)/dynvisc; if (k <= 100*laminarflow ) cf_x(k) = 0.664 /(sqrt(Re_x(k))); else % Turbulent cf_x(k) = 0.455 /(log(0.06 * Re_x(k))^2); end xcf_x(k)=x(k)*cf_x(k); end t = trapz(cf_x); Cf = 2 * delta_x * t /Chords(i,j); drag_section = q*Cf*findPerim(airfoil{i}, Chords(i,j), wingSpan).*dist; drag_wing=drag_wing+drag_section; else for k = 1:100 x(k)= k * Chords(i,j)/100; Re_x(k)= v_fps * x(k)/dynvisc; if (k <= 100*laminarflow -10) cf_x(k) = 0.664 /(sqrt(Re_x(k))); else % Turbulent cf_x(k) = 0.455 /(log(0.06 * Re_x(k))^2); end xcf_x(k)=x(k)*cf_x(k); end t = trapz(cf_x); Cf = 2 * delta_x * t /Chords(i,j); drag_section = q*Cf*findPerim(airfoil{i}, Chords(i,j), wingSpan).*dist; drag_fuse=drag_fuse+drag_section; end t = trapz(cf_x); Cf = 2 * delta_x * t /Chords(i,j); drag_section = q*Cf*findPerim(airfoil{i}, Chords(i,j), wingSpan).*dist; drag=drag+drag_section; end end if v_mach>.7 c_f=c_f*(1-.09*v_mach)^2; end % cd0engine=c_f*(s_engine/sref); D_lg=Sfront_lg*q*cd_lg; D_wing=drag_wing*2; D_fuse=drag_fuse*2; drag_tot=D_wing+D_engine+D_lg+D_fuse; cd0=drag_tot/(q*sref); cD_lg=D_lg/(q*sref); cD_engine=D_engine/(q*sref); cD_wing=D_wing/(q*sref); cD_fuse=D_fuse/(q*sref); end function finddragpolar(cd0, cl_alfa, cL0, e, ar, cd_correction, filename) alphadeg=[-12:10]; % alpha=(pi/180)*(alphadeg); alpha=alphadeg; cL_alfa=cl_alfa*(ar/(ar+(2*(ar+4)/(ar+2)))); for i=1:length(alpha) C_L(i)=cL_alfa*(alpha(i))+cL0; cd_i(i)=C_L(i)^2/(pi*ar*e); C_d(i)=cd0+cd_i(i)+cd_correction; end graphname='Drag Polar'; hold on grid on axisX0=linspace(0,max(cd)); axisY0=zeros(1,length(axisX0)); plot(axisX0, axisY0,'k--'); plot(C_d, C_L, 'k', 'LineWidth', 2); plot(C_d, C_L, 'k.'); title(graphname); xlabel('C_d', 'FontSize', 16); ylabel('C_L', 'FontSize', 16); axis([0, .02, -.3, .4]); set(gca, 'FontSize',16) % for i=1:length(alphadeg) % alphalabel=num2str(alphadeg(i)); % label=[' \alpha =' alphalabel]; % text(C_d(i), C_L(i), label) % end print ('-dtiff', filename) grid off close(1) end function q=finddynpress(alt, v_mach) if alt < 36152 temp=59-.00356*alt; press=2116*((temp+459.7)/518.6)^5.256; elseif alt < 82345 temp=-70; press=473.1*exp(1.73-.000048*alt); else print('Program Not Valid for given altitude') end %slpcuft: clugs ber cubic foot density_slpcuft=press/(1718*(temp+459.7)); v_fps=v_mach*sqrt(1.4*1718*(temp+459.7)); q=.5*density_slpcuft*v_fps^2; end function [q v_fps dynvisc]=findconstants(alt, v_mach) if alt < 36152 temp=59-.00356*alt; press=2116*((temp+459.7)/518.6)^5.256; elseif alt < 82345 temp=-70; press=473.1*exp(1.73-.000048*alt); else print('Program Not Valid for given altitude') end %slpcuft: clugs ber cubic foot density_slpcuft=press/(1718*(temp+459.7)); v_fps=v_mach*sqrt(1.4*1718*(temp+459.7)); q=finddynpress(alt, v_mach); S = 110.4; beta = 1.458e-6; temp_k=[(temp-32) / (1.8)] + 273 ; viscosity_pascalseconds = (beta*temp_k^(1.5))/(temp_k+S); %kgpmcu:kilograms per meter cubed %msqpse:meters squared per second %sqftpsec:square feet per second density_kgpmcu=density_slpcuft*515.4; dynamic_viscosity_msqpsec=viscosity_pascalseconds/density_kgpmcu; dynamic_viscosity_sqftpsec=dynamic_viscosity_msqpsec*10.76; dynvisc=dynamic_viscosity_sqftpsec; density=density_slpcuft; end function cLreq=findcLrequired(alt, s, v_mach, lift) q=finddynpress(alt, v_mach); cLreq=lift/(q*s); end function alpha=findalpharequired(cLalpha, cL0, cLrequired) alpha=(cLrequired-cL0)/cLalpha; end function cD=findcD(ar, e, cd0, C_L, cd_correction) % e=4.61*(1-.045*ar^.68)*cosd(sweep_LE)^.15-3.1; cd_i=C_L^2/(pi*ar*e); cd_i=C_L^2/(pi*ar*e); cD=cd0+cd_i+cd_correction; end function [pranformArea WetArea]=WhettedArea(spans, chords, offsets, airfoil, n) % fileName = 'AirData.mat'; wingSpan = 1; % n = 50; engineArea=0; Area = 0; pranformArea = 0; Chords = findChords(n, chords, spans); for i = 1:(size(chords)-1) for j = 1:n Chords(i,j); dist = wingSpan*(spans(i+1)-spans(i))./n; Area = Area+findPerim( airfoil{i}, Chords(i,j), wingSpan).*dist; end end for i = 1:(size(chords)-1) for j = 1:n dist = wingSpan*(spans(i+1)-spans(i))./n; pranformArea = pranformArea+Chords(i,j).*dist; % fprintf('Chord: %f pranformArea: %f \n', Chords(i,j), pranformArea) end end Area = Area *2+engineArea; WetArea = Area; pranformArea = pranformArea*2; end function Chords = findChords( n, Chord, X ) for( i = 1:(length(Chord)-1) ) Xtemp = linspace(X(i),X(i+1),n); Chords(i,:) = (Chord(i+1)-Chord(i))./(X(i+1)-X(i)).*(Xtemp-X(i))+Chord(i); end end function peri = findPerim(airFile, chord, span) peri = 0; Airfoil = load(airFile); Airfoil = Airfoil*chord*span; for i = 1:(size(Airfoil,1)-1) peri = peri + sqrt((Airfoil(i+1,1)-Airfoil(i,1))^2+(Airfoil(i+1,2)-Airfoil(i,2))^2); end end function [plotLE plotTE]=plotpranform(spans, chords, offsets, n) Chords = findChords(n, chords, spans); %plot planform for i = 1:(size(chords)-1) dist=0; for j = 1:n Chords(i,j); dist = (spans(i+1)-spans(i))./n; if i==1 xplotLE(i,j+1)=spans(1)+dist*(j); else xplotLE(i,1)=xplotLE(i-1,end); xplotLE(i,j+1)=xplotLE(i-1,end)+dist*(j); end fitX(1)=spans(i); fitX(2)=spans(i+1); fitY(1)=offsets(i); fitY(2)=offsets(i+1); a = polyfit(fitX, fitY,1); yplotLE(i,j) = a(1)*xplotLE(i,j)+a(2); end end plotLE=zeros(1,2); plotTE=zeros(1,2); for i=1:(size(chords)-1) for j=1:n plotLE=[plotLE; xplotLE(i,j) -yplotLE(i,j)]; plotTE=[plotTE; xplotLE(i,j) -(yplotLE(i,j)+Chords(i,j))]; end end plotTE(1,:)=[]; plotLE(1,:)=[]; end function [wakeplot]=findwake(airfoil, chords, spans, offsets, n, alt, v_mach, laminarflow, width) [q v_fps dynvisc]=findconstants(alt, v_mach); Chords = findChords(n, chords, spans); wakeplot=zeros(100,1); wakeplotX=zeros(100,1); wakeplotY=zeros(100,1); width=width/2; % for i = 1:(size(chords)-1) for i=1:1:(size(chords)-1) for j = 1:n % wakecoord(j,:)=linspace(0, Chords(i,j), n); dist = (spans(i+1)-spans(i))./n; y(j)=spans(i)+j*dist; if i==1 xplotLE(i,j+1)=spans(1)+dist*(j); else xplotLE(i,1)=xplotLE(i-1,end); xplotLE(i,j+1)=xplotLE(i-1,end)+dist*(j); end fitX(1)=spans(i); fitX(2)=spans(i+1); fitY(1)=offsets(i); fitY(2)=offsets(i+1); a = polyfit(fitX, fitY,1); yplotLE(i,j) = a(1)*xplotLE(i,j)+a(2); if spans(i)>width lam=laminarflow; else lam=laminarflow-.15; end x=linspace(0, Chords(i,j)); for k = 1:length(x) if (k <= 100*lam) delta_x=( 5.2 * x(k) )*(v_fps*x(k)/dynvisc)^-0.5; wake(k,j)= delta_x; else %if turbulent x_tr=x(100*lam); Re_xtr = v_fps * x_tr / dynvisc; delta_xtr=(5.2*x_tr)*(v_fps*x(k)/dynvisc)^-(0.5); delta_l = ((delta_xtr*((v_fps)^0.2))/(0.37*((dynvisc)^0.2)))^(1/0.8); lt = x(k) - x_tr + delta_l; Re_lt = v_fps * lt / dynvisc; delta = 0.37 * lt /(Re_lt)^0.2; wake(k,j) = delta; end wakeX(k,j)=x(k)+yplotLE(i,j); wakeY(k,j)=y(j); end for k = 1:100 wakeX(k,j)=x(k)+yplotLE(i,j); wakeY(k,j)=y(j); end wakeplot=[wakeplot, wake]; wakeplotX=[wakeplotX, wakeX]; wakeplotY=[wakeplotY, wakeY]; end end wakeplot(:,1)=[ ]; wakeplotX(:,1)=[ ]; wakeplotY(:,1)=[ ]; wakeplot(1,:)=[0]; flipud(wakeplot); end function thickness=findthickness(airfoil, Chord, xLoc) Airfoil = load(airfoil); Airfoil = Airfoil*Chord; % blah....goodnight AirfoilX=Airfoil(:,1); [min_diff, k]=min(abs(AirfoilX)); for j=1:k AirfoilYupper(j)=Airfoil(j,2); AirfoilXupper(j)=Airfoil(j,1); end for j=1:length(Airfoil)-k AirfoilYlower(j)=Airfoil(length(Airfoil)-(j-1), 2); AirfoilXlower(j)=Airfoil(length(Airfoil)-(j-1), 1); end AirfoilYupper=fliplr(AirfoilYupper); AirfoilXupper=fliplr(AirfoilXupper); AirfoilYlower=fliplr(AirfoilYlower); AirfoilXlower=fliplr(AirfoilXlower); AirfoilYlower=[0 AirfoilYlower]; AirfoilXlower=[0 AirfoilXlower]; [min_diff, k]=min(abs(AirfoilX)); [min_diff, kU]=min(abs(AirfoilXupper-xLoc)); [min_diff, kL]=min(abs(AirfoilXlower-xLoc)); if xLoc