Gliding Footprint Calculation Method for Aircraft with Thrust Failure Based on Six-Degree-of-Freedom Flight Envelope and Back-Propagation Artiﬁcial Neural Network

: In hostile environments, engine damage is of particular concern since the engine is the only component to generate thrust that affects survivability. For an aircraft suffering thrust failure, the forced landing sites should be identiﬁed within the gliding footprint, which is the reachable region on the ground. This paper proposes two calculation methods to obtain the gliding footprint by ﬁnding a series of boundary points with maximum gliding distance around the aircraft. Method 1 models the thrust-failed aircraft with six-degree-of-freedom (6-DOF) ﬂight dynamics and adopts a novel 6-DOF unpowered-ﬂight envelope to characterize its maneuvering capabilities. Given the initial altitude when thrust failure occurs, Method 1 determines all feasible gliding distances around the aircraft based on the constructed 6-DOF ﬂight envelope and selects the landing points of maximum gliding distances along different radial directions as the boundary points. Method 2 employs the Back-Propagation Artiﬁcial Neural Network (BP-ANN) to predict the boundary points. Using the well-trained BP-ANN, this method can estimate the maximum gliding distances with only the initial altitude and radial directions. Simulations are conducted to analyze these two methods. Compared with conventional methods using point-mass ﬂight dynamics, Method 1 considers more ﬂight constraints, and the gliding footprint area is reduced by 20.79%. These results are relatively conservative and can improve the safety threshold of forced landing sites. Method 2 can estimate the gliding footprints (encircled by the boundary points under the entire operational altitude and full radial direction) in real time, which reserves more response and action time for aircraft forced landing.


Introduction
Aircraft survivability is defined as the ability of an aircraft to avoid and/or withstand hostile environments [1].Attacks from air defense and terrorist weapons can bring in-flight damage to aircraft [2][3][4]].An attrition kill (the loss of aircraft) occurs when the aircraft loses one or more essential functions (i.e., structural integrity, lift, thrust, and control) for flight [1].To maintain such damaged aircraft in controlled flight, various technologies (including damage identification, control reconfiguration, envelope protection [5][6][7][8], etc.) have been developed.Subsequently, a successful landing [9,10] can prevent the attrition kill of a damaged aircraft, which actually improves its survivability.
Thrust failure is a common cause of kills/losses in military/civilian aircraft [11,12].In such cases, the determination of the forced landing sites is very significant since pilots have only one chance for landing [13].A well-known case is a passenger aircraft landing in the Hudson River with a twin-engine failure, where all 155 people onboard survived.In the accident report, the captain's decision to choose the Hudson River as the landing site rather than an airport was considered to provide the highest survival probability [14].Typically, the landing site for thrust-failed aircraft is determined in three steps [15].First, the gliding footprint, which is the reachable region on the ground, is estimated.Then, feasible landing sites are searched within it.These sites are finally prioritized by calculating a utility function to find the top-ranked one in terms of suitability as the target landing site.Knowledge of the gliding footprint can limit the search space for landing sites [16] and improve the efficiency of landing site determination.Therefore, the gliding footprint needs to be well developed.
To obtain the gliding footprint, a critical problem is finding the farthest gliding trajectories around the aircraft, because the footprint is encircled by corresponding landing points.Current methods can be summarized into two solutions.The first solution is based on the optimization theory.Maya et al. [17] employed a pseudo-spectral method to compute the farthest trajectories and gliding footprint, with a cost function built using the gliding distance.The control variables were lift coefficient and roll angle.Segal et al. [18] further calculated the footprint considering wind in the airspace, with velocity and bank angle set as control variables.Similar methods are also widely applied in the gliding footprint calculation of re-entry vehicles [19,20].To follow the optimal trajectories in [16,17], continuous maneuvers are required in a short period of time.This is a great challenge for unpowered aircraft with degraded performance, thus compromising landing safety.
The second solution employs a segmented landing trajectory that is relatively easy for aircraft to track.This trajectory is extended from a Dubins curve [21] and contains a steady straight segment with/without a steady turn segment.When calculating the gliding footprint under a given initial altitude, the focus is on finding the steady states that produce the farthest gliding trajectories over the entire range of maneuvering capabilities.Siegel et al. [22] adopted a segmented trajectory that included only a straight segment and assumed the gliding footprint to be a circle.The circle center was the vertical projection of the initial position, and the radius was the forward gliding distance with the straight segment at the best-glide angle (i.e., the maximum and also shallowest flight-path angle, negative value).This angle was chosen because the corresponding gliding distance is maximum for a given altitude loss [23].Some other researchers used a segmented trajectory consisting of a steady turn segment and a steady straight segment, and determined the gliding footprint based on the boundary points with maximum gliding distances around the aircraft.These studies all chose the best-glide state in the straight segment and differed mainly in the turn segment.Atkins et al. [15] predefined the turning states and obtained a circular gliding footprint based on three boundary points.Coombes et al. [24] also specified the turning states and obtained a gliding footprint encircled by a series of boundary points.Poissant et al. [25] established several groups of turning states and investigated the turning parameter effects on the gliding footprint.However, the determination of the best-turn states that corresponded to the farthest segmented trajectories was still unresolved.Seminal work was carried out by Di Donato et al. [26], who utilized the unpowered-flight envelope to characterize aircraft maneuvering capabilities.The best-turn states were computed by traversing the steady turning states in the envelope.
The flight envelope in [26] is a discrete library of attainable steady straight and steady turning states of the thrust-failed aircraft.An attainable steady state means that the aircraft can be trimmed at this state.Each steady state is characterized by a group of parameters (i.e., altitude, bank angle, and flight-path angle) which are set as known variables during the trim process.By setting proper ranges and step sizes for the characterization parameters, the flight envelope is computed and compiled into a three-dimensional matrix.Similar studies have also been conducted on powered aircraft [27][28][29][30].Altitude, velocity, turn rate, and flight-path angle are the characterization parameters, and the flight envelope is compiled into a four-dimensional matrix.The conventional velocity-altitude flight envelope is a subset of this matrix where the turn rate and flight-path angle are zero [31].The method in [26] requires extensive computations since every footprint boundary point is obtained after traversing all turning states in the flight envelope.An Artificial Neural Network (ANN) is an information-processing system inspired by biological neural networks [32].It can acquire knowledge by training and solve problems of the same class [33].The Back-Propagation Artificial Neural Network (BP-ANN) is one of the most widely applied neural networks.It has powerful computational advantages and can effectively map various nonlinear relationships [34].Current studies have applied ANN to reduce the computational time in computationally intensive situations.Norouzi et al. [27] utilized ANN on an aircraft with aileron failure and predicted the flight envelope in real time.Zhang et al. [35] further used BP-ANN to estimate the flight envelope of an aircraft with wing damage.These studies provide a good reference for gliding footprint calculation.
Although existing methods can obtain the gliding footprint well, there are still several challenges.First, all methods model the thrust-failed aircraft with point-mass flight dynamics.The aircraft is trimmed by balancing the forces only, without considering the moment equilibrium.This will lead to deviations in the flight parameters, especially the flight-path angle, which can accumulate into a considerable gliding distance error over the whole forced landing.Second, the unpowered-flight envelope contains only low-velocity steady states and does not match the aircraft's full maneuvering capabilities.Given a group of characterization parameters for steady states, two different velocities can be computed using force equilibrium equations.When compiling the flight envelope, only the state with a smaller velocity is selected.Third, the determination of best-turn states is computationally intensive and time-consuming, making the footprint results unsuitable for direct application to in-flight damage conditions that require the footprint in real time.
This paper proposes two methods to calculate the glide footprint of a fixed-wing aircraft with thrust failure.Method 1 extends the calculation process in [26], with the aircraft modeled using 6-DOF flight dynamics instead of point-mass dynamics.To characterize the maneuvering capabilities of the 6-DOF thrust-failed aircraft, a novel unpowered-flight envelope model is proposed.Method 2 utilizes BP-ANN to estimate the gliding footprint, with the network trained using the gliding footprint results obtained by Method 1.
In summary, this paper makes the following contributions: (1) We develop a gliding footprint calculation method based on the 6-DOF flight envelope.
Compared with conventional methods that employ point-mass dynamics, this 6-DOF method considers more constraints, including the moment equilibrium, control surface saturation, and aerodynamic contributions of the control surfaces.This paper is structured as follows.In Section 2, we describe the gliding footprint calculation problem and the research framework of this study.In Section 3, we detail the 6-DOF unpowered-flight envelope model, the flight envelope-based footprint calculation method, and the footprint results in comparison with existing methods.In Section 4, we present the BP-ANN model-based footprint estimation method and the footprint prediction results.Finally, we summarize our conclusions and future work in Section 5.

