Table of Contents s68

Table of Contents 1

Program Identification 2

Purpose of the Program 3

Mathematical Technique 3

Program Usage 5

Program Listing 5

Data Runs/Test Cases 6

Appendices 6

FLIGHT_ENVELOPE.m 7

NEWTON_RAPHSON.m 9

CDZERO.m 10

AK.m 10

F.m 11

ATMOS.m 12

Program Identification

Program Assignment: flight_envelope.m

Programmer: Gregory R. Whitney

Date: October 15, 2001

Requirements: Computer running MATLAB

Storage Required: 6.22KB

Special I/O Requirements: None

To plot the data, run flight_envelope.m. There is no input required.

Purpose of the Program

The purpose of the program is to plot the flight envelope for the Lear 23 at a weight of 11,000 pounds. The flight envelope is plotted for values of the altitude ranging from sea level to 45,000 ft. There are 5 curves generated by the program. The 5 curves are the stall velocity (Vstall), the maximum mach number (VM_Max), the maximum dynamic pressure velocity (VQ_max), the low speed solution (V_low), and the maximum speed solution (V_high).

Mathematical Technique

The mathematical technique of flight_envelope.m is based on the Newton-Raphson technique for solving for the roots of an equation. The equation that we need to find the roots of comes from the conditions of steady flight at full power setting. The equation is as follows:

(1)

where

T = Thrust

D = Drag

H = altitude (constant)

L = Lift (Constant)

Substituting the equation for thrrust and drag produces the next equation:

(2)

where

ρ* = Density at sea level

ρ = Density

M = Mach number

Cdo = Zero-lift drag coefficient

K = Induced drag factor

W = Weight (W=L)

V = velocity

Sw = Planform Area

The function f(V) can be solved for V numerically using the Newton-Raphson method.

The Newton Raphson Method yields

(3)

(4)

where

ΔV = 0.0001

By solving Eq.[2] for V, we get two solutions: a high speed solution and a low speed solution. As altitude changes, you have to recalculate Eq[2], Eq.[3], and Eq.[4]. There is quite a lot of computation involved therefore we have the computer do it. We can then plot the high velocity curve and the low velocity curve versus altitude.

The remaining curves for the flight envelope are the stall velocity, max mach number, and max dynamic pressure. The stall speed is given by

(5)

where

CLMax = The maximum lift coefficient (1.24)

The velocity at the maximum mach number is

(6)

where

Mmax = The maximum mach number for the aircraft (0.85)

a = Speed of Sound at the cruising altitude

Finally, the velocity at the maximum allowable dynamic pressure is

(7)

where

Qmax = The maximum dynamic pressure (300 psf)

If we plot Eq[2] and Eqs[5-7], then we will get the flight envelope for the Lear 23 at a weight of 11,000 lbs.

Program Usage

Input format: There is no input for this program

Output format: The output of the program is a plot of the flight envelope of the Lear 23 at 11,000 lbs from sea level to 45,000 ft.

Error Provisions: None

User Instructions:

1. Start MATLAB

2. Change directory to Z:\My Documents\167M\lab4

3. At the prompt, type “flight_envelope”

4. To save the figures, simply click “Save As” and choose a filename and directory

Program Listing

See the source files, starting on page 7, for the actual program codes.

Data Runs/Test Cases

The following plot shows the flight envelope for the Lear 23 (SBJ) at a weight of 11,000 pounds. This is the only output of flight_envelope.m

Appendices

This lab was done using a home PC and running Windows 2000 and MATLAB 6.1.

The remaining pages include the source code for flight_envelope.m and the other functions that were used.

FLIGHT_ENVELOPE.m

%Gregory R. Whitney

%ASE 167M 10/9/2001

%Computer Assignment 2

%This m-file constructs a level flight envelope

%There are a few routines in this program that

%Should be mentioned. The program makes use of

%atmos.m which calculates the temperature, pressure,

%density, and speed of sound of the atmosphere at a

%particular altitude. The other function to make

%note of is the newton-raphson routine. The N-R

