How to deduct obtained fuel value from initial weight

Hi, the following code is the initial calculation, I have obtained the fuel consumed to climb at each altitude (53 x 1 results) (the final code), I want to deduct this fuel consumption from the weight at every altitude (the weight is denoted is W) because weight reduced as the aircraft fly.

May I know how do I code it? I heard can use for-loop with different iteration? (because I have used a for-loop for other calculation before getting onto the coding below)

Please help. Very much appreciate :)

%% Coefficient of Lift
CDO = 0.023; % Zero-lift drag coefficient
K = 0.044; % Lift-induced drag factor
MTOW =  79010; % Maximum take-off mass [kg]
W = MTOW*9.81; % Weight in N

A= Thrust_22K/W;
B = sqrt((A.^2)+(12*CDO*K));

Cl = (6*CDO)./(A+B);

%% Drag Produced 
S = 125; % Wing area [m^2]
C = CDO + (K*(Cl.^2));
D = 0.5*Density_altitude.*(TAS.^2)*S.*C; % [N]

%% Power Available
P_available = TAS.*Thrust_22K;

%% Power Required
P_required = TAS.*D;

%% Excess Power 
P_excess = P_available - P_required; 

%% Rate of Climb (ROC)

ROC = (TAS.*(Thrust_22K - D))./W; % [m/s]

if any(ROC <=1.5)
    disp('Service Ceiling.')
else 
    disp ('Absolute Ceiling.')
end 

%% Climb Angle 
X = (Thrust_22K - D)./W;
Climb_angle = asind( X ); % [Degree]

%% Time Taken to Climb 
Time = Altitude./ROC; % [sec]

%% Mass Flow Rate
m_dot_f = TSFC_22K.*Thrust_22K; % [kg/s]

%% Fuel Consumed to Climb 
Fuel = m_dot_f.*Time;

[Dyuman Joshi]

Is the variable Altitude 53x1?

Do you want to iteratively calculate the weight as the plane ascends? Say, may be every x meters or for values in Altitude.

[Ke Yeun Yong]

The altitude is also 53x1. There is an initial weight which I used to calculate the fuel consumed (fuel results is 53x1). The weight is constant in the above coding. Now I have to go in loop to deduct the fuel consumed from every change of weight with altitude.

Yes, I want to iteratively calculate the weight as plane ascends in regards to the 53x1 altitude and fuel consumed

[William Rose]

@Ke Yeun Yong,

[edit: added TSFC_22K to list]

Please run your code in the window before you post it, to identify and fix simple errors. For example, @Dyuman Joshi noted that Altitude is not defined in your script. Other undefined variables are Thrust_22K, TAS, Density_altitude, TSFC_22K. When you add them, please include the units in the comments.

[William Rose]

As @Dyuman Joshi pointed out, it is not clear what quantities are vectors and what quatities are scalars.

[Dyuman Joshi]

It will be better if you can attach your whole code including values for variables, as has been mentioned before.

Answer 1

[Torsten]

Most probably something like this. I don't know where you use some of the defined variables - they seem to be superfluous.

%% Coefficient of Lift
CDO = 0.023; % Zero-lift drag coefficient
K = 0.044; % Lift-induced drag factor
W(1) = MTOW*9.81; % Weight in N;
Time(1) = 0.0;
Fuel(1) = 0.0;
for i = 1:numel(Altitude)-1
  A = Thrust_22K/W(i);
  B = sqrt((A.^2)+(12*CDO*K));
  Cl = (6*CDO)./(A+B);
  %% Drag Produced 
  S = 125; % Wing area [m^2]
  C = CDO + (K*(Cl.^2));
  D = 0.5*(Density_altitude(i)+Density_altitude(i+1))/2.*(TAS.^2)*S.*C; % [N]
  %% Power Available
  P_available = TAS.*Thrust_22K;
  %% Power Required
  P_required = TAS.*D;
  %% Excess Power 
  P_excess = P_available - P_required; 
  %% Rate of Climb (ROC)
  ROC = (TAS.*(Thrust_22K - D))./W(i); % [m/s]
  if any(ROC <=1.5)
    disp('Service Ceiling.')
  else 
    disp ('Absolute Ceiling.')
  end 
  %% Climb Angle 
  X = (Thrust_22K - D)./W(i);
  Climb_angle = asind( X ); % [Degree]
  %% Time Taken to Climb 
  Time(i+1) = (Altitude(i+1)-Altitude(i))./ROC; % [sec]
  %% Mass Flow Rate
  m_dot_f = TSFC_22K.*Thrust_22K; % [kg/s]
  %% Fuel Consumed to Climb 
  Fuel(i+1) = m_dot_f.*Time(i+1);
  W(i+1) = W(i) - Fuel(i+1);
end

[William Rose]

@Ke Yeun Yong,

@Torsten recognized that you start with a pre-defined vector of target altitudes before entering the for loop. And he recognized that time is not cumulative (as I was thinking) but rather is a vector of time increments from one altitude to the next. Fuel is also a vector, each element of which is the incremental fuel needed to climb from one altitude to the next. @Torsten is very good at figuring out these things.

