asc/senior-design
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
main
senior-design/ronald/Ronald.m
asc 5bacd30bad reorganize file structure
2010-04-17 12:00:00 -05:00

521 lines
18 KiB
Matlab
Executable File
Raw Permalink Blame History

%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)
%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/superGrimace Outputs';
%-----------------------------------------------------------------------]
n=50;
cd(paramLOCATION);
[variablenames, variablevalues]=textread(paramfile, '%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
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=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=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=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=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(pranformLOCATION)
filenameplan=regexprep(pranformfile, '.plan', 'ronaldplus.tiff');
filenameoutput=regexprep(pranformfile, '.plan', 'ronaldplus.txt');
plottitle=regexprep(pranformfile, '.plan', '');
fid = fopen(filenameoutput, 'w');
fprintf(plottitle);
fprintf('---------------------------------------------------\n');
fprintf(fid, plottitle);
fprintf(fid,'---------------------------------------------------\n');
fprintf('Wingspan: %f, AR: %f, Planform Area: %f, Wetted Area: %f \n', wingspan, AR, sref, swet);
fprintf(fid, 'Wingspan: %f, AR: %f, Planform Area: %f, Wetted Area: %f \n', wingspan, AR, sref, swet);
fprintf('\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----------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('\n----------Cruise---------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_cruise, cd0_cruise, cLreq_cruise, cD_cruise, D_cruise, cM_cruise, eloverdee_cruise);
fprintf(fid, '\n----------Cruise---------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_cruise, cd0_cruise, cLreq_cruise, cD_cruise, D_cruise, cM_cruise, eloverdee_cruise);
fprintf('\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----------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('\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, '\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);
fclose(fid);
cd(rootLOCATION)
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=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<AirfoilXupper(kU)
xFitU(1)=AirfoilXupper(kU-1);
xFitU(2)=AirfoilXupper(kU);
yFitU(1)=AirfoilYupper(kU-1);
yFitU(2)=AirfoilYupper(kU);
else
xFitU(1)=AirfoilXupper(kU);
xFitU(2)=AirfoilXupper(kU+1);
yFitU(1)=AirfoilYupper(kU);
yFitU(2)=AirfoilYupper(kU+1);
end
if xLoc<AirfoilXlower(kL)
xFitL(1)=AirfoilXlower(kL-1);
xFitL(2)=AirfoilXlower(kL);
yFitL(1)=AirfoilYlower(kL-1);
yFitL(2)=AirfoilYlower(kL);
else
xFitL(1)=AirfoilXlower(kL);
xFitL(2)=AirfoilXlower(kL+1);
yFitL(1)=AirfoilYlower(kL);
yFitL(2)=AirfoilYlower(kL+1);
end
aU = polyfit(xFitU, yFitU,1);
aL = polyfit(xFitL, yFitL,1);
% yplotLE(i,j) = a(1)*xplotLE(i,j)+a(2);
% thickness = AirfoilYupper(xLoc_indu)-AirfoilYlower(xLoc_indl);
thickness = (aU(1)*xLoc+aU(2))-(aL(1)*xLoc+aL(2));
end
Reference in New Issue View Git Blame Copy Permalink