Seaplane Camper Project - Airfoil & Geometry

function AirfoilandGeometry

%---------------------------------------------------------------------%

ClCRZ = 0.5; %assumed design lift

%coeff (Raymer, pg. 48)

ATMOS = xlsread('ATMOS.xlsx','A3:N85'); %std day atmosphere

ATMOSh = ATMOS(1:83,1); %std day altitude in feet

ATMOSP = ATMOS(1:83,6); %std day pressure inHg

hCRZ = 12500; %ft, design crz altitude

%to avoid prs-ized cabin

PCRZ = interp1(ATMOSh,ATMOSP,hCRZ); %interpolated cruise

%pressure in inches Hg

ATMOST = ATMOS(1:83,3); %std day temperature degR

TCRZ = interp1(ATMOSh,ATMOST,hCRZ); %interpoloated cruise

%temperature, degR

rhoCRZ = 0.041206*(PCRZ/TCRZ); %sl/ft^3, cruise density

%(P&W 79500, pg.71)

VCRZ = 220.0050; %design crz spd in fps

qCRZ = rhoCRZ * VCRZ^2/2; %lb/ft^2, cruise dyn prs

%(Raymer, pg.44)

WSCRZ = ClCRZ * qCRZ; %lb/ft^2 wing loading for

%design lift coefficient

%(Raymer, pg. 48)

%---------------------------------------------------------------------%

% Airfoil Selection Mission Profile: %

% Simple Cruise %

% Wing: NACA 23015 %

% Max t/c 15% Crz %

% Max L/D 90.9524 2-----3 %

% AOA 10.00 Clb / \ %

% / 4 Ltr %

% Tail: NACA 0006 TO / \ %

% Max t/c 6% 0----1 5 Lnd %

% Max L/D 56.6287 %

% AOA 4.50 (Raymer, pg. 19) %

%---------------------------------------------------------------------%

awRAW = xlsread('NACA23015.xlsx','A3:D20'); %wing airfoil coords

awUSTN = awRAW(1:18,1); %upper station coords

awUORD = awRAW(1:18,2); %upper ordinate coords

awLSTN = flipud( awRAW(1:18,3) ); %lower station coords

awLORD = flipud( awRAW(1:18,4) ); %lower ordinate coords

awdata = xlsread('NACA23015.xlsx','F12:L167'); %wing airfoil data

awdataAOA = awdata(1:156,1); %angle of attack

awdataCL = awdata(1:156,2); %lift coefficient

awdataCD = awdata(1:156,3); %drag coefficient

awdataCM = awdata(1:156,5); %moment coefficient

figure('name','NACA 23015','numbertitle','off') %wing airfoil plot

hold on

plot(awUSTN,awUORD,'color','blue')

plot(awLSTN,awLORD,'color','blue')

axis equal

grid on

figure('name','NACA 23015','numbertitle','off') %lift curve

hold on

plot(awdataAOA,awdataCL,'color','blue')

plot(awdataAOA,awdataCM,'color','red')

xlabel('AOA')

legend({'CL','CM'},'location','northwest')

grid on

figure('name','NACA 23015','numbertitle','off') %drag polar

hold on

plot(awdataCL,awdataCD,'color','blue')

plot(awdataCL,awdataCM,'color','red')

xlabel('CL')

legend({'CD','CM'},'location','southeast')

grid on

atRAW = xlsread('NACA0006.xlsx','A3:D20'); %tail airfoil coords

atUSTN = atRAW(1:18,1);

atUORD = atRAW(1:18,2);

atLSTN = flipud( atRAW(1:18,3) );

atLORD = flipud( atRAW(1:18,4) );

atdata = xlsread('NACA0006.xlsx','F12:L78'); %tail airfoil data

atdataAOA = atdata(1:67,1);

atdataCL = atdata(1:67,2);

atdataCD = atdata(1:67,3);

atdataCM = atdata(1:67,5);

figure('name','NACA 0006','numbertitle','off') %tail airfoil plot

hold on

plot(atUSTN,atUORD,'color','blue')

plot(atLSTN,atLORD,'color','blue')

axis equal

grid on

figure('name','NACA 0006','numbertitle','off') %lift curve

hold on

plot(atdataAOA,atdataCL,'color','blue')