%function iteratively calculates the roots of a function.

%Calculate the atmospheric properties using atmos.m

function flight_envelope

%Some Constants

S = 232.0; %[ft^2]

Qmax = 300; %[lbf/ft^2]

M_max = 0.85;

Cl_max = 1.24;

RHO_STAR = 7.0163E-4; %[slugs/ft^3]

W = 11000; %[lbs]

%Counter

i=1;

for h = 0:1000:45000

[TEMP,PRESS,RHO,SOS] = atmos(h);

%Calculate Stall Speed

Vstall(i) = sqrt((2*W)/(RHO*S*Cl_max));

%Calculate the Maximum Speed based on Dynamic Pressure

VQ_max(i) = sqrt((2*Qmax)/RHO);

%Calculate the maxium speed based on Max Mach Number

VM_max(i) = M_max * SOS;

%Try to get a good initial guess

if h < 25000

V_high = VQ_max(i)+150;

V_low = Vstall(i)-50;

elseif h < 35000

V_high = VQ_max(i);

V_low = Vstall(i)-50;

elseif h >=35000

V_high = VQ_max(i)+500;

V_low = Vstall(i)-50;

end

%Run the Newton Raphson Method to find the roots of Steady-Flight

V_low_root(i) = newton_raphson(V_low,SOS,RHO);

V_high_root(i) = newton_raphson(V_high,SOS,RHO);

alt(i) = h;

i = i + 1;

end

plot(V_low_root(:),alt(:),'b-')

hold

plot(V_high_root(:),alt(:),'g-')

plot(Vstall(:),alt(:),'r-')

plot(VM_max(:),alt(:),'m-')

plot(VQ_max(:),alt(:),'c-')

NEWTON_RAPHSON.m

%Gregory R. Whitney

%ASE 167M 10/9/2001

%Computer Assignment 2

%This function utilizes the Newton-Raphson

%Method of determining the roots of a function.

%This funciton already exists in MATLAB, but as

%an exercise, I have created my own function.

function [root] = newton_raphson(V,a,RHO)

dV = 10; %Just a large number so that the loop will start

delta_v = .0001;

i = 0;

while abs(dV) > delta_v & i < 100

Cdo = CDZERO(V,a);

K = AK(V,a);

f = F(V,a,RHO,Cdo,K);

f_prime = (F((V+delta_v),a,RHO,CDZERO(V+delta_v,a),AK(V+delta_v,a)) - f)/delta_v;

dV = (f/f_prime);

V = V - dV;

root = V;

i = i + 1;

end

CDZERO.m

%Gregory R. Whitney

%ASE 167M 10/9/2001

%Computer Assignment 2

%This m-file determines the value Cdo for the

%Lear jet 23. It recieves a Mach Number as

%Input and returns the value of Cdo

function [Cdo] = CDZERO(V,a)

Mach = V/a;

if Mach < 0.7

Cdo = 0.023;

elseif Mach > 0.7 & Mach < 0.9

Cdo = 0.0385 - 0.1*sqrt(0.024 - 0.6*(Mach - 0.7)^2);

elseif Mach > 0.9

Cdo = 0.0385;

End

AK.m

%Gregory R. Whitney

%ASE 167M 10/9/2001

%Computer Assignment 2

%This m-file determines the value K for the

%Lear jet 23. It recieves a Mach Number as

%Input and returns the value of K

function [K] = AK(V,a)

Mach = V/a;

if Mach < 0.7

K = 0.075;

elseif Mach > 0.7 & Mach < 0.95

K = 0.325 - sqrt(0.0625 - (Mach - 0.7)^2);

elseif Mach > 0.95

K = 0.325;

End

F.m

%Gregory R. Whitney

%ASE 167M 10/9/2001

%Computer Assignment 2

%This m-file creates the equation that will

%be used in conjucntion with newton_raphson.m

%and flight_envelope.m

%This routine defines the function used in determining

%The level flight envelope of the Lear 23 (SBJ)

function f = F(V,a,RHO,Cdo,K)

M = V/a;