Problem Definition
In this section, we first describe the gliding footprint calculation problem and clarify the research scope.Then, we present the research framework.

Problem Description
The gliding footprint of a fixed-wing aircraft with thrust failure is shown in Figure 1 [15,26].Optional landing sites are searched within it.To calculate the gliding footprint, this research focuses on finding the footprint outer boundary.For each boundary point, the radial angle ξ is used to characterize its radial azimuth relative to the initial heading direction, and its gliding distance d ξmax is the maximum value along angle ξ.The gliding footprint is finally encircled by a series of boundary points d ξmax , ξ around the aircraft.
The gliding footprint of a fixed-wing aircraft with thrust failure is shown in [15,26].Optional landing sites are searched within it.To calculate the gliding fo this research focuses on finding the footprint outer boundary.For each boundar the radial angle  is used to characterize its radial azimuth relative to the initial direction, and its gliding distance d max  is the maximum value along angle .The footprint is finally encircled by a series of boundary points ( ) To further clarify this issue, the following notes and assumptions are provid

•
The thrust-failed aircraft suffers no other failures.

•
The aircraft will not gain energy from the atmosphere during the forced land

•
The gliding trajectory includes a steady straight segment with/without a stea segment (as in Figure 1).It contains a straight segment when landing along th heading direction, while it includes a turn and a straight segment for other di

•
The transition between the turn and straight segments is accomplished ins ously.

•
There is no wind or obstacles in the airspace.

Research Framework
The research framework of this study is shown in Figure 2 and includes two Phase 1 calculates the gliding footprint based on the 6-DOF flight envelope, wh library of attainable steady straight and steady turning states.These states are co by a 6-DOF trim using the aircraft data at the initial altitude.The footprint outer bo is gained by extracting the maximum values of all feasible gliding distances aro aircraft, with the straight segment in the best-glide state and the turn segment tra each steady turning state in the flight envelope.To further clarify this issue, the following notes and assumptions are provided:

•
The thrust-failed aircraft suffers no other failures.

•
The aircraft will not gain energy from the atmosphere during the forced landing.

•
The gliding trajectory includes a steady straight segment with/without a steady turn segment (as in Figure 1).It contains a straight segment when landing along the initial heading direction, while it includes a turn and a straight segment for other directions.

•
The transition between the turn and straight segments is accomplished instantaneously.

•
There is no wind or obstacles in the airspace.

Research Framework
The research framework of this study is shown in Figure 2 and includes two phases.Phase 1 calculates the gliding footprint based on the 6-DOF flight envelope, which is a library of attainable steady straight and steady turning states.These states are computed by a 6-DOF trim using the aircraft data at the initial altitude.The footprint outer boundary is gained by extracting the maximum values of all feasible gliding distances around the aircraft, with the straight segment in the best-glide state and the turn segment traversing each steady turning state in the flight envelope.Details of Phase 1 are shown in Section 3.
The flight envelope-based method in Phase 1 is computationally intensive and cannot satisfy the real-time requirement of in-flight damage conditions.Therefore, we introduce the BP-ANN model in Phase 2, aiming to estimate the gliding footprint in real time.To train the BP-ANN, we establish a database of gliding footprint results at various initial altitudes computed using the flight envelope-based method.For a given initial altitude, the trained BP-ANN is utilized to predict the gliding footprint.The whole process of Phase 2 is described in Section 4. The flight envelope-based method in Phase 1 is computationally intensive and cannot satisfy the real-time requirement of in-flight damage conditions.Therefore, we introduce the BP-ANN model in Phase 2, aiming to estimate the gliding footprint in real time.To train the BP-ANN, we establish a database of gliding footprint results at various initial altitudes computed using the flight envelope-based method.For a given initial altitude, the trained BP-ANN is utilized to predict the gliding footprint.The whole process of Phase 2 is described in Section 4.

Footprint Calculation Based on 6-DOF Flight Envelope
In this section, we detail the gliding footprint calculation based on the 6-DOF unpowered-flight envelope.First, we describe how to achieve the maximum gliding distance along segmented trajectories.Then, we show the 6-DOF flight envelope model and the 6-DOF trim process that computes attainable steady states.After this, we present a flight envelope-based calculation method that obtains the footprint by finding its boundary points.Finally, we analyze the simulation results of the flight envelope and gliding footprint.

Maximum Gliding Distance Achievement
This study employs a segmented landing trajectory that contains a straight segment with/without a turn segment.Given an initial altitude, the achievement of maximum gliding distances is presented below.

Gliding Distance Maximization
Figure 3 shows a top view of feasible segmented landing trajectories given the initial altitude and radial angle ξ.To achieve the maximum gliding distance, the straight segments are all at the best-glide angle bg , str γ .For different steady turning states (TS1, TS2…TSn) in the turn segment (AB1, AB2…ABn), the landing points (C1, C2…Cn) can have different . By extracting the maximum gliding distance, we will find the boundary point ( )

Footprint Calculation Based on 6-DOF Flight Envelope
In this section, we detail the gliding footprint calculation based on the 6-DOF unpoweredflight envelope.First, we describe how to achieve the maximum gliding distance along segmented trajectories.Then, we show the 6-DOF flight envelope model and the 6-DOF trim process that computes attainable steady states.After this, we present a flight envelopebased calculation method that obtains the footprint by finding its boundary points.Finally, we analyze the simulation results of the flight envelope and gliding footprint.

Maximum Gliding Distance Achievement
This study employs a segmented landing trajectory that contains a straight segment with/without a turn segment.Given an initial altitude, the achievement of maximum gliding distances is presented below.

Gliding Distance Maximization
Figure 3 shows a top view of feasible segmented landing trajectories given the initial altitude and radial angle ξ.To achieve the maximum gliding distance, the straight segments are all at the best-glide angle γ bg,str .For different steady turning states (T S1 , T S2 . ..T Sn ) in the turn segment (AB 1 , AB 2 . ..AB n ), the landing points (C 1 , C 2 . ..C n ) can have different gliding distances (d ξ1 , d ξ2 . ..d ξn ).By extracting the maximum gliding distance, we will find the boundary point d ξmax , ξ along radial angle ξ.If we set a step size for ξ, then we can find a series of discrete boundary points around the aircraft.To achieve the maximum gliding distance along a radial angel ξ, we need b g , str γ and all feasible segmented trajectories along ξ.The value of b g , str γ can be extracted from the flight-path angles of all attainable steady straight states, and the feasible segmented trajectories need all attainable steady turning states.They both require the unpowered-flight envelope (as detailed in Section 3.2), which is a library of attainable steady turning and To achieve the maximum gliding distance along a radial angel ξ, we need γ bg,str and all feasible segmented trajectories along ξ.The value of γ bg,str can be extracted from the flight-path angles of all attainable steady straight states, and the feasible segmented trajectories need all attainable steady turning states.They both require the unpoweredflight envelope (as detailed in Section 3.2), which is a library of attainable steady turning and steady straight states of the unpowered aircraft.