plot(atdataAOA,atdataCM,'color','red')

xlabel('AOA')

legend({'CL','CM'},'location','northwest')

grid on

figure('name','NACA 0006','numbertitle','off') %drag polar

hold on

plot(atdataCL,atdataCD,'color','blue')

plot(atdataCL,atdataCM,'color','red')

xlabel('CL')

legend({'CD','CM'},'location','northwest')

grid on

%---------------------------------------------------------------------%

degTOrad = 0.017453; %degrees to radians

radTOdeg = 57.296; %radians to degrees

%---------------------------------------------------------------------%

W0 = 9738.5275; %takeoff weight

W1overW0 = 0.970; %est warmup and takeoff weight fraction

W2overW1 = 0.985; %est climb wt fraction (Raymer, pg. 20)

W2overW0 = W1overW0 * W2overW1; %weight fraction at beginning of cruise

W2 = W0 * W2overW0; %weight at beginning of cruise

S = W2 / WSCRZ; %wing area in ft^2

A = 8; %statistical aspect ratio

%(Raymer, pg. 59)

b = sqrt( A * S ); %wing span in ft

tpr = 0.4; %design taper ratio, ideal for unswept

%(Raymer, pg. 63)

Croot = 2*S / (b*(1+tpr)); %root chord length in ft

Ctip = tpr * Croot; %tip chrd len ft (Raymer, pg. 55)

swpC4deg = 0; %1/4 chrd swp deg (Raymer, pg. 60)

swpC4rad = swpC4deg * degTOrad; %rad

swpLErad = atan( tan(swpC4rad) + ... %leading edge sweep angle in

(1-tpr)/(A*(1+tpr)) ); %radians (Raymer, pg. 55)

swpLEdeg = swpLErad * radTOdeg; %degrees

cmac = 2/3 * Croot * ... %mean aerodynamic chord in ft

(1+tpr+tpr^2)/(1+tpr); %(Raymer, pg. 56)

Ybar = b/6 * (1+2*tpr)/(1+tpr); %loc of cmac from centerline in ft

cac = 0.25 * cmac; %subsonic aerodynamic center in ft

twistdeg = 0; %design wing twist, deg

%(Raymer, pg. 66)

twistrad = twistdeg * degTOrad; %rad

ideg = 2; %design wing incidence in degrees

%(Raymer, pg. 66)

irad = ideg * degTOrad; %rad

dihdeg = 1; %design wing dihedral in degrees

%(Raymer, pg. 68)

dihrad = dihdeg * degTOrad; %rad

%---------------------------------------------------------------------%

disp(' ');

disp('AIRFOIL & GEOMETRY:');

disp(' ');

disp(['Cruise Altitude = ',num2str(hCRZ),' ft']);

disp(['Cruise W/S = ',num2str(WSCRZ),' lb/ft^2']);

disp(' ');

disp(['Cruise Weight = ',num2str(W2),' lbs']);

disp(['Wing Area = ',num2str(S),' ft^2']);

disp(['Wing Span = ',num2str(b),' ft']);

disp(['Taper Ratio = ',num2str(tpr)]);

disp(['Root Chord Length = ',num2str(Croot),' ft']);

disp(['Tip Chord Length = ',num2str(Ctip),' ft']);

disp(' ');

disp(['Qrtr Chord Sweep = ',num2str(swpC4deg),' deg']);

disp(['Leading Edge Sweep = ',num2str(swpLEdeg),' deg']);

disp(['Mean Aero Chord = ',num2str(cmac),' ft']);

disp(['MAC Ybar = ',num2str(Ybar),' ft']);

disp(['Aerodynamic Center = ',num2str(cac),' ft']);

disp(' ');

disp(['Wing Twist = ',num2str(twistdeg),' deg']);

disp(['Wing Incidence = ',num2str(ideg),' deg']);

disp(['Wing Dihedral = ',num2str(dihdeg),' deg']);

disp(' ');

%---------------------------------------------------------------------%

% ----------->| STN _ %

% _ __-------_____ /__________________ ^ %

% ^ / 1 --> -----___ / | / | %

