Nomenclature
- a
-
step size scaling factor
- A
-
ratio of cos αcto cos αs
- AFM
-
Aircraft Flight Manual
- APM
-
Aero Propulsive Modelling
- BADA
-
Base of Aircraft Data
- CD
-
drag coefficient
- CD0
-
parasite drag coefficient
- CDi
-
induced drag coefficient
- CDc
-
drag coefficient at climb
- CD glide
-
drag coefficient at glide
- CDglide*
-
glide drag coefficient without induced drag
- CDs
-
sustained turn drag coefficient
- CDs-c
-
difference between sustained turn and climb drag coefficients
- CL
-
lift coefficient
- CLc
-
lift coefficient at climb
- CLglide
-
lift coefficient at climb
- CLs
-
lift coefficient at sustained turn
- CFD
-
Computational Fluid Dynamics
- CSA
-
Cuckoo Search Algorithm
- CTAS
-
Center/TRACON Automation System
- D
-
drag force
- DOF
-
Degree of Freedom
- E
-
Cost Function
- Fi
-
Fitness Values
- FDR
-
Flight Data Records
- FMS
-
Flight Management System
- g
-
gravitational acceleration
- GA
-
Genetic Algorithm
- GAME
-
Global Aircraft Modeling Environment
- h
-
flight altitude
- H(u)
-
heaviside function
- IP
-
Intellectual Property
- L(λ, s)
-
Lévy flight function
- M
-
mach number
- MSE
-
Mean Squared Error
- MAE
-
Mean Absolute Error
- MAPE
-
Mean Absolute Percentage Error
- n
-
number of nests
- nz
-
load factor
- pa
-
the rate of discovery of alien eggs
- PEP
-
Performance Engineer Program
- PSO
-
Particle Swarm Optimisation
- q
-
dynamic air pressure
- R
-
linear correlation coefficient
- RoC
-
Rate of Climb
- s
-
step size
- S
-
wing reference area
- t/c
-
time to climb
- T
-
Thrust
- u
-
coefficients of the drag polar
- V
-
True Airspeed
- W
-
Weight
Greek Symbols
- α
-
angle of attack
- αc
-
angle of attack at climb
- αs
-
angle of attack at sustained turn
- Γ
-
constant Lévy flight function
- ε
-
a random number drawn from a uniform distribution
- λ
-
Lévy flight parameter
1.0 Introduction
The jet aircraft flight dynamics simulations generally include a separate aerodynamic model and a propulsion model. Computational fluid dynamics (CFD), wind tunnel testing and/or flight testing either without or with propulsors operating at a constant throttle setting are the basis for the creation of an aerodynamic model; whereas utilising analytical techniques, CFD, or ground testing, the characterisation of propulsive forces and moments are achieved. The separate aerodynamic and propulsion models are used in conjunction in order to estimate the forces and moments acting on the aircraft [Reference Simmons, Gresham and Woolsey1].
The ability of executing a correct performance assessment with the precise prediction of the thrust and drag forces, increasing the flight efficiency and airspace capacity as well as conducting trajectory predictions in decision support tools require the development of an accurate aero-propulsive model (APM) for the jet aircraft. Development of a proper APM is composed of an aerodynamic drag polar model and a propulsive thrust model. Evidently, there are very few studies on APMs in the existing literature; whereas to the authors’ knowledge there is no study related to creating an APM model based on developing a drag polar model solely in the lack of thrust model/data provided by the engine manufacturers, as presented in this study.
Eurocontrol’s Base of Aircraft Data (BADA) model [2–Reference Nuic, Poinsot, Iagaru, Gallo, Navarro and Querejeta5] and the Global Aircraft Modelling Environment (GAME) model are the most well-known APMs in the current literature. The accuracy of both models is very limited in achieving optimal trajectory estimations and in conducting satisfactorily accurate aircraft performance assessments [Reference Baklacioglu and Cavcar6]. Although computer programs such as Boeing aircraft performance program (INFLT)/REPORT, Boeing’s performance software, and PEP, Airbus’s Performance Engineer Program use the equations of motion of the aircraft and are capable of generating climb, cruise and descent trajectory estimations, as mentioned by Suchkov et al., because of not only the dimensions of their databases that are necessary for calculating aero-propulsive forces and engine fuel flow rate but also the computational speed for achieving the flight profiles, these aircraft manufacturers’ models are not convenient to be adapted for flight management system (FMS) applications [Reference Baklacioglu and Cavcar6, Reference Ghazi, Mennequin and Botez7].
Moreover, the reluctance of engine manufacturers’ in providing their data to the aircraft industry for aero-propulsive modelling research, trajectory prediction studies and the development of air traffic management tools would cause problems related to intellectual property (IP) rights in case of the usage of the manufacturers’ data for modelling studies. Therefore, usage of aircraft flight manual (AFM) yields the most convenient resource for creating an APM, which can also be enhanced by the inclusion of flight test data if available.
Gong and Chan [Reference Gong and Chan8] proposed an APM to be utilised in the prediction of aircraft trajectories for the Center/TRACON Automation System (CTAS). Adopting INFLT in its development, this model was applied to estimate the climbing phase of Boeing 737-400 and Learjet 60. Using AFM data, Cavcar and Cavcar [Reference Cavcar and Cavcar9] derived an APM to find the time-to-climb values for Boeing 737-400. Although both of these models include the compressibility and wing camber effects, they are not incorporating the compressible drag rise effect above the critical Mach number. With regards to their work, Gong and Chan admitted that additional research on their aerodynamic model equation might introduce improvements in terms of model accuracy. Also, Gong and Chan used an engine scaling factor in their propulsive model rather than deriving an optimisation model to estimate the coefficients of an empirical formula for the engine thrust, taking into account flight altitude and Mach number effects [Reference Gong and Chan8]. Cavcar and Cavcar presented an empirical model for their propulsive model at various constant flight altitudes [Reference Cavcar and Cavcar9]. Baklacioglu and Cavcar [Reference Cavcar and Cavcar9] developed a genetic algorithm- (GA) based APM derived from Boeing 737-400 AFM data and capable of making climb and descent trajectory estimations. Their propulsive model took into account both flight altitude and Mach number as input parameters simultaneously [Reference Baklacioglu and Cavcar6]. But both propulsive models in [Reference Baklacioglu and Cavcar6, Reference Cavcar and Cavcar9] were developed utilising the available JT9D-7A turbofan engine data and these models were further normalised by using (a) the sea level maximum static take-off thrust with the Boeing 737-400’s CFM56-3-B1 turbofan engine’s sea level maximum static take-off thrust and (b) the thrust value at 0.8 Mach number at 35,000 ft flight altitude with the CFM56-3-B1’s corresponding value because of the lack of CFM56-3-B1 engine data. In a recent study, Simmons et al. [Reference Simmons, Gresham and Woolsey1] presented two methodologies to identify an integrated propulsion–airframe aerodynamic model and a decoupled propulsion model for fixed-wing aircraft with propellers using flight data. In their model, in order to generate data with high-quality information content for model identification, orthogonal phase-optimised multisine inputs were applied to both the control surfaces and propulsion system. Propulsion explanatory variables derived from propeller aerodynamics theory utilised in combination with the traditional aircraft modelling variables provided the basis of their aero-propulsive modelling. Oruc and Baklacioglu [Reference Oruc and Baklacioglu10] presented an aircraft performance model consisting of aerodynamic model, thrust and fuel flow rate models for B737-800 aircraft and obtained specific excess power contours implementing the energy manoeuverability method by using this model. They used flight data records (FDR) data to create their model, but their thrust model was based on a PW4056 turbofan engine model, and they normalised this model using the same procedure as in [Reference Baklacioglu and Cavcar6, Reference Cavcar and Cavcar9] for the CFM56-7B engines used in B737-800 aircraft.
In this study, for the aim of developing a drag polar model in the absence of a thrust model, the utilisation of a non-conventional optimisation technique was preferred because of the highly non-linear nature of the problem. In this sense, various drag polar estimation models were derived implementing cuckoo search algorithm (CSA) in order to optimise the model coefficients in the empirical formulation. The originality of this study lies in the fact that in the existing literature, this study constitutes the first attempt of creating a drag model: (i) utilising 6-DOF model data; (ii) solely by the use of sustained turn, climb and glide manoeuver data, which eliminates the need for a thrust model; and (iii) via a metaheuristic optimisation algorithm; namely, CSA.
This newly developed CSA model can be useful in real world applications in aircraft performance modelling and assessments; decision support applications like conflict detection, direct routing, arrival metering, precise planning of air traffic flow; minimisation of delays, congestion, as well as environmental effects; in combination with fuel flow rate prediction models; and performance of aero-vehicles in different flight phases [Reference Dinc, Caliskan and Ekici11–Reference Oruc, Baklacioglu and Turan15].
2.0 Methodology
In order to calculate the forces acting on an aircraft under specific flight conditions, an aero-propulsive model has to be used. This model includes two parts; namely, the first part is utilised so as to obtain the aerodynamic drag force, while the second provides the computation of the propulsive thrust force. Using the equations of motion and the difference between these two forces, at a specific weight, speed schedule and under given atmospheric conditions, the aircraft performance at climbing flight can be obtained. Inversely, knowing the climb performance of an aircraft enables the computation of the difference between these forces. In this regard, modelling that specific aircraft’s climb performance yields a complementary aero-propulsive model for this aircraft.
Using the kinetic equations of motion, rate of climb (ROC) of the aircraft may be stated as follows:
\begin{align}RoC = \frac{{dh}}{{dt}} = \frac{{\left( {T - D} \right)V}}{{W\left( {1 + \frac{V}{g}\frac{{dV}}{{dh}}} \right)}}\end{align}
where RoC abbreviates the rate of climb, T is thrust, D indicates drag, W shows aircraft weight, V denotes true airspeed, h represents flight altitude and g is the gravitational acceleration.
Drag and thrust are the two unknown parameters in Equation (1), so these are the variables to be modeled to accomplish an aero-propulsive model utilising this specific aircraft’s performance data; however, in the absence of a thrust model for an aircraft, new estimation modelling to obtain an aircraft drag model had to be proposed in this study. With regards to this, the thrust parameter in Equation (1) is eliminated by the following assumption:
Assuming that the aircraft flights with available excess power without accelerating, it performs a sustained turn manoeuver at a constant altitude where the corresponding equation of motion can be stated as:
\begin{align}\sum {F_s} = T - D = 0\end{align}
Then, the drag force for the manoeuver is:
\begin{align}{D_{Sustained\;Turn}} = {T_{Sustained\;Turn}}\end{align}
Considering the fact that the aircraft uses all of its excess power for the sustained turn manoeuver, the available thrust for the climb will be equal to available thrust for the sustained turn:
\begin{align}{T_{Sustained\;Turn}} = {T_{climb}}\end{align}
Hence, the available thrust force in Equation (1) can be replaced by the required thrust force for the sustained turn, which allows the elimination for the need of the thrust parameter in deriving the aircraft performance model. In this regard, Equation (1) can be rewritten as follows:
\begin{align}RoC = \frac{{\left( {{D_{Sustained\;Turn}} - {D_{climb}}} \right)V}}{{W\left( {1 + \frac{V}{g}\frac{{dV}}{{dh}}} \right)}}\end{align}
Consequently, modelling the drag polar for the aircraft would be adequate in estimating the aircraft performance utilising the above assumption eliminating the need for a thrust model via a reverse engineering process using Equation (5). Therefore, two different modelling approaches have been presented in this study: (i) implementing the 6-DOF model data for sustained turn and climb flight to accomplish an induced drag model and then incorporating the glide data to achieve the total drag polar model; and (ii) utilising the 6-DOF model data in order to include the C L -α dependency into the model.
3.0 Cuckoo search algorithm and analysis
Being one of the most recent nature-inspired metaheuristic algorithms and created by Xin-She Yang and Suash Deb in 2009 [Reference Oruc and Baklacioglu16–Reference Yang and Deb19], cuckoo search algorithm (CSA) is inspired by some cuckoo species having a brood parasitism; furthermore, rather than utilising simple isotropic random walks, this algorithm is improved by the Lévy flights [Reference Walton, Hassan, Morgan and Brown20]. In the current literature, CSA is reported to be potentially far more efficient compared to particle swarm optimisation (PSO) and GA [Reference Yang17–Reference Yang and Deb19]. Having an aggressive reproduction strategy, some species such as the Ani and Guira cuckoos lay their eggs in the communal nests so as to increase the probability of hatching their own eggs and removing other bird’s eggs without hesitation.
Describing the standard cuckoo search, there are three idealised rules applied as stated below:
-
• Each cuckoo lays one egg at a time, and drops it in a randomly selected nest.
-
• The best nests with high-quality eggs will be transferred over to the next generations.
-
• The number of available host nests is fixed, and the egg laid by a cuckoo is discovered by the host bird with a probability p a ε (0,1). There are two possibilities in that case; i.e. the host bird may either throw the egg out of its nest or leave the nest and build a new nest elsewhere [Reference Yang17, Reference Gunji, Deepak, Saraswathi and Mogili21, Reference Joshi, Kulkarni, Kakandikar and Nandedkar22].
The last assumption may be accomplished by replacing a fraction p a of the n host nests with new nests, which correspond to new random solutions. Each egg in a nest is supposed to represent a solution and only one egg (thus representing one solution) can be laid at a time by each cuckoo from the application point of view. The purpose of this process is to obtain a new and potentially more accurate solution in order to replace a not-so-good solution in the nests.
Although a more complicated version of the algorithm where each nest possesses multiple eggs indicating a set of solutions may be implemented, a simplified approach where each nest has only a single egg has been utilised in this study. In this regard, since each nest corresponds to one egg that also represents one cuckoo, no distinction between an egg, a nest or a cuckoo is made.
Adopting a switching parameter p a , a local random walk and the global explorative random walk have been combined in the developed algorithm. Here, the switching parameter p a represents the rate of discovery of alien eggs, which controls the combination of local and global random walk. The formulation for the local random walk can be given as:
\begin{align}x_i^{t + 1} = x_i^t + as \otimes H\left( {{p_a} - \in } \right) \otimes (x_j^t - x_k^t)\end{align}
where s indicates the step size and a is the step size scaling factor. x
j
t
and x
k
t
are two randomly chosen solutions. x
i
t+1
denotes a new solution and x
i
t
shows the current solution. H(u) is the heaviside function. ε indicates a random number drawn from a uniform distribution, and
$ \otimes $
is entry-wise product.
\begin{align}x_i^{t + 1} = x_i^t + aL\left( {s,\lambda } \right)\end{align}
\begin{align}L\left( {s,\lambda } \right) = \frac{{\lambda {\rm{\Gamma }}\!\left( \lambda \right){\rm{sin}}\left( {\frac{{\pi \lambda }}{2}} \right)}}{\pi }\frac{1}{{{s^{1 + \lambda }}}}\;,\;(s\, \gg {s_0} \gt 0)\end{align}
where the function Γ represents a constant for the Lévy flight parameter λ. In Equation (7), the first term on the right side indicates the current solution and the second term shows the probability of transition [Reference Yang17, Reference Yang23–Reference Oruc, Baklacioglu, Turan and Aydin25]. The flowchart for the CSA [Reference Oruc and Baklacioglu16] is shown in Fig. 1.
Flowchart of the CSA model [Reference Oruc and Baklacioglu16].
4.0 Drag polar model
Total drag of the aircraft can be represented as follows:
\begin{align}D = \frac{1}{2}\rho {V^2}S{C_D}\end{align}
As stated in the previous study, the following format of empirical formulation can be used for modelling the drag coefficient, where u 1–9 indicate the coefficients of the drag polar model:
\begin{align}{C_D} = ({u_1} + {u_2}M + {u_3}{M^2}){C_L} + ({u_4} + {u_5}M + {u_6}{M^2})C_L^2 + \;{u_7}{M^2} + {u_8}M + {u_9}\end{align}
The development of the CSA drag polar model was achieved in two stages in the absence of the thrust data. Using the sustained turn and climb data obtained from the 6-DOF model of the case study jet aircraft, firstly a relevant induced drag polar model has been created. As a second step, the remaining portion of the drag polar has been modelled by implementing the gliding flight data yielded by the 6-DOF model.
4.1 CSA drag polar model with small angle assumption for angle-of-attack
Induced Drag Polar Model (Level Sustained Turn-Climb Flight)
An induced drag polar model was developed based on the 6-DOF model data of the case study aircraft. Since the required thrust values can be initially defined for zero acceleration climb and level sustained turn utilising the excess power characteristics, it is possible to obtain the induced drag coefficient using these findings.
The rate of climb expression can be written as follows regarding the principle mentioned above:
\begin{align}RoC = \frac{{\left( {{D_{Sustained\;Turn}} - {D_{climb}}} \right)V}}{{W\left( {1 + \frac{V}{g}\frac{{dV}}{{dh}}} \right)}}\end{align}
From Equation (11), the drag coefficient for sustained turn-climb can be given as:
\begin{align}{C_{Ds - c}} = \frac{{ROC.W}}{{VqS}}\end{align}
where C Ds-c , q and S indicates the difference between the drag coefficients for the sustained turn and climb, dynamic air pressure and reference wing area, respectively.
The drag polar terms with the coefficients
${u_{7 - 9}}\;$
would be equal in sustained turn and climb phase, so these terms cancel each other resulting in the following empirical formulation for the induced drag coefficient:
\begin{align}{C_{{D_{S - C}}}} = ({u_1} + {u_2}M + {u_3}{M^2})\left( {{C_{Ls}} - {C_{Lc}}} \right) + \left( {{u_4} + {u_5}M + {u_6}{M^2}} \right)\left( {{C_{Ls}}^2 - {C_{Lc}}^2} \right)\end{align}
Eventually, an optimisation procedure, in which the coefficients u 1–6 of the following cost function was accomplished to achieve the optimal solution:
\begin{align}E = {C_{{D_{S - C}}}} - ({u_1} + {u_2}M + {u_3}{M^2})\left( {{C_{Ls}} - {C_{Lc}}} \right) - \left( {{u_4} + {u_5}M + {u_6}{M^2}} \right)\left( {{C_{Ls}}^2 - {C_{Lc}}^2} \right)\end{align}
where u 1–6 correspond to the coefficients of the induced drag polar CSA model; M indicates Mach number; C Ds-c, C Ls , and C Lc shows the induced drag coefficient obtained from sustained turn-climb flight, lift coefficient for sustained turn and lift coefficient for climb flight, respectively.
The CSA has been utilised to develop the induced drag polar model and the parameters of the algorithm were chosen as the number of nests n = 25, discovery rate of alien eggs p a = 0.25, number of iterations 100,000. The CSA model was obtained in a piecewise approach for Mach number values both less than and greater than 0.75. For M < 0.75 and M > 0.75, the model coefficients of the induced drag polar empirical formula were obtained as shown in Table 1.
The coefficients of the induced drag polar estimation model
Achieving an error evaluation, error values corresponding to each induced drag coefficient value versus Mach number were found as shown in Fig. 2(a) (M < 0.75) and Fig. 2(b) (M > 0.75) below. Likewise, the comparison of the estimated induced drag values achieved from the CSA model and the six DOF model induced drag values are depicted in Fig. 3(a) and 3(b) for M < 0.75, and in Fig. 3(c) and 3(d) for M > 0.75.
Error values corresponding to the induced drag polar model for the sustained turn-climb flight.
Comparison of [C Ds – C Dc ) induced ] six DOF model and [C Ds – C Dc ) induced ]CSA model for the sustained turn-climb (i) M < 0.75 ((a), (b)); (ii) M > 0.75 ((c), (d)).
The variation of lift coefficient values versus induced drag values were shown in Fig. 4 for (i) M < 0.75 (a) climb, (b) sustained-turn flight; and (ii) M > 0.75 (c) climb, (d) sustained-turn flight.
Lift coefficient values achieved from six DOF model versus calculated and predicted induced drag coefficient values for (i) M < 0.75 (a) climb, (b) sustained-turn flight; and (ii) M > 0.75 (c) climb, (d) sustained-turn.
It should be noted that in order to create the CSA model, the C D values calculated using Equation (12) were taken into consideration; whereas, the trimmed C D values taken from the 6-DOF model were compared with the CSA drag polar values later in Fig. 13. Several error types were achieved with respect to the calculated and predicted induced drag coefficients including minimum error, maximum error, mean squared error (MSE), mean absolute error (MAE), mean absolute percentage error (MAPE) and linear correlation coefficient (R), as shown in Table 2.
Error evaluation for calculated and predicted induced drag coefficients
The linear correlation coefficient values R = 0.999 and 0.9971 indicate that the developed CSA model has been able to predict the induced drag coefficient values accurately. It should be noted that for a model in which M < 0.75 and M > 0.75 distinction was not made, R value was found as 0.9982.
Modelling the Remaining Portion (including CD0) of the Drag Polar (Gliding Flight)
Feeding glide 6-DOF model data into Equation (10) and utilising the model coefficients u 1–6 obtained from the sustained turn-climb flight, the induced drag values for the gliding flight were calculated from Equation (16), which provided the formula given in Equation (17). Subtracting these values from the total drag values achieved from the 6-DOF model data for the gliding flight, the values regarding the remaining portion of the total drag polar were calculated.
The coefficients of the remaining portion of the total drag polar model corresponding to parasite drag whose equation format was given in Equation (18) were obtained as shown in Table 3.
\begin{align}{C_{{D_{total,\;glide}}}} = \frac{{{C_{{L_{glide}}}}}}{{{{\left( {\frac{L}{D}} \right)}_{glide}}}}\end{align}
\begin{align}{C_{{D_{induced\;\;glide}}}} = \left( {{u_1} + {u_2}M + {u_3}{M^2}} \right){C_L} + \left( {{u_4} + {u_5}M + {u_6}{M^2}} \right)C_L^2\end{align}
\begin{align}C_{{D_{glide}}}^* = \frac{{{C_{{L_{glide}}}}}}{{{{\left( {\frac{L}{D}} \right)}_{glide}}}} - \left( {{u_1} + {u_2}M + {u_3}{M^2}} \right){C_L} + \left( {{u_4} + {u_5}M + {u_6}{M^2}} \right)C_L^2\end{align}
\begin{align}C_{{D_{glide}}}^* = {u_7}{M^2} + {u_8}M + {u_9}\end{align}
The model coefficient u9 in Equation (18) refers to the constant parasite drag coefficient; whereas, the other coefficients u7 and u8 may be thought of as corresponding to the remaining skin friction drag, pressure drag and wave drag effects. Comparing the calculated and forecasted values for this remaining portion of the drag polar, the calculated delta error is shown in Fig. 5 while Fig. 6 indicates the calculated versus predicted
$C_{D}{}^{\ast}$
values.
The coefficients of the remaining portion of the CSA drag polar model
Delta error values between the calculated and predicted values for the remaining portion of the drag polar.
Comparison of calculated versus estimated
$C_{D}{}^{\ast}$
drag coefficients.
Error assessment results corresponding to the calculated versus predicted
$C_{D}{}^{\ast}$
drag coefficients for the glide flight are given in Table 4.
The CSA model created to predict the C D * values using the glide 6-DOF model data was able to converge highly satisfactorily taking into account the linear correlation coefficient R = 0.9994 given in Table 4. In Fig. 7, a direct comparison between the drag coefficient values taken from the 6-DOF model were made with the drag polar values predicted by the CSA model having the small angle-of-attack assumption for the sustained turn, climb and glide phases. It should be noted that these C D values of 6-DOF model were trimmed values of C D taken from the aerodynamic database based on CFD and wind tunnel tests incorporated into the 6-DOF model.
4.2 CSA drag polar model incorporating the effect of angle-of-attack
In order to increase the generalisation capability of the CSA drag polar model developed, the objective functions for sustained turn-climb and glide flight were modified incorporating the effect of angle-of-attack into the model. The C L -α dependency relationship was achieved utilising an aero-database based on wind tunnel test and CFD analysis results.
Error evaluation for calculated versus estimated CD* drag coefficients
Comparison of trimmed and predicted drag coefficients for the sustained turn, climb and glide phases utilising the CSA model with small angle assumption for angle-of-attack (for M < 0.8).
With regards to this, so as to include the angle-of-attack effect in the new drag polar formulation, a simple linear optimisation model for C L (α) was obtained by using the aero-database of case study jet aircraft. The angle-of-attack values achieved from the aero-database data were given in Equation (19). This C L -α relationship have been further incorporated into the sustained turn, climb and glide drag polar models:
\begin{align}\alpha = \frac{{{C_L} - 0.0151}}{{0.0657}}\end{align}
The equations of motion are modified to include the effect of angle-of-attack into the empirical formulation for the drag polar as given below:
For the climb phase:
\begin{align}\frac{{ROC.W}}{V} = Tcos{\alpha _c} - {D_c}\end{align}
For the sustained turn phase:
\begin{align}Tcos{\alpha _s} = {D_s}\end{align}
From the above relations, it is possible to obtain the following expressions:
\begin{align}\frac{{ROC.W}}{V} = {D_s}\frac{{cos{\alpha _c}}}{{cos{\alpha _s}}} - {D_c}\end{align}
\begin{align}\frac{{ROC.W}}{{VqS}} = {C_D}_s\frac{{cos{\alpha _c}}}{{cos{\alpha _s}}} - {C_D}_c\end{align}
Equation (22b) can be rewritten as in the below if the term cos α c /cos α s is denoted by A:
\begin{align}\frac{{ROC.W}}{{VqS}} = A{C_D}_s - {C_D}_c\end{align}
The second degree drag polar empirical formulation may be written as follows:
\begin{align}{C_D} = {C_{{D_0}}} + {K_1}{C_L} + {K_2}C_L^2\end{align}
Indicating the parasite drag and induced drag coefficients as a function of Mach number, it is possible to write:
\begin{align}{C_{{D_0}}} = {u_1} + {u_2}M + {u_3}{M^2}\end{align}
\begin{align}{K_1} = {u_4} + {u_5}M + {u_6}{M^2}\end{align}
\begin{align}{K_2} = {u_7} + {u_8}M + {u_9}{M^2}\end{align}
Therefore, Equation (23) can be rewritten as:
\begin{align} \frac{{ROC.W}}{{VqS}} &= ({u_1} + {u_2}M + {u_3}{M^2})\left( {A - 1} \right) + \left( {{u_4} + {u_5}M + {u_6}{M^2}} \right)\left( {A{C_{{L_s}}} - {C_L}_c} \right)\nonumber\\&\quad + \left( {{u_7} + {u_8}M + {u_9}{M^2}} \right)\left( {AC_{{L_s}}^2 - C_{{L_c}}^2} \right)\end{align}
It should be noted that a small angle-of-attack assumption should be still valid so as to be able to calculate lift coefficients in the sustained turn and climb phase.
\begin{align}{C_{{L_c}}} = W/qS\end{align}
\begin{align}{C_{{L_s}}} = {C_{{L_c}}}.nz\end{align}
where nz refers to the load factor. Receiving the values of ROC, W, and V from the 6-DOF model, these C L values were used to calculate the corresponding difference of C Ds and C Dc values.
Similarly, an expression for the glide can be given as follows:
\begin{align}{C_{{D_{glide}}}} &= ({u_1} + {u_2}M + {u_3}{M^2}) + ({u_4} + {u_5}M + {u_6}{M^2}){C_{{L_{glide}}}}\nonumber \\ &\quad+ \left( {{u_7} + {u_8}M + {u_9}{M^2}} \right)C_{{L_{glide}}}^2\end{align}
Obtaining the (L/D) glide values from glide 6-DOF model data and substituting the C Lglide values, C Dglide values were also calculated.
Eventually, an optimisation procedure, in which the coefficients u1–9 of the following two cost functions were optimised simultaneously, was accomplished to achieve the optimal solution:
\begin{align}{E_1} &= \frac{{ROC.W}}{{VqS}} - [({u_1} + {u_2}M + {u_3}{M^2})\left( {A - 1} \right) + \left( {{u_4} + {u_5}M + {u_6}{M^2}} \right)\left( {A{C_{{L_s}}} - {C_L}_c} \right)\nonumber \\& \quad + \left( {{u_7} + {u_8}M + {u_9}{M^2}} \right)\left( {AC_{{L_s}}^2 - C_{{L_c}}^2} \right)]\end{align}
\begin{align}{E_2} = \frac{{{C_{{L_{glide}}}}}}{{{{\left( {\frac{L}{D}} \right)}_{glide}}}} - [({u_1} + {u_2}M + {u_3}{M^2}) + ({u_4} + {u_5}M + {u_6}{M^2}){C_{{L_{glide}}}} + \left( {{u_7} + {u_8}M + {u_9}{M^2}} \right)C_{{L_{glide}}}^2]\end{align}
where E 1 and E 2 indicates the error values between the calculated and estimated drag polar values for sustained turn-climb and glide flight, respectively.
Finetuning the coefficients of the cost functions given above, the errors E 1 and E 2 were minimised in a Matlab script.
Different formats of drag polar formulations were tried in order to obtain a good convergence between the predicted drag coefficient values by the CSA model and the data provided by the 6-DOF model. Different variations including (u i +u i+1 M+u i+2 M 2), (u i +u i+1 M 2), (u i +u i+1 M+u i+2 M 2+u i+3 M 4), (u i +u i+1 M+u i+2 M 4), (u i +u i+1 M 4) terms as well as 2nd, 3rd, and 4th degree lift coefficient formats of the induced drag were experimented. Consequently, it was deduced that the best accuracy was achieved using the following drag polar expression when the effect of angle-of-attack was incorporated:
For the sustained turn-climb flight:
\begin{align}\frac{{ROC.W}}{{VqS}} &= ({u_1} + {u_2}M + {u_3}{M^2})\left( {A - 1} \right) + ({u_4} + {u_5}M + {u_6}{M^2})\left( {A{C_{{L_s}}} - {C_L}_c} \right)\nonumber \\ &\quad + ({u_7} + {u_8}M + {u_9}{M^2})\left( {AC_{{L_s}}^2 - C_{{L_c}}^2} \right)\end{align}
For the gliding flight:
\begin{align} {C_{{D_{glide}}}} &= ({u_1} + {u_2}M + {u_3}{M^2}) + ({u_4} + {u_5}M + {u_6}{M^2}){C_{{L_{glide}}}} \nonumber\\& \quad + \left( {{u_7} + {u_8}M + {u_9}{M^2}} \right)C_{{L_{glide}}}^2\end{align}
The results shown in Table 5 were obtained when a simultaneous optimisation of the drag polar coefficients was conducted.
The coefficients of the drag polar model including the effect of angle-of-attack
For the sustained turn-climb, comparing the predicted and calculated drag coefficients, the delta error was obtained as shown in Fig. 8. It should be noted that this time, the difference between the C D terms for the sustained turn and climb gives the equation for the total drag but not the induced drag portion as in the case of small angle assumption because of the parasite drag term having a multiplier of (A−1); i.e. a difference term dependent on the cosine ratios of climb and sustained turn angle-of-attack. Eventually, taking into account the effect of angle-of-attack in the drag polar equation provided the total drag polar formulation; furthermore, simultaneous optimisation of the drag polar coefficients achieved the optimal drag polar formula.
Error values for the drag polar model including angle-of-attack effect for the sustained turn-climb.
For the sustained turn-climb flight, the estimated CSA drag model values and calculated drag coefficient values using the six DOF model were compared as shown in Fig. 9(a) and 9(b). The lift coefficient values taken from the six DOF model versus calculated and predicted drag polar values for: (a) climb and (b) sustained turn flight were depicted in Fig. 10. A simultaneous optimisation of the drag polar coefficients was accomplished for the sustained turn-climb and glide data from the six DOF model and the error evaluation comparing the calculated and predicted drag coefficient values for the glide is shown in Fig. 11 while Fig. 12 indicates the calculated 6-DOF model drag coefficients versus predicted values by the CSA model processing the glide data.
Comparison of (C Ds – C Dc ) six DOF model and (C Ds – C Dc )CSA model including α for the sustained turn-climb.
Lift coefficient values obtained from the six DOF model versus calculated and predicted drag polar values from the CSA model including the effect of angle-of-attack for: (a) climb and (b) sustained-turn.
Error values between the calculated and estimated drag coefficient values for the gliding flight utilising the CSA model including angle-of-attack.
Comparison of calculated and predicted drag coefficients using gliding flight 6-DOF model data utilising CSA model including angle-of-attack.
Different types of errors were computed corresponding to the calculated versus estimated C D drag coefficients by the CSA model incorporating the angle-of-attack for the sustained turn-climb flight are shown in Table 6.
As shown in Table 6, the linear correlation coefficient R was found as 0.999 for the CSA model, which means that C D values were estimated highly accurately via the CSA including the α effect for the sustained turn-climb flight.
Similarly, for the glide flight, the error values corresponding to the CSA model including α effect are shown in Table 7.
The value of the linear correlation coefficient R = 0.9824 as seen in Table 7 shows that the CSA model incorporating α effect managed to predict C D values for the gliding flight pretty accurately. Taking into account sustained turn-climb and glide 6-DOF model data together, the overall R of the simultaneous CSA optimisation model is equal to 0.9998, which is very close to 1 showing the derived model’s highly accurate converging ability.
Finally, a direct comparison between the drag coefficient values taken from the 6-DOF model were made with the drag polar values predicted by the CSA model (including the effect of angle-of-attack) for the sustained turn, climb and glide phases. It should be again noted that these CD values of 6-DOF model were trimmed values of C D taken from the aerodynamic database based on CFD and wind tunnel tests incorporated into the 6-DOF model. This comparison is given in Fig. 13. Comparing the results in Fig. 7 versus Fig. 13, it was seen that the convergence of the CSA model was improved by incorporating the effect of angle-of-attack in the drag polar empirical formula. The created CSA model covers the whole Mach number interval and provides a complete drag polar model, but the developed methodology based on CSA is able to predict the trimmed C D values satisfactorily for subsonic flow regime. For the transonic region, an improvement on the drag polar empirical model incorporating the wave drag effects is planned to be implemented as a future study.
Error evaluation for calculated versus estimated C D drag coefficients (CSA model including α effect) for the sustained turn-climb
Error evaluation for calculated versus estimated CD drag coefficients (CSA model including α effect) for the glide
Comparison of trimmed and predicted drag coefficients for the sustained turn, climb and glide phases utilising the CSA model including angle-of-attack (for M < 0.8).
5.0 Conclusions
A reverse engineering methodology to achieve a drag polar estimation for a case-study jet aircraft was implemented throughout this study without having access to the thrust data of the aircraft’s engine. Two approaches: namely, (i) implementation of 6-DOF model data for sustained-turn and climb-flight to achieve an induced drag model, and then incorporating the glide data to accomplish the total drag polar model; and (ii) utilisation of 6-DOF model data together with the incorporation of the effect of CL-α dependency into the model were utilised to achieve an optimum and generic modelling approach using the CSA optimisation technique to optimise the model coefficients of the proposed empirical drag polar formulation. The following results were achieved:
-
• The originality of the study lies in the facts that to the authors’ knowledge, (i) this methodology is the first attempt utilising a stochastic metaheuristic optimisation algorithm such as CSA for drag polar modelling in the existing literature, and (ii) for the first time in the literature, in the case of no access to the thrust data (or absence of an individual thrust model), a drag polar modelling methodology enabling the elimination for the need of the thrust parameter in deriving an aero-propulsive model implementing equations of motion was proposed here.
-
• Including the variation of angle- of-attack in the drag polar formulation in the developed CSA model was seen to increase the converging ability of the model. It should be added that a small angle assumption should be still made in order to achieve the C L -α relationship to be inserted in the model.
-
• Within different trials, it was seen that only the 1st and 2nd degree lift coefficient terms in the drag polar empirical formulation should be preserved while incorporating the 3rd and 4th degree lift coefficient terms tends to decrease the convergence of the model. On the other hand, including 1st and 2nd degree Mach terms (drag polar format of u 1+u 2M+u 3 M 2) provided the best converging results for the CSA model.
-
• Compatibility of the estimated drag polar values obtained from the model with the trimmed C D values of 6-DOF model were found to be satisfactory within the high subsonic and transonic interval of C L values. This highlights the limitations of this drag polar estimation model since it has been derived utilising the climb, sustained turn and gliding flight data within the range of specified C L values. Including data from other stages of flight would probably improve the model coefficient results to provide better accuracy and be more applicable to other flight phases. Nonetheless, the methodology proposed in this study is the first of its kind in the literature in the estimation of drag polar without using engine thrust input data, about which manufacturers are reluctant to provide information; furthermore within the limitations of climb-sustained turn and glide data, the developed methodology was proven to work satisfactorily for the climb phase and seen to be accurate enough and practically usable within a certain C L range of the case study jet’s climbing flight.
-
• The CSA drag polar model developed in this study initiates a basis for creating similar accurate drag polar prediction models utilising AFM data as well as flight test data.
-
• Integrating the inverse of Equation (1) with respect to flight altitude provides the time-to-climb given as:
(34)
\begin{align}t/c = \smallint \frac{{W\left( {1 + \frac{{VdV}}{{gdh}}} \right)}}{{\left( {T - D} \right)V}}dh\end{align}
This equation is the foundation for the derivation of an aero-propulsive model, since it relates the climb performance of an aircraft to aero-propulsive forces. With regards to this, utilising the drag polar model derived in this study, a climb performance estimation model is planned to be created as another study.
Competing interests
None.