Gliding Distance Formula Derivation
Given the initial altitude, the best-glide angle γ bg,str , and a steady state in the turn segment, the gliding distance d ξ of a feasible segmented trajectory landing along the radial angle ξ is derived as follows.
Figure 4 shows a top view of the segmented trajectory landing on the right-half gliding footprint along radial angle ξ, with necessary angles and auxiliary lines added.The coordinate origin refers to the initial position when thrust failure occurs.The forward direction is consistent with the initial heading direction, while the lateral direction is perpendicular to it.The range of ξ is [−180 • , 180 • ], with the value positive if the landing point is on the right side of the forward direction and negative if it is on the left [26].The angle ∆ψ is the angular change in the heading direction during the turn segment and is also the central angle.The lengths of the arc and straight line in the horizontal plane are denoted by d trn and d str .The gliding distance d ξ is the distance between the landing point and the coordinate origin.To achieve the maximum gliding distance along a radial angel ξ, we need b g , str γ and all feasible segmented trajectories along ξ.The value of b g , str γ can be extracted from the flight-path angles of all attainable steady straight states, and the feasible segmented trajectories need all attainable steady turning states.They both require the unpowered-flight envelope (as detailed in Section 3.2), which is a library of attainable steady turning and steady straight states of the unpowered aircraft.

Gliding Distance Formula Derivation
Given the initial altitude, the best-glide angle b g , str γ , and a steady state in the turn segment, the gliding distance d ξ of a feasible segmented trajectory landing along the radial angle ξ is derived as follows.
Figure 4 shows a top view of the segmented trajectory landing on the right-half gliding footprint along radial angle ξ, with necessary angles and auxiliary lines added.The coordinate origin refers to the initial position when thrust failure occurs.The forward direction is consistent with the initial heading direction, while the lateral direction is perpendicular to it.The range of ξ is [−180°, 180°], with the value positive if the landing point is on the right side of the forward direction and negative if it is on the left [26].The angle Δψ is the angular change in the heading direction during the turn segment and is also the central angle.The lengths of the arc and straight line in the horizontal plane are denoted by d trn and d str .The gliding distance d ξ is the distance between the landing point and the coordinate origin.In the turn and straight segments, the altitude loss is calculated as C 1 ( ) In the turn and straight segments, the altitude loss is calculated as where γ trn is the flight-path angle in the turn segment, and R is the horizontal turn radius.
Based on the geometric relationships in Figure 4a, given an initial altitude H 0 , we obtain [26] where ∆ψ satisfies ξ < ∆ψ < 2ξ.In this study, Equations ( 3) and ( 4) are verified to also be applicable for 90 • ≤ ξ ≤ 180 • as in Figure 4b.In addition, if the turn segment is conducted at a left turning (−180 • ≤ ξ < 0 • ), the value of ξ in Equations ( 3) and (4) should be replaced by its absolute value.
If the parameters ξ, H 0 , γ bg,str , γ trn , and R are known, ∆ψ and d ξ can be calculated by Equations ( 3) and (4), respectively.Among them, the values of ξ and H 0 are known, and γ bg,str is extracted from the trim solutions of all attainable steady straight states.For a given steady turning state, γ trn is within the trim solution and R is computed as follows.
In the flight-path frame, the unpowered aircraft satisfies the following dynamic equation of motion: mV .
where m is the aircraft mass; V is the velocity; χ is the flight-path azimuth angle; the over-dot indicates the time derivative; L and Y are the lift and side forces; and µ is the bank angle.
This study considers the difference between the bank angle µ and roll angle φ.Their values are not equal, and µ is derived from other aerodynamic and flight-path angles [36]: where α and β are the angles of attack and sideslip, and θ is the pitch angle.
Then, the values of µ and .
χ can be calculated with Equations ( 6) and ( 5), respectively, where α, β, φ, θ, L, and Y are all obtained with the trim solution of the given turning state.
Therefore, the radius R is calculated as

Unpowered-Flight Envelope Modeling
Using point-mass dynamics, [26] studied the unpowered-flight envelope.But the envelope contained only the low-velocity steady states.As an extension, we propose a novel model of a 6-DOF unpowered-flight envelope consisting of all attainable steady states.

Six-Degree-of-Freedom Flight Envelope Model
The 6-DOF nonlinear dynamic equations of motion describing the translational and rotational accelerations of an aircraft with thrust failure are given as where f is a vector of nonlinear functions; x and u are the state and control vector; p, q, and r are the roll, pitch, and yaw rates; δ e , δ a , and δ r are the angular deflections in the elevator, aileron, and rudder.A steady state is the condition in which all forces and moments are constant or zero.This requires all linear and angular velocity rates and aerodynamic angles rates of the thrust-failed aircraft to be zero.Then, the dynamic state vector x d satisfies .
with the following additional constraints according to the flight condition: Steady straight state : . φ, where h is the flight altitude; ψ is the yaw angle; and the asterisk superscript denotes the desired trimmed conditions.Among the parameters ), we find the trim condition through the 6-DOF trim.A general method for deriving the trim condition is to solve the nonlinear constrained optimization problem that minimizes the cost function.A typical choice of the cost function is: where Q is a positive definite weighting matrix describing the contribution of each variable in .
x d to the cost function J.For a trimmed condition, J satisfies To trim the thrust-failed aircraft given the set (h * , V * , .ψ * ), ten parameters in the trim vector (x * , u * ) should be calculated.They are α, β, p, q, r, φ, θ, δ e , δ a , and δ r .As with existing studies on the gliding footprint, only the coordinated flight states are considered to land in this paper.So, the value of β is zero.Meanwhile, four of the aforementioned parameters (i.e., p, q, r, and θ) can be derived and computed for the others: where γ is introduced in the trim process and Then, six parameters, including α, γ, φ, δ e , δ a , and δ r , are left to be calculated.They are set as the trim control parameters, which are directly solved by the numerical optimization method.Meanwhile, Equation (11) corresponds to six dynamic equations of motion.The number of equations equals the number of unknown trim control parameters.So, the trim problem has a unique solution if the solution exists [37].Moreover, the inequality constraints of the trim control parameters are (21) This study employs the Sequential Quadratic Programming (SQP) technique [27][28][29] to solve the nonlinear optimization problem.A condition is considered to be a feasible trim condition if the cost function J < 10 −7 .After obtaining the trim solution, the value of .h * in Equations ( 12) and ( 13) can be computed as .
During the trim process, J is calculated using the following equations of motion [38]: . .
where D is the drag force; L, M, and N are the roll, pitch, and yaw moments; g is the gravitational acceleration; I i is the moment of inertia, and I ij is the product of inertia, with i, j = x, y, z.When calculating the aerodynamic forces and moments, this study considers the effects of both the angle of attack α and the control surface deflections δ e , δ a , and δ r .