% | ( _____\ | / 30deg | %

% | `---________-------- /______________|/ _ _ z %

% ORD <-- 2 / ------->| y %

% %

% Airfoil Coords: Hoerner Wing tip: %

% STN - Station Lower surface is canted %

% ORD - Ordinate approx 30 deg to horiz. %

% %

% z = m*y + b %

% min(awtipsZraw) = m*(b/2) + bint %

% bint = min(awtipsZraw) - m*(b/2) %

% awtipsYhrn = (awtipsZraw - bint)/m %

%---------------------------------------------------------------------%

awrootXraw = cat(1,awUSTN,awLSTN) * Croot; %raw wing root af x coords

awrootZraw = cat(1,awUORD,awLORD) * Croot; %raw wing root af z coords

awrootYraw = zeros( length(awUSTN) + ...

length(awLSTN), 1); %raw wing root af y coords

awtipsXraw = cat(1,awUSTN,awLSTN) * Ctip + ...

b/2 * tan(swpLErad); %raw wng tips x coords

awtipsZraw = cat(1,awUORD,awLORD) * Ctip + ...

b/2 * tan(dihrad); %raw wng tips z coords

awtipsYraw = zeros( length(awUSTN) + ...

length(awLSTN), 1) + b/2; %raw wng tips y coords

HOERNdeg = 30; %Hoerner (Raymer, pg. 72)

HOERNrad = HOERNdeg * degTOrad; %radians

m = tan(HOERNrad); %slope of Hoerner

bint = min(awtipsZraw) - m*(b/2); %Hoerner origin intercept

awtipsYhrn = (awtipsZraw - bint)/m; %Hoerner y coords

awrootXrot1 = awrootXraw - 0.25*Croot; %translate a/c to origin

awrootXrot2 = awrootXrot1*cos(irad) + ...

awrootZraw*sin(irad); %rotate airfoil

awrootXrot = awrootXrot2/cos(irad) + ...

0.25*Croot; %translate a/c back

awrootZrot1 = awrootXraw - 0.25*Croot; %translate a/c to origin

awrootZrot2 = -awrootZrot1*sin(irad) + ...

awrootZraw*cos(irad); %rotate airfoil

awrootZrot = awrootZrot2/cos(irad); %translate a/c back

awtipsXrot1 = awtipsXraw - ...

min(awtipsXraw) - 0.25*Ctip; %trans a/c to orgn

awtipsXrot2 = awtipsXrot1 * cos(twistrad+irad) + ...

(awtipsZraw - b/2*tan(dihrad)) * ...

sin(twistrad+irad); %rotate airfoil

awtipsXrot3 = awtipsXrot2 / cos(twistrad+irad); %unshrink from rot

awtipsXrot = awtipsXrot3 + ...

min(awtipsXraw) + 0.25*Ctip; %trans a/c back

awtipsZrot1 = awtipsXraw - ...

min(awtipsXraw) - 0.25*Ctip; %trans a/c to org

awtipsZrot2 = -awtipsZrot1 * sin(twistrad+irad) + ...

(awtipsZraw - b/2*tan(dihrad)) * ...

cos(twistrad+irad); %rotate airfoil

awtipsZrot3 = awtipsZrot2 / cos(twistrad+irad); %unshrink frm rot

awtipsZrot = awtipsZrot3 + b/2*tan(dihrad); %trans a/c back

%---------------------------------------------------------------------%

wingX = zeros(10*length(awUSTN),1); %initialize plot coords size

wingY = zeros(10*length(awUSTN),1);

wingZ = zeros(10*length(awUSTN),1);

i = 1; %wing coordinate position number

n = 1; %airfoil data position number

while n < length(awrootXraw)

wingX(i+0) = awtipsXrot(n+0);

wingX(i+1) = awtipsXrot(n+0);

wingX(i+2) = awrootXrot(n+0);

wingX(i+3) = awtipsXrot(n+0);

wingX(i+4) = awtipsXrot(n+0);

wingX(i+5) = awtipsXrot(n+1);

wingX(i+6) = awtipsXrot(n+1);

wingX(i+7) = awrootXrot(n+1);

wingX(i+8) = awtipsXrot(n+1);

wingX(i+9) = awtipsXrot(n+1);

wingY(i+0) = -awtipsYhrn(n+0);

wingY(i+1) = -awtipsYraw(n+0);

wingY(i+2) = awrootYraw(n+0);

wingY(i+3) = awtipsYraw(n+0);

wingY(i+4) = awtipsYhrn(n+0);

wingY(i+5) = awtipsYhrn(n+1);

wingY(i+6) = awtipsYraw(n+1);

wingY(i+7) = awrootYraw(n+1);

wingY(i+8) = -awtipsYraw(n+1);

wingY(i+9) = -awtipsYhrn(n+1);

wingZ(i+0) = awtipsZrot(n+0);

wingZ(i+1) = awtipsZrot(n+0);

wingZ(i+2) = awrootZrot(n+0);

wingZ(i+3) = awtipsZrot(n+0);

wingZ(i+4) = awtipsZrot(n+0);

wingZ(i+5) = awtipsZrot(n+1);

wingZ(i+6) = awtipsZrot(n+1);

wingZ(i+7) = awrootZrot(n+1);

wingZ(i+8) = awtipsZrot(n+1);

wingZ(i+9) = awtipsZrot(n+1);

i = i + 10; %increment count for next pass

n = n + 2; %increment count for airfoil data

end

figure('name','Perspective','numbertitle','off')

view(-45,15.264)

camproj('perspective')

hold on

plot3(awtipsXrot,-awtipsYhrn,awtipsZrot,'color','blue')

plot3(awtipsXrot,-awtipsYraw,awtipsZrot,'color','blue')

plot3(awrootXrot,awrootYraw,awrootZrot,'color','blue')

plot3(awtipsXrot,awtipsYraw,awtipsZrot,'color','blue')

plot3(awtipsXrot,awtipsYhrn,awtipsZrot,'color','blue')

plot3(wingX,wingY,wingZ,'color','blue')

axis equal

axis off

figure('name','Three-View','numbertitle','off')

subplot(2,2,1)

view(0,0) %side view

hold on

plot3(awtipsXrot,-awtipsYhrn,awtipsZrot,'color','blue')

plot3(awtipsXrot,-awtipsYraw,awtipsZrot,'color','blue')

plot3(awrootXrot,awrootYraw,awrootZrot,'color','blue')

plot3(awtipsXrot,awtipsYraw,awtipsZrot,'color','blue')

plot3(awtipsXrot,awtipsYhrn,awtipsZrot,'color','blue')

plot3(wingX,wingY,wingZ,'color','blue')

axis equal

grid on

subplot(2,2,2)

view(-90,0) %head-on view

hold on

plot3(awtipsXrot,-awtipsYhrn,awtipsZrot,'color','blue')

plot3(awtipsXrot,-awtipsYraw,awtipsZrot,'color','blue')

plot3(awrootXrot,awrootYraw,awrootZrot,'color','blue')

plot3(awtipsXrot,awtipsYraw,awtipsZrot,'color','blue')

plot3(awtipsXrot,awtipsYhrn,awtipsZrot,'color','blue')

plot3(wingX,wingY,wingZ,'color','blue')

axis equal

grid on

subplot(2,2,3)

view(0,90) %top view

hold on

plot3(awtipsXrot,-awtipsYhrn,awtipsZrot,'color','blue')

plot3(awtipsXrot,-awtipsYraw,awtipsZrot,'color','blue')

plot3(awrootXrot,awrootYraw,awrootZrot,'color','blue')

plot3(awtipsXrot,awtipsYraw,awtipsZrot,'color','blue')

plot3(awtipsXrot,awtipsYhrn,awtipsZrot,'color','blue')

plot3(wingX,wingY,wingZ,'color','blue')

axis equal

grid on

subplot(2,2,4)

view(-45,35.264) %isometric view

hold on

plot3(awtipsXrot,-awtipsYhrn,awtipsZrot,'color','blue')

plot3(awtipsXrot,-awtipsYraw,awtipsZrot,'color','blue')

plot3(awrootXrot,awrootYraw,awrootZrot,'color','blue')

plot3(awtipsXrot,awtipsYraw,awtipsZrot,'color','blue')

plot3(awtipsXrot,awtipsYhrn,awtipsZrot,'color','blue')

plot3(wingX,wingY,wingZ,'color','blue')

axis equal

grid on

%---------------------------------------------------------------------%

OUTPUT