RHO_STAR = 7.0163E-4; %[slugs/ft^3]

S = 232; %[ft^2]

W = 11000; %[lbf]

f = 2000*(RHO/RHO_STAR)*(0.89*(M-0.4)^2 + 0.64) - (0.5*(Cdo)*RHO*S*V^2) - (2*K*(W^2))/(RHO*S*V^2); %[lbf]

ATMOS.m

% The function atmos was created by Gregory Whitney 9/26/2001 ASE 167

% atmos requires an input of altitude (H) in ft,

% and returns four parameters

% 1. TEMP = Temperature (Rankine)

% 2. PRESS = Pressure (Lbs/ft^2)

% 3. DENS = Density (slugs/ft^3)

% 4. SOS = Speed of Sound (ft/sec)

% All of the equations for the four parameters are modeled after

% the 1962 standard atmosphere

function [TEMP,PRESS,DENS,SOS] = atmos(H)

% Definitions for Sea-Level values of Temperature, density, and pressure

RHOSL = 0.002376899341; % slugs/ft^3

TEMPSL = 518.69; % Rankine

PASL = 2116.199414; % lbs/ft^2

% Layer 1: For units, see definitions above

if H >= 0 & H < 36089.2

TEMP = TEMPSL - 0.00356616*H;

PRESS = (1.137193514E-11)*(TEMP)^5.2560613;

DENS = (RHOSL * TEMPSL / PASL)*(PRESS / TEMP);

SOS = 49.0214 * (TEMP)^.5;

% Layer 2

elseif H > 36089.2 & H < 65616.8

TEMP = 389.9901385;

PRESS = 2678.38681*exp((-4.80624E-5)*H);

DENS = (RHOSL * TEMPSL / PASL) * PRESS / TEMP;

SOS = 49.0214 * (TEMP)^.5;

% Layer 3

elseif H > 65616.8 & H < 104986.9

TEMP = 353.9901374 + 0.00054864 * H;

PRESS = (3.802075778E+90)*(TEMP)^(-34.1643987);

DENS = (RHOSL * TEMPSL / PASL)*(PRESS / TEMP);

SOS = 49.0214 * (TEMP)^.5;

% Layer 4

elseif H > 104986.9 & H < 154199.5

TEMP = 250.31033243 + 0.00153619 * H;

PRESS = (1.44217649E+33)*(TEMP)^(-12.20158686);

DENS = (RHOSL * TEMPSL / PASL) * PRESS / TEMP;

SOS = 49.0214 * (TEMP)^.5;

% Layer 5

elseif H > 154199.5 & H < 170603.7

TEMP = 487.1900543;

PRESS = (873.6400826)*exp((-3.84736009E-5)*H);

DENS = (RHOSL * TEMPSL / PASL) * PRESS / TEMP;

SOS = 49.0214 * (TEMP)^.5;

% Layer 6

elseif H > 170603.7 & H < 200131.2

TEMP = 674.3900821 - 0.00109728*H;

PRESS = (1.509860827E-46)*(TEMP)^17.08219937;

DENS = (RHOSL * TEMPSL / PASL) * PRESS / TEMP;

SOS = 49.0214 * (TEMP)^.5;

% Layer 7

elseif H > 200131.2 & H < 259186.0

TEMP = 893.9900454 - (2.19456e-3)*H;

PRESS = (7.577930948E-24)*(TEMP)^8.541099698;

DENS = (RHOSL * TEMPSL / PASL) * PRESS / TEMP;

SOS = 49.0214 * (TEMP)^.5;

% Layer 8

elseif H > 259186.0 & H <= 2000000

TEMP = 325.1908172;

PRESS = (6.669501806E+4)*exp((-5.763986782E-5)*H);

DENS = (RHOSL * TEMPSL / PASL) * PRESS / TEMP;

SOS = 49.0214 * (TEMP)^.5;

else

disp('That value of H is unacceptable');

TEMP = 'Not Defined';

PRESS = 'Not Defined';

DENS = 'Not Defined';

SOS = 'Not Defined';

end