Gliding Footprint Calculation Method
Conventional methods simplify the aircraft as a point mass.In this study, we develop a gliding footprint calculation method with the aircraft modeled using 6-DOF flight dynamics that considers more flight constraints.
The gliding footprint is encircled by a series of boundary points with maximum gliding distances.In this study, we define the flight parameters in steady states that produce the segmented trajectories to footprint boundary points as the best-footprint parameters, including γ bf , φ bf , V bf , ∆ψ b f , etc.Meanwhile, when the initial altitude H 0 is insufficient, there may exist certain unreachable areas within the footprint outer boundary.The simple connectivity of the gliding footprint is judged by the following equation [26,39]: where R min is the minimum turn radius, and l sc is the minimum horizontal projection length of the segmented trajectories arriving at the footprint outer boundary.
In [26], the footprint boundary points are calculated based on the unpowered-flight envelope characterized by (h * , γ * , µ * ).This study follows this idea, and the extension lies in two aspects: (1) the aircraft is modeled with 6-DOF dynamics and the steady states are computed by the 6-DOF trim; (2) a novel 6-DOF flight envelope model is developed with the steady states characterized by (h * , V * , .

ψ *
).The gliding footprint calculation process of a fixed-wing aircraft is illustrated in Figure 5, including the following steps.
ξ within the range [0°, 180°].During the computation, the straight segment is at b g , str and the turn segment is at every right turning state in the flight envelope.
Step 4. Extract the farthest landing point along every discrete ξ from the computed feasible gliding distances.Corresponding landing points then make up the right-half footprint outer boundary and the best-turn states are also obtained.The left-half outer boundary is determined via the symmetry along the initial heading direction.
Step 5. Identify the simple connectivity of the gliding footprint with Equation ( 27).If the equation is satisfied, the gliding footprint is simply connected and consists of all the points within the outer boundary.Otherwise, the inner boundary of the gliding footprint needs further calculation, which is beyond the scope of this paper.Step 1. Calculate the 6-DOF unpowered-flight envelope.The envelope is discretized by setting proper ranges and step sizes for V and .ψ at given H 0 .Every steady state in the flight envelope corresponds to a set (H 0 , V, .ψ) and is computed by the 6-DOF trim using the aircraft data.
Step 2. Extract the maximum flight-path angle γ bg,str from the trim solutions of all steady straight states in the flight envelope, and calculate the turn radius of all steady turning states using Equations ( 5)- (7).The minimum radius R min also needs to be extracted from the above turn radii.
Step 3. Compute all feasible gliding distances along different radial directions with Equations ( 3) and (4).The radial directions are discretized by setting proper step size for ξ within the range [0 • , 180 • ].During the computation, the straight segment is at γ bg,str and the turn segment is at every right turning state in the flight envelope.
Step 4. Extract the farthest landing point along every discrete ξ from the computed feasible gliding distances.Corresponding landing points then make up the right-half footprint outer boundary and the best-turn states are also obtained.The left-half outer boundary is determined via the symmetry along the initial heading direction.
Step 5. Identify the simple connectivity of the gliding footprint with Equation ( 27).If the equation is satisfied, the gliding footprint is simply connected and consists of all the points within the outer boundary.Otherwise, the inner boundary of the gliding footprint needs further calculation, which is beyond the scope of this paper.

Simulation Results and Discussion
According to above footprint calculation method, simulation results including the 6-DOF unpowered-flight envelope and gliding footprint were generated and are discussed below.

Simulation Setup
In this study, simulations were conducted on a Learjet-24-like subsonic business jet (SBJ) under flaps-up configuration.We performed the simulation using a workstation with 3.70 GHz Intel Xeon eight-core processor under the Windows 7 operating system.The data of the SBJ and the simulation settings were as follows.
(1) Data of the SBJ The aerodynamic data of the SBJ were calculated using the CFD method, with flight conditions set with altitudes of 0, 5, and 10 km and Ma = 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, and 0.8, respectively.For other flight conditions, the aerodynamic data were obtained through interpolation.In [40], the experimental data of lift coefficient C L and drag coefficient C D are shown with Reynolds number Re = 8.6 × 10 6 .To verify the credibility of the CFD simulated results, the experimental data were chosen for comparison, as shown in Figure 6.Between the two groups of results, the lift curve slope error was less than 5% and the drag coefficient data were close to each other.Thus, the CFD results were considered reliable.
According to above footprint calculation method, simulation results including the 6-DOF unpowered-flight envelope and gliding footprint were generated and are discussed below.

Simulation Setup
In this study, simulations were conducted on a Learjet-24-like subsonic business jet (SBJ) under flaps-up configuration.We performed the simulation using a workstation with 3.70 GHz Intel Xeon eight-core processor under the Windows 7 operating system.The data of the SBJ and the simulation settings were as follows.
(1) Data of the SBJ The aerodynamic data of the SBJ were calculated using the CFD method, with flight conditions set with altitudes of 0, 5, and 10 km and Ma = 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, and 0.8, respectively.For other flight conditions, the aerodynamic data were obtained through interpolation.In [40], the experimental data of lift coefficient L C and drag coefficient D C are shown with Reynolds number Re = 8.6 ×10 6 .To verify the credibility of the CFD simulated results, the experimental data were chosen for comparison, as shown in Figure 6.Between the two groups of results, the lift curve slope error was less than 5% and the drag coefficient data were close to each other.Thus, the CFD results were considered reliable.The geometric and inertial data of the SBJ are shown in Table 1.The geometric and inertial data of the SBJ are shown in Table 1.The SBJ also satisfies some physical constraints as in Equation (21).The ranges of control surface deflections are given in [43], satisfying |δ e | ≤ 15 • , |δ a | ≤ 18 • , and |δ r | ≤ 30 • .The upper limit of α is defined as the maximum α in the linear segment of the lift coefficient to avoid aerodynamic stalls [25,28], being approximately 12 • , as shown in Figure 6.The lower limit of α is mainly used to trim the SBJ and set it at −12 • , which is much smaller than the trimmed solutions.Considering the structural limit represented by n max (Table 1), this paper also assumes that |φ| ≤ 60 • as in [26].Moreover, the range of γ is assumed to be [−40

Parameters
(2) Simulation settings During simulations on the SBJ, the initial altitude was 0.5 km.The other settings of the flight envelope and gliding footprint are shown in Table 2. Unless otherwise specified, the following settings will apply to all simulations.) was computed.For a specified initial altitude H 0 , the flight envelope is represented by a two-dimensional matrix (V * , .

ψ *
). to avoid aerodynamic stalls [25,28], being approximately 12°, as shown in Figure 6.The lower limit of α is mainly used to trim the SBJ and set it at −12°, which is much smaller than the trimmed solutions.Considering the structural limit represented by nmax (Table 1), this paper also assumes that φ ≤ 60° as in [26].Moreover, the range of γ is assumed to be (2) Simulation settings During simulations on the SBJ, the initial altitude was 0.5 km.The other settings of the flight envelope and gliding footprint are shown in Table 2. Unless otherwise specified, the following settings will apply to all simulations.

Flight envelope setting
Range of velocity (Ma) 0~0.8 Step size of velocity (Ma) 0.001 Range of turn rate ψ& (°/s)

