Cannot plot the full interval
1 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
采纳的回答
Simon Chan
2021-11-28
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]
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Guidance, Navigation, and Control (GNC) 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!