Obviously Density_altitude is a vector of air densities [kg/m^3] at the various altitudes. I would rename that variable, because the term "density altitude" has a well established meaning in aviaiton, and it is an altitude (units of length), not a density.

[Torsten]

@William Rose

I took the freedom to make "Time(i+1)" and "Fuel(i+1)" incremental values (i.e. from Altitude(i) to Altitude(i+1)). Seemed reasonable in order to deduce W(i+1) from W(i).

sum(Fuel) should then give the total fuel consumption from Altitude(1) to Altitude(end) and sum(Time) should give the total time needed to reach Altitude(end).

[Ke Yeun Yong]

@Torsten@William Rose the following is the whole code.

% Clear all previous variables created in the workspace

clear all

% Clear all previous figures whose handles are not hidden

close all

% Clears all previous text from the Command Window

clc

%% Engine Data

load data.txt

Altitude = data(:,1); % m

Thrust_22K = (data(:,2)).*2; % N

TSFC_22K_hour = data(:,3); % per hour

TSFC_22K = TSFC_22K_hour/3600; % per sec

%% Density

To = 288.15; % Temperature in ISA [K]

% pre-allocate Ta to be the same size as Altitude (53x1), with all elements 216.65

Ta = 216.65*ones(size(Altitude));

% calculate the first 45 elements of Ta:

Ta(1:45) = To-0.0065*Altitude(1:45); % Temperature at altitude [K]

for i = 1:53

Po = 101325; % Pressure in ISA [Pa]

Pa_R(i) = Po*(1-2.2558E-5*Altitude(i))^5.25864; % Pressure in altitude (row format) [Pa]

Pa = Pa_R.';

R = 287.05; % [J/kgK]

Density_altitude_R(i) = Pa(i)/(R*Ta(i)); % Density at each altitude (row format) [kg/m^3]

Density_altitude = Density_altitude_R.';

%% Relative Density

Relative_density_R(i) = Density_altitude(i)/1.225; % (row format)

Relative_density = Relative_density_R.';

%% True Airspeed

IAS = 250/1.944; % [m/s]

CCF = 0.996; % Compressibility correction factor

TAS_R(i) = (IAS*CCF)*sqrt((1/Relative_density(i))); % (row format) [m/s]

TAS = TAS_R.';

end

%% Coefficient of Lift

CDO = 0.023; % Zero-lift drag coefficient

K = 0.044; % Lift-induced drag factor

MTOW = 79010; % Maximum take-off mass [kg]

W = MTOW*9.81; % Weight in N

A= Thrust_22K/W;

B = sqrt((A.^2)+(12*CDO*K));

Cl = (6*CDO)./(A+B);

%% Drag Produced

S = 125; % Wing area [m^2]

C = CDO + (K*(Cl.^2));

D = 0.5*Density_altitude.*(TAS.^2)*S.*C; % [N]

%% Power Available

P_available = TAS.*Thrust_22K;

%% Power Required

P_required = TAS.*D;

%% Excess Power

P_excess = P_available - P_required;

%% Rate of Climb (ROC)

ROC = (TAS.*(Thrust_22K - D))./W; % [m/s]

if any(ROC <=1.5)

disp('Service Ceiling.')

else

disp ('Absolute Ceiling.')

end

%% Climb Angle

X = (Thrust_22K - D)./W;

Climb_angle = asind( X ); % [Degree]

%% Time Taken to Climb

Time = Altitude./ROC; % [sec]

%% Mass Flow Rate

m_dot_f = TSFC_22K.*Thrust_22K; % [kg/s]

%% Fuel Consumed to Climb

Fuel = m_dot_f.*Time; % [kg]

[Torsten]

You should be able to adapt it with what I posted.

Make a loop over the number of altitudes, compute the time needed to go from altitude i to altitude i+1 as time(i+1) and deduce the fuel consumed from altitude i to altitude i+1 as Fuel(i+1) = time(i+1)*m_dot_fuel. At the end of the loop, compute the actual weight at altitude i+1 as the weight at altitude i minus the fuel consumed from altitude i to altitude i+1, i.e. Fuel(i+1).

[Ke Yeun Yong]

@Torsten

I applied your coding and i got the following error in command window;

Unable to perform assignment because the left and right sides have a different number of elements.

Error in Copy_of_Calculation_22K (line 86)

Time(i+1) = (Altitude(i+1)-Altitude(i))./ROC; % [sec]

Furthermore, the time and fuel increases with altitude so i think i dont have to do sum for them. As the calculation shows the time and fuel at each altitude and not the difference between two altitude. Please correct me if I am wrong.

Thank you.

[Torsten]

ROC = (TAS.*(Thrust_22K - D))./W(i); % [m/s]

Use the local value for all arrays that are part of this expression.

Maybe

ROC = TAS(i)*(Thrust_22K(i)-D)/W(i);

?

We don't know. You must use the local values (at Altitude(i+1)) everywhere where needed.