Script DragPolar6.m
% Script to determine drag polar from rate of climb data.
% For the Baron 58 we can use rate of climb data from the
% Pilot's Operating Handbook (POH) page 5-30
close all
clear
disp(' '); disp('Start Here'); disp(' ');
load Baron58RCdata2a.txt % Pilot's Operating Handbook page 5-30
Hp2=Baron58RCdata2a(:,1) % Pressure altitude (ft)
OAT2=Baron58RCdata2a(:,2) % degC
hdot2=Baron58RCdata2a(:,3) % rate of climb (ft/min)
Wt=Baron58RCdata2a(:,4) % Weight (lbf)
ifig=0;
index=1:length(Hp2);
ifig=ifig+1; figure(ifig)
subplot(411);plot(index,Hp2,'x')
title('Baron 58 2-Engine Climb Data From POH page 5-30, 105 knots')
ylabel('pressure altitude (ft)')
subplot(412);plot(index,OAT2,'x')
ylabel('OAT (degC)')
subplot(413); plot(index,hdot2,'x')
ylabel('climb rate (ft/min)')
subplot(414); plot(index,Wt,'x')
ylabel('Weight (lbf)')
flapVector=0*ones(size(Hp2)); % No Flaps
% Vcalibrated=AS1Baron58(Vindicated,flap)
% function [Vequivalent,qcOverp0,qcOverp]=AS2(Vcalibrated,p)
% function [temp,press,rho,Hgeopvector]=atmosphere4(Hvector,GeometricFlag)
% function Hcalibrated=ALTBaron58(Hindicated,Vindicated,flap)
% function [VtrueKnots,M,aknots,rho,sigma]=AS3(Vequivalent,TdegR,Hpcalibrated)
Vi=105; % knots
ViVector=Vi*ones(size(Hp2)); % knots
VcVector=AS1Baron58(ViVector,flapVector); % knots
HcVector=ALTBaron58(Hp2,ViVector,flapVector); % a vector same length as Hp2
[tVector,pVector,rho,Hgeopvector]=atmosphere4(HcVector,0); % altitude is geopotential, p is a vector
R=1716.55; %ft^2/(sec^2degR)
TempDegR=9*OAT2/5+32+459.67;
rho2=pVector./(R*TempDegR); % Non-standard atmosphere density. slugs/ft^3
rho; % Use standatd atmosphere density. slugs/ft^3
%rho=rho2 % Alternatively, use rho2. slugs/ft^3
[VeVector,qcOverp0,qcOverp]=AS2(VcVector,pVector); % knots
[VtrueKnots,M,aknots,rho3,sigma]=AS3(VeVector,tVector,HcVector);
VtFtperSec=VtrueKnots*1.687808; % Convert knots to ft/sec
HdotFtperSec=hdot2/60; % convet to ft/sec
singamma=HdotFtperSec./VtFtperSec; % sine(gamma)
gamma=asin(singamma); % gamma in radians
cosgamma=cos(gamma); % cos(gamma)
% Get per engine performance
BHP=zeros(size(Hp2)); % Set up matrices of the correct size.
thrust=BHP; % Set up matrices of the correct size.
Mtip=BHP; % Set up matrices of the correct size.
BHPmax=300; % See POH page 1-11
nblades=3; % See POH page 1-11
Din=76; % We will stick to the diameter for which we have prop data for.
RPM=2700;
S=199.2; % wing area (ft^2)
B=37+10/12; % Wing span (ft)
AR=B*B/S % Aspect ratio
% Compute thrust from tabulated propeller data and gngine graphical data
for i=1:length(Hp2)
VtrueKnots(i);
HcVector(i);
RPM;
tVector(i);
nblades;
BHPmax;
Din;
[BHP(i),thrust(i),Mtip(i)]=Gthrust2(VtrueKnots(i),HcVector(i),RPM,tVector(i),nblades,BHPmax,Din);
end
% Compute thrust from propeller polar
%C=.16; m=1.40; b=-0.032 % good agreement with Gthrust2
%C=.12; m=1.38; b=-0.032 % best agreement with drag polar
%C=.12; m=1.38; b=-0.042 % changing b changes CD0
C=.12; m=1.38; b=-0.032; % Changing m spreads out the drag data that should be coincident
%m=1.76210205056895; b=-0.0884525793118286; % 3 blades, 76 in diameter
[BHP2,thrust2]=Sthrust4(VtrueKnots,HcVector,RPM,tVector,nblades,BHPmax,Din,C,m,b);
Thrust2engine=2*thrust2; % Use thrust from propeller polar
% Plot thrust comparison
ifig=ifig+1; figure(ifig)
plot(thrust,'rd'); hold on; plot(thrust2,'gx'); hold off
xlabel('Data point number');ylabel('One Engine Thrust (lbf)')
legend('Engine and Prop Tabulated data','Propeller polar')
str2=['Propeller Polar: Cp/J*J= ',num2str(m),'*Cp/J*J + ',num2str(b),', C= ',num2str(C)];
text2(.1,.1,str2)
ifig=ifig+1; figure(ifig)
plot(BHP,'rd'); hold on; plot(BHP2,'gx'); hold off
legend('Engine Tabulated data','Gagg-Ferrar rule (function PHI)')
str2=['Propeller Polar: Cp/J*J= ',num2str(m),'*Cp/J*J + ',num2str(b),', C= ',num2str(C)];
text2(.1,.1,str2)
xlabel('Data point number');ylabel('One Engine BHP (hp)')
% Assume alpha is non-zero, and iterate to find it.
Cl0=.353; Clalpharad=5.62; % for PA-28
% Iterate on alpha and gamma
alpharad=zeros(size(Hp2)); % initialize and get matrix sizes correct
singamma=(hdot2/60)./VtFtperSec;
gammarad=asin(singamma);
for i=1:20
Lift=Wt.*cos(gammarad)-Thrust2engine.*sin(alpharad);
CL=2*Lift./(rho.*VtFtperSec.*VtFtperSec*S);
dVdHpersec=dVdH(ViVector,flapVector,Hp2);
% See equation 12.11 page 122 of Kimberlin.
multiplier=(1+VtFtperSec.*dVdHpersec/32.17);
Drag=Thrust2engine.*cos(alpharad)-Wt.*singamma.*multiplier;
CD=Drag./(.5*rho.*VtFtperSec.*VtFtperSec*S);
alpharad=(CL-Cl0)/Clalpharad;
alphadeg=alpharad*57.3;
end
ifig=ifig+1; figure(ifig)
plot(CL.*CL,CD,'rx')
xlabel('CL^2')
ylabel('CD')
title('Drag Polar, climb')
CD0=.037 % Paracite drag coefficient
e=.72 % Airplane efficiency factor
CL2=0:.1:1;
k=1/(pi*AR*e)
CD2=CD0+k*CL2.*CL2;
hold on; plot(CL2.*CL2,CD2); hold off
str1=['Line is from drag polar with CD= ',num2str(CD0),' + ',num2str(k),'*CL^2'];
str2=['Propeller Polar: Cp/J*J= ',num2str(m),'*Cp/J*J + ',num2str(b)];
text2(.15,.1,str2)
text2(.15,.2,str1)
text2(.05,.9,'If the propeller polar is known, we can compute the drag polar using the red x data.')
text2(.05,.85,'Unfortunately, from this method we cannot separate CD0 from b of the propeller polar.')
text2(.05,.8,'We must use glide data to determine CD0 and perhaps k.')
text2(.05,.75,'This plot can be used to determinme m of the propeller polar and k of the drag polar.')
text2(.05,.70,'Adjust m until the CD data lies on a single straight line.')
qbar1=.5*0.00237691267925741*VeVector.*VeVector*1.687808*1.687808;
qbar2=.5*rho.*VtFtperSec.*VtFtperSec; % same as qbar1
ifig=ifig+1; figure(ifig)
subplot(411)
plot(index,rho,'x')
ylabel('density slug/ft^3')
subplot(412)
plot(index,qbar1,'x')
ylabel('qbar1')
subplot(413)
plot(index,CD,'rx')
hold on; plot(CD0+k*CL.*CL,'gx'); hold off
legend('POH','Polar')
ylabel('CD')
subplot(414)
plot(index,CL,'x')
ylabel('CL')
function [BHP,thrust]=Sthrust4(Vtrueknots,Hpress,RPM,Temp,Nblades,BHPmax,Din,C,m,b)
% function [BHP,thrust]=Sthrust4(Vtrueknots,Hpress,RPM,Temp,Nblades,BHPmax,Din,C,m,b)
% Function to compute PER ENGINE brake horsepower and thrust.
%
% Vtrueknots, Hpress, RPM, Temp, BHP, thrust can be vectors of the same size
% Vtrueknots, knots
% Hpress, ft
% RPM is engine revolutions per minute of the propeller
% Temp is the Outside air temperature (degR) of the engine.
% BHP is brake horsepower PER ENGINE. <
% thrust in lbf. per engine.
% Total thrust for two engines is twice the returned value.
% This version uses the altitude drop-off factor with C=.12
% and either 2-bladed or 3-bladed generic propellers.
% BHPmax=max BHP per engine
% Typically BHPmax=300 hp
% Din in propeller diameter in inches
% Typically for Nblades=2 for 2-bladed prop, Din=78 inches)
% for Nblades=3 for 3-bladed prop, Din=76 inches)
% C in an input (typically C=.12)
% m is the slope of the prop[eller polar
% b is the intercept of the propeller polar.
% Typical numbers
% if Nblades==2
% m=1.34710448406077; b=-0.0424066524285002; % 2 bladed propeller polar 2700 RPM
% %m=1.39222; b=-.043472; % 2 bladed propeller polar 2750RPM
% elseif Nblades==3
% m=1.76210205056895; b=-0.0884525793118286; % 3 bladed propeller polar 2700RPM
% else
% m=1.34710448406077; b=-0.0424066524285002; % default for other blade numbers
% end
% Get atmospheric pressure for the given pressure altitude.
[temp,press,rho,Hgeopvector]=atmosphere4(Hpress,0);
% Compute density for given temperature (Temp) and pressure (press).
rho=press./(1716.55*Temp);
sigma=rho/0.00237691267925741;
%C=.12; % C in an input
phi=(sigma-C)/(1-C); % altitude drop-off factor
BHP=BHPmax*phi; % Compute hp from altitude drop-off factor
Dft=Din/12;
Vtruefps=Vtrueknots*1.687808; % Convert knots to ft/sec
n=RPM/60;
Cpi=550*BHP./(rho.*n.^3*Dft^5);
J=Vtruefps/(n*Dft);
CtoverJ2=m*(Cpi./(J.*J))+b; % Use propeller polar
Cti=CtoverJ2.*J.*J;
thrust=Cti.*rho.*n.^2*Dft^4;
RESULTS OF RUNNING DragPolar6.m
Propeller polar determined from Tabulated propeller data.
RESULTS OF RUNNING DragPolar6.m
Propeller polar determined by adjusting m (and b).