Cannot plot the full interval

7 views (last 30 days)
function [] = FuelJetPlot
H_climb = 0:500:35000; % Altitude [ft]
for k1 = 1:length(H_climb)
% CLIMB PHASE
% Thrust calculation
C_Tc_1 = .75979E+06;
C_Tc_2 = .52423E+05;
C_Tc_3 = .40968E-10;
Thr_jet_climb_ISA (k1)= C_Tc_1* (1-(H_climb(k1)/C_Tc_2)+C_Tc_3*H_climb(k1)^2); % Maximum climb and take-off thrust for jet engine [N]
Thr (k1)= Thr_jet_climb_ISA(k1)/1000; % [N] -> [kN]
% Vtas calculation
Vcl_1 = 335; % Standard calibrated airspeed [kt]
Vcl_2 = 172.3; % Standard calibrated airspeed [kt] -> [m/s] (To find the Mach transition altitude)
Vcl_2_in_knots = 335; % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cl = 0.86; % Standard calibrated airspeed [kt]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
g0 = 9.80665; % Gravitational acceleration [m/s2]
a0= 340.294; % Speed of Sound [m/s]
CV_min = 1.3; % Minimum speed coefficient
Vstall_TO = .14200E+03; % Stall speed at take-off [KCAS]
Vd_CL1 = 5; % Climb speed increment below 1500 ft (jet)
Vd_CL2 = 10; % Climb speed increment below 3000 ft (jet)
Vd_CL3 = 30; % Climb speed increment below 4000 ft (jet)
Vd_CL4 = 60; % Climb speed increment below 5000 ft (jet)
Vd_CL5 = 80; % Climb speed increment below 6000 ft (jet)
% CLIMB SPEEDS
% 1) Jet aircraft
CAS_climb = Vcl_2;
Mach_climb = M_cl;
delta_trans = (((1+((K-1)/2)*(CAS_climb/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_climb^2)^(K/(K-1)))-1); % Pressure ratio at the transition altitude
teta_trans = delta_trans ^ (-Bt*R/g0); % Temperature ratio at the transition altitude
H_p_trans_climb = (1000/0.348/6.5)*(T0*(1-teta_trans)); % Transition altitude for climb [ft]
if (0<=H_climb(k1)&&H_climb(k1)<=1499)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;
elseif (1500<=H_climb(k1)&&H_climb(k1)<=2999)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
elseif (3000<=H_climb(k1)&&H_climb(k1)<=3999)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL3;
elseif (4000<=H_climb(k1)&&H_climb(k1)<=4999)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL4;
elseif (5000<=H_climb(k1)&&H_climb(k1)<=5999)
Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL5;
elseif (6000<=H_climb(k1)&&H_climb(k1)<=9999)
Vnom_climb_jet(k1) = min(Vcl_1,250);
elseif (10000<=H_climb(k1)&&H_climb(k1)<=H_p_trans_climb)
Vnom_climb_jet(k1) = Vcl_2_in_knots;
elseif (H_p_trans_climb<H_climb(k1))
Vnom_climb_jet (k1)= M_cl;
end
Vcas(k1) = Vnom_climb_jet(k1)* 0.514; % [kn] -> [m/s]
H_climb (k1)= H_climb(k1) * 0.3048; % [feet] -> [m]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0; % Value of the real temperature T in ISA conditions [K]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
P0 = 101325; % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665; % Gravitational acceleration [m/s2]
p0 = 1.225; % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;
T(k1) = T0 + deltaT + Bt * H_climb(k1); % Temperature [K]
P (k1)= P0*((T(k1)-deltaT)/T0).^((-g0)/(Bt*R)); % Pressure [Pa]
p (k1)= P(k1) ./ (R*T(k1)); % Density [kg/m^3]
Vtas(k1) = (2*P(k1)/visc/p(k1)*((1 + P0/P(k1)*((1 + visc*p0*Vcas(k1)*Vcas(k1)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
Vtas (k1)= Vtas(k1) * 1.94; % [m/s] -> [knots]
% Fuel consumption calculation
C_f1 = .58259E+00;
C_f2 = .12839E+04;
n_jet(k1) = C_f1*(1+Vtas(k1)/C_f2); % Thrust specific fuel consumption for jet engine [kg/(min·kN)]
f_nom_jet(k1) = n_jet(k1) * Thr(k1); % Nominal fuel flow for jet engine [kg/min], can be used in climb phase
end
H_cruise = 0:500:35000; % Altitude[ft]
for k2 = 1:length(H_cruise)
% CRUISE PHASE
% Speed Calculation
% CRUISE SPEEDS
Vcr_1 = 320; % Standard calibrated airspeed [kt]
Vcr_2 = 164.62; % Standard calibrated airspeed [kt] -> [m/s] (To find the Mach transition altitude)
Vcr_2_in_knots = 320; % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cr = 0.84; % Standard calibrated airspeed [kt]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
%deltaT = 0; % Value of the real temperature T in ISA conditions [K]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
g0 = 9.80665; % Gravitational acceleration [m/s2]
a0= 340.294; % Speed of Sound [m/s]
% 1) Jet Aircraft
CAS_cruise = Vcr_2;
Mach_cruise = M_cr;
delta_trans = (((1+((K-1)/2)*(CAS_cruise/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_cruise^2)^(K/(K-1)))-1); % Pressure ratio at the transition altitude
teta_trans = delta_trans ^ (-Bt*R/g0); % Temperature ratio at the transition altitude
H_p_trans_cruise = (1000/0.348/6.5)*(T0*(1-teta_trans)); % Transition altitude for cruise [ft]
if (0<=H_cruise(k2)&&H_cruise(k2)<=2999)
Vnom_cruise_jet(k2) = min(Vcr_1,170);
elseif (3000<=H_cruise(k2)&&H_cruise(k2)<=5999)
Vnom_cruise_jet(k2) = min(Vcr_1,220);
elseif (6000<=H_cruise(k2)&&H_cruise(k2)<=13999)
Vnom_cruise_jet(k2) = min(Vcr_1,250);
elseif (14000<=H_cruise(k2)&&H_cruise(k2)<=H_p_trans_cruise)
Vnom_cruise_jet(k2) = Vcr_2_in_knots;
elseif (H_p_trans_cruise < H_cruise(k2))
Vnom_cruise_jet(k2) = M_cr;
end
Vcas (k2)= Vnom_cruise_jet(k2) * 0.514; %[knots] -> [m/s]
% Thrust Calculation
L = .44225E+06; % .44225E+03 tons = .44225E+06 kg and W = L.
S = .51097E+03; % Surface Area [m^2]
H_cruise(k2) = H_cruise(k2) * 0.3048; % [feet] -> [m]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0; % Value of the real temperature T in ISA conditions [K]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
P0 = 101325; % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665; % Gravitational acceleration [m/s2]
p0 = 1.225; % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;
T (k2)= T0 + deltaT + Bt * H_cruise(k2); % Temperature [K]
P (k2)= P0*((T(k2)-deltaT)/T0).^((-g0)/(Bt*R)); % Pressure [Pa]
p (k2)= P(k2) ./ (R*T(k2)); % Density [kg/m^3]
Vtas(k2) = (2*P(k2)/visc/p(k2)*((1 + P0/P(k2)*((1 + visc*p0*Vcas(k2)*Vcas(k2)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
Cl(k2) = 2*L / (p(k2)*Vtas(k2)*Vtas(k2)*S); % Lift coefficient
C_D0cr = .25669E-01;
C_D2cr = .39082E-01;
C_D(k2) = C_D0cr + C_D2cr * (Cl(k2)).^2; % Drag Coefficient
D (k2)= 1/2 * p(k2) * Vtas(k2)* Vtas(k2) * S * C_D(k2); % Drag Force [N]
Thr(k2) = D(k2)/1000; % Cruise drag (Thrust) [N] -> [kN]
Vtas(k2) = Vtas(k2) * 1.94; % [m/s] -> [knots]
% Fuel Consumption Calculation
C_fcr = .89625E+00;
C_f1 = .58259E+00;
C_f2 = .12839E+04;
n_jet(k2) = C_f1*(1+Vtas(k2)/C_f2); % Thrust specific fuel consumption for jet engine [kg/(min·kN)]
f_cr_jet(k2) = n_jet(k2) * Thr(k2) * C_fcr; % Cruise fuel flow [kg/min]
end
H_descent = 0:500:35000; % Altitude [ft]
for k3 = 1:length(H_descent)
% DESCENT PHASE
% The minimum fuel flow, fmin [kg/min], corresponds to idle thrust descent conditions for both jet and turboprop engines
% so Vtas and Thrust calculations are unnecessary.
% Fuel Consumption Calculation
C_f3 = .45380E+02;
C_f4 = .72929E+05;
f_min_jet(k3) = C_f3 * (1-H_descent(k3)/C_f4); % Minimum fuel flow [kg/min]
end
figure(1)
plot(H_climb,f_nom_jet);
xlabel('Altitude [ft]');
ylabel('Nominal fuel flow for jet engine [kg/min]');
title('H vs f for climb phase');
figure(2)
plot(H_cruise,f_cr_jet);
xlabel('Altitude [ft]');
ylabel('Cruise fuel flow for jet engine [kg/min]');
title('H vs f for cruise phase');
figure(3)
plot(H_descent,f_min_jet);
xlabel('Altitude [ft]');
ylabel('Minimum fuel flow for jet engine [kg/min]');
title('H vs f for descent phase');
end
Hi, the code plots fuel consumption for climb, cruise and descent phases of a jet engine aircraft with respect to altitude. Whenever I run the ccode, Fig1 and Fig2's altitude stops at H=x=10668 ft but I want them to end at H=x=35,000 ft. How to fix this problem?
And if you see any mistakes in cruise phase part of the code, please inform me because its plot looks wrong.
Thanks.

Accepted Answer

Simon Chan
Simon Chan on 28 Nov 2021
You converted H_climb & H_cruise from [feet] to [m] by multiplying the value by 0.3048.
However, there is no such conversion for variable H_descent.
You can just remove the conversion and get your required figures.
%H_climb (k1)= H_climb(k1) * 0.3048; % [feet] -> [m]
%
%H_cruise(k2) = H_cruise(k2) * 0.3048; % [feet] -> [m]
  1 Comment
Turgut Ataseven
Turgut Ataseven on 28 Nov 2021
Thanks for pointing the mistake out!
Those conversions were necessary for climb and cruise phases, so I redefined altitudes at the last row of climb and cruise phases as
H_climb = 0:500:35000; % Redefining altitude[ft] because it was converted to [m] above.
H_cruise = 0:500:35000; % Redefining altitude[ft] because it was converted to [m] above.
Hope it won't change the calculations.

Sign in to comment.

More Answers (0)

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!