Gliding footprint setting
Range of radial angle ξ (°) Step size of radial angle ξ (°) Figure 7 shows the 6-DOF unpowered-flight envelope results after sparse processing.The envelope contains a total of 752 thousand discrete steady states and is symmetrical along the zero ψ& line, showing equivalent maneuvering capability to turn right ( 0 and left ( 0 The maximum ψ & is 12.39 °/s at a velocity of 85.11 m/s (0.251 Ma), where the minimum turn radius R min 382.26m = is also obtained.The velocity range of the flight envelope is 59.34~261.10m/s (0.175~0.770Ma).As the velocity increases, the number of corresponding steady states shows a trend of first increasing, and then, decreasing.During the 6-DOF trim process, steady states should satisfy the inequality constraints of α, γ, φ, δ e , δ a , and δ r , as in Equation (21). Figure 8 illustrates some trimmed parameters of the steady states, trying to find the parameters constraining the envelope boundary.The boundary is divided into four sections, namely ABC, CD, DE, and AE.The leftmost section, ABC, is the minimum-velocity boundary at different .ψ and is constrained by the upper limit of α (in Figure 8a).The maximum lift force coefficients at different .ψ are obtained (in Figure 8d) in this section.For flight states with smaller velocity, the dynamic pressure will not be able to produce adequate lift force to trim the aircraft.The rightmost section DE is the maximum-velocity boundary and is limited by the lower limit of γ (the steepest flight path in Figure 8b).The maximum drag forces at different .ψ are obtained (in Figure 8e) in this section.For flight states with larger velocity, the drag force will be too great for the gravitational component along the flight-path angle to balance.The conditions of the above two sections are consistent with the conventional velocity-altitude flight envelope of powered aircraft.The remaining two sections, CD and AE, are restricted by the lower and upper limits of φ, respectively (in Figure 8c).The maximum lift forces at different velocities are gained (in Figure 8f).For flight states with bigger .ψ , the lift force component along φ will be insufficient to provide centripetal force during the turn.
upper limit of α (in Figure 8a).The maximum lift force coefficients at different ψ & are ob- tained (in Figure 8d) in this section.For flight states with smaller velocity, the dynamic pressure will not be able to produce adequate lift force to trim the aircraft.The rightmost section DE is the maximum-velocity boundary and is limited by the lower limit of γ (the steepest flight path in Figure 8b).The maximum drag forces at different ψ & are obtained (in Figure 8e) in this section.For flight states with larger velocity, the drag force will be too great for the gravitational component along the flight-path angle to balance.The conditions of the above two sections are consistent with the conventional velocity-altitude flight envelope of powered aircraft.The remaining two sections, CD and AE, are restricted by the lower and upper limits of φ, respectively (in Figure 8c).The maximum lift forces at different velocities are gained (in Figure 8f).For flight states with bigger ψ& , the lift force component along φ will be insufficient to provide centripetal force during the turn.The altitude loss in the turn segment affects the total gliding distance along the segmented trajectory and is determined by R, trn γ , and Δψ, as in Equation (1). Figure 9 shows the altitude loss over a complete Δψ = 360° turn for steady states in the flight envelope.For a constant velocity, the altitude loss decreases gradually with the increase in the turn rate ψ& , when the value of φ also increases (in Figure 8c).For a constant value of ψ & , such a trend also arises as the velocity decreases.The minimum altitude loss is 509.07 m at the steady states with maximum ψ& , when R min is also obtained simultaneously.When ψ& is near zero, the turning descent is close to straight gliding, and magnitude of the altitude loss is rather huge.The altitude loss in the turn segment affects the total gliding distance along the segmented trajectory and is determined by R, γ trn , and ∆ψ, as in Equation (1). Figure 9 shows the altitude loss over a complete ∆ψ = 360 • turn for steady states in the flight envelope.For a constant velocity, the altitude loss decreases gradually with the increase in the turn rate .ψ , when the value of |φ| also increases (in Figure 8c).For a constant value of .ψ, such a trend also arises as the velocity decreases.The minimum altitude loss is 509.07 m at the steady states with maximum .ψ , when R min is also obtained simultaneously.When .ψ is near zero, the turning descent is close to straight gliding, and magnitude of the altitude loss is rather huge.

Gliding Footprint
Figure 10 illustrates the segmented trajectories arriving at the footprint outer boundary, with a step size of 30 • for ξ.The gliding footprint is the area within the outer boundary that is discretized into 73 boundary points.For the right-half boundary, the turn segments are conducted with right turning, and for the left half, with left turning.The flight trajectory corresponding to zero ξ contains only a straight segment, while the 180 • ξ has two equivalent flight trajectories with the turn segment at left and right turnings, respectively.In addition, the values of the gliding distance and altitude loss in the straight segment are larger than those of the turn segment.To better illustrate the above segmented trajectories, several best-footprint parame ters were also analyzed.Figure 11 Δh (km) V (m/s)   To better illustrate the above segmented trajectories, several best-footprint parameters were also analyzed.Figure 11 shows the results of V bf , ∆ψ bf , µ bf , and R bf regarding the radial angle.In this figure, the values at ξ = 0 • are parameters of the best-glide state, and those at nonzero ξ are of the best-turn states.The range of V bf is 68.15~87.15m/s.Combined with the flight envelope results in Figure 7 where the velocity range is 59.34~261.10m/s, both the best-glide and best-turn states are found to be located in the low-velocity region of the envelope.Take the right-half footprint (positive ξ), for example; V bf first decreases, and then, increases with an increase in ξ, and the values of ∆ψ bf increase linearly.The maximum V bf is obtained at zero ξ.When ξ = 5  When ξ = 5°, 90°, and 180°, the results of R bf are 12.76, 0.55, and 0.49 km.The minimum R bf is obtained when ξ = 180° and is still larger than R min 382.26m = .To investigate the effects of 6-DOF flight dynamics on the gliding footprint, we compared the gliding footprint results between the proposed method and the method in [26].
As a critical factor influencing the gliding distance, the flight-path angle was first analyzed.Figure 12 shows the comparative results of flight-path angles in attainable steady straight states, which were computed using two trim methods.The point-mass-based trim computes str γ by balancing the forces, while the 6-DOF trim considers both force and moment equilibrium.The values of str γ in the 6-DOF trim at different velocities are all smaller than those of the former method, implying a steeper descent path.The best-glide angle b g , str γ of the former method is −4.386° at a velocity of 82.64 m/s, which is larger than the value −4.738° of the latter method at a velocity of 87.15 m/s.The absolute and relative errors of b g , str γ are 0.352° and 8.03%, respectively.Then, according to Equation (2), the best-glide distance in the straight segment (i.e., the maximum straight gliding distance) of the two trim methods are 13.04 and 12.07 times the altitude loss.The relative error of the maximum straight gliding distance is 7.44%.Such an error can lead to significant deviations in the gliding footprint, as demonstrated below.To investigate the effects of 6-DOF flight dynamics on the gliding footprint, we compared the gliding footprint results between the proposed method and the method in [26].
As a critical factor influencing the gliding distance, the flight-path angle was first analyzed.Figure 12 shows the comparative results of flight-path angles in attainable steady straight states, which were computed using two trim methods.The point-mass-based trim computes γ str by balancing the forces, while the 6-DOF trim considers both force and moment equilibrium.The values of γ str in the 6-DOF trim at different velocities are all smaller than those of the former method, implying a steeper descent path.The best-glide angle γ bg,str of the former method is −4.386 • at a velocity of 82.64 m/s, which is larger than the value −4.738 • of the latter method at a velocity of 87.15 m/s.The absolute and relative errors of γ bg,str are 0.352 • and 8.03%, respectively.Then, according to Equation (2), the best-glide distance in the straight segment (i.e., the maximum straight gliding distance) of the two trim methods are 13.04 and 12.07 times the altitude loss.The relative error of the maximum straight gliding distance is 7.44%.Such an error can lead to significant deviations in the gliding footprint, as demonstrated below.The comparative results of the gliding footprint are shown in Figure 13 and Table 3.
According to Table 3, the straight segment in the proposed method is steeper ( b g , str The comparative results of the gliding footprint are shown in Figure 13 and Table 3.According to Table 3, the straight segment in the proposed method is steeper (γ bg,str = −4.738• ) than that of the method in [26] (γ bg,str = −4.386• ).When ξ = 0 • , the maximum gliding distances d ξmax are 6.52 and 6.03 km, respectively.The relative error is 7.44% and is consistent with the above results in Figure 12.With the increase in ξ, the two groups of d ξmax both decrease obviously, which can also be seen in Figure 13.When ξ = 0 • , 45 • , 90 • , 135 • , and 180 • , the absolute errors of d ξmax between the two methods are 0.49, 0.52, 0.61, 0.73, and 0.90 km.These errors present an increasing trend, which is the result of the combined effects of the gliding distance deviations in both the straight and turn segments.The gliding footprint areas of the two methods are 93.48 and 74.05 km 2 , respectively, with a relative error of 20.79%.Therefore, there are significant deviations in both the maximum gliding distances and gliding footprints between the two methods, which cannot be ignored.The comparative results of the gliding footprint are shown in Figure 13 and Table 3.
According to Table 3, the straight segment in the proposed method is steeper ( b g , str γ = −4.738°)than that of the method in [26] ( b g , str γ = −4.386°).When ξ = 0°, the maximum gliding distances max ξ d are 6.52 and 6.03 km, respectively.The relative error is 7.44% and is consistent with the above results in Figure 12.With the increase in ξ, the two groups of max ξ d both decrease obviously, which can also be seen in Figure 13.When ξ = 0°, 45°, 90°, 135°, and 180°, the absolute errors of max ξ d between the two methods are 0.49, 0.52, 0.61, 0.73, and 0.90 km.These errors present an increasing trend, which is the result of the combined effects of the gliding distance deviations in both the straight and turn segments.The gliding footprint areas of the two methods are 93.48 and 74.05 km 2 , respectively, with a relative error of 20.79%.Therefore, there are significant deviations in both the maximum gliding distances and gliding footprints between the two methods, which cannot be ignored.Lateral distance (km) Method in [26] The proposed method Figure 13.Gliding footprint comparison of the method in [26] and the proposed method.To further analyze how the trim method affects the gliding footprint, we also compared the flight-path angles in the best-turn states.The comparative results regarding the radial angle ξ are shown in Figure 14.With the increase in ξ, the two groups of γ bf show a consistent trend with each other.When ξ = 180 • , the absolute error of γ bf between the two methods reaches 0.908 • , which is considerably larger than that of the best-glide state (0.352 • ).The values of γ bf in the proposed method are smaller over the whole range of ξ.Therefore, the farthest segmented trajectories arriving at the footprint outer boundary of the proposed method are steeper in both the straight and turn segments than those of the method in [26], which directly leads to deviations in the gliding footprint.
consistent trend with each other.When ξ = 180°, the absolute error of bf γ between the two methods reaches 0.908°, which is considerably larger than that of the best-glide state (0.352°).The values of bf γ in the proposed method are smaller over the whole range of ξ.
Therefore, the farthest segmented trajectories arriving at the footprint outer boundary of the proposed method are steeper in both the straight and turn segments than those of the method in [26], which directly leads to deviations in the gliding footprint.Combining the analysis of the results in Figures 12-14, it is proven that the 6-DOF flight dynamics used to model the unpowered aircraft does have a significant effect on the gliding footprint.The reason is that the 6-DOF trim considers more constraints (including the moment equilibrium, control surface saturation, and aerodynamic contributions of the control surfaces) when computing the gliding flight parameters, resulting in steeper segmented trajectories to the footprint outer boundary and a smaller gliding footprint under a given initial altitude.However, it takes about 6 h to calculate the discrete boundary points encircling the gliding footprint.Thus, the proposed method is unsuitable for direct application to in-flight damage conditions.

Footprint Estimation Based on BP-ANN
According to the simulation results in Section 3, the relationships between initial altitude H 0 and footprint boundary points ( , d max ξ ξ ) are complex and difficult to build directly.It costs a lot of computation to obtain even a single boundary point.To satisfy the real-time requirement of in-flight damage conditions, this section further develops a footprint estimation method based on BP-ANN.

Network Architecture
In this study, the BP-ANN model was adopted instead of simulation to construct the relationship between the initial altitude and footprint boundary points.A typical BP-ANN model is shown in Figure 15, including the input layer, the hidden layer, and the output layer.The Mean Square Error (MSE) and correlation coefficient (R) were used to evaluate Method in [26] The proposed method Combining the analysis of the results in Figures 12-14, it is proven that the 6-DOF flight dynamics used to model the unpowered aircraft does have a significant effect on the gliding footprint.The reason is that the 6-DOF trim considers more constraints (including the moment equilibrium, control surface saturation, and aerodynamic contributions of the control surfaces) when computing the gliding flight parameters, resulting in steeper segmented trajectories to the footprint outer boundary and a smaller gliding footprint under a given initial altitude.However, it takes about 6 h to calculate the discrete boundary points encircling the gliding footprint.Thus, the proposed method is unsuitable for direct application to in-flight damage conditions.

Footprint Estimation Based on BP-ANN
According to the simulation results in Section 3, the relationships between initial altitude H 0 and footprint boundary points (d ξmax , ξ) are complex and difficult to build directly.It costs a lot of computation to obtain even a single boundary point.To satisfy the real-time requirement of in-flight damage conditions, this section further develops a footprint estimation method based on BP-ANN.

Network Architecture
In this study, the BP-ANN model was adopted instead of simulation to construct the relationship between the initial altitude and footprint boundary points.A typical BP-ANN model is shown in Figure 15, including the input layer, the hidden layer, and the output layer.The Mean Square Error (MSE) and correlation coefficient (R) were used to evaluate the performance of the network [27].The three-layer BP-ANN has been proved to approximate any complex rational function [34].It has strong ability to fit the non-linear mapping and can predict the maximum gliding distance d ξmax under different pairs of (H 0 , ξ).In the network of this paper, there are two input units and an output unit.The input parameters are the initial altitude H 0 and radial angle ξ.The output parameter is the maximum gliding distance along radial angle ξ.To determine the number of hidden units that are very significant to the performance of the network, a trial-and-error process [32] was utilized.The number was first set to 2 during the network training, and then, adjusted In the network of this paper, there are two input units and an output unit.The input parameters are the initial altitude H 0 and radial angle ξ.The output parameter is the maximum gliding distance along radial angle ξ.To determine the number of hidden units that are very significant to the performance of the network, a trial-and-error process [32] was utilized.The number was first set to 2 during the network training, and then, adjusted continuously.When it was adjusted to 4, the predicted results were excellent.Therefore, the architecture of the neural network in this paper is determined to be 2-4-1.

Network Training
Using the method in Section 3, we computed the gliding footprint with H 0 = 0.5, 1, 2, and 3 km and employed them to train the network.The step of ξ was 5 • , and there were 292 groups of sample points (H 0 , ξ, d ξmax ) in total.These sample data were randomly divided into the training, validation, and test sets with proportions of 70%, 15%, and 15%, respectively.The training set was used to train the neural network.The validation set was used to monitor the generalization ability of the network and stop the iteration training when the ability stopped improving.The test set could evaluate the performance of the trained network.The active function of the hidden layer was f(x) = (1 + e −x ) −1 and that of the output layer was linear.The weights were initialized between [0,1].The initial learning rate was 0.001 and the maximum epoch number was 1000.The Levebverg-Marquardt method was selected as the training method.Figure 16  We also conducted a linear regression analysis on the predicted and target v d max ξ .When R is close to 1, the two groups of values will be regarded as strongl lated.Figure 17 shows the results for the training, validation, and test.The R va all greater than 0.99, and there is almost a linear relationship between the BP-AN mated results and targets.We also conducted a linear regression analysis on the predicted and target values of d ξmax .When R is close to 1, the two groups of values will be regarded as strongly correlated.Figure 17 shows the results for the training, validation, and test.The R values are all greater than 0.99, and there is almost a linear relationship between the BP-ANN estimated results and targets.
To test the prediction effect of the constructed BP-ANN model, six scenarios were set with the initial altitude H 0 selected randomly within 0.  We also conducted a linear regression analysis on the predicted and target values of . When R is close to 1, the two groups of values will be regarded as strongly correlated.Figure 17 shows the results for the training, validation, and test.The R values are all greater than 0.99, and there is almost a linear relationship between the BP-ANN estimated results and targets.The range of relative errors for the footprint area is 0.01%~0.62%,which is also very small.Figure 18 shows the gliding footprint results of the three above scenarios, with H 0 selected to be 0.955, 1.579, and 2.777 km.The gliding footprints under each initial altitude are quite close to each other.There is little difference between the simulated and BP-ANN predicted gliding footprints.The range of relative errors for the footprint area is 0.01%~0.62%,which is also very small.Figure 18 shows the gliding footprint results of the three above scenarios, with H 0 selected to be 0.955, 1.579, and 2.777 km.The gliding footprints under each initial altitude are quite close to each other.There is little difference between the simulated and BP-ANN predicted gliding footprints.

Estimation Results with BP-ANN
Using the trained network, we estimated the maximum gliding distance

Estimation Results with BP-ANN
Using the trained network, we estimated the maximum gliding distance d ξmax for different pairs (H 0 , ξ).Then, we also investigated the initial altitude effects on the footprint.
Figure 19 shows the estimation results of d ξmax regarding radial angle and initial altitude.The ranges of ξ and H 0 are −180~180 • and 0.5~3 km, and the step sizes are 5 • and 0.1 km, respectively.There are 1898 groups of boundary points (H 0 , d ξmax , ξ) in total.The minimum d ξmax is 2.62 km, obtained when H 0 = 0.5 km and ξ = ±180 • .The maximum value is 36.01km when H 0 = 3 km and ξ = 0 • .The estimation costs a total of 0.072 s and the average time for each initial altitude is 0.003 s.Compared with the flight envelope-based method that requires 6 h to obtain a gliding footprint, this method significantly improves the computational efficiency and can obtain the footprint in real time.To investigate the initial altitude effects, we took 0.6, 1.2, 1.8, and 2.4 km as the initial altitudes.The data of the boundary points are extracted from Figure 19, and the gliding footprint results are shown in Figure 20.Take the right half footprint, for example; the maximum gliding distance decreases as ξ increases for a given initial altitude, and the decreasing trend gradually weakens as the initial altitude increases.This weakening can be seen by comparing the maximum gliding distances along different radial angles.When H 0 = 0.6 km, the maximum gliding distance along ξ = 180° is 53.42% of that along ξ = 0, while the ratio increases to 86.92% when H 0 = 2.4 km.The reason is that the higher the initial altitude, the smaller the ratios of the gliding distance and corresponding altitude loss of the turn segment relative to the whole landing trajectory.Additionally, with the increase in the initial altitude, the shape of footprint outer boundary gradually becomes closer to a circle.When H 0 = 0.6 km, the distance between the boundary points along ξ = 0° and 180° is 88.84% of the footprint projection length in the lateral direction, while this ratio increases to 97.40% when H 0 = 2.4 km and the gliding footprint is rather near a circle.To investigate the initial altitude effects, we took 0.6, 1.2, 1.8, and 2.4 km as the initial altitudes.The data of the boundary points are extracted from Figure 19, and the gliding footprint results are shown in Figure 20.Take the right half footprint, for example; the maximum gliding distance decreases as ξ increases for a given initial altitude, and the decreasing trend gradually weakens as the initial altitude increases.This weakening can be seen by comparing the maximum gliding distances along different radial angles.When H 0 = 0.6 km, the maximum gliding distance along ξ = 180 • is 53.42% of that along ξ = 0, while the ratio increases to 86.92% when H 0 = 2.4 km.The reason is that the higher the initial altitude, the smaller the ratios of the gliding distance and corresponding altitude loss of the turn segment relative to the whole landing trajectory.Additionally, with the increase in the initial altitude, the shape of footprint outer boundary gradually becomes closer to a circle.When H 0 = 0.6 km, the distance between the boundary points along ξ = 0 • and 180 • is 88.84% of the footprint projection length in the lateral direction, while this ratio increases to 97.40% when H 0 = 2.4 km and the gliding footprint is rather near a circle.
while the ratio increases to 86.92% when 0 = 2.4 km.The reason is that the higher the initial altitude, the smaller the ratios of the gliding distance and corresponding altitude loss of the turn segment relative to the whole landing trajectory.Additionally, with the increase in the initial altitude, the shape of footprint outer boundary gradually becomes closer to a circle.When H 0 = 0.6 km, the distance between the boundary points along ξ = 0° and 180° is 88.84% of the footprint projection length in the lateral direction, while this ratio increases to 97.40% when H 0 = 2.4 km and the gliding footprint is rather near a circle.

Conclusions
In this paper, we propose two methods to calculate the gliding footprint of an aircraft with thrust failure.Method 1 models the aircraft with 6-DOF flight dynamics and determines the footprint based on the best-glide and best-turn states that are decided within the 6-DOF unpowered-flight envelope.Method 2 is based on the BP-ANN model, which is trained by a footprint database computed using Method 1.
The simulations show that 6-DOF flight dynamics has significant effects on the gliding footprint.Compared with the conventional point-mass dynamics method, Method 1 produces steeper trajectories to the footprint outer boundary and a smaller footprint area.The deviation in footprint area is 20.79%, and that in the maximum straight gliding distance is 7.44%.The results of Method 1 are relatively conservative, which improves the safety threshold of the landing sites.Both the best-glide and best-turn states are located in the low-velocity region of the flight envelope.
Using the BP-ANN model, Method 2 can predict the gliding footprint under the entire operational altitude and full radial direction while satisfying the real-time requirement of in-flight damage conditions.Given an initial altitude, the maximum gliding distance keeps decreasing as the radial angle increases.With the increase in the initial altitude, the decreasing trend weakens and the gliding footprint gradually becomes closer to a circle.
The proposed methods are significant for improving the forced landing capability of a thrust-failed aircraft, and hence, its survivability.Future work will include experimental flights to validate the proposed methods, gliding footprint calculation methods under additional constraints such as wind and obstacles, and effects of downward flap deflection on the gliding footprint.

( 2 )
We propose a 6-DOF unpowered-flight envelope model by defining three novel characterization parameters of the steady states.Each group of characterization parameters corresponds to a unique steady state.The flight envelope is compiled into a three-dimensional matrix that represents the full maneuvering capabilities of the aircraft.(3)We introduce BP-ANN to estimate the gliding footprint.The trained BP-ANN effectively constructs the relationship between the initial altitude and footprint boundary points, and predicts the gliding footprint in real time.

Figure 1 .
Figure 1.Gliding footprint of a fixed-wing aircraft with thrust failure.

Figure 1 .
Figure 1.Gliding footprint of a fixed-wing aircraft with thrust failure.

Figure 2 .
Figure 2. The research framework of this study.

ξ
along radial angle ξ.If we set a step size for ξ, then we can find a series of discrete boundary points around the aircraft.

Figure 2 .
Figure 2. The research framework of this study.

Aerospace 2024 ,of 24 Figure 3 .
Figure 3. Top view of feasible segmented trajectories landing along radial angle ξ.

Figure 3 .
Figure 3. Top view of feasible segmented trajectories landing along radial angle ξ.

Figure 3 .
Figure 3. Top view of feasible segmented trajectories landing along radial angle ξ.

Figure 4 .
Figure 4. Top view of the segmented trajectory landing on the right-half footprint along radial angel in x * and u * , we choose h * , V * , and .ψ * as the characterization parameters of the steady states.Each group of (h * , V * , .ψ * ) corresponds to a unique flight state (as derived in Section 3.2.2).Then, .h * , x * , and u * are actually functions of (h * , V * , .ψ * ), which can represent the trim state families.By setting proper ranges and step sizes for them, a three-dimensional matrix will be constructed.The 6-DOF trim is used to judge whether a state (h * , V * , .ψ * ) in the matrix is attainable.Then, the unpowered-flight envelope is obtained and consists of all attainable steady straight ( states.This flight envelope model differs from that of powered aircraft, which is characterized by (h * , V * , γ * , .ψ * ) [27-30].The reason is that the flight-path angle γ is used to trim the aircraft due to the loss of engine thrust control input and is no longer independent from h * , V * , and .ψ * .3.2.2.Six-Degree-of-Freedom Trim Procedure For a specified set (h * , V * , .ψ *

Figure 5 .
Figure 5. Footprint calculation process for fixed-wing aircraft based on the 6-DOF flight envelope.

Figure 7
shows the 6-DOF unpowered-flight envelope results after sparse processing.The envelope contains a total of 752 thousand discrete steady states and is symmetrical along the zero .ψ line, showing equivalent maneuvering capability to turn right ( .ψ > 0) and left ( .ψ < 0).The maximum .ψ is 12.39 • /s at a velocity of 85.11 m/s (0.251 Ma), where the minimum turn radius R min = 382.26m is also obtained.The velocity range of the flight envelope is 59.34-261.10m/s (0.175~0.770Ma).As the velocity increases, the number of corresponding steady states shows a trend of first increasing, and then, decreasing.

5 3.4. 2 .(&
Unpowered-Flight EnvelopeAccording to the footprint calculation process in Section 3.3, the 6-DOF unpoweredflight envelope characterized by h V was computed.For a specified initial altitude H 0 , the flight envelope is represented by a two-dimensional matrix V

Figure 8 .
Figure 8.Some trimmed parameters of steady states in the flight envelope: (a) angle of attack; (b) flight-path angle; (c) roll angle; (d) lift force coefficient; (e) drag force; (f) lift force.

Figure 8 .
Figure 8.Some trimmed parameters of steady states in the flight envelope: (a) angle of attack; (b) flight-path angle; (c) roll angle; (d) lift force coefficient; (e) drag force; (f) lift force.

Figure 9 .
Figure 9. Altitude loss over a complete Δψ = 360° turn for steady states in the flight envelope.
Figure 10 illustrates the segmented trajectories arriving at the footprint outer bound ary, with a step size of 30° for ξ.The gliding footprint is the area within the outer boundary that is discretized into 73 boundary points.For the right-half boundary, the turn segment are conducted with right turning, and for the left half, with left turning.The flight trajec tory corresponding to zero ξ contains only a straight segment, while the 180° ξ has two equivalent flight trajectories with the turn segment at left and right turnings, respectively In addition, the values of the gliding distance and altitude loss in the straight segment ar larger than those of the turn segment.

Figure 10 .
Figure 10.Segmented trajectories to the footprint outer boundary along different radial directions.
shows the results of V bf , bf ψ Δ , bf μ , and R bf regarding th radial angle.In this figure, the values at ξ = 0° are parameters of the best-glide state, and those at nonzero ξ are of the best-turn states.The range of V bf is 68.15~87.15m/s.Com bined with the flight envelope results in Figure 7 where the velocity range is 59.34~261.10m/s, both the best-glide and best-turn states are found to be located in the low-velocity region of the envelope.Take the right-half footprint (positive ξ), for example; V bf first de creases, and then, increases with an increase in ξ, and the values of bf ψ Δ increase linearlyThe maximum V bf is obtained at zero ξ.When ξ = 5°, 90°, and 180°, the results of bf

Figure 9 .
Figure 9. Altitude loss over a complete ∆ψ = 360 • turn for steady states in the flight envelope.

Figure 9 .
Figure 9. Altitude loss over a complete Δψ = 360° turn for steady states in the flight envelope.

Figure 10 .Figure 10 .
Figure 10.Segmented trajectories to the footprint outer boundary along different radial directions , 90 • , and 180 • , the results of ∆ψ bf are 5.57 • , 97.23 • , and 200.99 • , respectively.The values of ∆ψ bf are larger than the corresponding radial angles ξ.As the positive ξ increases, µ bf becomes larger, and R bf becomes smaller.When ξ = 5 • , 90 • , and 180 • , the results of R bf are 12.76, 0.55, and 0.49 km.The minimum R bf is obtained when ξ = 180 • and is still larger than R min = 382.26m. 5.57°, 97.23°, and 200.99°, respectively.The values of bf ψ Δ are larger than the corresponding radial angles ξ.As the positive ξ increases, bf μ becomes larger, and R bf becomes smaller.

Aerospace 2024 , 24 Figure 12 .
Figure 12.Flight-path angle comparison of the two trim methods in attainable steady straight states.

Figure 12 .
Figure 12.Flight-path angle comparison of the two trim methods in attainable steady straight states.

Figure 12 .
Figure 12.Flight-path angle comparison of the two trim methods in attainable steady straight states.

Figure 13 .
Figure13.Gliding footprint comparison of the method in[26] and the proposed method.

Figure 14 .
Figure 14.Comparison of the best-footprint flight-path angle regarding the radial angle.

Figure 14 .
Figure 14.Comparison of the best-footprint flight-path angle regarding the radial angle.

Aerospace 2024 ,
11, x FOR PEER REVIEW 18 of 24 the performance of the network [27].The three-layer BP-ANN has been proved to approximate any complex rational function [34].It has strong ability to fit the non-linear mapping and can predict the maximum gliding distance d max ξ under different pairs of ( , H 0 ξ ).

Figure 15 .
Figure 15.Architecture of the BP-ANN model.

Figure 15 .
Figure 15.Architecture of the BP-ANN model.
Figure 16.MSE of the BP-ANN model during the training process.

Figure 16 .
Figure 16.MSE of the BP-ANN model during the training process.
5~3.0 km.Using the flight envelope-based and BP-ANN model-based methods, six groups of gliding footprints were computed.The step size of ξ was 5 • , and each gliding footprint contained 73 discrete boundary points (d ξmax , ξ).The maximum relative error for the gliding distance d ξmax of all discrete boundary points under each initial altitude, and the relative error for the footprint area, were computed and are shown in Table 4.The range of maximum relative errors of d ξmax under different H 0 is 0.12%~1.55%.Considering the small relative errors, the results of d ξmax in the two methods match very well.

Figure 16 .
Figure 16.MSE of the BP-ANN model during the training process.

Figure 17 .Figure 17 .
Figure 17.Linear regression plots for training, validation, and test.To test the prediction effect of the constructed BP-ANN model, six scenarios were set with the initial altitude H 0 selected randomly within 0.5~3.0km.Using the flight envelope-based and BP-ANN model-based methods, six groups of gliding footprints were computed.The step size of ξ was 5°, and each gliding footprint contained 73 discrete boundary points ( max ,ξ Aerospace 2024, 11, x FOR PEER REVIEW 20 of 24 all discrete boundary points under each initial altitude, and the relative error for the footprint area, were computed and are shown in Table 4.The range of maximum relative errors of d max ξ under different H 0 is 0.12%~1.55%.Considering the small relative errors, the results of d max ξ in the two methods match very well.

Aerospace 2024 , 24 Figure 19 .
Figure 19.The maximum gliding distance regarding radial angle and initial altitude.

Figure 19 .
Figure 19.The maximum gliding distance regarding radial angle and initial altitude.

Table 2 .
Settings of flight envelope and gliding footprint.

Table 2 .
Settings of flight envelope and gliding footprint.

Table 3 .
Gliding footprint parameters of the two methods.

Table 4 .
Relative errors for gliding distance and footprint area between the two methods.

Table 4 .
Relative errors for gliding distance and footprint area between the two methods.